CN101598543A - 一种实用的遥感影像大气校正方法 - Google Patents

Abstract

一种实用的遥感影像大气校正方法，包含三个步骤：第一步骤是基于浓密植被算法借助于6S模型查找表从MODIS数据反演大气气溶胶光学厚度。第二个步骤是利用MODIS数据波段比值方法计算总大气水蒸汽含量，然后计算大气水蒸汽光学厚度。第三个步骤是将反演得到的大气气溶胶光学厚度和大气水蒸气光学厚度与DOS模型相结合反演地表反射率。本发明计算得到的地表反射率精度高，该方法实用性强，适用于各种陆地卫星影像，而且MODIS数据每天都可以免费获取，因此本方法可以作为一种业务化的遥感影像大气校正方法。

Description

一种实用的遥感影像大气校正方法

技术领域

本发明涉及一种对遥感影像进行大气校正的方法，能够应用在农业、林业、气象、生态环境以及国防军事等遥感部门。

背景技术

遥感影像大气校正是地表反射率反演前的一个重要的预处理步骤。在遥感大气校正研究中，精确的大气校正需要卫星过境时刻同步测量的大气探空数据，而这种数据的获得往往很困难，理想的大气校正方法是仅需要图像本身的信息而不需要卫星过境时刻的同步测量数据，暗目标减法(Dark Object Subtraction，DOS)就是为了满足这种需要而产生的。该方法的基本原理就是在地表均匀、朗伯面反射、大气性质均一，忽略大气多次散射的前提下，假定待校正的遥感图像上存在反射率为0的黑暗像元区域，由于大气的影响，使得这些像元的反射率实际不为0，可以认为增加的这部分辐射是由于大气程辐射而产生的。利用暗目标像元值计算出程辐射，借助于DOS模型即可实现大气校正。DOS大气校正方法自出现以来，由于它的简单实用，获得了广泛的应用。大气透过率是决定DOS大气校正方法校正精度的关键因素。在影响大气透过率的四种因素(大气分子、臭氧含量、气溶胶和总大气水蒸汽含量)中，大气分子和臭氧含量比较稳定，在空间上差异不大，而气溶胶和总大气水蒸汽含量在空间上变动较大，因此必需提高气溶胶光学厚度和大气水蒸汽光学厚度的精度才能提高大气透过率的精度，从而提高最终大气校正的精度。高精度地获取气溶胶光学厚度和大气水蒸汽光学厚度的空间分布非常困难，通常的做法是不考虑这两个因素或者利用气象点数据来近似代替面数据，这样就严重影响了DOS大气校正方法的精度。

因此，为了提高DOS方法大气校正的精度，必须提高气溶胶光学厚度和大气水蒸汽光学厚度的计算精度。本发明针对以往做法的不足，提出利用MODIS数据来获取大面积区域的气溶胶光学厚度和大气水蒸汽光学厚度，然后结合DOS模型进行大气校正。由于MODIS数据每天都可以免费获得，本方法操作起来实用性很强。

发明内容

本发明的目的在于提供一种对遥感影像进行大气校正的方法，以克服现有的DOS模型的缺点。利用本方法进行大气校正能够获得高精度的地表反射率，从而促进相关领域遥感应用的发展。

为实现上述目的，本发明提出的遥感影像大气校正方法为：

第一步、反演大气气溶胶光学厚度

1-1)对MODIS数据进行几何校正

1-2)根据所获影像的特征选择浓密植被区域；

1-3)根据浓密植被区域MODIS数据中红外波段(2.1μm)反射率与红(0.66μm)、蓝(0.47μm)波段反射率之间的关系，由中红外波段的反射率得到红、蓝波段的反射率；

1-4)利用6S模型查找表求得蓝、红波段的气溶胶光学厚度；

1-5)由蓝、红波段的气溶胶光学厚度根据

公式得到各个波段的气溶胶光学厚度；

第二步、反演大气水蒸汽光学厚度

2-1)利用MODIS波段比值方法计算总大气水蒸汽含量；

2-2)由总大气水蒸汽含量计算大气水蒸汽光学厚度；

第三步、大气校正和地表反射率反演

3-1)对陆地卫星影像进行几何校正；

3-2)将陆地卫星影像的亮度值转换为星上辐射亮度；

3-3)计算大气程辐射；

3-4)计算天空光漫射到地表面的光谱辐照度；

3-5)计算臭氧吸收光学厚度和瑞利散射光学厚度；

3-6)计算大气透过率；

3-7)利用DOS模型进行大气校正，得到地表反射率；

3-8)对地表反射率反演结果进行精度验证。

