CN103954957A - 一种基于分离大气mtf的遥感影像大气邻近效应校正方法 - Google Patents

Abstract

一种基于分离大气MTF的遥感影像大气邻近效应校正方法，首先分别获取遥感器在轨成像MTF、遥感器静态MTF和卫星平台MTF，其中遥感器在轨MTF利用刃边法从遥感影像上获取，遥感器静态MTF通过遥感器实验室测试获得，卫星平台MTF通过卫星平台在轨测试参数进行仿真获得。然后根据在轨遥感影像退化机理，由遥感器在轨MTF、遥感器静态MTF和卫星平台MTF计算出大气MTF。随后，依据有限长数字滤波器原理利用大气MTF构建消除大气邻近效应的二维滤波卷积核。最后，针对绝对辐射校正后的光学遥感影像，利用MODTRAN大气辐射传输模型进行大气校正，再利用构建的二维滤波卷积核进行卷积，最终获得消除大气邻近效应的遥感影像。

Description

一种基于分离大气MTF的遥感影像大气邻近效应校正方法

技术领域

本发明属于航天遥感技术领域，涉及一种遥感影像的大气散射校正方法，可用于光学遥感图像的大气校正和定量遥感反演。

背景技术

随着高分辨率光学遥感卫星的发展和国外高分辨率遥感影像的竞争，用户对高分辨率光学遥感影像的质量要求也越来越高。目前，高分辨率影像常常出现在轨MTF较低和图像模糊的现象，这其中大气的影响不容忽视。

目前，消除大气影响的图像复原方法很多，但多基于对图像的滤波处理，与遥感定量应用缺少有效、深入的结合，一方面导致图像复原比较盲目，另一方面这些图像复原技术将会影响影像的定量应用，以致无法处理多光谱影像。为了更好地消除大气影响，确保卫星在轨成像质量和影像的定量应用，基于大气散射机理的邻近效应校正研究已凸显其重要性，特别是实时对大气邻近效应的在轨测量，将更有利于实时遥感影像的大气校正。

针对光学遥感影像邻近效应校正的研究，国内外研究人员已经开展了大量工作。目前，遥感影像邻近效应校正主要通过大气点扩展函数来实现。2011年49卷7期《IEEE Transactions On Geoscience And Remote Sensing》上，Semenov等通过从影像估算归一化大气点扩展函数，从而构建了一个空间滤波器进行邻近效应校正，以达到影像复原和提高图像质量的效果。这种直接从影像上估算归一化大气点扩展函数的方法，没有借助遥感器的相关先验知识，属于一种盲提取，同时估算过程也较复杂；同时文章中也指出它的校正处理是一种盲复原，并没有从大气辐射传输对地表反射能量的定量影响进行分析。

发明内容

本发明解决的技术问题是：克服现有技术的不足，提供了一种基于分离大气MTF的遥感影像大气邻近效应校正方法，可以解决现有遥感影像大气影响消除方法中存在的盲复原、未考虑定量应用和大气点扩展函数估算不准确的问题。

本发明的技术解决方案是：一种基于分离大气MTF的遥感影像大气邻近效应校正方法，步骤如下：

(1)获取遥感影像的在轨MTF，记为MTF_on-orbit；

(2)获取遥感器的静态MTF，记为MTF_sensor；

(3)获取卫星平台的MTF，记为MTF_sat；

(4)根据步骤(1)、(2)、(3)的结果，根据公式计算得到大气的MTF，记为MTF_atm，MTF_on-orbit＝MTF_atm×MTF_sensor×MTF_sat；

(5)利用获取的MTF_atm，计算得到一维卷积系数h(k)，|H(e^jw)|＝A(w,k)*h(k)，其中h(k)是中心对称的，h(k)的长度为M；w为频率，|H(e^jw)|为MTF_atm的大小，k取1,2,…,

(6)利用得到的一维卷积系数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}}{\mathrm{\Σ}}}\underset{t_y=0}{\overset{t_x}{\mathrm{\Σ}}}x(t_x,t_y)\×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$

式中 $\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}}{\mathrm{\Σ}}}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)为卷积系数的重复次数，f_x,f_y在零频到奈奎斯特频率范围内选取；

(7)实时获取遥感影像，结合MODTRAN大气辐射传输模型对遥感影像进行大气校正，从中获得未消除邻近效应的地物平均反射率ρ；

(8)将ρ表示为目标反射率ρ_t和背景平均反射率ρ_e两者之间与大气点扩展函数PSF_atm(x,y)之间的卷积关系：

