CN101699509A - 一种利用气象数据进行大气模糊遥感影像恢复的方法 - Google Patents

Abstract

本发明提供了一种利用气象数据进行大气模糊遥感影像恢复的方法，利用遥感影像成像时刻当地的局部气象数据，先对大气湍流调制传输函数和气溶胶调制传输函数进行估计，然后通过这两种调制传输函数的组合模型得到整个大气效应的调制传输函数，进而利用改进的直接解卷积算法实现遥感影像的去大气模糊操作，完成对受到大气湍流模糊的遥感影像的高清晰恢复与重建。

Description

一种利用气象数据进行大气模糊遥感影像恢复的方法

技术领域

本发明属于遥感影像的大气纠正和数字信号处理领域，涉及对高分辨遥感影像的去大气模糊恢复与重建，特别是涉及对已知成像时刻局部区域气象数据的大气模糊遥感影像的高分辨率恢复技术领域。

背景技术

从上世纪六七十年代至今，我国的遥感信息获取与处理技术已经处于国际先进行列，成功研制并发射了一系列遥感卫星。但是，随着现代遥感卫星地面分辨率的不断提高，大气湍流和气溶胶成分对卫星遥感影像的影响也越来越严重。研究表明，当大气湍流比较弱时，它对卫星光学遥感系统地面分辨力的影响只有1cm左右。然而，湍流较强时，这一影响可以增大到10cm的水平，对于现代高分辨率的光学遥感系统，这种影响不容忽视，所以，需要对卫星遥感数据进行后处理，消除大气湍流对影像所造成的影响。大气对遥感影像的退化作用主要是由湍流、气溶胶和大气微粒的吸收和散射造成。大气湍流造成长曝光成像波前相位的倾斜，进而造成像平面影像的扭曲；气溶胶影响变化比较缓慢，主要造成光线的扩散和散射；大角度的光线散射造成光强度的衰减，而小角度的光线散射则造成影像的模糊。

对模糊遥感影像进行恢复，通常可以把观测影像模型化为原始影像和大气调制传输函数(MTF：Modulation Transform Function)的卷积过程，如果能够准确估计大气MTF值和确定各种噪声的统计特性，就可以利用解卷积技术来消除大气效应对遥感影像的影响。

对受到大气模糊的遥感影像进行恢复与重建是遥感影像处理中的一个难点，人们常用的算法主要有带约束的复原算法、正则化算法、蒙特-卡罗算法、最大熵算法等，这些算法需要传感器平台信息，光学成像系统以及成像时刻大气状况的参数信息，然而，这些参数通常是很难得到精确满足的。遥感影像的获取与传输过程要经过大气、光学系统、CCD等一系列环节，各个环节均可能对影像质量产生退化作用，引起影像质量的下降，给遥感影像的后续应用造成严重的障碍。而造成影像模糊的因素通常可以用信息获取与传输过程中各个环节的点扩散函数或者调制传输函数来描述。

对于一个线性位移不变的成像系统来说，成像过程可以用式子(1)来描述：

g＝f*h+n    (1)

其中，f和g分别为理想图像和模糊图像，h为退化函数，也即点扩散函数，n为加性噪声。

式子(1)的傅立叶变换形式为：

G＝FH+N     (2)

对式子(2)进行变换则可以得到对理想图像的一个估计：

$F=G-\frac{N}{H}---\left(3\right)$

这就是最简单的直接逆滤波方法。但是，由于噪声的存在，利用(3)式并不能准确得到理想图像f。可采用一种基于G类点扩散函数的直接解卷积算法来得到更准确的图像估计，下面首先介绍G类点扩散函数的定义。

一个二维对称Le′vyStable概率密度函数h(x，y)通常用其傅立叶(Fourier)变换进行定义：

$\hat{h}(\mathrm{\&epsiv;},\mathrm{\&eta;})={\&Integral;}_{R^2}h(x,y)e^{-2\mathrm{\&pi;i}({\mathrm{\&epsiv;}}^2+{\mathrm{\&eta;}}^2)}\mathrm{dxdy}=e^{-\mathrm{\&alpha;}{({\mathrm{\&epsiv;}}^2+{\mathrm{\&eta;}}^2)}^{\mathrm{\&beta;}}},(\mathrm{\&alpha;}>\mathrm{0,0}<\mathrm{\&beta;}\&le;1)---\left(4\right)$