附图说明

图1是本发明的主流程示意图

具体实施方式

遥感地表反射率反演是以辐射传输方程为基础的，在假定地表均匀、朗伯面反射、大气性质均一，忽略大气多次散射的前提下，星上辐射亮度和地表反射率的关系如式1所示[Kaufman Y J，Sendra C.Algorithm for automatic atmospheric correction to visible and near-IRsatellite imagery，International Journal of Remote Sensing，1988，9(8)：1357-1381.]：

$L_{\mathrm{sat}}=L_p+\frac{\mathrm{\ρ}{F}_d{T}_v}{\mathrm{\π}(1-\mathrm{s\ρ})d^2}$ (式1)

式1中：L_sat是星上辐射亮度，L_p是程辐射，Fd_是地表接收到的辐照度，T_v是传感器观测方向的大气透过率，s是大气下界的半球反射率，ρ是地表反射率，d是日地距离。其中F_d＝E_b+E_down，E_down是由天空光漫射到地表面的光谱辐照度，E_b是太阳直射辐照度，E_b＝E₀cos(θ_z)T_z，E₀是大气层外相应波长的太阳光谱辐照度，可由探测器响应函数计算得到。T_z是太阳照射方向上的大气透过率。θ_z是太阳天顶角，由影像头文件获得。因为s值很小，通常可以忽略，所以由式1可以得到地表反射率的计算公式(式2)：

$\mathrm{\ρ}=\frac{\mathrm{\π}(L_{\mathrm{sat}}-L_p)d^2}{T_v(E_0\cos \left({\mathrm{\θ}}_z\right)T_z+E_{\mathrm{down}})}$ (式2)

星上辐射亮度L_sat可以由像元亮度值经辐射定标得到(式3)。

$L_{\mathrm{sat}}=\frac{L_{\max }-L_{\min }}{{\mathrm{QCAL}}_{\max }-{\mathrm{QCAL}}_{\min }}(\mathrm{QCAL}-{\mathrm{QCAL}}_{\min })+L_{\min }$ (式3)

其中，L_sat为星上辐射亮度(w/(m²μmsr))，QCAL为某一像元的亮度值，QCAL_max和QCAL_min分别为像元可以取得的最大和最小亮度值，L_max和L_min分别为QCAL＝QCAL_max和QCAL＝QCAL_min时的光谱辐射亮度值，QCAL_max、QCAL_min、L_max和L_min可以由影像头文件得到。

程辐射L_p由式4计算。

$L_p=\frac{L_{\max }-L_{\min }}{{\mathrm{QCAL}}_{\max }-{\mathrm{QCAL}}_{\min }}{\mathrm{QCAL}}_{\mathrm{dark}}+L_{\min }$ (式4)

QCAL_dark是影像中暗目标的亮度值，暗目标可以选择水体或者阴影。

日地距离d(天文单位)根据式5和式6计算[Spencer J W.Fourier series representation of theposition of the Sun.Search，1971，2(5)：172-201.]

1/d²＝1.000110+0.034221cosΓ+0.001280sinΓ+0.000719cos2Γ+0.000077sin2Γ(式5)

Γ＝2π(d_n-1)/365                                                       (式6)

d_n为儒略日(即遥感影像获取日期距离1月1日的天数)，如果遇到闰年，则用366代替式6中的365。

太阳照射方向和传感器观测方向的大气透过率T_z和T_v是大气光学厚度τ的函数。

T_z＝exp(-τ/cosθ_z)＝exp{(-τ_r-τ_a-τ_o-τ_w)/cosθ_z}    (式7)

T_v＝exp(-τ/cosθ_v)＝exp{(-τ_r-τ_a-τ_o-τ_w)/cosθ_v}    (式8)

其中，τ_r、τ_a、τ_o和τ_w分别是瑞利散射光学厚度、气溶胶光学厚度、臭氧吸收光学厚度和大气水蒸汽光学厚度。θ_v是传感器观测天顶角，由影像头文件获得。

瑞利散射光学厚度τ_r和臭氧吸收光学厚度τ_o相对稳定，可以表达为波长λ和高程h的函数[Stephens G L.Remote Sensing of the Lower Atmosphere.Oxford University Press，1994，pp.523.]、[Guzzi R，Rizzi R and Zibordi G.Atmospheric correction of data measured by aflying platform over the sea：elements of a model and its experimental validation.Applied Optics，1987，26：3043-3051.]：

τ_r＝[0.0088λ^{(-4.15+0.2λ)}][exp(-0.1188h-0.00116h²)]  (式9)

