CN105427267A - 获取刃边法测量光学遥感载荷在轨mtf精度的方法 - Google Patents

Abstract

本发明提供了一种获取基于刃边法测量光学遥感载荷在轨MTF的测量精度的方法。该方法基于实际测试的MTF值设置MTF参考值，针对刃边法MTF测试特点，基于实际测试时载荷获取的刃边图像提取关键参数生成参考图像，利用刃边法对参考图像进行计算得到MTF计算值，然后将计算值与MTF参考值对比，获得绝对误差和相对误差，由此得到MTF的测量精度(或不确定度)，从而解决了现有精度分析方法无法获得或估计刃边法实际测量光学遥感载荷在轨MTF测量精度的问题。

Description

获取刃边法测量光学遥感载荷在轨MTF精度的方法

技术领域

本发明涉及对地观测遥感技术领域，特别涉及一种获取基于刃边法测量光学遥感载荷在轨MTF的测量精度的方法。

背景技术

调制传递函数MTF(ModulationTransferFunction)是评价光学遥感载荷成像性能的一个重要指标，它表示了各个空间频率信号经过光学成像系统后调制度损失的百分比，反映成像系统在对目标成像过程中信号的扩散与削弱程度，是目前国际上普遍使用的评定光学遥感载荷成像性能以及空间分辨率的指标之一。

根据选用靶标的不同，目前测量光学遥感载荷在轨MTF的方法主要有刃边法、三线靶标法、点源/点阵法、脉冲法、辐射状靶标法。在这些方法中，刃边法因其能够获得不同空间频率的MTF曲线、能够更加完整地刻画光学遥感载荷空间响应特性、在实际应用中靶标布设与选取条件相对宽松(既可以是人工布设的靶标，也可以选取合乎要求的刃边类型的地物目标)，因而是评估中高分辨率光学遥感载荷在轨MTF最普遍使用的方法。诸如SPOT5/6/7、IKONOS、QUICKBIRD、GEOEYE、ZY02/03-HR、GF1/2等高分辨率卫星的光学遥感载荷均采用刃边法进行在轨MTF评测。

为了提高刃边法评测载荷在轨MTF的精度，各卫星运行机构相继布设高质量的刃边靶标，基于刃边法的MTF在轨评测技术也日趋成熟。图1为国内外典型的高分辨率刃边靶标。其中，A为美国德克萨斯州BigSpring靶标；B为台湾澎湖靶标；法国SalondeProvince靶标；D为加拿大Mt.Albert，Ontario靶标。

图2为现有技术采用刃边法评测光学遥感载荷在轨MTF的数据处理流程图。具体处理流程如下：

步骤S101，提取刃边子图像：从经过相对辐射校正的含有刃边靶标或刃边地物(以下统称为刃边靶标)的图像数据中提取不小于15像元×15像元的有效计算区域(刃边子图像)；特别说明一下，根据刃边子图像中刃边的方向不同，可以分别测量得到载荷阵列方向(对应于所获取图像的x方向)的MTF和载荷运动方向(对应于所获取图像的y方向)的MTF，如图3所示，为了描述简单，后续描述如不特别说明，均针对载荷阵列方向的MTF，载荷运动方向的MTF计算与精度测量按照同样方法参照实现。

步骤S102，边缘探测：提取刃边子图像的每一行的灰度数据，采用一定的模型进行差值，利用微分法、质心法、Hough变换法、Fermi函数法等方法，确定该行数据灰度变换廓线的边缘位置到亚像元精度(刃边边缘的中心位置)。对每一行数据采用同样方法，获得刃边子图像每一行数据的亚像元精度的刃边边缘中心位置，然后通过直线拟合对探测到的边缘中心位置进行进一步调整，使探测到的各行边缘中心位置点在一条直线上，计算该拟合直线的斜率，可以得到刃边的倾角θ；

步骤S103，提取边扩散函数ESF(EdgeSpreadFunction)：利用步骤S102探测到的每一行的刃边边缘中心位置，对刃边子图像的各行数据进行配准、合并、插值处理，得到亚像元精度的ESF；

步骤S104，求取线扩散函数LSF(LineSpreadFunction)：对步骤S103得到的ESF进行差分，并采用窗函数进行截断处理，得到宽度为-5像元～+5像元(或-10像元～+10像元)亚像元精度的LSF；

步骤S105，MTF求取：对步骤S104得到的LSF进行傅里叶变换，对变换结果取模值并做归一化处理，即可得到MTF序列曲线和Nyquist频率下的MTF值MTF_Nyquist(对宽度为-5像元～+5像元的LSF，其MTF_Nyquist为归一化MTF序列中的第6个值；对宽度为-10像元～+10像元的LSF，其MTF_Nyquist为归一化MTF序列中的第11个值)。

作为一种普遍使用的光学遥感载荷在轨MTF测量方法，现有的文献更多地涉及到MTF评测方法的改进和优化；也有部分文献采用仿真方法分析了亮暗区域对比度、靶标倾角、计算区域尺寸、噪声等因素对于MTF算法评测精度的影响，但是上述分析并没有与载荷实际成像条件结合起来。比如，文献中通常假设刃边靶标的亮暗区域的灰度值分别为200、100，假设噪声的标准差为0.7、1、1.5，靶标倾角为5°；此外，文献中通常假设刃边亮暗区域的噪声水平是相同的，并且往往没有考虑随机相位的影响。这些假设往往与载荷实际成像条件及获取图像并不一致，因而无法真实反映光学遥感载荷在轨实际成像条件下刃边法评测MTF的精度。

发明内容

(一)要解决的技术问题

本发明提出了一种基于载荷实际成像条件获取刃边法测量光学遥感载荷在轨MTF测量精度的方法，以解决现有精度分析方法无法获得或估计刃边法实际测量光学遥感载荷在轨MTF测量精度的问题。

(二)技术方案

本发明获取基于刃边法测量光学遥感载荷在轨MTF在第一方向的测量精度的方法包括：

步骤A：在光学遥感载荷获取的实测刃边靶标图像中截取尺寸为m行n列的多个刃边子图像，由该多个刃边子图像获取光学遥感载荷在第一方向的在轨MTF参考值和刃边倾角参考值；

步骤B：由光学遥感载荷在第一方向的在轨MTF参考值，获取光学遥感载荷在该第一方向的扩散尺度，构建亚像元尺度下的参考图像退化模板psf(x，y)；

步骤C：由光学遥感载荷获取的实测刃边靶标图像提取被刃边分割的两个均匀区域的灰度均值和随机噪声标准差；

步骤D：根据两个均匀区域的随机噪声标准差，获得用于获取后续步骤中亚像元水平下刃边噪声参考图像的噪声参数：刃边噪声参考图像中两个区域的噪声标准差；