考虑到成像系统可能由多个满足(4)式的元件组成，因此总的PSF(点扩散函数，Point Spread Function)是由单个PSF通过卷积得到：

$\hat{h}(\mathrm{\&epsiv;},\mathrm{\&eta;})=e^{-{\mathrm{\&Sigma;}}_{i=1}^J{\mathrm{\&alpha;}}_i{({\mathrm{\&epsiv;}}^2+{\mathrm{\&eta;}}^2)}^{{\mathrm{\&beta;}}_i}},({\mathrm{\&alpha;}}_i\&GreaterEqual;\mathrm{0,0}\&le;{\mathrm{\&beta;}}_i\&le;1)---\left(5\right)$

所有傅立叶(Fourier)变换满足(5)式的PSFh(x，y)定义为G类函数。

对于G类点扩散函数，定义卷积积分算子如下

K^tf＝F^-1{H^t(ε，η)F(ε，η)}，0≤t≤1      (6)

G类点扩散函数与传导过程有着密切的联系，根据定理，线性滤波器经局部化和迭代作用的极限状态是传导方程初值问题的解。因此若点扩散函数是高斯函数，那么模糊图像g(x，y)可以看作是t＝1时刻传导方程的解，而所期望的清晰图像f(x，y)则是该传导问题的初值问题

u^t＝λΔu，0≤t ≤1

                                            (7)

u(1)≈g(x，y)

实际上，u(x，y，t)＝K^tf，0≤t≤1，K表示高斯卷积。相应的，对高斯模糊图像进行去模糊在数学上等价于求解(7)的初值。更一般的，对于G类点扩散函数，上面的传导方程就变为一个具有小数指数Laplacian(拉普拉斯算子)算子的广义线性传导方程：

$u_i=-\underset{i=1}{\overset{J}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_i{(-\mathrm{\&Delta;})}^{{\mathrm{\&beta;}}_i}u,{\mathrm{\&lambda;}}_i={\mathrm{\&alpha;}}_i{\left(4{\mathrm{\&pi;}}^2\right)}^{-{\mathrm{\&beta;}}_i},0\&le;t\&le;1---\left(8\right)$

u(x，y，1)≈g(x，y)

因此解方程

$\mathrm{Kf}=\underset{-\&infin;}{\overset{+\&infin;}{\&Integral;}}\underset{-\&infin;}{\overset{+\&infin;}{\&Integral;}}h(x-x^{\&prime;},y-y^{\&prime;})f(x^{\&prime;},y^{\&prime;})\mathrm{dxdy}=h*f=g(x,y)---\left(9\right)$

在数学上等价于对广义传导方程求解其初值u(x，y，0)＝f(x，y)。若f(x，y)已知，u(x，y，t)＝K^tf则是式(8)在时刻t的解。直接解卷积DD算法就是一种来解决式(8)病态问题的正则化方法，并且考虑了t＝1时刻模糊图像数据g(x，y)中存在噪声的情况。u(t)为式(8)的解，ξ，M为已知的正常数，

||u(0)||₂＝||f||₂≤M||u(1)-g||₂＝||n||₂≤ξξ＜＜M    (10)

这里|| ||₂代表L²范数。对于任意常数k＞0且k＜＜M/ξ，定义s^*为

s^*＝log{M/(M-kξ)}/log(M/ξ)         (11)

慢演化约束应用于解病态方程(8)时，就是要求存在一个已知的小常数k(k＞0)和一个已知固定的小常数s(s＞0，s/s^*＞＞1)，满足：

||u(s)-u(0)||₂≤kξ                  (12)

正则化参数k和s表示的是关于病态方程(8)解的一些先验信息。给定k和s，病态方程(8)式逆向问题的SECB(Slow Evolution Constraint-Based，基于约束的慢演化方法)解定义为求u⁺(0)的初值问题，相当于在L²里所有可选择的初值u(0)中选择最合适的u(0)来对下式进行最小化：

${\left|\right|u\left(1\right)-g\left|\right|}_2^2+(\mathrm{\&xi;}/M){\left|\right|u\left(0\right)\left|\right|}^2+k^{-2}{\left|\right|u\left(s\right)-u\left(0\right)\left|\right|}_2^2---\left(13\right)$

直接解卷积DD算法得到的去模糊图像f⁺(x，y)≡u⁺(0)，用其傅立叶(Fourier)形式表达为

$F^{+}(\mathrm{\&epsiv;},\mathrm{\&eta;})=\frac{H^{*}(\mathrm{\&epsiv;},\mathrm{\&eta;})G(\mathrm{\&epsiv;},\mathrm{\&eta;})}{{\left|H(\mathrm{\&epsiv;},\mathrm{\&eta;})\right|}^2+{(\mathrm{\&xi;}/M)}^2+k^{-2}{|1-H(\mathrm{\&epsiv;},\mathrm{\&eta;})|}^2}---\left(14\right)$

其中H^*表示H的共轭。实际上IFFT(逆傅立叶变换)算法得到的f⁺(x，y)同样会出现个别像素灰度为负或超过8位图像最大值255的情况。同样，所有负值置零，所有大于255的值设置为255。

从式子(14)中可以发现，H(ε，η)作为系统的点扩散函数是仅基于观测模糊图像数据经参数估计得到的，未包含气象数据信息，不能准确的描述系统在成像过程中所受到的模糊。

发明内容

本发明的目的是提供一种利用气象数据进行大气模糊遥感影像恢复的方法，更符合实际成像物理过程，大气模糊遥感影像恢复结果更加准确。

一种利用气象数据进行大气模糊遥感影像恢复的方法，其中：包括如下步骤：

1)、进行系统初始化；之后，进入步骤2)；