${\mathrm{\τ}}_o=0.03(1-\frac{1.0183}{1+0.0183\exp (h/5)})\exp [-277{(\mathrm{\λ}-0.6)}^2]$ (式10)

其中，λ是影像各波段的中心波长(μm)，h是海拔高度(km)。

本发明的关键在于利用MODIS数据来反演大气气溶胶光学厚度和大气水蒸汽光学厚度，然后和DOS模型结合进行大气校正。

利用MODIS数据反演大气气溶胶光学厚度使用的是浓密植被算法。它的主要思路是基于较低反射率地表(浓密植被区域)在中红外波段(2.1μm)反射率与红(0.66μm)、蓝(0.47μm)可见光波段反射率之间的统计关系[Kaufman Y J，Wald A E，Remer L A，et al.The MODIS 2.1μm channel-correlation with visible reflectance for use in remote sensing of aerosol.IEEETransactions on Geoscience and Remote Sensing，1997，35(5)：1286-1298.](式11和式12)。

ρ_0.47＝ρ_2.1/4     (式11)

ρ_0.66＝ρ_2.1/2     (式12)

中红外波段受大气影响很小，其星上反射率就近似等于地表反射率，由中红外波段地表反射率利用式11和式12确定蓝、红波段地表反射率。然后利用6S模型查找表求得蓝、红波段的气溶胶光学厚度。

[

A.The parameters of atmospheric turbidity.Tellus，1964，16：64-75.]提出大气气溶胶光学厚度随波长的变化可写为：

τ_a(λ)＝βλ^-α                 (式13)

其中，λ为波长(μm)，τ_a(λ)是该波长的气溶胶光学厚度，β为

大气浑浊度系数，它是波长为1μm处大气气溶胶光学厚度，β＜0.1，表示大气非常清洁；β＞0.2，则表示大气是浑浊的。α为波长指数，提供气溶胶粒子尺度信息，α越大，粒子尺度越小，反之亦然。利用两个波长(0.47μm和0.66μm)的气溶胶光学厚度可以按式14和式15得到α和β值，然后利用式13得到各个波长的气溶胶光学厚度。

$\mathrm{\α}=-\frac{{\ln \mathrm{\τ}}_a\left({\mathrm{\λ}}_1\right)-\ln {\mathrm{\τ}}_a\left({\mathrm{\λ}}_2\right)}{{\ln \mathrm{\λ}}_1-{\ln \mathrm{\λ}}_2}$ (式14)

β＝τ_a(λ₁)λ₁ ^α＝τ_a(λ₂)λ₂ ^α  (式15)

水蒸汽光学厚度τ_w可表达为[Leckner B.The spectral distribution of solar radiation at theearth′s surface-elements ofa model.Solar Energy，1978，20：143-150.]：

${\mathrm{\τ}}_w=\frac{0.2385a_{\mathrm{w\λ}}\mathrm{wM}}{{(1+20.07a_{\mathrm{w\λ}}\mathrm{wM})}^{0.45}}$ (式16)

其中w是可降水汽(cm)，a_wλ是水汽吸收系数，各波长的a_wλ值可以查文献得到[Bird R E，Riordan C.Simple solar spectral model for direct and diffuse irradiance on horizontal and tiltedplanes at the earth′s surface for cloudless atmospheres.Journal of climate and applied meteorology，1986，25：87-97.]。相对大气量M由下式获得[Kasten F.A new table and approximate formula forrelative optical air mass.Archives of Meteorology Geophysics and Bioclimatology SeriesB，1966，14(2)：206-223.]

M＝[cosθ_z+0.15(93.885-θ_z)^-1.253]^-1    (式17)

θ_z是太阳天顶角，由影像头文件得到。

其中可降水汽w(cm)与总大气水蒸汽含量W(g/cm²)在数值上相等，可见计算水蒸汽光学厚度τ_w的关键在于获取总大气水蒸汽含量W。

总大气水蒸汽含量可以由MODIS的波段比值方法计算得到[Sobrino J A，Kharraz J EL.Surface temperature and water vapor retrieval from MODIS data[J].International Journal ofRemote Sensing，2003，24(24)：5161-5182.]，具体计算过程如下：

首先按式18-式20计算水汽吸收通道与大气窗口通道的辐射亮度比值G₁₇，G₁₈和G₁₉：

$G_{17}=\frac{L_{17}}{L_2}$ (式18)

$G_{18}=\frac{L_{18}}{L_2}$ (式19)

$G_{19}=\frac{L_{19}}{L_2}$ (式20)

L₁₇、L₁₈和L₁₉分别是MODIS 17、18和19波段传感器接收到的辐射亮度。辐射亮度(L)按式21计算：