步骤E：利用刃边靶标参数，构建亚像元水平的第一方向的理想刃边参考图像I(x，y)，其中，所述刃边靶标参数包括：刃边子图像尺寸、靶标图像被刃边分割的两个均匀区域的灰度均值、光学遥感载荷在第一方向的刃边倾角参考值；

步骤F：将所述理想刃边参考图像I(x，y)与参考图像退化模板psf(x，y)进行二维卷积，截取卷积结果矩阵中有效数据，得到亚像元水平的第一方向的无噪声刃边退化图像I_psf(x，y)；

步骤G：对所述刃边退化图像I_psf(x，y)中的数据进行数据聚合，生成整像元水平的第一方向的无噪声刃边退化参考图像g(x，y)；

步骤H：利用所述刃边噪声参考图像中两个区域的噪声标准差，构建K幅亚像元水平的第一方向的带噪声刃边参考图像I_k(x，y)，k＝1，2，…K；

步骤I：将带噪声刃边参考图像I_k(x，y)k＝1，2，…K，分别与参考图像退化模板psf(x，y)进行二维卷积，截取卷积结果矩阵中的有效数据得到K幅亚像元水平的第一方向的带噪声刃边退化图像I_psf-k(x，y)，k＝1，2，…K；

步骤J：对所述带噪声刃边退化图像I_psf-k(x，y)，k＝1，2，…K，进行数据聚合处理，生成K幅整像元水平的阵列方向带噪声退化参考图像g_k(x，y)，k＝1，2，…K；以及

步骤K：利用刃边法MTF测试算法，对所述整像元水平阵列无噪声刃边退化图像g(x，y)，以及所述K幅整像元水平的阵列方向带噪声刃边退化参考图像g_k(x，y)，k＝1，2，…K，分别进行第一方向的MTF计算，将Nyquist频率处的归一化MTF计算值与第一方向的在轨MTF参考值进行对比，得到光学遥感载荷在轨MTF精度；

其中，所述第一方向为载荷阵列方向和载荷运动方向其中之一。

(三)有益效果

从上述技术方案可以看出，本发明法获取刃边法测量光学遥感载荷在轨MTF测量精度的方法具有以下有益效果：

(1)基于载荷实际成像条件下所获取的图像，提取刃边法在轨MTF评测中的关键参数(刃边亮暗区域灰度、刃边靶标倾角、MTF计算有效区域尺寸、噪声等)，并考虑刃边两侧区域随机噪声水平不同以及随机采样相位的影响，得到刃边法测量光学遥感载荷在轨MTF的测量精度，该测量精度与MTF在轨实际测量密切相关，可为MTF在轨实际测量结果提供更可靠的定量化精度和置信度量；

(2)基于实测条件获得各种MTF评测方法的测量精度，便于更有效对比不同MTF评测方法的优劣。

附图说明

图1为国内外典型的高分辨率刃边靶标；

图2为现有技术采用刃边法评测光学遥感载荷在轨MTF的数据处理流程图；

图3为现有技术从刃边靶标图像中提取的刃边子图像(MTF计算有效区域)的示意图；

图4为根据本发明实施例获取刃边法评测光学遥感载荷在轨MTF精度的方法流程图；

图5A为测量载荷阵列方向MTF的刃边子图像的刃边及倾角示意图。

图5B为测量载荷飞行方向MTF的刃边子图像的刃边及倾角示意图；

图6为按照步骤B生成的退化模板示意图(σ₁＝σ₂＝0.61244)；

图7为实施步骤E时参考图像P、原点O、倾角为θ_x-ref的刃边直线的示意图；

图8为从某星载光学载荷获取的阵列方向刃边靶标图像(方框内区域为用于阵列方向MTF计算的刃边子图像)；

图9为基于图8进行MTF计算及所获得图像参数，生成用于获取其MTF精度所需要的刃边参考图像的过程示意图。

具体实施方式

本发明基于实际测试的MTF值设置MTF参考值，针对刃边法MTF测试特点，基于实际测试时载荷获取的刃边图像提取关键参数生成参考图像，利用刃边法对参考图像进行计算得到MTF计算值，然后将计算值与MTF参考值对比，获得绝对误差和相对误差，由此得到MTF的测量精度。

在本发明的一个示例性实施例中，提供了一种获取刃边法测量光学遥感载荷在轨MTF在阵列方向的测量精度的方法。需要说明的是，本实施例是以载荷阵列方向的MTF，即MTF_x为例进行说明。关于载荷运动方向的MTF，即MTF_y可按照同样方法参照实现。

图4为根据本发明实施例获取刃边法评测光学遥感载荷在轨MTF精度的方法流程图。如图4所示，本实施例获取光学遥感载荷在轨MTF精度的方法包括：

步骤A：由光学遥感载荷获取实测刃边靶标图像，在该实测刃边靶标图像中截取的尺寸为m×n的多个刃边子图像在归一化频率下的MTF值MTF_x-Nyquist和刃边倾角θ_x经平均获取光学遥感载荷在阵列方向的在轨MTF参考值MTF_x-ref和刃边倾角θ_x-ref；

MTF_x-ref采用多次实际测试过程中测量得到的MTF值的平均值：即对测试中载荷获取的刃边靶标图像，采用刃边法计算载荷的归一化Nyqiust频率(f＝0.5)下的MTF值(MTF_x-Nyquist)作为MTF_x-ref。由于受噪声和计算区域等不确定因素影响，每次MTF_x-Nyquist计算的结果会有差别，MTF_x-ref取多次计算的平均值，步骤A又可以分为以下子步骤：

子步骤A1：在光学遥感载荷获取的实测刃边靶标图像中截取满足MTF计算要求的多个刃边子图像D；

其中，该多个刃边子图像D具有相同尺寸m行×n列(或接近相同尺寸)，后续将基于该尺寸m行×n列构建参考图像。多个刃边子图像D在刃边靶标图像的位置并没有特别的限定。优选地，在选取刃边子图像D时保证刃边位于刃边子图像D中间位置效果会更好一些。

子步骤A2：对于多个刃边子图像中的每一个，采用刃边法计算其在载荷阵列方向的归一化Nyquist频率下的MTF，即MTF_x-Nyquist，并在计算过程中获得刃边倾角θ_x；

其中，θ_x为刃边子图像中的刃边与图像列方向的夹角，范围为-90°～90°。图5A显示了测量载荷阵列方向MTF的4种刃边子图像中的刃边以及相应的刃边倾角θ_x。其中，图(a)和图(c)中由亮暗区域所形成的刃边向图像列方向偏左的方向倾斜，θ_x的取值范围为-90°～0°；图(b)和图(d)中由亮暗区域所形成的刃边向图像列方向偏右的方向倾斜，θ_x的取值范围为0°～90°。

