CN102339460B - 卫星图像自适应复原方法 - Google Patents

Abstract

本发明为一种卫星图像自适应复原方法，包括步骤：S1，利用周期加平滑图像分解技术，将模糊卫星图像Y分解成为周期信号P和平滑信号S两部分；S2，对卫星图像退化参数进行估计；S3，利用所述参数估计结果对所述周期信号P进行基于复数小波包域的自适应复原，得到周期信号P的复原结果X_P；以及S4，将所述周期信号P的复原结果X_P加上所述平滑信号S，作为卫星图像最终的复原结果X。本发明复原方法可以避免边界振铃现象，实现在去模糊的同时抑制噪声(特别是彩色噪声)的增益，以及保持强方向性纹理区域，实现对高频信息丰富的卫星图像的有效复原。

Description

卫星图像自适应复原方法

技术领域

本发明涉及遥感数字图像处理领域，特别涉及一种基于周期加平滑图像分解和复数小波包域的卫星图像自适应复原方法。

背景技术

在过去的几十年中，遥感图像的复原一直是国内外科研机构的热点问题，大致可分为单尺度和多尺度两种方法。比较经典的单尺度方法包括Wiener反卷积方法，Andre Jalobeanu的RHEA方法，Richrdson-Lucy方法，Tikhonov正则化方法等等。所有这一类单尺度去卷积方法的问题在于，对于强方向性或者细节性的局部纹理，不能起到很好的增强作用，反而会出现一些虚假瑕疵；另外由于优化过程需要迭代，所以单尺度方法的速度往往较慢。以小波变换为基础的多尺度方法是目前研究的重点，常见的有实数小波、小波包等等。多尺度方法的优势在于：对图像的局部纹理可以有非常好的增强效果，而同时抑制平坦区域的噪声；如果不涉及迭代，通常速度较快。但是多尺度方法通常需要在分解域中修改系数并反变换回空间域，因此很容易产生振铃效应。

由于卫星图像的特殊性以及图像复原本身固有的病态性、复杂性，卫星图像复原还存在着一些难以解决的问题：(1)边界振铃。卫星图像尺寸往往非常巨大，所以在实际的工程应用中一般会先将图像进行分割，并行处理后再进行图像拼接。但是图像的分割会造成边界截断，现有方法在进行复原处理时都会产生或多或少的边界振铃，当复原后的子图像拼接为最终结果时，在接边处这种人工瑕疵视觉上非常明显，并且这种造假一旦引入就很难消除。(2)噪声增益。通常情况下反卷积的卷积核都为低通滤波器，在高频区域大部分值都接近于零，所以直接求解反卷积就会过于放大高频域的噪声。如果卷积核存在估计误差，就会加剧这种病态性。(3)方向性纹理丰富。高分辨率的卫星图像，包含大量的方向性和细节性的局部纹理，单尺度方法对这类纹理不能起很好的增强作用，而传统小波工具对方向的选择特性较弱，对这类强方向性纹理的保持也具有一定局限性。

综上，现有遥感图像复原技术还没有能够同时解决边界振铃、噪声增益和保持强方向性纹理的有效方法。

发明内容

(一)要解决的技术问题

本发明要解决的技术问题包括：如何解决现有复原方法难以避免的边界振铃现象，如何实现在去模糊的同时抑制噪声(特别是彩色噪声)的增益，以及如何保持强方向性纹理区域，实现对高频信息丰富的卫星图像的有效复原。

(二)技术方案

为了解决上述技术问题，本发明提供一种卫星图像自适应复原方法，包括步骤：S1，利用周期加平滑图像分解技术，将模糊卫星图像Y分解成为周期信号P和平滑信号S两部分；S2，对卫星图像退化参数进行估计；S3，利用所述参数估计结果对所述周期信号P进行基于复数小波包域的自适应复原，得到周期信号P的复原结果X_P；以及S4，将所述周期信号P的复原结果X_P加上所述平滑信号S，作为卫星图像最终的复原结果X。

优选地，所述步骤S2中图像退化参数估计包括MTF参数估计和噪声参数估计。

优选地，所述步骤S3中基于复数小波包域的自适应复原包括：利用复数小波包工具，对周期信号P先去卷积再多频段建模分析降噪，进行去模糊和降噪的复原处理。

(三)有益效果

本发明的技术方案具有如下有益效果：不仅较好地解决了现有复原方法难以避免的边界振铃现象，还通过合理加强先验模型的假设和限制，提高了退化参数自适应估计的准确度，并且充分发挥了复数小波包工具的近似平移不变性和更丰富的方向选择性(6个方向响应)的优点，在去模糊的同时实现对噪声(特别是彩色噪声)的抑制，出色地保持了强方向性纹理区域，实现对高频信息丰富的卫星图像的有效复原。实验证明，除了视觉效果的明显优势以外，该方法在信噪比(SNR)以及结构相似度(SSIM)指标上都优于其他方法，可以支持大规模图像数据的并行计算和业务化运行，具有实际的工程应用价值。