2)、获取实验数据，实验数据包括多帧遥感影像数据及其相对应的成像时刻的局部区域气象数据，并将获取的数据放入指定的系统数据存储区域；之后，同时进入步骤3)、步骤4)；

3)、读取系统数据存储区域中的遥感影像数据，同时进行对遥感影像数据进行主观质量评价和客观质量评价，其中客观质量评价是利用数字影像的数学参数进行评价，包括信噪比、信息熵、方差、均值、灰度平均梯度中的一种或任意组合的评价，其中，灰度平均梯度对影像进行客观评价的计算公式如下所述：

$\mathrm{GMG}=\frac{1}{(M-1)(N-1)}\underset{i=1}{\overset{M-1}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N-1}{\mathrm{\&Sigma;}}}\sqrt{\frac{{[f(i,j+1)-f(i,j)]}^2+{[f(i+1,j)-f(i,j)]}^2}{2}}$

其中，M，N为影像的尺寸，f(i，j)表示影像中位于坐标(i，j)处的像素值，f(i，j+1)，f(i+1，j)表示该点f(i，j)相邻的像素值；

将当前图像的灰度平均梯度GMG值与系统预设的灰度平均梯度阈值进行比较，若小于，则认为图像质量不达标，直接进入步骤5)；若大于或等于，则认为图像质量达标，进入步骤6)；

4)、读取系统数据存储区域中的气象数据，计算大气调制传输函数，包括湍流调制传输函数和气溶胶调制传输函数，同时进入步骤41)、42)：

41)、计算短曝光湍流调制传输函数为：

${\mathrm{MTF}}_{\mathrm{se}}=\exp \{-57.3v^{5/3}{C}_n^2{\mathrm{\&lambda;}}^{-1/3}R[1-\mathrm{\&mu;}{\left(\frac{\mathrm{\&lambda;\&upsi;}}{D}\right)}^{1/3}\left]\right\}$

其中，v为角空间频率，C_n ²为折射率结构常数，R为光程，单位为m，λ为辐射波长，单位为m，D为光圈直径，单位为m，μ为整个成像系统系数，无量纲，近场时μ＝1，远场时μ＝0.5，其中，基于相对时间概念的计算折射率结构常数C_n ²的公式如下：

$C_n^2=5.9\&times;{10}^{-15}W+1.6\&times;{10}^{-15}T-3.7\&times;{10}^{-15}\mathrm{RH}+6.7\&times;{10}^{-17}{\mathrm{RH}}^2-3.9\&times;{10}^{-19}{\mathrm{RH}}^3$

