一种基于空域的自适应MTFC遥感图像复原方法

摘要

一种基于空域的自适应MTFC遥感图像复原方法，首先根据实验室实测的MTF值以及用户的需求值确定遥感器系统从零频到奈奎斯特频率范围内不同频率处MTF的补偿值。然后依据有限长滤波器原理，计算获取一维卷积系数，并在此基础上，构造二维有限长滤波器而不是传统的傅里叶变换来计算二维去卷积系数，大大简化了计算过程。同时，本发明方法根据实验室实测的不同成像参数下不同灰度的SNR和不同的景物内容进行自适应抑噪，可以在提升图像MTF的同时保证图像的SNR下降很小。本发明方法能有效提升图像的MTF，提高图像的清晰度，同时不会放大噪声，改善图像质量。

说明书

技术领域

本发明属于图像处理技术领域，具体涉及一种遥感影像的复原方法。

背景技术

遥感成像过程受大气和成像系统中光学、平台、探测器、电子器件等的影 响，会出现成像退化的情况，造成图像模糊，降低图像的质量。如果充分利用 地面MTFC（调制传递函数补偿），则可以一定程度上改善图像质量，提高图像 信息的解译能力。

MTFC是一种由光学遥感器系统MTF引申出来的一种数学概念，最早由 Leachtenauer在其著作中首次提出。其思想被美国政府主导开发的通用图像质 量方程（GIQE）所采纳，并在IKONOS、QuickBird-2、OrbView-3和GeoEye-1 中均进行了成功应用，使得图像质量得到显著提高。但由于技术封锁，很少有 公开发表的文献透露其具体的技术细节。

从1997年起，法国航天局（CNES）组织开展了针对卫星图像不同处理方 法的研究，主要有滤波、线性规整、概率、小波、小波包等。研究结果表明， 利用线性规整的方法会带来振铃现象，产生虚假响应，引起一定的信息损失。 小波方法可以很好地去除白噪声，对彩色噪声无能为力。而小波包的方法可以 避免这些缺点。法国航天局将复数小波包方法成功应用于Pleiades图像的地面 处理中，但是经过前期试验发现复数小波包方法同样也会带来振铃现象。

国内方面，顾行发等进行了CBERS-02卫星CCD相机MTF在轨测量及 图像MTF的补偿研究，针对CBERS-02卫星的CCD相机，采用迭代法、维 纳滤波法及修正的逆滤波器法（MIF）对CCD图像进行MTF补偿，根据实测 45度方向0.5频率处的MTF值插值建立二维MTF矩阵，然后依据在轨测试的 MTF进行补差。但是由于脉冲法和刃边法受噪声影响较大，因此测得的MTF 精度有限，恢复效果并不好。

发明内容

本发明的技术解决问题是：克服现有技术的不足，针对空间光学遥感器光 学系统、探测器与电子元器件等在遥感成像过程中造成的图像质量退化，提出 了一种基于空域反卷积的自适应图像复原方法，可以在不放大噪声的情况下提 高图像的清晰度、改善图像质量，同时处理的速度和效果满足工程实践需求。

本发明的技术解决方案是：一种基于空域的自适应MTFC遥感图像复原方 法，步骤如下：

（1）在地面通过实测得到遥感器系统从零频到奈奎斯特频率范围内的 MTF值；

（2）根据成像需求确定理想的MTF提升曲线，并与步骤（1）的结果进 行比较，得到遥感器系统从零频到奈奎斯特频率范围内不同频率处MTF的补偿 值|H(e^jw)|；

（3）利用步骤（2）的结果，计算得到一维卷积系数h(k)，|H(e^jw)|=A(w,k)*h(k)， 其中h(k)是中心对称的，h(k)的长度为M；w为频率， k取

（4）利用步骤（3）得到的一维卷积系数h(k)构造二维卷积核h(i,j)，h(i,j) 为一个M*M的正方形矩阵且满足中心对称和轴对称，其中i,j的取值范围分别 为h(i,j)中不同的卷积系数的个数为 象限的二维频域响应函数表示为：

$\left|H(f_x,f_y)\right|=\underset{t_x=0}{\overset{\frac{M-1}{2}}{\Sigma }}\underset{t_y=0}{\overset{t_x}{\Sigma }}x(t_x,t_y)\times k(t_x,t_y)[\cos \left(f_x{\times t}_x\right)\times \cos (f_y\times t_y)+\cos (f_y\times t_x)\times \cos (f_x\times t_y)]/2$

式中$\left|H(f_x,f_y)\right|=H\left(\sqrt{{f_x}^2+{f_y}^2}\right)=\underset{k=0}{\overset{\frac{M-1}{2}}{\Sigma }}h\left(k\right)·2\cos (\frac{M-1}{2}-k)\sqrt{{f_x}^2+{f_y}^2},$x(t_x,t_y)为一个长度 的函数，其中t_y∈[0,t_x]，|H(f_x,f_y)|为二维频率处的 MTF补偿值，k(t_x,t_y)为卷积系数的重复次数；