ρ＝ρ_t·∫∫_tPSF_atm(x,y)dxdy+ρ_e·∫∫_ePSF_atm(x,y)dxdy

利用步骤(6)获取的二维卷积核h(i,j)对地物平均反射率ρ进行卷积，即可得到消除邻近效应的最终目标反射率影像：

ρ_t＝ρ*h(i,j)

其中，∫∫_tPSF_atm(x,y)dxdy表示大气点扩展函数对目标像元的影响权重，∫∫_ePSF_atm(x,y)dxdy)表示大气点扩展函数引起邻近像元对目标像元的影响权重，x和y分别表示像元在遥感影像水平和垂直方向上离目标像元的距离，下标t和e分别表示目标和背景环境，且∫∫_tPSF_atm(x,y)dxdy+∫∫_ePSF_atm(x,y)dxdy)＝1。

所述步骤(1)中获取遥感影像的在轨MTF的方法为刃边法或者脉冲法。

所述步骤(2)中获取遥感器的静态MTF的方法为矩形靶标法或者低频靶标法。

所述步骤(3)中获取卫星平台的MTF的方法为：

MTF_sat＝MTF_linear×MTF_random

其中MTF_sat为卫星平台的MTF，MTF_linear为卫星线性运动引起的MTF，MTF_random为随机颤振引起的MTF。

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

(1)现有的直接从影像上估算归一化大气点扩展函数的方法，没有借助遥感器的相关先验知识，属于一种盲提取，同时估算过程也较复杂。本发明方法通过从影像上在轨MTF的获取和遥感器静态MTF的测试来得到实时大气的MTF，并在此基础上进行临近效应校正，方法简单且保证了结果的可靠性和准确性；

(2)现有的邻近效应校正处理方法是基于大气点扩展函数生成的空间滤波器实现的，是一种盲复原方法，没有从大气辐射传输对地表反射能量的定量影响进行分析。本发明方法从大气辐射传输理论出发，考虑在大气散射作用下邻近像元对目标像元的贡献权重，从而利用获取的大气MTF构建卷积核，进行邻近效应校正，不仅可以提升图像质量，也有利于提高遥感影像的定量化反演精度。

附图说明

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

具体实施方式

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

(1)获取遥感影像的在轨MTF。

利用目前较为常用的刃边法或脉冲法，可以获取遥感影像的在轨MTF，表示为MTF_on-orbit。

刃边法的原理及具体步骤可以参考《航天返回与遥感》2009年第30卷第2期上的《基于刃边法的航天光学遥感器在轨MTF测试研究》的文章。

脉冲法的原理及具体步骤可以参考Stennis Space Center2008年8月的报告《Quickbird Satellite in-orbit Modulation Transfer Function(MTF)Measurement Using Edge,Pulse and Impulse Methods for Summer2003》。

(2)获取遥感器的MTF。

在实验室利用高对比度矩形靶标或者低频靶标进行测量，可以获得遥感器的MTF，表示为MTF_sensor。采用高对比度矩形靶标或者低频靶标测量遥感器MTF的方法可参考《航天返回与遥感》2006年第27卷第4期上的《航天光学遥感器MTF测试技术研究》的文章。

(3)获取卫星平台的MTF。

通过卫星平台在轨参数可以计算卫星平台的MTF。卫星平台对影像的影响主要表现为积分时间内卫星线性运动产生的模糊和随机颤振引起的模糊。线性运动的模糊可以认为是一个矩形波函数，主要与卫星的速度和积分时间有关。随机颤振可认为是一个高斯形点扩展函数，主要与积分时间的颤振量有关，通过这些参数的获取即可计算卫星平台的MTF。

MTF_sat＝MTF_linear×MTF_random

其中MTF_sat为卫星平台的MTF，MTF_linear为卫星线性运动引起的MTF，MTF_random为随机颤振引起的MTF。

卫星平台线性运动引起的MTF为MTF_linear(f)：

MTF_linear(f)＝sinc(π·f·v·t_int)

其中，f为空间频率，sinc()表示sinc函数，v表示卫星运行速度，t_int为遥感器的积分时间。

卫星颤振的MTF为MTF_random(f)：

${\mathrm{MTF}}_{\mathrm{random}}\left(f\right)=e^{-2{[\mathrm{\π}\·\mathrm{jitter}\·p\·f]}^2}$

其中，p为遥感器的探测器像元尺寸，jitter为卫星的随机振动rms值，由卫星测试数据统计获得。