$-3.7\&times;{10}^{-15}\mathrm{WS}+1.3\&times;{10}^{-15}{\mathrm{WS}}^2-8.2\&times;{10}^{-17}{\mathrm{WS}}^3+2.8\&times;{10}^{-14}\mathrm{SF}-1.8\&times;{10}^{-14}\mathrm{TCSA}$

$+1.4\&times;{10}^{-14}{\mathrm{TCSA}}^2-3.9\&times;{10}^{-13}$

其中，TCSA为中间系数：

TCSA＝9.69×10^-4RH-2.75×10^-5RH²+4.86×10^-7RH³-4.48×10^-9RH⁴+1.66×10^-11RH⁵-6.26×10^-3lnRH-1.34×10^-5SF⁴+7.30×10^-3

其中，RH为相对湿度，SF为太阳辐射通量，单位kw/m²，W为时间权值，无量纲，WS为风速，单位m/s，T为温度，单位为摄氏度℃；

42)计算气溶胶调制传输函数为：

${\mathrm{MTF}}_a=\left\{\begin{array}{c}\exp [-A_{a}R-S_{a}R{(v/v_c)}^2],v\&le;v_c\\ \exp [-(A_a+S_a)R],v>v_c\end{array}\right.$

其中，S_a、A_a分别为大气散射和吸收系数，无量纲，它们之和为衰减系数Q_e，一般情况下，A_a远远小于S_a，将可见光波段的衰减系数近似为Q_e＝3.912/R_m，其中R_m为能见距离，单位为m；

43)由步骤41)、42)得到的湍流调制传输函数和气溶胶调制传输函数，计算大气调制传输函数为：MTF_atmos表示大气调制传输函数，MTF_turb表示湍流调制传输函数，MTF_aero表示气溶胶调制传输函数；

MTF_atmos＝MTF_turb×MTF_aero

5)、对步骤4)中得到的大气调制传输函数通过改进的直接解卷积方法对大气模糊遥感影像进行影像恢复操作，具体步骤如下：

51)、对待恢复的图像f(x，y)进行傅立叶变换，并进行中心化，得到G(ε，η)；进入步骤52)；

52)、对G(ε，η)进行采样，得到G(ε，0)的采样值；进入步骤53)；

53)、利用步骤52)中得到的G(ε，0)的采样值，对式子-α|ε|^2β-A进行最小二乘拟合，得到参数α，β的值，根据G类点扩散函数的定义，系统的点扩散函数即光学传递函数

可表示为：

$\hat{H}(\mathrm{\&epsiv;},\mathrm{\&eta;})=e^{-{\mathrm{\&Sigma;}}_{i=1}^J{\mathrm{\&alpha;}}_i{({\mathrm{\&epsiv;}}^2+{\mathrm{\&eta;}}^2)}^{{\mathrm{\&beta;}}_i}}$

54)、利用式子

对光学传输函数OTF进行纠正；

最后，根据改进的光学传输函数OTF利用下式进行计算：

$F(\mathrm{\&epsiv;},\mathrm{\&eta;})=\frac{H_{\mathrm{new}}^{*}(\mathrm{\&epsiv;},\mathrm{\&eta;})G(\mathrm{\&epsiv;},\mathrm{\&eta;})}{{\left|H_{\mathrm{new}}(\mathrm{\&epsiv;},\mathrm{\&eta;})\right|}^2+{(\mathrm{\&xi;}/M)}^2+k^{-2}{|1-H_{\mathrm{new}}(\mathrm{\&epsiv;},\mathrm{\&eta;})|}^2}$

得到经过处理后的遥感影像，其中，参数ξ，k，M为系统参数，需要在计算过程中进行人工调节，通常设s＝ξ/M，并把初值定为：s＝0.01，k＝100；返回步骤3)；

6)、把最终质量达标的结果数据存入系统指定的图像存储区域，完成恢复后图像的存储。

本发明采用上述技术方案将达到如下的技术效果：

本发明的利用气象数据进行大气模糊遥感影像恢复的方法，是利用估计的大气调制传输函数来对大气模糊遥感影像进行纠正，得到更符合实际成像物理过的光学传递函数，这样，改进本发明方法的影像恢复结果则更加准确。

附图说明

图1为本发明的利用气象数据进行大气模糊遥感影像恢复的方法的流程图；

图2为一个区域一特定时刻的遥感影像成像时刻的气象数据参数和经纬度参数表；