（5）实测遥感器系统在不同成像参数下不同灰度下的SNR，利用拟合的 方式得到遥感器系统在不同成像参数下不同灰度下的SNR曲线；

（6）根据步骤（5）得到的SNR曲线构造抑噪查找表系数M(i)， M(i)=10^((y(i)-y_max)/20)，i在[0,2^N-1]内取整，N为量化位数，i为某一确定成像 参数下的某一DN值，y(i)为i与对应的SNR，y_max为同一个确定的成像参数下 的最大SNR；

（7）利用步骤（4）得到的二维卷积核h(i,j)对待复原的图像进行卷积，然 后进行噪声抑制后输出复原后的图像；所述复原后的图像中的每一个像元的灰 度值i₂满足条件i₂=i₀+k(i₁-i₀)，k=K||i₁-i₀||M([i₀])，其中i₀为像元的原始灰度值， i₁为为像元经过二维卷积后的灰度值，[]表示取整函数，M([i₀])通过查找步骤（6） 得到的抑噪查找表系数获得，K为大于0的整数，是可变的抑噪阈值参数，根 据当前像素补偿前的SNR而定，K的取值应保证补偿后该像素DN值的均匀景 物SNR下降不超过设定值。

本发明与现有技术相比的优点在于：

（1）本发明方法根据实验室测得的遥感器系统静态MTF以及理想MTF 提升曲线确定不同频率点处的MTF补偿值，在获取一维卷积系数的基础上通过 构造二维有限长滤波器计算二维去卷积系数。由于在空域中直接进行图像复原， 避免了在频域中进行MTF提升引起的相位发生变化引起的图像异常现象；

（2）本发明方法根据遥感器系统在不同成像参数下的SNR也不一致，结 合实测的不同参数下的SNR曲线对去卷积后的图像进行不同程度的抑噪，在 提升MTF的同时有效保证图像的信噪比下降很小；

（3）本发明方法在去卷积的过程中，根据景物内容自适应的进行抑噪处理。 对于平坦区域，高频分量中噪声所占比例较大，对于灰度层次丰富的区域，高 频分量中噪声所占比例相对较小。因此本发明方法根据图像灰度层次丰富程度 来确定反卷积复原时的抑噪参数。通过合理的约束，能够很好的避免高频噪声 放大过大。

附图说明

图1为本发明方法的流程框图。

具体实施方式

如图1所示，为本发明方法的流程图，主要步骤如下：

（1）根据实验室测得的遥感器系统不同频率处的MTF值，精确计算一维 卷积系数；

（11）通过实验室实测得到遥感器系统从零频到奈奎斯特频率范围内的 MTF值；

（12）利用步骤（11）的结果，同时结合用户需求值确定复原后理想 MTF曲线，两者进行比较计算得到遥感器系统从零频到奈奎斯特频率范围内不 同频率处MTF的补偿值（补偿的频率点数受处理速度及资源的限制，补偿的频 率点越多需要的处理空间越大，速度也会越慢）；

（13）根据步骤（12）得到的不同频率点的MTF补偿值构造有限长数 字滤波器，计算一维卷积系数。

有限长数字滤波器的构造方法如下：

根据有限长数字滤波器的基本原理，设某频率处的MTF补偿值为|H(e^jw)|， 有限长数字滤波器为h(k)，长度为M，那么根据点扩散函数的中心对称性，它 的频率响应特性如下：

