[发明专利]气象数据修正的增量维纳滤波大气模糊遥感影像恢复方法

权利要求书

1.一种气象数据修正的增量维纳滤波大气模糊遥感影像恢复方法，其特征在于：包括如下步骤：

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

2)、读取系统数据存储区域中的气象数据，计算短曝光湍流调制传输函数MTF_se、气溶胶调制传输函数MTF_a和大气调制传输函数MTF_atmos；

3)、对步骤2)中得到的大气调制传输函数MTF_atmos通过修正的增量维纳滤波方法对大气模糊遥感影像进行影像恢复操作，并把结果数据存入系统指定的图像存储区域，完成恢复后图像数据的存储。

2.如权利要求1所述的气象数据修正的增量维纳滤波大气模糊遥感影像恢复方法，其特征在于：步骤2)中所述的计算短曝光湍流调制传输函数MTF_se和气溶胶调制传输函数MTF_a、大气调制传输函数MTF_atmos，具体步骤如下，同时进入步骤21)、22)：

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

MTFse=exp{-57.3v5/3Cn2λ-1/3R[1-μ(λvD)1/3]}]]>

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

Cn2=5.9×10-15W+1.6×10-15T-3.7×10-15RH+6.7×10-17RH2-3.9×10-19RH3]]>

-3.7×10-15WS+1.3×10-15WS2-8.2×10-17WS3+2.8×10-14SF-1.8×10-14TCSA]]>

+1.4×10-14TCSA2-3.9×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为温度，单位为摄氏度℃；

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

MTFa=exp[-AaR-SaR(v/vc)2],v≤vcexp[-(Aa+Sa)R],v>vc]]>

其中，v为角空间频率，R为光程，v_c为气溶胶截止频率；S_a、A_a分别为大气散射系数和大气吸收系数，无量纲，S_a、A_a之和为衰减系数Q_e，一般情况下，A_a远远小于S_a，将可见光波段的衰减系数近似为Q_e＝3.912/R_m，其中R_m为能见距离，单位为米；

23)、由步骤21)、22)得到的短曝光湍流调制传输函数MTF_se和气溶胶调制传输函数MTF_a，计算大气调制传输函数MTF_atmos为：

MTF_atmos＝MTF_se×MTF_a。

3.如权利要求1所述的气象数据修正的增量维纳滤波大气模糊遥感影像恢复方法，其特征在于：步骤3)中所述的对大气调制传输函数MTF_atmos通过修正的增量维纳滤波方法对大气模糊遥感影像进行影像恢复操作，并把结果数据存入系统指定的图像存储区域，完成恢复后图像数据的存储，具体步骤如下：

31)、输入观测得到的原始图像数据f(x，y)，作为初始估计量其中，初始估计量中的k＝0，k表示循环迭代的次数，每完成一个迭代过程，k的值都要增加1；x、y分别表示原始图像中像素的横向、纵向的坐标；

32)、判断当前代表循环迭代次数的k是等于0还是大于0：若k＝0，对点扩散函数h_k(x，y)进行初值设置，令初始的点扩散函数h_k(x，y)中的k＝0，把点扩散函数h_k(x，y)(k＝0)的初始值设置为满足高斯分布且和为1的矩阵，且令点扩散函数h_k(x，y)(k＝0)的尺寸大小等于步骤31)中输入的所观测得到的原始图像的尺寸大小；若k＞0，则令点扩散函数h_k+1(x，y)＝h_k(x，y)，k+1表示当前循环迭代的次数k的值增加1；

33)、对当前估计量当前点扩散函数h_k(x，y)分别进行傅里叶变换，得到当前估计量的傅立叶变换矩阵F_k(p，q)和当前点扩散函数h_k(x，y)的傅立叶变换矩阵即光学传递函数H_k(p，q)；其中，p、q分别表示傅里叶域中的元素坐标；

34)、利用大气调制传输函数MTF_atmos对光学传递函数H_k(p，q)进行修正，其修正模型为：

|Hknew(p,q)|=|Hk(p,q)|+MTFatmos2]]>

H_k^new(p，q)为修正的光学传递函数；

35)、根据式子

Fk+1(p,q)=Fk(p,q)+Hk*(p,q)S1(p,q)|Hk(p,q)|2+Snn(p,q)/Sxx(p,q)]]>

k+1表示当前循环迭代的次数k的值增加1，F_k+1(p，q)表示F_k(p，q)的迭代更新值；H^*_k(p，q)表示光学传递函数H_k(p，q)的复共轭；S_l(p，q)＝F(p，q)-F_k(p，q)H_k(p，q)，F(p，q)为步骤31)中输入的原始图像数据f(x，y)的傅立叶变换形式，S_nn(u，v)和S_xx(u，v)分别表示加性噪声与图像的功率谱；

36)、对F_k+1(p，q)进行逆傅里叶变换，得到对初始估计量(k＝0)的改进值(k＞0)，令改进值(k＞0)的参数k＝k+1，k+1表示当前循环迭代的次数k的值增加1，然后，判断当前循环迭代的次数k是否等于N，N为预先设定的最大迭代次数：若是，终止迭代，得到重建的图像数据f’(x、y)；若否，开始下一次的迭代：将改进值(k＞0)作为估计量，计算点扩散函数FFT^-1表示逆傅立叶变换，作为当前的点扩散函数，返回步骤33)；

37)、把重建的图像数据f’(x、y)存入系统指定的图像存储区域，完成恢复后图像数据的存储。