图3为相对时刻权值的对照查询表。

具体实施方式

本发明提供了一种利用气象数据进行大气模糊遥感影像恢复的方法，具体步骤如下：

1)、进行系统初始化；之后，进入步骤2)；

2)、获取实验数据，实验数据包括多帧遥感影像数据及其相对应的成像时刻的局部区域气象数据，并将获取的数据放入指定的系统数据存储区域；之后，同时进入步骤3)、步骤4)；

其中，遥感影像主要是从卫星应用地面站获得的地面遥感影像，数据为*.jpg格式的数字图像，图像大小为4096×4096像素，另外需要获取数据成像波段的波长；而气象数据则需要根据具体遥感影像的成像区域和成像时刻到当地气象观测站获取，或者在获取遥感影像的同时进行气象数据的同步观测与记录，主要包括：遥感影像成像当天的日出时刻、日落时刻、成像时间、日平均温度、成像时刻的风速、成像时刻的相对湿度、成像时刻的太阳辐射照度、成像时刻的能见距离等，如图2表中所示的各数据；

3)、读取系统数据存储区域中的影像数据，首先对其质量进行评价，包括主观质量评价和客观质量评价，其中主观质量评价主要是利用人眼对影像的可观测性进行评价，而客观质量评价是利用数字影像的数学参数进行评价，包括信噪比、信息熵、方差、均值、灰度平均梯度，可采用前述数学参数评价中的一种或多种组合，本发明中主要利用影像的灰度平均梯度来对影像进行客观评价，下面对参数的具体含义与计算公式进行详细介绍：

灰度平均梯度(GMG)能较好的反映图像清晰度和图像纹理的变化特征，特别对于图像的边缘信息，GMG值越大表示图像对应像素点的边缘越清晰，图像的品质越好，相反，图像的品质越差，该灰度平均梯度(GMG)的计算公式为：

$\mathrm{GMG}=\frac{1}{(M-1)(N-1)}\underset{i=1}{\overset{M-1}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N-1}{\mathrm{\&Sigma;}}}\sqrt{\frac{{[f(i,j+1)-f(i,j)]}^2+{[f(i+1,j)-f(i,j)]}^2}{2}}$

其中，M，N为影像的大小，f(i，j)为影像中位于坐标(i，j)点处的像素值，f(i，j+1)，f(i+1，j)表示该点f(i，j)相邻的像素值，图像的边缘也清晰，则每一像素附近的灰度值变化越大，即GMG值越大，

通常，取灰度平均梯度(GMG)的阈值为40，当图像的灰度平均梯度值(GMG)小于该阈值40时，认为图像质量不达标，进入步骤5)，对该遥感影像进行恢复处理；否则，图像质量达标，进入步骤6)，即把结果数据存入输出寄存器单元，准备进行数据的输出与存储；

4)、利用步骤2)中获取的气象数据，计算大气调制传输函数MTF，包括湍流调制传输函数和气溶胶调制传输函数，具体计算步骤如下所示：

大气调制传输函数MTF包括两个分量，即湍流调制传输函数和气溶胶调制传输函数：

湍流调制传输函数根据曝光时间的长短可以分为长曝光湍流调制传输函数和短曝光湍流调制传输函数，在遥感对地观测成像中，曝光时间一般都不超过几毫秒，因此，本发明中的湍流调制传输函数可以视为短曝光湍流调制传输函数，短曝光湍流调制传输函数可以概括为：

${\mathrm{MTF}}_{\mathrm{se}}=\exp \{-57.3v^{5/3}{C}_n^2{\mathrm{\&lambda;}}^{-1/3}R[1-\mathrm{\&mu;}{\left(\frac{\mathrm{\&lambda;\&upsi;}}{D}\right)}^{1/3}\left]\right\}$

其中，v为角空间频率，C_n ²为折射率结构常数，R为光程，单位为m，λ为辐射波长，单位为m，D为光圈直径，单位为m，μ为整个成像系统系数，无量纲，近场时μ＝1，远场时μ＝0.5，其中，光程、辐射波长、光圈直径都是可以获取的，因此对于湍流调制传输函数MTF的估算主要集中在对折射率结构常数C_n ²的估计，首先介绍相对时间的概念，将日出和日落之间时间的1/12作为相对时间的一个小时，称作标准相对时间；而计算相对时间的方法是：用当前时刻减去日出时刻再除以标准相对时间，则得到当前时刻的相对时间；每一个相对时间都对应于一个时间权值，具体的权值可查询图3中的表，而基于相对时间概念的计算折射率结构常数C_n ²的方法：

