抗差估计在倾斜刃边法MTF检测中的应用分析

第37卷第1期2017年1月海洋测绘HYDROGRAPHIC SURVEYINGAND CHARTINGVol. 37,No. 1Jan. ,2017D 01:10.3969/j .issn .1671-3044.2017.01.018抗差估计在倾斜刃边法M TF 检测中的应用分析曲梦雅，张永生(郑州测绘学院，河南郑州450052)摘要：传统倾斜刃边法没有充分考虑噪声的影响，在进行刃边直线拟合和ESF (edge spread function )拟合时均 米用最小二乘拟合估计，当影像噪声较大时会使得调制传递函数MTF (modulation transfer function )重建结果失真。

采用抗差估计代替最小二乘估计，优化了传统倾斜刃边法中刃边直线拟合和ESF 拟合过程，并采用模拟仿真的方 法对倾斜刃边法的边缘直线拟合精度及MTF 重建精度进行了评估与分析。

实验结果表明，基于抗差估计的倾斜 刃边法可有效抑制噪声影响，边缘直线拟合结果及MTF 重建结果均与真实值吻合较好，且稳定性强。

关键词：遥感影像;调制传递函数;倾斜刃边法;抗差估计;ESF 拟合;图像噪声中图分类号：P 237文献标志码：B文章编号：1671-3044(2017)01-0071-041引言调制传递函数（modulation transfer function, M TF)是评价光学系统成像质量的一个重要指标[1]，它反映了成像系统在不同空间频率下的响应 度。

利用卫星的遥感图像数据信息来检测成像系统 的MTF ，是对卫星成像系统进行在轨检测的重要手 段，对于卫星遥感图像的质量评估与图像复原具有重要的意义。

李小英等[2]基于MTF 对CBERS-02B卫星W FI 图像进行补偿，部分恢复了 W F I 图像的 细部信息。

张瑛等[3]通过MTF 滤波对北京一号小 卫星遥感影像进行融合，效果良好。

Kohm [4]详细描 述了 MTF 估计方法，并对0RBVIEW -3卫星高分辨 率遥感影像质量进行了评估。

陈博洋[5]基于MTF 进行图像复原，提高了风云二号（F Y -2)红外波段遥感图像的空间分辨率。

MTF 检测的方法主要有：点源法、刃边法和脉冲法等[6]，其中倾斜刃边法应用最为广泛。

在倾斜刃边法中，刃边边缘直线的确定和ESF (edge spreadfunction)的拟合重采样是其关键步骤。

现有算法在确定边缘点亚像素位置后，常采用最小二乘直线拟 合来确定刃边边缘直线。

ESF 拟合重采样的常用方 法有:三次样条插值法、高次多项式拟合法及费米 (Ferm i)函数拟合法[7]。

其中三次样条插值法拟合 精度较差;高次多项式拟合法由于没有固定形状，在 局部范围内拟合精度较高，但容易受到噪声等粗差的影响，拟合精度也不高;Ferm i 函数拟合法由于其 形状与E S F 相似，可在一定程度上抑制噪声的影响,拟合精度较高[7]。

实际中常用Ferm i 函数拟合法,例如国际上的IS012233倾斜刃边法[8]就采用Ferm i 函数进行E SF 拟合重采样。

采用Ferm i 函数进行ESF 拟合解算时，常用算法有模拟退火算法[9]和迭代最小二乘算法[10]，其中模拟退火算法由于收 敛速度较慢,算法性能与参数初始值相关，拟合精度 和计算效率均不高,而迭代最小二乘算法较优。

传统的倾斜刃边法没有考虑噪声的影响，在进 行刃边直线拟合和E SF 拟合时均采用最小二乘拟 合估计。

但最小二乘估计抗粗差能力较差，崩溃污 染率较低[11]，当成像噪声大时，会严重影响刃边直 线拟合精度和ESF 拟合精度，导致最终MTF 检测精 度较差。

本文采用抗差估计[12-14]代替最小二乘估 计，进而改进算法中刃边直线拟合和E SF 拟合过 程，提升倾斜刃边法的抗粗差性能。

最后采用模拟 仿真的方法评估了本文算法在不同噪声条件下的刃 边直线拟合精度和MTF 重建精度，并与现有方法进 行了比较。

2基于抗差估计的倾斜刃边法2.1抗差最小二乘原理可得刃边直线拟合和E SF 拟合的拟合模型可 统一表示为[11]:V = A X - L(1)收稿日期：2016-10-09;修回日期：2016-12-14 基金项目：国家自然科学基金（41001262)。

作者简介：曲梦雅（1991-)，女，山东烟台人，硕士研究生，主要从事遥感影像MTF 与图像复原的研究。

72海洋测绘第37卷式中为模型待估参数，对应文中的刃边直线拟合 参数和E S F 拟合函数待估参数;4为系数矩阵;L 为 自由项;F 为拟合残差。

在观测值相互独立的条件 下设其权矩阵最小二乘算法的目标函数和相应解算结果分 别为：nX 户夕2 = min(2)i = 1为粗差,权值赋为零。

2.2基于抗差估计的倾斜刃边法步骤基于抗差最小二乘估计的倾斜刃边法具体步骤可总结为：(1)刃边区域的选择。

选择满足要求即具有刃边边缘的图像数据，尽量保证选择的刃边图像两侧 区域的灰度值有一定差异并且每一个区域内的灰度比较均勻。

X = (a t P A ) ~1 A t P L(3)式中,n 为拟合样本的个数。