(4)根据在轨遥感影像退化机理，遥感影像的在轨MTF由大气MTF(表示为MTF_atm)、遥感器的MTF和卫星平台的MTF构成，因此由前三步得到的结果可计算出大气MTF。

MTF_on-orbit＝MTF_atm×MTF_sensor×MTF_sat

(5)利用大气MTF构建用于消除大气邻近效应的二维卷积核。具体步骤如下：

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

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

其中，w为选取的频率值(从零频到奈奎斯特频率范围内选取)，由此可以根据|H(e^jw)|和A(w,k)计算得到序列h(k)。

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

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

$\begin{array}{c}\left|H(f_x,f_y)\right|=|\underset{n1=0}{\overset{\frac{M-1}{2}}{\mathrm{\Σ}}}\underset{n2=0}{\overset{\frac{M-1}{2}}{\mathrm{\Σ}}}h(n1,n2)\×4\cos (f_x\×n1)\×\cos (f_y\×n2)|\\ =\underset{t_x=0}{\overset{\frac{M-1}{2}}{\mathrm{\Σ}}}\underset{t_y=0}{\overset{t_x}{\mathrm{\Σ}}}x(t_x,t_y)\×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\\ =A(f_x,f_y,t_x,t_y)*x(t_x,t_y)\end{array}$

其中，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。

iii)根据得到的一维卷积系数计算|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}}{\mathrm{\Σ}}}h\left(k\right)\·2\cos (\frac{M-1}{2}-k)\sqrt{{f_x}^2+{f_y}^2}$

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

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

(6)结合MODTRAN大气辐射传输模型对光学遥感影像进行大气校正，通过步骤(5)获得的二维卷积核消除由大气漫散射引起的邻近效应，获得消除邻近效应的遥感影像。

下面具体分三步进行操作：

i)计算遥感影像的入瞳辐亮度和表观反射率。

首先通过遥感影像定标系数将影像的DN值转换成入瞳辐射亮度L，从而根据下式计算影像表观反射率ρ_TOA：

${\mathrm{\ρ}}_{\mathrm{TOA}}=\frac{\mathrm{\π}{\mathrm{Ld}}^2}{E_s\cos {\mathrm{\θ}}_s}$

其中，d为日地距离(JD为卫星过境的Julian日)，E_s为太阳平均辐照度，θ_s为太阳天顶角。

ii)进行未考虑邻近效应的大气校正。

由大气辐射传输理论，可用如下形式表示入瞳反射率，即大气层顶的表观反射率：

${\mathrm{\ρ}}_{\mathrm{TOA}}={\mathrm{\ρ}}_p+{\mathrm{\ρ}}_t\frac{T\left({\mathrm{\θ}}_s\right)e^{-\mathrm{\τ}/{\mathrm{\μ}}_v}}{1-{\mathrm{\ρ}}_{e}S}+{\mathrm{\ρ}}_e\frac{T\left({\mathrm{\θ}}_s\right)t_d\left({\mathrm{\θ}}_v\right)}{1-{\mathrm{\ρ}}_{e}S}$

其中，ρ_TOA为大气层顶的表观反射率，ρ_p为大气自身反射率，ρ_t为目标反射率，ρ_e为背景平均反射率，T(θ_s)为太阳天顶角方向的大气透过率，τ为气溶胶光学厚度，θ_v为遥感器观测天顶角，μ_v为θ_v的余弦，S为大气半球反照率，t_d(θ_v)为θ_v方向的大气上行散射透过率。

由于在实际大气校正中，相邻像元反射特性在大气邻近效应影响下，对目标反射的影响不容易估算，常常是将受大气邻近效应影响的地物平均反射率ρ作为校正对象，从而将上式进行简化，即有如下关系：

${\mathrm{\ρ}}_{\mathrm{TOA}}={\mathrm{\ρ}}_p+\frac{\mathrm{\ρ}}{1-\mathrm{\ρS}}T\left({\mathrm{\θ}}_s\right)T\left({\mathrm{\θ}}_v\right)$

其中，T(θ_v)为遥感器观测天顶角方向的大气透过率。

利用大气辐射传输模型MODTRAN可获得ρ_p、S、T(θ_s)、T(θ_v)等大气参数，从而可求得受邻近效应影响的地物平均反射率ρ：