$C_n^2=5.9\&times;{10}^{-15}W+1.6\&times;{10}^{-15}T-3.7\&times;{10}^{-15}\mathrm{RH}+6.7\&times;{10}^{-17}{\mathrm{RH}}^2-3.9\&times;{10}^{-19}{\mathrm{RH}}^3$

$-3.7\&times;{10}^{-15}\mathrm{WS}+1.3\&times;{10}^{-15}{\mathrm{WS}}^2-8.2\&times;{10}^{-17}{\mathrm{WS}}^3+2.8\&times;{10}^{-14}\mathrm{SF}-1.8\&times;{10}^{-14}\mathrm{TCSA}$

$+1.4\&times;{10}^{-14}{\mathrm{TCSA}}^2-3.9\&times;{10}^{-13}$

其中，TCSA为中间系数：

TCSA＝9.69×10^-4RH-2.75×10^-5RH²+4.86×10^-7RH³-4.48×10^-9RH⁴+1.66×10^-11RH⁵-6.26×10^-3lnRH-1.34×10^-5SF⁴+7.30×10^-3

其中，RH为相对湿度，SF为太阳辐射通量，单位kw/m²，W为时间权值，无量纲，WS为风速，单位m/s，T为温度，单位为摄氏度℃，得到了折射率结构常数C_n ²就可以预测出湍流调制传输函数；

在高空间频率时，气溶胶调制传输函数可以近似为常数，则实际成像系统的气溶胶调制传输函数影响可以用估计的气溶胶调制传输函数来近似表示，其表达式可概括为：

${\mathrm{MTF}}_a=\left\{\begin{array}{c}\exp [-A_{a}R-S_{a}R{(v/v_c)}^2],v\&le;v_c\\ \exp [-(A_a+S_a)R],v>v_c\end{array}\right.$

其中，S_a、A_a分别为大气调制传输函数MTF散射和吸收系数，无量纲，它们之和为衰减系数Q_e，一般情况下，A_a远远小于S_a，气溶胶调制传输函数的衰减系数计算起来相当复杂，为了简化计算，可以将可见光波段的衰减系数近似为Q_e＝3.912/R_m，其中R_m为能见距离，单位m，在得到湍流调制传输函数和气溶胶调制传输函数后，就可以通过下列基于气象数据估计大气调制传输函数MTF的公式来得到大气调制传输函数MTF：

MTF_atmos＝MTF_turb×MTF_aero

MTF_atmos表示大气调制传输函数，MTF_turb表示湍流调制传输函数，MTF_aero表示气溶胶调制传输函数；

5)、利用大气调制传输函数MTF对直接解卷积方法进行改进，并利用改进后的直接解卷积方法对大气调制传输函数MTF模糊遥感影像进行影像恢复操作，具体的实现步骤如下：

首先，对于待恢复的图像数f(x，y)据进行傅立叶变换，并进行中心化，得到G(ε，η)；(傅立叶变换和中心化公式为现有常用公式，这里不再赘述)

其次，对G(ε，η)进行采样，得到G(ε，0)的采样值；

然后，利用G(ε，0)的采样值，对式子-α|ε|^2β-A进行最小二乘拟合，得到参数α，β的值，设点扩散函数为G类点扩散函数，则可估计系统的光学传递函数

为：

$\hat{H}(\mathrm{\&epsiv;},\mathrm{\&eta;})=e^{-{\mathrm{\&Sigma;}}_{i=1}^J{\mathrm{\&alpha;}}_i{({\mathrm{\&epsiv;}}^2+{\mathrm{\&eta;}}^2)}^{{\mathrm{\&beta;}}_i}}$

然后，利用式子

对光学传输函数OTF进行纠正(注：纠正即更改，利用新的表达式替代原来的表达式)；

最后，根据改进的光学传输函数OTF H_new(ε，η)，利用下式进行计算即得到经过处理后的遥感影像：(注，把获取的遥感数据代入下面的式子，计算得到的结果即为恢复后的遥感影像)