公式(2)表明最小二乘法目标函数是残差平方 和的函数，若样本中存在粗差,会对残差平方和造成 很大影响，故最小二乘估计不具有抗差性，为使估计 稳健化,要控制样本粗差对平差结果的影响,需寻找 一个增长慢的函数p (〜）来代替^，抗差最小二乘估 计的目标函数为[11]:X P i P ( v i ) = min⑷i =1dp ( Vi )-取加t =-------，等价权 P = diagj 户，t } = diagV i dV i 丨-i 1,则抗差估计解的结果可表示为：X = (A t P A )_1At P L(5)式(5)表明抗差最小二乘保留了最小二乘算法 的优点，并通过选择适当的等价权函数来增强算法 的抗差性,其等价权函数是残差的函数，在实际解算 时需迭代求解，本文以最小二乘结果作为迭代的 初值。

不同的等价权函数具有不同的抗差性,常用的 等价权函数有H u b er 权函数、IGG 1权函数以及 IGG 3权函数等[12]。

本文采用IGG 3等价权函数作 为抗差估计的等价权函数，该等价权函数将观测样 本分为三类:有用信息、可用信息和无用信息。

通过 比较标准化残差和等价权函数中临界值来赋予三类 信息不同的权重,从而达到较好的抗粗差效果。

其 数学模型为[12]:1, v ^ *〇IGG 3( v )=h - A 0k(2)，k 〇 < v 矣 k,(6)〇，v > k1于1〇〇3权函数一般取&0 = 1.5和&1=2.5。

式(7)表明，IGG3权函数分为三段:①当v 矣k 〇时, 样本为有用信息，采用最小二乘法，样本权值不变;②当k 〇< v ^k ,时，样本为可用信息，对其进行降 权处理;③当v >k ,时，样本为无用信息，即将其视(2)刃边边缘直线抗差最小二乘拟合。

对于选定的刃边区域，采用三次多项式二阶导为零的方 法[12]确定刃边边缘点的亚像素位置,并进行抗差最 小二乘直线拟合得到刃边边缘直线。

(3)依据亚像素边缘位置进行E S F 采样。

在获取刃边边缘拟合直线后，将各采样点沿刃边边缘直 线方向投影到刃边边缘梯度方向上获得初始ESF 样本。

(4)基于抗差估计的F erm i 函数E S F 拟合重采样。

E S F 初始样本受刃边边缘定位误差等影响，往 往存在粗差，并且在局部区域内是不规则的,需对其 进行拟合重采样。

实验表明，采用3个F erm i 函数 的线性组合即可获得足够的拟合精度[10]:E SF (X )= X—— + d (7)〜〔〒卜式中,-i A和^为待估拟合参数，采用抗差估计解算得到Ferm i 函数参数，从而得到规则的E S F 曲线。

(5)L S F 米样。

对于离散样本,L SF(line spreadfunction )可通过E S F 样本直接差分得到：E SF t - E S F t 1LSFi =—i i 1 (8)在获得L S F 样本后,首先对L S F 进行尾部截取 处理以减弱L S F 尾部抖动影响。

L S F 的形状近似为 高斯型[7]，但并不完全符合高斯表达式，为了减少 P S F 估计的误差，对得到的L S F 曲线进行一维高斯 函数拟合。

(6)对归一化的L S F 进行离散傅里叶变换，取各分量的模,以零频率处的M T F 值作为基准进行归 一化重建得到M T F 曲线。

3实验分析3.1刃边直线拟合精度分析为评估基于抗差估计的刃边直线拟合精度，采 用计算机仿真不同噪声条件下的退化边缘图像，并 采用抗差估计对其进行刃边直线拟合，将所得结果 与模拟真值比较。

退化图像的生成过程为:①采用第1期曲梦雅，等：抗差估计在倾斜刃边法MTF检测中的应用分析73计算机模拟生成无退化的刃边图像,并对其添加高斯白噪声;②采用高斯函数构造PSF;③利用PSF对含噪声的刃边图像进行卷积处理得到模拟的退化刃边图像。

模拟退化边缘图像参数设置见表1,其中0为刃边倾角'为高斯函数参数。

表1退化刃边图像参数设置高DN值低DN值0〇■刃边图像大小2404015° 1.0 100像素x100像素对图像添加不同的高斯白噪声，噪声参数分别设置为0. 1、〇. 2、0. 3、0. 4、0. 5、0. 6。

为评价算法的稳定性，每个噪声条件下的退化图像分别随机生成100幅，采用抗差估计对其进行刃边直线拟合,将拟合直线倾角与真值作差，得到拟合误差序列，并计算误差序列RMS值,将所得结果与常规最小二乘算法比较。

不同噪声条件下的刃边倾角拟合误差RMS值见图1,其刃边倾角拟合误差序列见图2,其中LS表示最小二乘估计,RLS表示抗差估计。

图1刃边倾角拟合偏差RMS值示意图0.3 ^0-2 〇0.1 ^ 〇.3-0.1 ^-0.2 ■0.3sigma=0.1—LS■命RLS60 80 100number (n)60 80 100number (n)'0.53-0.5 -1.01si g m a=6.5+L S■-R L S60 80 100number (n)2' sigma=6.21-1-2夺LS♦R L SJ20 40 60 80 100number (n)u 'sigma: d).4-*-LS I♦R L Sw20 40 60 80 100number (n)u ' sigma=6.65"♦L S-R L Snumber (n)图2刃边倾角拟合误差序列示意图(2)当噪声增加时，传统拟合算法受噪声影响严重，拟合误差RMS值明显增大，刃边直线拟合精度严重下降，且算法稳定性较差，拟合误差序列波动较大,特别当噪声参数较大时尤为明显。