对于测量载荷飞行方向MTF的刃边子图像，刃边倾角θ_y为子图像中的刃边与图像行方向的夹角，范围为-90°～90°。图5B显示了测量载荷飞行方向MTF的4种刃边子图像中的刃边以及相应的刃边倾角θ_y。其中，图(a)和图(c)中由亮暗区域所形成的刃边向图像行方向偏上的方向倾斜，θ_y的取值范围为0°～90°；图(b)和图(d)中由亮暗区域所形成的刃边向图像列方向偏下的方向倾斜，θ_y的取值范围为-90°～0°。

子步骤A3：对多个刃边子图像在载荷阵列方向的归一化Nyquist频率下的MTF，即MTF_x-Nyquist和刃边倾角θ_x分别求平均，获取光学遥感载荷在载荷阵列方向的在轨MTF参考值MTF_x-ref和刃边倾角参考值θ_x-ref。

步骤B：由光学遥感载荷在阵列方向的在轨MTF参考值MTF_x-ref，获取光学遥感载荷在该方向的扩散尺度σ_x，构建亚像元尺度(水平)下的参考图像退化模板psf(x，y)；

光学遥感载荷系统响应函数PSF通常采用高斯模型表示，结合光电系统数字采样效应，表示为：

$\begin{array}{c}PSF\left(x,y\right)=\frac{1}{2{\mathrm{\π\σ}}_x{\mathrm{\σ}}_y}{e}^{-\frac{x^2}{2{\mathrm{\σ}}_x^2}}{e}^{-\frac{y^2}{2{\mathrm{\σ}}_y^2}}*rect\left(x\right)*rect\left(y\right)\\ =\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_x}{e}^{-\frac{x^2}{2{\mathrm{\σ}}_x^2}}*rect\left(x\right)\×\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_y}{e}^{-\frac{y^2}{2{\mathrm{\σ}}_y^2}}*rect\left(y\right)\\ =PSF\left(x\right)\×PSF\left(y\right)\end{array}---\left(1-1\right)$

其中：

$PSF\left(x\right)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_x}{e}^{-\frac{x^2}{2{\mathrm{\σ}}_x^2}}*rect\left(x\right),PSF\left(y\right)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_y}{e}^{-\frac{y^2}{2{\mathrm{\σ}}_y^2}}*rect\left(y\right)---(1-2)$

其中，x，y分别表示光学遥感载荷在阵列方向和运动方向的坐标，σ_x、σ_y分别表示载荷系统函数PSF在阵列方向和运动方向的扩散尺度。则MTF参考值MTF_ref与该光学遥感载荷系统函数PSF的扩散尺度σ之间存在如下关系：

${\mathrm{MTF}}_{x-ref}=e^{-2{\mathrm{\π}}^2{{\mathrm{\σ}}_x}^2/4}\×\frac{2}{\mathrm{\π}}---(2-1)$

${\mathrm{MTF}}_{y-ref}=e^{-2{\mathrm{\π}}^2{{\mathrm{\σ}}_y}^2/4}\×\frac{2}{\mathrm{\π}}---(2-2)$

利用公式(2)，已知MTF参考值，即可得到光学遥感载荷系统函数PSF的关键参数——扩散尺度。比如，当MTF_x-ref＝0.1时，经由式2-1计算，对应的系统函数PSF的扩散尺度σ_x为0.61244。

由公式(1-1)可知，光学遥感载荷的系统响应函数PSF在阵列方向和运动方向是可以分离的，当光学遥感载荷在运动方向的在轨MTF参考值MTF_y-ref可以测量得到时，可以按照以上同样方法获得载荷运动方向的扩散尺度σ_y；如果无法测量MTF_y-ref时，可以令MIF_y-ref＝MTF_x-ref，从而有σ_y＝σ_x，将不会影响本发明在阵列方向MTF的测量精度。

基于系统函数PSF的扩散尺度σ_x、σ_y，可以生成亚像元水平的参考图像归一化退化模板为：

$psf(x,y)=e^{-\frac{x^2}{2{\mathrm{\σ}}_x^2}}{e}^{\frac{y^2}{2{\mathrm{\σ}}_y^2}}---\left(3\right)$

本实施例中，参考图像退化模板的尺寸选为整像元水平下的10像元×10像元(-5像元≤x≤4.9像元，-5像元≤y≤4.9像元)，两数据点之间的间隔为0.1像元。图6为按照步骤B生成的亚像元水平参考图像退化模板psf(x，y)的示意图，图中，psf(x，y)的尺寸选为10像元×10像元(-5像元≤x≤4.9像元，-5像元≤y≤4.9像元)，扩散尺度σ_x＝σ_y＝0.61244，该模板对应的MTF_x-rey＝MTF_y-ref＝0.1。

需要说明的是，本发明并不以上述实施例为限，本领域技术人员可以合理设置参考图像退化模板的尺寸以及数据点之间的间隔。但需要注意的是，参考图像退化模板psf(x，y)中数据点之间的间隔应当介于0.01～0.1像元之间。

步骤C：由光学遥感载荷获取的实测刃边靶标图像提取被刃边分割的两个均匀区域的灰度均值和随机噪声标准差，其中，设该刃边靶标图像中刃边左侧均匀区域Z_L的平均灰度为μ_L，随机噪声标准差σ_noise-L；刃边右侧均匀区域Z_R的平均灰度为μ_R，Z_R的随机噪声标准差σ_noise-R；

本步骤中，上述实测刃边靶标图像中被刃边分割的两个均匀区域分为以下两种情况：当测量载荷阵列方向的MTF时，实测刃边靶标图像被刃边分割成左右两个均匀区域，刃边左侧为亮区域，右侧为暗区域，或者刃边左侧为暗区域，右侧为亮区域；当测量飞行方向的MTF时，图像被刃边分割成上下两个均匀区域，刃边上部为亮区域，下部为暗区域，或者刃边上部为暗区域，下部为亮区域。

本实施例中，按照阵列方向MTF测量进行描述。具体的方法是：基于光学遥感载荷获取的实测刃边靶标图像，在靶标图像中的刃边两侧的均匀区域中分别选取一个尽可能大的矩形区域(Z_L和Z_R)，为避免邻近区域的干扰，选取的区域应当满足：与图像中的刃边中心位置距离大于3个像元，与靶标的各个边缘的距离大于2个像元。分别计算这两个所选取矩形区域的灰度均值和标准差，得到该图像中刃边左右两侧均匀区域的灰度均值μ_L、μ_R和区域Z_L、Z_R的随机噪声标准差σ_noise-L、σ_noise-R。该过程对于本领域技术人员是公知的，在相关的技术文件中也可以得到，此处不再进行赘述。