附图说明

图1为本发明卫星图像自适应复原方法的流程图；

图2A-2F为将一幅卫星图像分别利用Wiener方法、Richardson-Lucy方法、[S.Cho et a1.2009]方法、实数小波方法、复数小波包方法以及本发明复原方法进行复原后的效果示例图；

图3A-3F分别为另一幅卫星图像的原始清晰图像、模糊加噪声的图像、利用Wiener方法、[S.Cho et al.2009]方法、实数小波方法以及本发明复原方法进行复原后的效果示例图。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不是限制本发明的范围。

如图1所示，本发明所述的卫星图像自适应复原方法包括步骤：S1，利用周期加平滑图像分解技术，将模糊卫星图像Y分解成为周期信号P和平滑信号S两部分；S2，对卫星图像退化参数进行估计；S3，利用所述参数估计结果对所述周期信号P进行基于复数小波包域的自适应复原，得到周期信号P的复原结果X_P；以及S4，将所述周期信号P的复原结果X_P加上所述平滑信号S，作为卫星图像最终的复原结果X。

所述步骤S1，首先利用周期加平滑图像分解技术，将模糊图像分解成为周期信号P和平滑信号S两部分，只对其中的周期信号P进行反卷积处理，通过这种方法来消除频域处理造成的边界振铃。

对于一个M×N的离散图像U∈R^Ω，Ω＝{0，..M-1}×{0，..N-1}，其周期元素P的计算如下公式所示：

$\∀(q,r)\∈\mathrm{\Ω}\left(\mathrm{0,0}\right),$

$\widehat{\mathrm{per}\left(u\right)}(q,r)=\hat{u}(q,r)-\hat{v}(q,r)/(4-2\cos \left(\frac{2\mathrm{\πq}}{N}\right)-2\cos \left(\frac{2\mathrm{\πr}}{M}\right)),$

并且 $\widehat{\mathrm{per}\left(u\right)}\left(\mathrm{0,0}\right)=\hat{u}\left(\mathrm{0,0}\right),$

其中v＝v₁+v₂