L＝radiance_scales(SI-radiance_offsets)    (式21)

其中，radiance_scales和radiance_offsets分别是MODIS数据各波段的辐射率缩放比和辐射率截距，由头文件得到，SI为MODIS影像中的亮度值。

然后，计算从17、18和19波段得到的总大气水蒸汽含量：

$W_{17}=26.314-54.434G_{17}+28.449G_{17}^2$ (式22)

$W_{18}=5.012-23.017G_{18}+27.884G_{18}^2$ (式23)

$W_{19}=9.446-26.887G_{19}+19.914G_{19}^2$ (式24)

在同样的大气情况下，大气水汽在三个通道的吸收情况并不相同，需要为各通道赋予一个权重，这样最终总大气水蒸汽含量(W)按式25确定：

W＝0.192W₁₇+0.453W₁₈+0.355W₁₉    (式25)

如果不考虑地表和大气间的多次反射，天空光漫射到地表面的光谱辐照度E_down可以表示为：

E_down＝E_r+E_a        (式26)

其中，E_r为瑞利散射分量，E_a为气溶胶散射分量。

$E_r=E_0/d^2\cos \left({\mathrm{\θ}}_z\right)T_o{T}_w{T}_{\mathrm{aa}}(1-T_r^{095})0.5$ (式27)

$E_a=E_0/d^2\cos \left({\mathrm{\θ}}_z\right)T_o{T}_w{T}_{\mathrm{aa}}{T}_r^{1.5}(1-T_{\mathrm{as}})F_s$ (式28)

其中T_o、T_w和T_r分别是臭氧吸收透过率、水汽吸收透过率和瑞利散射透过率，T_aa和T_as分别是气溶胶吸收透过率和气溶胶散射透过率，两者之积为气溶胶透过率，F_s是气溶胶散射向下分量的比例，各参数的具体计算参见文献[Bird R E，Riordan C.Simple solar spectral model fordirect and diffuse irradiance on horizontal and tilted planes at the earth′s surface for cloudlessatmospheres.Journal ofclimate and applied meteorology，1986，25：87-97.]。

利用本发明的方法对CBERS 02 CCD陆地卫星影像进行大气校正，利用卫星过境时实地测量的地表反射率对本发明的地表反射率反演结果进行了精度验证，验证结果如表1所示。

表1大气校正前后地表反射率与实测值的对比

ρ^insitu表示实测地表反射率；ρ^TOA表示大气校正前的地表反射率；ρ^TOC表示利用本发明进行大气校正后的地表反射率

从表1可知，经过本发明大气校正方法校正后，各个波段的精度都有所提高。对于测点1，大气校正前的均方根误差(root mean squared error，RMSE)为0.083，大气校正后降为0.022；对于测点2，大气校正前的均方根误差(root mean squared error，RMSE)为0.092，大气校正后降为0.037，表明了本发明大气校正方法的有效性。

Claims (5)

1、一种实用的遥感影像大气校正方法，其步骤为：

第一步、反演大气气溶胶光学厚度

第二步、反演大气水蒸汽光学厚度

第三步、大气校正和地表反射率反演。

2、如权利要求1所述的方法，其中，第一步包括以下五个步骤：

2-1)对MODIS数据进行几何校正

2-2)根据所获影像的特征选择浓密植被区域；

2-3)根据浓密植被区域MODIS数据中红外波段反射率与红、蓝波段反射率之间的关系，由中红外波段的反射率得到红、蓝波段的反射率；

2-4)利用大气辐射传输模型查找表求得蓝、红波段的气溶胶光学厚度；

2-5)由蓝、红波段的气溶胶光学厚度根据

公式得到各个波段的气溶胶光学厚度。

3、如权利要求2所述的方法，其中，第二步包括以下两个步骤：

3-1)利用MODIS波段比值方法计算总大气水蒸汽含量；

3-2)由总大气水蒸汽含量计算大气水蒸汽光学厚度。

4、如权利要求3所述的方法，其中，第三步包括以下八个步骤：

4-1)对陆地卫星影像进行几何校正；

4-2)将陆地卫星影像像元亮度值转换为星上辐射亮度；

4-3)计算大气程辐射；

4-4)计算天空光漫射到地表面的光谱辐照度；

4-5)计算臭氧吸收光学厚度和瑞利散射光学厚度；

4-6)计算大气透过率；

4-7)利用DOS模型进行大气校正，得到地表反射率；

4-8)对地表反射率反演结果进行精度验证。

5、如权利要求4所述的方法，其中大气辐射传输模型为6S模型。

Images