步骤D：根据刃边左侧均匀区域Z_L(m_ZL行×n_ZL列)的随机噪声标准差σ_noise-L和刃边右侧均匀区域Z_R(m_ZR行×n_ZR列)的随机噪声标准差σ_noise-R，获得用于构建后续步骤中亚像元水平下刃边噪声参考图像的噪声参数：刃边噪声参考图像中刃边左右两个区域的噪声标准差σ_nL和σ_nR；

其中，获取刃边噪声参考图像刃边左区域的噪声标准差σ_nL可以分为以下子步骤：

子步骤D1：构建一个均值为0、标准差σ_INZL的初始值为σ_noise-L、尺寸为[(m_ZL+10)×10行]×[(n_ZL+10)×10列]的亚像元水平的噪声矩阵IN_ZL(x，y)，该噪声矩阵IN_ZL(x，y)对应的数据点间距为0.1像元；

子步骤D2：令噪声矩阵IN_ZL(x，y)与步骤B得到的参考图像退化模板Psf(x，y)进行卷积：IN_ZL(x，y)*psf(x，y)，截取卷积结果矩阵中最中间的有效数据((m_ZL×10+1)行×(n_ZL×10+1)列)的数据(在Matlab中，使用卷积函数conv2的‘valid’参数)IPN_ZL；

子步骤D3：对矩阵IPN_ZL中的数据每10×10个数据为一组进行聚合，得到一个m_ZL行×n_ZL列的矩阵PN_ZL，PN_ZL数据点间距为1像元；

子步骤D4：计算PN_ZL的噪声标准差σ_PNZL；

子步骤D5：对比σ_PNZL和σ_noise-L，如果σ_PNZL不等于σ_noise-L，则调整子步骤D1中噪声矩阵IN_ZL(x，y)的标准差σ_INZL，然后重复子步骤D2～D5，直到σ_PNZL＝σ_noise-L。

则此时步骤D1中噪声矩阵IN_ZL(x，y)的标准差σ_INZL，即为后续构建亚像元水平下刃边参考图像刃边左区域噪声标准差σ_nL。

采用同样的方法，获取刃边参考图像刃边右区域的噪声标准差σ_nR。

步骤E：利用步骤A、B、C、D中提取的刃边靶标参数，包括：刃边子图像尺寸，即行数m和列数n；靶标图像中被刃边分割的左右两侧均匀区域的灰度均值μ_L、μ_R；光学遥感载荷在载荷阵列方向的刃边倾角参考值θ_x-ref，构建亚像元水平的阵列方向的理想刃边参考图像I(x，y)；

步骤E又可以分为以下子步骤：

子步骤E1：基于步骤A提取的刃边子图像D的尺寸(m行n列，数据点间隔为1像元，构建一个(m+1+10)×10行，(n+1+10)×10列的亚像元水平的参考图像P，参考图像P的数据点(像素点)间距为0.1像元。

子步骤E2，在参考图像P中，令与图像P的列方向(如图7所示的y轴方向)成θ_x-ref角度的一条直线l，从参考图像P中第 像素的左上角点(以下简称O点)斜穿而过，则根据该直线l穿越P中某一像素时，分割该像素为左右两个部分，设参考图像P中各个像素的面积为1，根据P中各个像素被直线l分割的左右两部分的面积为权重，对参考图像P进行灰度赋值：

(1)对于P中没有被直线l穿过的像素p(x，y)，如果p(x，y)位于直线l左侧，则p的灰度I(x，y)赋值为：

I(x，y)＝μ_L

(2)对于参考图像P中没有被直线l穿过的像素p(x，y)，如果p(x，y)位于直线l右侧，则该p(x，y)的灰度I(x，y)赋值为：

I(x，y)＝μ_R

(3)对于被直线l穿过的像素p(x，y)，则计算直线l穿过p时将p分割成为左右两部分，所占面积分别为S₁、S₂(S₁+S₂＝1)，则该像素p(x，y)的灰度I(x，y)赋值为：

I(x，y)＝S₁·μ_L+S₂·μ_R

图7显示了实施步骤E时构建的参考图像P、原点O、倾角为θ_x-ref的刃边直线的示意图。图中，O点四周的四个像素分别为O₁、O₂、O₃、O₄，O₁没有被直线l穿过，且位于l左侧，所以O₁的灰度值为μ_L；O₂被直线l穿过，设O₂面积为1，直线l分割O₂左右两部分的面积分别为：tgθ_x-ref和1-tgθ_x-ref，则O₂的灰度值为μ_L·tgθ_x-ref+μ_R·(1-tgθ_x-ref)；O₃同样被直线l穿过，直线l分割O₃左右两部分的面积分别为：1-tgθ_x-ref和tgθ_x-ref，则O₃的灰度值为μ_L·(1-tgθ_x-ref)+μ_R·tgθ_x-ref；O₄没有被直线l穿过，且位于l右侧，所以O₄的灰度值为μ_R。

由此，构建出一幅亚像元水平的载荷阵列方向理想刃边参考图像I(x，y)。

步骤F：将步骤E得到的亚像元水平的阵列方向的理想刃边参考图像I(x，y)与步骤B得到的参考图像退化模板psf(x，y)进行二维卷积，截取卷积结果矩阵中有效数据(在Matlab中，使用卷积函数conv2的‘valid’参数)，得到(m+1)×10+1行、(n+1)×10+1列的亚像元水平的阵列方向的无噪声刃边退化图像I_psf(x，y)；

步骤G：对步骤F得到的亚像元水平的阵列方向的刃边退化图像I_psf(x，y)，以I_psf(x，y)中位于行、列的像素为起始位置，在图像I_psf(x，y)中上下左右每10×10个像素为1组进行数据聚合，生成整像元水平的阵列方向的无噪声刃边退化参考图像g(x，y)；

步骤H：利用步骤D得到的亚像元水平下刃边参考图像刃边左右两侧区域的随机噪声标准差σ_nL、σ_nL，构建K幅(K通常取100)亚像元水平的阵列方向的带噪声刃边参考图像I_k(x，y)，k＝1，2，…K；

步骤H又可分为以下子步骤：

子步骤H1：按照子步骤E1所述方法，再次构建K幅(m+1+10)×10行、In+1+10)×10列的亚像元水平的参考图像P_k，k＝1，2，…K，参考图像P_k的数据点(像素点)间距为0.1像元。

子步骤H2：分别构建K个均值为0、标准差σ_nL、尺寸为(m+1+10)×10行、(n+1+10)×10列的噪声矩阵NPL_k，k＝1，2，…K，和K个均值为0、标准差σ_nR、尺寸为(m+1+10)×10行、(n+1+10)×10列的噪声矩阵NPR_k，k＝1，2，…K，噪声矩阵NPL_k和NPR_k的数据点间距为0.1像元；