$\mathrm{\ρ}=\frac{{\mathrm{\ρ}}_{\mathrm{TOA}}-{\mathrm{\ρ}}_p}{({\mathrm{\ρ}}_{\mathrm{TOA}}-{\mathrm{\ρ}}_p)S+T\left({\mathrm{\θ}}_s\right)T\left({\mathrm{\θ}}_v\right)}$

iii)消除大气邻近效应。

在实际辐射传输中，ρ与目标反射率ρ_t和背景平均反射率ρ_e之间有着如下关系：

ρ＝ρ_t·∫∫_tPSF_atm(x,y)dxdy+ρ_e·∫∫_ePSF_atm(x,y)dxdy

其中，∫∫_tPSF_atm(x,y)dxdy代表了大气点扩展函数对目标像元影响的权重，∫∫_ePSF_atm(x,y)dxdy)代表了大气点扩展函数引起的邻近像元对目标像元影响的权重，x和y分别表示像元在影像水平和垂直方向上离目标像元的距离，下标t和e分别表示目标和背景环境，且∫∫_tPSF_atm(x,y)dxdy+∫∫_ePSF_atm(x,y)dxdy)＝1。

通过步骤(5)建立的二维滤波核h(i,j)对地物平均反射率影像进行卷积(*表示卷积)，即可得到消除邻近效应的最终目标反射率影像ρ_t：

ρ_t＝ρ*h(i,j)

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

Claims (4)

1.一种基于分离大气MTF的遥感影像大气邻近效应校正方法，其特征在于步骤如下：

(1)获取遥感影像的在轨MTF，记为MTF_on-orbit；

(2)获取遥感器的静态MTF，记为MTF_sensor；

(3)获取卫星平台的MTF，记为MTF_sat；

(4)根据步骤(1)、(2)、(3)的结果，根据公式计算得到大气的MTF，记为MTF_atm，MTF_on-orbit＝MTF_atm×MTF_sensor×MTF_sat；

(5)利用获取的MTF_atm，计算得到一维卷积系数h(k)，|H(e^jw)|＝A(w,k)*h(k)，其中h(k)是中心对称的，h(k)的长度为M；w为频率，|H(e^jw)|为MTF_atm的大小，k取1,2,…,

(6)利用得到的一维卷积系数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}}{\mathrm{\Σ}}}\underset{t_y=0}{\overset{t_x}{\mathrm{\Σ}}}x(t_x,t_y)\×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$

式中 $\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}}{\mathrm{\Σ}}}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)为卷积系数的重复次数，f_x,f_y在零频到奈奎斯特频率范围内选取；

(7)实时获取遥感影像，结合MODTRAN大气辐射传输模型对遥感影像进行大气校正，从中获得未消除邻近效应的地物平均反射率ρ；

(8)将ρ表示为目标反射率ρ_t和背景平均反射率ρ_e两者之间与大气点扩展函数PSF_atm(x,y)之间的卷积关系：

ρ＝ρ_t·∫∫_tPSF_atm(x,y)dxdy+ρ_e·∫∫_ePSF_atm(x,y)dxdy

利用步骤(6)获取的二维卷积核h(i,j)对地物平均反射率ρ进行卷积，即可得到消除邻近效应的最终目标反射率影像：

ρ_t＝ρ*h(i,j)

其中，∫∫_tPSF_atm(x,y)dxdy表示大气点扩展函数对目标像元的影响权重，∫∫_ePSF_atm(x,y)dxdy)表示大气点扩展函数引起邻近像元对目标像元的影响权重，x和y分别表示像元在遥感影像水平和垂直方向上离目标像元的距离，下标t和e分别表示目标和背景环境，且∫∫_tPSF_atm(x,y)dxdy+∫∫_ePSF_atm(x,y)dxdy)＝1。

2.根据权利要求1所述的一种基于分离大气MTF的遥感影像大气邻近效应校正方法，其特征在于：所述步骤(1)中获取遥感影像的在轨MTF的方法为刃边法或者脉冲法。

3.根据权利要求1或2所述的一种基于分离大气MTF的遥感影像大气邻近效应校正方法，其特征在于：所述步骤(2)中获取遥感器的静态MTF的方法为矩形靶标法或者低频靶标法。

4.根据权利要求1或2所述的一种基于分离大气MTF的遥感影像大气邻近效应校正方法，其特征在于：所述步骤(3)中获取卫星平台的MTF的方法为：

MTF_sat＝MTF_linear×MTF_random

其中MTF_sat为卫星平台的MTF，MTF_linear为卫星线性运动引起的MTF，MTF_random为随机颤振引起的MTF。