$\∀(x,y)\∈\mathrm{\Ω},v_1(x,y)=\left\{\begin{array}{cc}u(x,y)-u(M-1-x,y)& \mathrm{if\; x}=0\mathrm{orx}=M-1,\\ 0& \mathrm{else},\end{array}\right.$

$\∀(x,y)\∈\mathrm{\Ω},v_2(x,y)=\left\{\begin{array}{cc}u(x,y)-u(x,N-1-y)& \mathrm{if\; y}=0\mathrm{ory}=N-1,\\ 0& \mathrm{else},\end{array}\right.$

通过上述公式，最终可得到图像U分解后的周期元素为P＝per(u)，光滑元素为S＝U-P。

步骤S2，预估计卫星图像退化参数调制传递函数MTF和噪声方差σ²。

对于卫星图像，导致模糊的因素主要包括大气、光学系统、相对运动、以及传感器这些方面，总的MTF可以写成：

MTF＝MTF_int·F_u，v    (1)

这里MTF_int是已知的MTF初始值，

F_u，v＝exp(-α_uu²-α_vv²)    (2)

因此我们需要估计的MTF参数，相当于是(2)中F_u，v的α_u和α_v的值。

对于清晰的卫星地形图像，图像在频域上的能量谱近似满足概率分布：

$F{\left[X\right]}_{u,v}=G_{u,v}{w}_0{r}_{u,v}^{-q}---\left(3\right)$

其中w₀和q是关于地形图像的未知参数。

根据此模型，我们可以使用最大似然方法估计出(2)中参数α_u和α_v的值。

对于卫星图像退化参数噪声的估计，根据Andre Jalobeanu的方法，我们将信号变换到频域：

F(Y)＝F(H)F(X)+F(N)    (4)

选取高频段[r₀，r_max]，假设在这个频段内F(H)F(X)的能量与F(N)相比可以忽略；而作为高斯白噪声，F(N)的能量平均地分布到频域，其密度为σ²。我们计算F(Y)的平均能量，作为噪声方差σ²的估计值。

步骤S3，基于复数小波包域的自适应复原，利用复数小波包工具，对周期信号P先去卷积再多频段建模分析降噪，进行去模糊和降噪的复原处理。

首先我们按照公式X＝(Y-N)/H来计算Y/N，这里H需要做强制的下界限制，我们选取经验值0.05。由于高频噪声被放大很多，因此Y/N包含的噪声是彩色的，我们使用复数小波包分解工具，分频带对彩色噪声进行建模和处理，假设每个频带K中，Y/N包含噪声方差为

的高斯白噪声，

可以由参数H和σ唯一决定，按如下公式进行计算：

${\mathrm{\σ}}_K^2={\mathrm{\σ}}^2\×{\mathrm{\Σ}}_{i,j}{|F{\left(W_K\right)}_{i,j}/F{\left(H\right)}_{i,j}|}^2---\left(5\right)$

其中W_K为第K个频段的冲击响应。

基于以上各频带的噪声模型，我们由观测数据x₀＝ξ+u可以写出概率分布：

$P\left(x_0\right|\mathrm{\ξ})=\exp (-\frac{{|x_0-\mathrm{\ξ}|}^2}{2{\mathrm{\σ}}_K^2})/2\mathrm{\π}{\mathrm{\σ}}_K^2---\left(6\right)$

这里x₀为各频带的信号，ξ为对原始图像信号的一个估计值，而ξ的先验分布可以通过一个粗糙的去卷积结果给出。假设一个粗糙的去卷积方法所得到的结果是s，ξ的先验分布可以表示成如下形式：

$P\left(\mathrm{\ξ}\right|s)=\exp (-\frac{{\left|\mathrm{\ξ}\right|}^2}{2s^2})/2{\mathrm{\πs}}^2---\left(7\right)$

结合(6)和(7)，基于最大似然法则，我们估计ξ的值使得乘积P(x₀|ξ)×P(ξ|s)最大。ξ的最理想估计值如(8)所示：

$\mathrm{\ξ}=x_0\×s^2/(s^2+{\mathrm{\σ}}_K^2)---\left(8\right)$

所以要对结果ξ进行计算，还需要得到一个对s²的较精确的估计。这可以先通过一个单尺度的快速估计方法(如RHEA，或Richardson-Lucy方法，或TOG09的方法)获得估计值ξ′，再加上软阈值滤波修正获得：

$s^2={\left|{\mathrm{\ξ}}^{\′}\right|}^2/4-{\mathrm{\σ}}_K^{\′2}$ 若 $\left|{\mathrm{\ξ}}^{\′}\right|\≥4{\mathrm{\σ}}_K^{\′2}$

s²＝0  若 $\left|{\mathrm{\ξ}}^{\′}\right|\≥4{\mathrm{\σ}}_K^{\′2}---\left(9\right)$

这里

是关于ξ的估计方差，也可以使用

替代。

结合(8)和(9)我们可以分频带获得所有ξ的值，最后再通过复数小波包逆变换还原到空间域上，即为去卷积的结果X_P。

步骤S4，信号合成。周期信号P的复原结果X_P与光滑元素S相加，得到最终的复原结果X＝X_P+S。

如图2A-2F、3A-3F所示，为现有方法及本发明基于周期加平滑分解和复数小波包域的卫星图像自适应复原方法的效果示例图。图2A-2F中，首先将一幅真实的卫星图像分割成4部分，然后分别采用不同的方法进行复原处理，最后再合成一幅图像。从结果中可以清楚地看出：Wiener方法和[S.Cho et al.2009]方法在复原的过程中引入了大量的噪声和振铃，Richardson-Lucy方法、实数小波方法、复数小波包的方法在图像边缘处有明显的振铃，这在拼接合成的时候尤为明显，而本发明方法不仅很好地恢复了图像的高频细节部分，同时抑制了噪声的放大，而且确实有效解决了边界振铃这一难题。图3A-3F中，从复原结果可见，Wiener方法噪声增益很大，[S.Cho et al.2009]的方法振铃现象明显，与实数小波方法相比，本发明方法对强方向性的局部纹理保持得更好。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和替换，这些改进和替换也应视为本发明的保护范围。

Claims (3)

1.一种卫星图像自适应复原方法，其特征在于，包括步骤：

S1，利用周期加平滑图像分解技术，将模糊卫星图像Y分解成为周期信号P和平滑信号S两部分，对所述周期信号P进行反卷积处理；

S2，对卫星图像退化参数进行估计；

S3，利用所述参数估计结果对所述周期信号P进行基于复数小波包域的自适应复原，得到周期信号P的复原结果X_P；以及

S4，将所述周期信号P的复原结果X_P加上所述平滑信号S，作为卫星图像最终的复原结果X。

2.如权利要求1所述的卫星图像自适应复原方法，其特征在于，所述步骤S2中图像退化参数估计包括MTF参数估计和噪声参数估计。

3.如权利要求1或2所述的卫星图像自适应复原方法，其特征在于，所述步骤S3中基于复数小波包域的自适应复原包括：利用复数小波包工具，对周期信号P先去卷积再多频段建模分析降噪，进行去模糊和降噪的复原处理。

Images