子步骤H3，按照类似子步骤E2的方法，在参考图像P_k，k＝1，2，…K中，令与参考图像P_k的列方向(如图7所示的y轴方向)成θ_x-ref角度的一条直线l从参考图像P_k中第像素的左上角点(以下简称O点)斜穿而过，则根据该直线l穿越P_k中某一像素时，分割该像素为左右两个部分，设参考图像P_k中各个像素的面积为1，根据参考图像P_k中各个像素被直线l分割的左右两部分的面积为权重，对参考图像P_k进行灰度赋值：

(1)对于P_k中没有被直线l穿过的像素p_k(x，y)，如果p_k(x，y)位于直线l左侧，则p_k(x，y)的灰度I_k(x，y)赋值为：

I_k(x，y)＝μ_L+NPL_k(x，y)；

(2)如果p_k(x，y)位于直线l右侧，则该P_k(x，y)的灰度I_k(x，y)赋值为：

I_k(x，y)＝μ_R+NPR_k(x，y)。

(3)对于被直线l穿过的像素p_k(x，y)，则计算直线l穿过p_k时将p_k分割成为左右两部分，所占面积分别为S_1k、S_2k(S_1k+S_2k＝1)，则该像素p_k(x，y)的灰度I_k(x，y)赋值为：

I_k(x，y)＝S_1k·[μ_L+NPL_k(x，y)]+S_2k·[μ_R+NPR_k(x，y)]

由此，构建出K幅亚像元水平的阵列方向带噪声刃边参考图像I_k(x，y)，k＝1，2，…K。

步骤I：将步骤H得到的亚像元水平的阵列方向带噪声刃边参考图像I_k(x，y)，k＝1，2，…K，分别与步骤B得到的参考图像退化模板psf(x，y)进行二维卷积，截取卷积结果矩阵中最中间的有效数据(在Matlab中，使用卷积函数conv2的‘valid’参数)，得到K幅(m+1)×10+1行、(n+1)×10+1列的亚像元水平的阵列方向带噪声刃边退化图像I_psf-k(x，y)，k＝1，2，…K；

步骤J：对步骤I得到的亚像元水平的阵列方向的带噪声刃边退化图像I_psf-k(x，y)，k＝1，2，…K，进行数据聚合处理，生成K幅整像元水平的阵列方向带噪声退化参考图像g_k(x，y)，k＝1，2，…K；

需要注意的是，采样相位是MTF测量中另一个必须考虑的随机因素。载荷在对地物成像并以像元为单位进行数字化采样时，在不同的采样相位下，获取的图像灰度值会不同，并最终影响MTF的计算结果，因而在步骤J生成整像元水平的带噪声退化参考图像g_k(x，y)k＝1，2，…K时，需要将随机采样相位因素引入。步骤K又可以分为以下子步骤：

子步骤J1：在(-5、-4、-3、-2、-1、0、1、2、3、4、5)这11个数中每次随机选出2个数，共随机选出K对数值(S_xk，S_yk)k＝1，2，…K，作为随机相位的偏移值；

子步骤J2：对第k幅亚像元水平的阵列方向带噪声刃边退化图像I_psf-k(x，y)，k＝1，2，…K，以I_psf-k(x，y)中第行、列的像素为聚合起始位置，在图像I_psf-k(x，y)中上下左右每10×10个像素为1组进行聚合，即对该10×10个像素组内所有像素数据点的灰度数值累加求平均并四舍五入取整，最终得到K幅带有随机噪声和随机相位的、间距为1像元的整像元水平的阵列方向带噪声退化刃边参考图像g_k(x，y)，k＝1，2，…K。

步骤G得到的整像元水平的阵列方向无噪声刃边退化图像，没有包含随机噪声和随机采样相位造成的随机误差，可以用于分析光学遥感载荷在轨MTF阵列方向的MTF系统误差。

步骤J得到的整像元水平的阵列方向带噪声刃边退化参考图像g_k(x，y)，k＝1，2，…K，综合了MTF实测过程中刃边靶标亮暗反射区域不同辐亮度反射、刃边靶标倾角、随机噪声、CCD相机随机采样相位、CCD相机成像时的模糊退化效应、MTF计算图像区域尺寸等多重关键影响因素，因而使用MTF刃边算法针对g_k(x，y)，k＝1，2，…K，计算得到的MTF值，最能够反映实际测量的情况。

图8和图9显示了按照上述步骤A～步骤J生成用于获取实测某星载光学载荷阵列方向MTF的测量精度时，所生成的一幅带噪声刃边参考图像的过程。图8显示了某星载光学载荷获取的刃边靶标图像，图中方框内的区域为用于测量阵列方向MTF的刃边子图像。通过该图像，按照步骤A得到MTF_x-ref＝0.070，θ＝8.75°，子图像尺寸m＝40，n＝40；按照步骤C，得到μ_L＝804.25，μ_R＝191.68，σ_noise-A＝8.79，σ_{noise_B}＝6.16；按照步骤D，得到σ_nL＝130，σ_nR＝190。图9为基于图8进行MTF测量及所获得图像参数，生成用于获取其MTF精度所需要的刃边参考图像的过程示意图。其中，图9a为按照步骤E基于上述参数(μ_L＝804.25，μ_R＝191.68，θ_x-ref＝8.75°，m＝40，n＝40)构建的亚像元水平阵列方向理想刃边参考图像I(x，y)，图像尺寸为510×510；图9b为按照步骤H得到的1幅亚像元水平下有噪声参考图像I₁(x，y)，图像尺寸为510×510；图9c为按照步骤I得到的1幅亚像元水平下有噪声退化参考图像，图像尺寸为411×411；图9d为按照步骤J，经过数据聚合并考虑随机相位得到的整像元水平下有噪声退化参考图像，图像尺寸为40×40。

载荷飞行方向的无噪声和带噪声刃边退化参考图像可以参照同样的方法得到。

步骤K：利用刃边法MTF测试算法，对步骤G得到的整像元水平阵列无噪声刃边退化图像g(x，y)，以及步骤J生成的K幅整像元水平的阵列方向带噪声刃边退化参考图像g_k(x，y)，k＝1，2，…K，分别进行阵列方向的MTF计算，然后将Nyquist频率处的MTF计算值MTF_cal与步骤A中的阵列方向的MTF参考值MTF_x-ref进行对比，得到光学遥感载荷在轨MTF精度。

本步骤中，由于载荷阵列方向和飞行方向精度计算方法相同，为描述方便，以下将MTF计算值(MTF_x-cal、MTF_y-cal)统一表示为MTF_cal，MTF参考值(MTF_x-ref、MTF_y-ref)统一表示为MTF_ref。

该获取MTF测量精度的过程又可以分为以下两种情况：

一、不考虑随机因素的情况