$\left(\begin{array}{c}\left|H\left(e^{\mathrm{jw}}\right)\right|=|\underset{k=0}{\overset{M-1}{\Sigma }}h\left(k\right)·e^{-\mathrm{jwk}}|=|\underset{k=0}{\overset{\frac{M-1}{2}}{\Sigma }}h\left(k\right)·[e^{-\mathrm{jwk}}+e^{-\mathrm{jw}(M-1-k)}|\\ =\underset{k=0}{\overset{\frac{M-1}{2}}{\Sigma }}h\left(k\right)·2\cos (\frac{M-1}{2}-k)w=A(w,k)*h\left(k\right)\end{array}\right)$

其中，w为选取的频率值（从零频到奈奎斯特频率范围 内选取），$k=\mathrm{1,2},...,\frac{M-1}{2}.$

由此可以根据|H(e^jw)|和A(w,k)计算得到序列h(k)。

（2）利用步骤（1）计算得到的一维卷积系数构造二维卷积核；

假设遥感器系统的退化函数在二维频域中是关于中心点(0,0)对称的，那么 二维频率点(f_x,f_y)处的提升倍数等价于一维频率处的提升倍数。因此， 可以由一维的卷积系数构造出一组二维的卷积核。具体如下：

（21）根据步骤（13）中的描述，由于一维卷积系数的长度为M，可 以构造出二维卷积核的大小为M*M的一个正方形矩阵h(i,j)，其中取值范围 由二维卷积核的中心对称性可知，只要求出 一个象限内的卷积系数即可，又由二维卷积核的轴对称性，即h(i,j)=h(j,i)可知 不同卷积系数的个数为

（22）先设定一个长度为的函数x(t_x,t_y)，其中 t_y∈[0,t_x]，即取值范围为象限。设某二维频率处的MTF补偿值为 |H(f_x,f_y)|，由步骤（21）可知二维卷积核的对称性可以得到每个系数重复的次数 为k(t_x,t_y)，则二维频域响应函数可以构造成：

$\left(\begin{array}{c}\left|H(f_x,f_y)\right|=|\underset{n1=0}{\overset{\frac{M-1}{2}}{\Sigma }}\underset{n2=0}{\overset{\frac{M-1}{2}}{\Sigma }}h(n1,n2)\times 4\cos (f_x\times n1)\times \cos (f_y\times n2)|\\ =\underset{t_x=0}{\overset{\frac{M-1}{2}}{\Sigma }}\underset{t_y=0}{\overset{t_x}{\Sigma }}x(t_x,t_y)\times k(t_x,t_y)[\cos (f_x\times t_x)\times \cos (f_y\times t_y)+\cos (f_y\times t_x)\times \cos (f_x\times t_y)]/2\\ =A(f_x,f_y,t_x,t_y)*x(t_x,t_y)\end{array}\right)$

其中，A(f_x,f_y,t)=k(t_x,t_y)[cos(f_x×t_x)×cos(f_y×t_y)+cos(f_y×t_x)×cos(f_x×t_y)]/2。

（23）根据步骤（13）得到的一维卷积系数计算|H(f_x,f_y)|

$\left|H(f_x,f_y)\right|=H\left(\sqrt{{f_x}^2+{f_y}^2}\right)=\underset{k=0}{\overset{\frac{M-1}{2}}{\Sigma }}h\left(k\right)·2\cos (\frac{M-1}{2}-k)\sqrt{{f_x}^2+{f_y}^2}$

注：f_x,f_y在零频到奈奎斯特频率范围内选取。

（24）由步骤（22）和步骤（23）得到象限处的卷积系数x(t_x,t_y)后， 由于二维卷积核的轴对称性，即h(i,j)=h(j,i)，可以构造出一个象限内的卷积系 数，同时根据二维卷积核的中心对称性，构造出整个二维卷积核h(i,j)。

（3）根据遥感器系统成像时所设参数确定的不同灰度下的SNR（信噪比） 信息，计算抑噪的查找表系数；

根据遥感器成像时所设参数（包括积分级数、增益），实验室实测的不同灰 度值下的SNR。一般是低灰度值下的信噪比较低，高灰度值下的信噪比较高。

在逐点进行空域去卷积复原时，根据当前像素的DN值查表得到探测器的 SNR。

（31）根据不同成像参数下多个辐亮度处（对应多档DN值）的实测SNR， 拟合得到整个DN值动态范围内的SNR曲线；

（32）根据SNR曲线构造抑噪查找表系数，构造的前提是保证该亮度 下的SNR信息，构建方法如下：

M(i)=10^((y(i)-y_max)/20)

上式中i在[0,2^N-1]（N为量化位数）内取整，y(i)为某一DN值处的SNR，y_max为该相机参数状态下最大的SNR。

（4）用步骤（2）得到的二维卷积核对待复原的图像进行卷积，卷积过程 中根据抑噪参数以及当前图像内容进行自适应的抑噪。

对于平坦区域，高频部分噪声所占比例较大，对于灰度层次丰富的区域， 高频部分噪声所占比例相对较小。

本发明方法根据图像灰度层次的丰富程度来确定去卷积复原的约束参数。

假设某一像元原始灰度值为i₀，卷积运算后结果为i₁，取i₂=i₀+k(i₁-i₀)为最 终滤波输出值。当k=1时，则滤波后MTF为全额补偿，若k=0.5，则MTF为半 额补偿，若k=0，则为零补偿。

在此，取k=K||i₁-i₀||M([i₀])，其中K为大于0的整数，是可变的抑噪阈值参 数，K的取值要保证补偿后该像素DN值的均匀景物SNR下降不超过一定值， 具体大小根据当前DN值处补偿前的SNR而定，例如要求800DN（10bit量化） 值处的SNR不小于48dB，而补偿前800DN值处的SNR达到了51dB,那么可 以设定补偿后SNR下降2dB，从而可以通过补偿后的SNR反推得到K值。因 此通过调节K对SNR下降情况进行了不同程度的控制。M([i₀])是步骤（3）中 得到的抑噪查找表系数，它是关于i₀的函数，[]在此处表示取整函数。

本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。