$F(\mathrm{\&epsiv;},\mathrm{\&eta;})=\frac{H_{\mathrm{new}}^{*}(\mathrm{\&epsiv;},\mathrm{\&eta;})G(\mathrm{\&epsiv;},\mathrm{\&eta;})}{{\left|H_{\mathrm{new}}(\mathrm{\&epsiv;},\mathrm{\&eta;})\right|}^2+{(\mathrm{\&xi;}/M)}^2+k^{-2}{|1-H_{\mathrm{new}}(\mathrm{\&epsiv;},\mathrm{\&eta;})|}^2}$

其中，参数ξ，k，M为系统参数，需要在计算过程中进行人工调节(注：因为这几个系统参数不可能精确的地确定，所以需要在算法迭代过程中进行调节，可以根据给定的初值为参考进行调节；)，通常设s＝ξ/M，并把初值定为：s＝0.01，k＝100；

之后，返回步骤3)对进行上述恢复后的卫星遥感影像进行主观质量评价和客观质量评价；

6)、把最终质量达标的结果数据存入系统指定的图像存储区域，完成恢复后图像的存储。

本发明是在总结已有恢复算法缺陷的基础上，针对受到大气模糊的中巴地球资源卫星二号星(CBERS-02)的遥感影像数据，提出一种利用卫星成像时刻气象参数估计大气调制传输函数的遥感影像恢复的实用系统，利用遥感影像成像时刻当地的局部气象数据，先对大气湍流调制传输函数和气溶胶调制传输函数进行估计，然后通过这两种调制传输函数的组合模型得到整个大气效应的调制传输函数，进而利用改进的直接解卷积算法实现遥感影像的去大气模糊操作，完成对受到大气湍流模糊的遥感影像的高清晰恢复与重建；经实际应用，恢复效果良好。

Claims (1)

1.一种利用气象数据进行大气模糊遥感影像恢复的方法，其特征在于：包括如下步骤：

1)、进行系统初始化；之后，进入步骤2)；

2)、获取实验数据，实验数据包括多帧遥感影像数据及其相对应的成像时刻的局部区域气象数据，并将获取的数据放入指定的系统数据存储区域；之后，同时进入步骤3)、步骤4)；

3)、读取系统数据存储区域中的遥感影像数据，同时进行对遥感影像数据进行主观质量评价和客观质量评价，其中客观质量评价是利用数字影像的数学参数进行评价，包括信噪比、信息熵、方差、均值、灰度平均梯度中的一种或任意组合的评价，其中，灰度平均梯度对影像进行客观评价的计算公式如下所述：

$\mathrm{GMG}=\frac{1}{(M-1)(N-1)}\underset{i=1}{\overset{M-1}{\mathrm{\&Sigma;}}}\stackrel{N-1}{\underset{j=1}{\mathrm{\&Sigma;}}}\sqrt{\frac{[f(i,j+1)-f(i,j){]}^2+{[f(i+1,j)-f(i,j)]}^2}{2}}$

其中，M，N为影像的尺寸，f(i，j)表示影像中位于坐标.(i，j)处的像素值，f(i，j+1)，f(i+1，j)表示该点f(i，j)相邻的像素值；

将当前图像的灰度平均梯度GMG值与系统预设的灰度平均梯度阈值进行比较，若小于，则认为图像质量不达标，直接进入步骤5)；若大于或等于，则认为图像质量达标，进入步骤6)；

4)、读取系统数据存储区域中的气象数据，计算大气调制传输函数，包括湍流调制传输函数和气溶胶调制传输函数，同时进入步骤41)、42)：

41)、计算短曝光湍流调制传输函数为：

${\mathrm{MTF}}_{\mathrm{se}}=\exp \{-57.3v^{5/3}{C}_n^2{\mathrm{\&lambda;}}^{-1/3}R[1-\mathrm{\&mu;}{\left(\frac{\mathrm{\&lambda;\&upsi;}}{D}\right)}^{1/3}\left]\right\}$

其中，v为角空间频率，C_n ²为折射率结构常数，R为光程，单位为m，λ为辐射波长，单位为m，D为光圈直径，单位为m，μ为整个成像系统系数，无量纲，近场时μ＝1，远场时μ＝0.5，其中，基于相对时间概念的计算折射率结构常数C_n ²的公式如下：