如果不考虑测量中的随机因素(噪声和采样相位影响)，基于步骤G得到的无噪声退化参考图像g(x，y)，和步骤A得到MTF参考值MTF_ref，得到MTF测量的系统误差。

即基于步骤G得到的退化参考图像g(x，y)，按照图1所示的刃边法MTF数据处理流程计算MTF值得到MTF_cal，然后以步骤A得到的MTF参考值MTF_ref为基准，按照公式(13)计算相对误差，即为系统误差。

$\mathrm{\Δ}=\frac{{\mathrm{MTF}}_{cal}-{\mathrm{MTF}}_{ref}}{{\mathrm{MTF}}_{ref}}\×100\%---\left(13\right)$

二、考虑随机因素的情况

考虑测量中的随机噪声和采样相位等随机因素，基于步骤J得到的K幅带噪声刃边参考图像g_k(x，y)，和步骤A得到MTF参考值MTF_ref，采用统计方法计算MTF整体测量精度。

即对步骤J生成的K幅考虑随机噪声和采样相位的刃边随机参考图像g_k(x，y)，k＝1，2，…K，分别利用图1所示的刃边法计算MTF，并提取出Nyquist频率下的MTF值(MTF_k-cal，k＝1，2，…，K，K常取100或1000)，与步骤A的MTF_ref进行比较，按照下式得到：相对误差、平均相对误差和误差标准差。

相对误差：

${\mathrm{\Δ}}_k=\frac{{\mathrm{MTF}}_{k-cal}-{\mathrm{MTF}}_{ref}}{{\mathrm{MTF}}_{ref}}\×100\%,(k-1,2,...,K)---\left(13\right)$

平均相对误差：

$\begin{array}{c}\stackrel{\‾}{\mathrm{\Δ}}=avg(\frac{{\mathrm{MTF}}_{k-cal}-{\mathrm{MTF}}_{ref}}{{\mathrm{MTF}}_{ref}}\×100\%)\\ =\frac{1}{N}\underset{k=1}{\overset{K}{\Σ}}\left(\frac{{\mathrm{MTF}}_{k-cal}-{\mathrm{MTF}}_{ref}}{{\mathrm{MTF}}_{ref}}\right)\×100\%\end{array}---\left(14\right)$

均方根误差：

$\begin{array}{c}{\mathrm{\σ}}_{\mathrm{\Δ}}=std(\frac{{\mathrm{MTF}}_{k-cal}-{\mathrm{MTF}}_{ref}}{{\mathrm{MTF}}_{ref}}\×100\%)\\ =\sqrt{\frac{\underset{k=1}{\overset{K}{\Σ}}{({\mathrm{MTF}}_{k-cal}-\stackrel{\‾}{\mathrm{\Δ}})}^2}{K-1}}\×100\%\end{array}$

通常情况下，采用为在当前MTF测量下的平均测量误差，以σ_Δ作为MTF测量精度的1σ精度。

由于在MTF测量中存在随机因素，通常还会按照下述方式来表示MTF测量的置信度：以σ_Δ作为MTF测量精度的1σ精度，则MTF测量精度在区间的置信度为66.7％；以2σ_Δ作为MTF测量精度的2σ精度，则表示测量精度在区间的置信度为95.45％。

另外，由于通常有指标要求MTF测量的精度要满足优于10％，针对这种情况，可计算相对误差的绝对值|Δ_n|，统计|Δ_n|小于10％的个数N₁，计算则可认为MTF测量精度优于10％的置信度为

以上方法除了可以获取光学遥感载荷在轨MTF测量精度，还可以通过在步骤E调整MTF测量过程中的关键参数，如靶标对比度A、B、靶标尺寸、靶标倾角等，以及在步骤H和步骤J调整噪声标准差和随机相位，然后计算系统误差和整体精度，进行关键误差因素的敏感性分析，以实现算法的优化和靶标设计与布局的优化。

至此，已经结合附图对本实施例进行了详细描述。依据以上描述，本领域技术人员应当对本发明获取刃边法测量光学遥感载荷在轨MTF精度的方法有了清楚的认识。

需要说明的是，在附图或说明书正文中，未绘示或描述的实现方式，均为所属技术领域中普通技术人员所知的形式，并未进行详细说明。此外，上述对各元件和方法的定义并不仅限于实施例中提到的各种方式，本领域普通技术人员可对其进行简单地更改或替换，例如：

(1)带噪声刃边参考图像I_k(x，y)的数目可以根据需要进行调整，可以取50、100或者200，本领域技术人员可以根据需要进行选择；

(2)本文可提供包含特定值的参数的示范，但这些参数无需确切等于相应的值，而是可在可接受的误差容限或设计约束内近似于相应值；

(3)实施例中提到的方向用语，例如“上”、“下”、“前”、“后”、“左”、“右”等，仅是参考附图的方向，并非用来限制本发明的保护范围。

综上所述，本发明提供一种获取刃边法测量光学遥感载荷在轨MTF精度的方法。该方法针对现有刃边法评估光学MTF方法的技术特点和实际观测参数，提出了系统性获取其MTF测量精度(或置信度)的方法，该方法还可以用于关键技术参数的敏感性分析，从而可为刃边的靶标设计与布设以及实施刃边法的MTF在轨评估提供理论指导和依据。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (14)

1.一种获取基于刃边法测量光学遥感载荷在轨MTF在第一方向的测量精度的方法，其特征在于，包括：

步骤A：在光学遥感载荷获取的实测刃边靶标图像中截取尺寸为m行n列的多个刃边子图像，由该多个刃边子图像获取光学遥感载荷在第一方向的在轨MTF参考值和刃边倾角参考值；

步骤B：由光学遥感载荷在第一方向的在轨MTF参考值，获取光学遥感载荷在该第一方向的扩散尺度，构建亚像元尺度下的参考图像退化模板psf(x，y)；

步骤C：由光学遥感载荷获取的实测刃边靶标图像提取被刃边分割的两个均匀区域的灰度均值和随机噪声标准差；

步骤D：根据两个均匀区域的随机噪声标准差，获得用于获取后续步骤中亚像元水平下刃边噪声参考图像的噪声参数：刃边噪声参考图像中两个区域的噪声标准差；

步骤E：利用刃边靶标参数，构建亚像元水平的第一方向的理想刃边参考图像I(x，y)，其中，所述刃边靶标参数包括：刃边子图像尺寸、靶标图像被刃边分割的两个均匀区域的灰度均值、光学遥感载荷在第一方向的刃边倾角参考值；

步骤F：将所述理想刃边参考图像I(x，y)与参考图像退化模板psf(x，y)进行二维卷积，截取卷积结果矩阵中有效数据，得到亚像元水平的第一方向的无噪声刃边退化图像I_psf(x，y)；

步骤G：对所述刃边退化图像I_psf(x，y)中的数据进行数据聚合，生成整像元水平的第一方向的无噪声刃边退化参考图像g(x，y)；

步骤H：利用所述刃边噪声参考图像中两个区域的噪声标准差，构建K幅亚像元水平的第一方向的带噪声刃边参考图像I_k(x，y)，k＝1，2，…K；

步骤I：将带噪声刃边参考图像I_k(x，y)k＝1，2，…K，分别与参考图像退化模板psf(x，y)进行二维卷积，截取卷积结果矩阵中的有效数据得到K幅亚像元水平的第一方向的带噪声刃边退化图像I_psf-k(x，y)，k＝1，2，…K；

步骤J：对所述带噪声刃边退化图像I_psf-k(x，y)，k＝1，2，…K，进行数据聚合处理，生成K幅整像元水平的阵列方向带噪声退化参考图像g_k(x，y)，k＝1，2，…K；以及

步骤K：利用刃边法MTF测试算法，对所述整像元水平阵列无噪声刃边退化图像g(x，y)，以及所述K幅整像元水平的阵列方向带噪声刃边退化参考图像g_k(x，y)，k＝1，2，…K，分别进行第一方向的MTF计算，将Nyquist频率处的归一化MTF计算值与第一方向的在轨MTF参考值进行对比，得到光学遥感载荷在轨MTF精度；

其中，所述第一方向为载荷阵列方向和载荷运动方向其中之一。

2.根据权利要求1所述的方法，其特征在于，所述步骤A包括：

子步骤A1：在光学遥感载荷获取的实测刃边靶标图像中截取m行×n列的多个刃边子图像D，其中，在选取刃边子图像D时刃边位于刃边子图像D中间位置；

子步骤A2：对于多个刃边子图像中的每一个，采用刃边法计算其在第一方向的归一化Nyquist频率下的MTF和刃边倾角；

子步骤A3：对多个刃边子图像的在第一方向的Nyquist频率下的MTF和刃边倾角分别求平均，获取光学遥感载荷在第一方向的在轨MTF参考值和刃边倾角参考值。

3.根据权利要求1所述的方法，其特征在于，所述步骤B中，所述参考图像归一化退化模板为：

$psf(x,y)=e^{-\frac{x^2}{2{\mathrm{\σ}}_x^2}}{e}^{-\frac{y^2}{2{\mathrm{\σ}}_y^2}}$

其中，x，y分别表示光学遥感载荷在阵列方向和运动方向的坐标，σ_x、σ_y分别表示载荷系统函数PSF在阵列方向和运动方向的扩散尺度，其满足：

${\mathrm{MTF}}_{x-ref}=e^{-2{\mathrm{\π}}^2{{\mathrm{\σ}}_x}^2/4}\×\frac{2}{\mathrm{\π}}$

${\mathrm{MTF}}_{y-ref}=e^{-2{\mathrm{\π}}^2{{\mathrm{\σ}}_y}^2/4}\×\frac{2}{\mathrm{\π}}$

其中，MTF_x-ref为光学遥感载荷在载荷阵列方向的在轨MTF参考值，MTF_y-ref为光学遥感载荷在载荷运动方向的在轨MTF参考值；

其中，所述MTF_x-ref和MTF_y-ref的其中之一通过步骤A得到，其中另一通过步骤A得到或令其与所述其中之一相等。

4.根据权利要求1所述的方法，其特征在于，所述步骤C中：

当第一方向为载荷阵列方向时，实测刃边靶标图像被刃边分割成左右两个均匀区域；

当第一方向为载荷运动方向时，实测刃边靶标图像被刃边分割成上下两个均匀区域。

5.根据权利要求4所述的方法，其特征在于，所述第一方向为载荷阵列方向，所述步骤D包括：

子步骤D1：构建一个均值为0、标准差σ_INZL的初始值为σ_noise-L、尺寸为[(m_ZL+10)×10行]×[(n_ZL+10)×10列]的亚像元水平的噪声矩阵IN_ZL(x，y)，该噪声矩阵IN_ZL(x，y)对应的数据点间距为0.1像元；

其中，σ_noise-L为刃边靶标图像中刃边左侧均匀区域Z_L的随机噪声标准差；m_ZL和n_ZL分别为刃边左侧均匀区域Z_L的行数和列数；

子步骤D2：令噪声矩阵IN_ZL(x，y)与参考图像退化模板psf(x，y)进行卷积，截取卷积结果矩阵中的有效数据((m_ZL×10+1)行×(n_ZL×10+1)列)的数据IPN_ZL；

子步骤D3：对矩阵IPN_ZL中的数据每10×10个数据为一组进行聚合，得到一个m_ZL行×n_ZL列的矩阵PN_ZL，PN_ZL数据点间距为1像元；

其中，m_ZL和n_ZL分别为刃边左侧均匀区域Z_L的行数和列数；

子步骤D4：计算PN_ZL的噪声标准差σ_PNZL；

子步骤D5：对比σ_PNZL和σ_noise-L，如果σ_PNZL不等于σ_noise-L，则调整子步骤D1中噪声矩阵IN_ZL(x，y)的标准差σ_INZL，然后重复子步骤D2～D5，直到σ_PNZL＝σ_noise-L，则将此时噪声矩阵IN_ZL(x，y)的标准差σ_INZL作为后续构建亚像元水平下刃边参考图像刃边左区域噪声标准差σ_nL；

其中，按照子步骤D1～D5得到噪声矩阵IN_ZR(x，y)的标准差σ_INZR将作为后续构建亚像元水平下刃边参考图像刃边右区域噪声标准差σ_nR。

6.根据权利要求4所述的方法，其特征在于，所述第一方向为载荷阵列方向，所述步骤E包括：

子步骤E1：基于刃边子图像的尺寸，构建一(m+1+10)×10行，(n+1+10)×10列的亚像元水平的参考图像P，参考图像P的数据点间距为0.1像元；

其中，刃边子图像的数据点间隔为1像元；

子步骤E2，在参考图像P中，令与图像P的列方向成θ_x-ref角度的一条直线l，从参考图像P中第像素的左上角点O点斜穿而过，则该直线l穿越P中某一像素时，分割该像素为左右两个部分，设参考图像P中各个像素的面积为1，根据P中各个像素被直线l分割的左右两部分的面积为权重，对参考图像P进行灰度赋值，灰度赋值后的参考图像P即为亚像元水平的载荷阵列方向的理想刃边参考图像I(x，y)；

其中，θ_x-ref为光学遥感载荷在载荷阵列方向的刃边倾角参考值。

7.根据权利要求6所述的方法，其特征在于，所述子步骤E中对参考图像P进行灰度赋值包括：