$C_n^2=5.9\&times;{10}^{-15}W+1.6\&times;{10}^{-15}T-3.7\&times;{10}^{-15}\mathrm{RH}+6.7\&times;{10}^{-17}R{H}^2-3.9\&times;{10}^{-19}R{H}^3$

$-3.7\&times;{10}^{-15}\mathrm{WS}+1.3\&times;{10}^{-15}W{S}^2-8.2\&times;{10}^{-17}W{S}^3+2.8\&times;{10}^{-14}\mathrm{SF}-1.8\&times;{10}^{-14}\mathrm{TCSA}$

$+1.4\&times;{10}^{-14}\mathrm{TCS}{A}^2-3.9\&times;{10}^{-13}$

其中，TCSA为中间系数：

TCSA＝9.69×10^-4RH-2.75×10^-5RH²+4.86×10^-7RH³-4.48×10^-9RH⁴+1.66×10^-11RH⁵-6.26×10^-3lnRH-1.34×10^-5SF⁴+7.30×10^-3

其中，RH为相对湿度，SF为太阳辐射通量，单位kw/m²，W为时间权值，无量纲，WS为风速，单位m/s，T为温度，单位为摄氏度℃；

42)计算气溶胶调制传输函数为：

${\mathrm{MTF}}_a=\left\{\begin{array}{c}\exp [-A_{a}R-S_{a}R{(v/v_c)}^2],v\&le;v_c\\ \exp [-(A_a+S_a)R],v>v_c\end{array}\right.$

其中，S_a、A_a分别为大气散射和吸收系数，无量纲，它们之和为衰减系数Q_e，一般情况下，A_a远远小于S_a，将可见光波段的衰减系数近似为Q_e＝3.912/R_m，其中R_m为能见距离，单位为m；

43)由步骤41)、42)得到的湍流调制传输函数和气溶胶调制传输函数，计算大气调制传输函数为：MTF_atmos表示大气调制传输函数，MTF_turb表示湍流调制传输函数，MTF_aero表示气溶胶调制传输函数；

MTF_atmos＝MTF_turb×MTF_aero

5)、对步骤4)中得到的大气调制传输函数通过改进的直接解卷积方法对大气模糊遥感影像进行影像恢复操作，具体步骤如下：

51)、对待恢复的图像f(x，y)进行傅立叶变换，并进行中心化，得到G(ε，η)；进入步骤52)；

52)、对G(ε，η)进行采样，得到G(ε，0)的采样值；进入步骤53)；

53)、利用步骤52)中得到的G(ε，0)的采样值，对式子-α|ε|^2β-A进行最小二乘拟合，得到参数α，β的值，根据G类点扩散函数的定义，系统的点扩散函数即光学传递函数

可表示为：

$\hat{H}(\mathrm{\&epsiv;},\mathrm{\&eta;})=e^{-{\mathrm{\&Sigma;}}_{i=1}^J{\mathrm{\&alpha;}}_i{({\mathrm{\&epsiv;}}^2+{\mathrm{\&eta;}}^2)}^{{\mathrm{\&beta;}}_i}}$

54)、利用式子

对光学传输函数OTF进行纠正；最后，根据改进的光学传输函数OTF利用下式进行计算：

$F(\mathrm{\&epsiv;},\mathrm{\&eta;})=\frac{H_{\mathrm{new}}^{*}(\mathrm{\&epsiv;},\mathrm{\&eta;})G(\mathrm{\&epsiv;},\mathrm{\&eta;})}{{\left|H_{\mathrm{new}}(\mathrm{\&epsiv;},\mathrm{\&eta;})\right|}^2+{(\mathrm{\&xi;}/M)}^2+k^{-2}{|1-H_{\mathrm{new}}(\mathrm{\&epsiv;},\mathrm{\&eta;})|}^2}$

得到经过处理后的遥感影像，其中，参数ξ，k，M为系统参数，需要在计算过程中进行人工调节，通常设s＝ξ/M，并把初值定为：s＝0.01，k＝100；返回步骤3)；

6)、把最终质量达标的结果数据存入系统指定的图像存储区域，完成恢复后图像的存储。

Images