(1)对于参考图像P中没有被直线l穿过的像素p(x，y)，如果p(x，y)位于直线l左侧，则p的灰度I(x，y)赋值为：

I(x，y)＝μ_L

(2)对于参考图像P中没有被直线l穿过的像素p(x，y)，如果p(x，y)位于直线l右侧，则该p(x，y)的灰度I(x，y)赋值为：

I(x，y)＝μ_R

(3)对于被直线l穿过的像素p(x，y)，则计算直线l穿过p时将p分割成为左右两部分，所占面积分别为S₁、S₂，则该像素p(x，y)的灰度I(x，y)赋值为：

I(x，y)＝S₁·μ_L+S₂·μ_R

其中，μ_L、μ_R分别为靶标图像中被刃边分割的左右两侧均匀区域的灰度均值。

8.根据权利要求1所述的方法，其特征在于，所述步骤F中无噪声刃边退化图像I_psf(x，y)的尺寸为：(m+1)×10+1行、(n+1)×10+1列。

9.根据权利要求1所述的方法，其特征在于，所述步骤G中，对刃边退化图像I_psf(x，y)，以I_psf(x，y)中位于行、列的像素为起始位置，在图像I_psf(x，y)中上下左右每10×10个像素为1组进行数据聚合，生成整像元水平的第一方向的无噪声刃边退化参考图像g(x，y)。

10.根据权利要求1所述的方法，其特征在于，所述第一方向为载荷阵列方向，所述步骤H包括：

子步骤H1：构建K幅(m+1+10)×10行、(n+1+10)×10列的亚像元水平的参考图像P_k，k＝1，2，…K，参考图像P_k的数据点间距为0.1像元；

子步骤H2：分别构建K个均值为0、标准差σ_nL、尺寸为(m+1+10)×10行、(n+1+10)×10列的噪声矩阵NPL_k，k＝1，2，…K，和K个均值为0、标准差σ_nR、尺寸为(m+1+10)×10行、(n+1+10)×10列的噪声矩阵NPR_k，k＝1，2，…K，噪声矩阵NPL_k和NPR_k，k＝1，2，…K的数据点间距为0.1像元；

其中，σ_nL和σ_nR分别为亚像元水平下刃边参考图像刃边左区域和右区域的噪声标准差；

子步骤H3，在参考图像P_k，k＝1，2，…K中，令与参考图像Pk的列方向成θ_x-ref角度的一条直线l从参考图像P_k中第像素的左上角点O点斜穿而过，则根据该直线l穿越P_k中某一像素时，分割该像素为左右两个部分，设参考图像P_k中各个像素的面积为1，根据P_k中各个像素被直线l分割的左右两部分的面积为权重，对参考图像P_k进行灰度赋值，生成K幅亚像元水平的阵列方向带随机噪声的刃边参考图像I_k(x，y)，k＝1，2，…K；

其中，所述θ_x-ref为为光学遥感载荷在载荷阵列方向的刃边倾角参考值。

11.根据权利要求10所述的方法，其特征在于，所述子步骤H2中，对参考图像P_k，k＝1，2，…K，进行灰度赋值包括：

(1)对于P_k中没有被直线l穿过的像素p_k(x，y)，如果p_k(x，y)位于直线l左侧，则p_k(x，y)的灰度I_k(x，y)赋值为：

I_k(x，y)＝μ_L+NPL_k(x，y)；

(2)如果p_k(x，y)位于直线l右侧，则该p_k(x，y)的灰度I_k(x，y)赋值为：

I_k(x，y)＝μ_R+NPR_k(x，y)；

(3)对于被直线l穿过的像素p_k(x，y)，则计算直线l穿过p_k时将p_k分割成为左右两部分，所占面积分别为S_1k、S_2k，则该像素p_k(x，y)的灰度I_k(x，y)赋值为：

I_k(x，y)＝S_1k·[μ_L+NPL_k(x，y)]+S_2k·[μ_R+NPR_k(x，y)]。

12.根据权利要求1所述的方法，其特征在于，所述步骤J包括：

子步骤J1：在(-5、-4、-3、-2、-1、0、1、2、3、4、5)这11个数中每次随机选出2个数，共随机选出K对数值(S_xk，S_yk)k＝1，2，…K，作为随机相位的偏移值；

子步骤J2：对第k幅亚像元水平的阵列方向带噪声刃边退化图像I_psf-k(x，y)，k＝1，2，…K，以I_psf-k(x，y)中第行、列的像素为聚合起始位置，在图像I_psf-k(x，y)中上下左右每10×10个像素为1组进行聚合，得到K幅带有随机噪声和随机相位的、间距为1像元的整像元水平的阵列方向带噪声退化刃边参考图像g_k(x，y)，k＝1，2，…K。

13.根据权利要求1至12中任一项所述的方法，其特征在于，在考虑随机因素时，以MTF的均方根误差作为光学遥感载荷在轨MTF的1σ精度，以MTF测量精度在区间作为光学遥感载荷在轨MTF的1σ置信度，其中：

系统的平均相对误差：

$\begin{array}{c}\stackrel{\‾}{\mathrm{\Δ}}=avg(\frac{{\mathrm{MTF}}_{k-cal}-{\mathrm{MTF}}_{ref}}{{\mathrm{MTF}}_{ref}}\×100\%)\\ =\frac{1}{N}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\left(\frac{{\mathrm{MTF}}_{k-cal}-{\mathrm{MTF}}_{ref}}{{\mathrm{MTF}}_{ref}}\right)\×100\%\end{array}$

均方根误差：

$\begin{array}{c}{\mathrm{\σ}}_{\mathrm{\Δ}}=std(\frac{{\mathrm{MTF}}_{k-cal}-{\mathrm{MTF}}_{ref}}{{\mathrm{MTF}}_{ref}}\×100\%)\\ =\sqrt{\frac{\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{({\mathrm{MTF}}_{k-cal}-\stackrel{\‾}{\mathrm{\Δ}})}^2}{K-1}}\×100\%\end{array}$

其中，MTF_k-cal为刃边随机参考图像g_k(x，y)利用刃边法计算得到的MTF值，并提取的Nyquist频率下的MTF值，MTF_ref为光学遥感载荷在第一方向的在轨MTF参考值。

14.根据权利要求1至12中任一项所述的方法，所述步骤K中，在不考虑随机因素时，系统误差为：

$\mathrm{\Δ}=\frac{{\mathrm{MTF}}_{cal}-{\mathrm{MTF}}_{ref}}{{\mathrm{MTF}}_{ref}}\×100\%,(k=1,2,...,K)$

其中，MTF_cal为基于无噪声退化参考图像g(x，y)利用刃边法计算得到的MTF值，MTF_ref为光学遥感载荷在第一方向的在轨MTF参考值。