CN102288956A - 一种遥感卫星多光谱数据的大气订正方法 - Google Patents

Abstract

一种遥感卫星多光谱数据的大气订正方法，(1)读取数据，将数据转换为表观辐亮度数据和表观反射率数据；(2)依据表观反射率数据计算决策树法中的决策因子：比值植被指数、土壤可调植被指数和归一化水体指数，通过决策树法实现暗目标自动提取；(3)根据表观反射率数据获取上述暗目标区域的表观反射率数据，并计算大气路径反射率；(4)根据大气路径反射率和瑞利反射率比值与气溶胶光学厚度的关系，确定多光谱各波段的气溶胶光学厚度；并依据多光谱各波段的气溶胶光学厚度和瑞利散射光学厚度，确定出550nm处总的光学厚度；(5)依据上述550nm处总的光学厚度，利用查找表获取反演参数，通过辐射传输方程求解地表反射率，完成大气订正。

Description

一种遥感卫星多光谱数据的大气订正方法

技术领域

本发明涉及一种遥感卫星的大气订正方法，特别涉及一种与环境减灾卫星传感器波段设置相类似卫星的多光谱数据大气订正方法。

背景技术

环境减灾卫星自发射以来，在汶川地震连续检测、印度洋海啸灾情监测、澳大利亚火灾监测、日本地震和海啸灾情监测等灾害监测中发挥了独特的作用，并广泛应用于农业、林业、水利、土地利用、地况调查、城市规划等众多领域，取得了显著的效果。为了实现对不同地点、不同时间和不同的传感器获取的遥感数据进行比较和应用，及时为国家环境保护和灾害预报、抗灾、救灾提供准确、科学的决策依据，这必然要求遥感数据定量化。

对于遥感定量化，就要详细考虑大气对于成像波段传输过程的定量影响，对其影响做出准确的订正，实现地物目标辐射特性真实再现，达到提高遥感信息质量及定量化应用水平目的。

目前国外大气订正方面的成果很多，像LANDSAT、SPOT、MODIS、ASTER等其他国外遥感卫星都针对各自卫星的特点开展了相应的大气订正；国内学者也在CBERS-01卫星、CBERS-02、02B卫星等国产遥感卫星数据的大气订正研究产生了一些相关算法，但由于大气状况瞬息万变，当时、当地的大气对遥感数据影响程度和重要性不同，再加上各卫星的传感器波段设置不同，特性不同，研究目的和要求也不同，即使相同地区，大气状况也是复杂多变的，每景图像成像时的大气状况也存在明显差异，很难有一种算法是普遍适用的。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提供一种普适性强的遥感卫星多光谱数据大气订正方法。

本发明的技术解决方案是：一种遥感卫星多光谱数据的大气订正方法，步骤如下：

(1)读取遥感卫星多光谱数据，对该多光谱数据进行预处理，将多光谱数据转换为表观辐亮度数据和表观反射率数据；

(2)依据表观反射率数据计算决策树法中的决策因子：比值植被指数RVI、土壤可调植被指数SAVI和归一化水体指数NDWI，通过决策树法实现DDV暗目标自动提取；

(3)根据步骤(1)中的表观反射率数据获取上述提取的DDV暗目标区域的DDV表观反射率数据；并根据该DDV表观反射率数据结合DDV地表反射率计算Path反射率；

(4)根据Path反射率和Rayleigh反射率比值与气溶胶光学厚度的关系，确定多光谱各波段的气溶胶光学厚度；并依据多光谱各波段的气溶胶光学厚度和Rayleigh散射光学厚度，确定出550nm处总的光学厚度τ；

(5)依据上述550nm处总的光学厚度，利用查找表获取反演参数ρ_path(λ)、T(λ)、S(λ)，通过辐射传输方程求解地表反射率，完成大气订正。

所述步骤(3)中的DDV地表反射率的计算公式如下：

ρ_DDV(λ)＝S_vegρ_veg(λ)+(1-S_veg)ρ_soil(λ)

其中：ρ_DDV(λ)对应多光谱波段λ的DDV地表反射率；

ρ_veg(λ)、ρ_soil(λ)是实测多光谱波段λ的植被反射率和土壤反射率；

S_veg是植被调整因子，

所述步骤(5)中查找表的索引τ分别与ρ_path(λ)、T(λ)、S(λ)之间的关系满足下式：

ρ_Path(λ)＝c₀(λ)+c₁(λ)τ+c₂(λ)τ²

S(λ)＝d₀(λ)+d₁(λ)τ+d₂(λ)τ²

T(λ)＝e₀(λ)+e₁(λ)τ+e₂(λ)τ²

其中，c₀(λ)、c₁(λ)、c₂(λ)、d₀(λ)、d₁(λ)、d₂(λ)、e₀(λ)、e₁(λ)、e₂(λ)对应多光谱波段λ在确定的卫星成像几何参数太阳天顶角θ_S、卫星观测天顶角θ_V、太阳与卫星的方位角下的路径反射率ρ_path(λ)、大气透过率T(λ)、半球反射率S(λ)与总光学厚度τ的系数。

本发明与现有技术相比有益效果为：

(1)本发明基于浓密植被暗目标区(DDV-density dark vegetation)反射率与土壤可调植被指数及波谱库实测植被和土壤地表反射率的相对稳定关系，利用辐射传输方程确定出路径反射率，然后利用为环境减灾卫星多光谱传感器建立的路径反射率和瑞利散射反射率比值与气溶胶光学厚度关系的查找表确立气溶胶光学厚，为此实现与环境减灾卫星传感器波段设置相似类卫星的大气订正算法，方法完整、合理可行、精度较高。

(2)本发明针对遥感卫星多光谱波段设置特点，依据不同地物特性，选取大量遥感数据，分析多种植被指数和水体指数，甄选出比值植被指数RVI、土壤可调植被指数SAVI和归一化水体指数NDWI作为决策树法的决策因子，实现了暗目标自动提取，使此类卫星多光谱数据的大气订正更具有普适性。

(3)本发明利用MODTRAN辐射传输模型，分析遥感卫星在不同几何状态下的各种情况，针对此类卫星多光谱数据波段设置特点，为大气订正各反演参数与主变量间建立二次线性查找表，有效地缩减了大气订正产品生产时间，使得业务化大气订正得以实现。

(4)本发明铜鼓大量实测地物波谱数据，结合遥感卫星地表反演量，针对多光谱数据的DDV暗目标区，分析实测地物波谱与多种植被指数的线性关系，确立了用土壤可调植被指数SAVI作为植被实测植被和土壤地表反射率的调整因子，使DDV暗目标区的地表反射率计算更加合理，确保了方法的合理可行性。

附图说明

图1为本发明总体流程图；

图2为本发明具体流程图；

图3为本发明大气订正结果对比曲线。

具体实施方式

如图1、2所示，本发明一种遥感卫星多光谱数据的大气订正方法，步骤如下：

(1)读取遥感卫星多光谱数据，对该多光谱数据进行预处理，将多光谱数据转换为表观辐亮度数据和表观反射率数据；

(1.1)表观辐亮度

读取多光谱CCD数据，依据数据附带文件中的绝对定标系数，按照公式①计算多光谱传感器入瞳处辐亮度即表观辐亮度：

$L\left(\mathrm{\λ}\right)=\frac{\mathrm{DN}\left(\mathrm{\λ}\right)}{A\left(\mathrm{\λ}\right)}+L_0\left(\mathrm{\λ}\right)$ ①

其中：λ是多光谱的波段；

L(λ)对应多光谱波段λ的入瞳处辐亮度；

DN(λ)对应多光谱波段λ的数字DN值；

A(λ)、L₀(λ)对应多光谱波段λ的定标系数。

(1.2)表观反射率数据

将上述计算得到的传感器入瞳处辐射亮度，按公式②计算得到传感器入瞳处的表观反射率数据：

${\mathrm{\ρ}}_{\mathrm{TOA}}\left(\mathrm{\λ}\right)=d^2\×\frac{\mathrm{\πL}\left(\mathrm{\λ}\right)}{E_S\left(\mathrm{\λ}\right)\cos {\mathrm{\θ}}_S}$ ②

其中：ρ_TOA(λ)对应多光谱波段λ的表观反射率；

d²是日-地距离校正因子可通过公式③计算

$d^2=\frac{1}{1-0.01673\cos \left[0.9856\frac{\mathrm{\π}(J-4)}{180}\right]}$ ③

J是儒略历天数；

E_S(λ)对应多光谱波段λ的大气外太阳光谱辐照度；

θ_S是卫星成像时的太阳天顶角，可从数据附带文件中成像几何参数中获取；

上述数据附带文件是指遥感卫星成像时的已知状态参数，根据卫星下传数据得到，对本领域技术人员来讲为已知量。

(2)依据表观反射率数据计算决策树法中的决策因子：比值植被指数RVI、土壤可调植被指数SAVI和归一化水体指数NDWI，通过决策树法实现DDV暗目标自动提取；

$\mathrm{RVI}=\frac{{\mathrm{\ρ}}_{\mathrm{TOA}}\left({\mathrm{\λ}}_{\mathrm{nir}}\right)}{{\mathrm{\ρ}}_{\mathrm{TOA}}\left({\mathrm{\λ}}_{\mathrm{red}}\right)}$ ④

$\mathrm{SAVI}=(1+C)\frac{{\mathrm{\ρ}}_{\mathrm{TOA}}\left({\mathrm{\λ}}_{\mathrm{nir}}\right)-{\mathrm{\ρ}}_{\mathrm{TOA}}\left({\mathrm{\λ}}_{\mathrm{red}}\right)}{{\mathrm{\ρ}}_{\mathrm{TOA}}\left({\mathrm{\λ}}_{\mathrm{nir}}\right)+{\mathrm{\ρ}}_{\mathrm{TOA}}\left({\mathrm{\λ}}_{\mathrm{red}}\right)+C}$ ⑤

$\mathrm{NDWI}=\frac{{\mathrm{\ρ}}_{\mathrm{TOA}}\left({\mathrm{\λ}}_{\mathrm{green}}\right)-{\mathrm{\ρ}}_{\mathrm{TOA}}\left({\mathrm{\λ}}_{\mathrm{red}}\right)}{{\mathrm{\ρ}}_{\mathrm{TOA}}\left({\mathrm{\λ}}_{\mathrm{green}}\right)+{\mathrm{\ρ}}_{\mathrm{TOA}}\left({\mathrm{\λ}}_{\mathrm{red}}\right)}$ ⑥

其中：ρ_TOA(λ_green)是多光谱数据绿光波段的表观反射率；

ρ_TOA(λ_red)是多光谱数据红光波段的表观反射率；

ρ_TOA(λ_nir)是多光谱数据近红外波段的表观反射率；

C是土壤调整因子，C∈[0，1]。

(3)针对此类卫星多光谱传感器的响应特性，确立大气组成中吸收、散射等衰减对多光谱数据的主要影响因子(Rayleigh散射衰减和气溶胶衰减)，依此根据步骤(1)中的表观反射率数据获取上述提取的DDV暗目标区域的DDV表观反射率数据ρ_Path-DDV(λ)；并根据该DDV表观反射率数据结合DDV地表反射率计算Path反射率；

(3.1)瑞利(Rayleigh)散射反射率

利用公式⑦，计算DDV暗目标选取区相对应DEM高程数据处多光谱波段λ的Rayleigh散射反射率：

ρ_Ray-DDV(λ)＝a₀(λ)+a₁(λ)H+a₂(λ)H²⑦

其中：H是DEM高程数据获得海拔高度；

a₀(λ)、a₁(λ)、a₂(λ)对应多光谱波段λ的Rayleigh散射反射率与海拔高度H关系的系数；

(3.2)DDV地表反射率

利用大量实测地物波谱库，建立DDV地表反射率与土壤可调植被指数SAVI的拟合关系，由公式⑧确定多光谱波段λ的DDV地表反射率：

ρ_DDV(λ)＝S_vegρ_veg(λ)+(1-S_veg)ρ_soil(λ)⑧

其中：ρ_DDV(λ)对应多光谱波段λ的DDV地表反射率；

ρ_veg(λ)、ρ_soil(λ)是实测地物波谱库中对应多光谱波段λ的植被反射率和土壤反射率；

S_veg是植被调整因子，

⑨

上述实测地物波谱库是由不同时间、不同地区、不同类型的地物(植被、水体、土壤等)实地测量组成。

(3.3)Path反射率

在DDV暗目标选取区，利用公式⑩计算多光谱波段λ的大气路径反射率：

ρ_Path-DDV(λ)＝ρ_TOA-DDV(λ)-ρ_DDV(λ)⑩

(4)根据Path反射率和Rayleigh反射率比值与气溶胶光学厚度的关系，确定多光谱各波段的气溶胶光学厚度；并依据多光谱各波段的气溶胶光学厚度和Rayleigh散射光学厚度，确定出550nm处总的光学厚度τ；

(4.1)多光谱波段λ的气溶胶光学厚度

利用DDV暗目标选取区Path路径反射率和Rayleigh散射反射率的比值与气溶胶光学厚度关系，通过公式

求出对应多光谱波段λ的气溶胶光学厚度：

$\frac{{\mathrm{\ρ}}_{\mathrm{Path}-\mathrm{DDV}}\left(\mathrm{\λ}\right)}{{\mathrm{\ρ}}_{\mathrm{Ray}-\mathrm{DDV}}\left(\mathrm{\λ}\right)}=b_0\left(\mathrm{\λ}\right)+b_1\left(\mathrm{\λ}\right){\mathrm{\τ}}_{\mathrm{Aer}}\left(\mathrm{\λ}\right)+b_2\left(\mathrm{\λ}\right){\mathrm{\τ}}_{\mathrm{Aer}}{\left(\mathrm{\λ}\right)}^2$

其中：τ_Aer(λ)对应多光谱波段λ的气溶胶光学厚度；

b₀(λ)、b₁(λ)、b₂(λ)对应多光谱波段λ的Path路径反射率和Rayleigh散射反射率的比值与气溶胶光学厚度关系的系数。

(4.2)550nm处的气溶胶光学厚度τ₅₅₀

依据气溶胶光学厚度与波长之间的

公式

将各波段到的气溶胶光学厚度进行最小二乘法拟合，以确定

公式系数α、β：

τ_Aer(λ)＝αλ^-β

由上式确定出550nm处的气溶胶光学厚度τ₅₅₀：

τ_Aer(550)＝α×550^-β

(4.3)550nm处的Rayleigh散射光学厚度τ_Ray(550)

依据DEM高程数据，通过公式

可求550nm处的Rayleigh散射光学厚度τ_Ray(550)：

τ_Ray(550)＝f₀(550)+f₁(550)H+f₂(550)H²

其中：f₀(550)、f₁(550)、f₂(550)对应550nm处Rayleigh散射光学厚度与海拔高度H大气关系的系数。

(4.4)大气总的光学厚度τ

由上可知550nm处大气总的光学厚度τ为：

τ＝τ_Aer (550)+τ_Ray(550)

(5)依据上述550nm处总的光学厚度，利用查找表获取反演参数ρ_path(λ)、T(λ)、S(λ)，通过辐射传输方程求解地表反射率，完成大气订正。

由步骤(4)确定的大气总的光学厚度τ，再依据确定的卫星成像几何参数，通过查找表获取大气订正方程中的参数ρ_path(λ)、T(λ)、S(λ)，然后由公式

求出地表反射率：

$\mathrm{\ρ}\left(\mathrm{\λ}\right)=\frac{1}{S\left(\mathrm{\λ}\right)+\frac{T\left(\mathrm{\λ}\right)}{{\mathrm{\ρ}}_{\mathrm{TOA}}\left(\mathrm{\λ}\right)-{\mathrm{\ρ}}_{\mathrm{Path}}\left(\mathrm{\λ}\right)}}$

上述查找表的建立为利用MODTRAN辐射传输模型，依据几何参数、大气模式、气溶胶模式等设置建立的；查找表的索引τ分别与ρ_path(λ)、T(λ)、S(λ)之间的关系满足下式：

ρ_Path(λ)＝c₀(λ)+c₁(λ)τ+c₂(λ)τ²

S(λ)＝d₀(λ)+d₁(λ)τ+d₂(λ)τ²

T(λ)＝e₀(λ)+e₁(λ)τ+e₂(λ)τ²

c₀(λ)、c₁(λ)、c₂(λ)、d₀(λ)、d₁(λ)、d₂(λ)、e₀(λ)、e₁(λ)、e₂(λ)对应多光谱波段λ在确定的卫星成像几何参数太阳天顶角θ_S、卫星观测天顶角θ_V、太阳与卫星的方位角

下的大气路径反射率ρ_path(λ)、大气透过率T(λ)、半球反射率S(λ)与总光学厚τ的系数。

(6)对步骤(5)中输出的大气订正产品进行精度检验

利用与其他同类卫星反演、其他方法反演和实地测量等多种方法进行对比验证，以确保大气订正的精度，如不符合精度，就从步骤(2)开始重新进行，符合精度输出大气订正产品。

附图3是植被的表观反射率、订正反射率与实测地表反射率曲线对比，从图比较可知，表观反射率虽然经过了绝对辐射校正，但它是地面反射率和大气反射率的总和，不能反映地物的特性，而经过本发明大气订正的反射率与实测地表反射率非常接近误差非常小，确保了订正地表反射率真实反映地物的辐射特性，这使大气订正后的图像更加适宜在地物识别、定量反演和模拟中应用。

本发明未详细说明部分属于本领域技术人员公知常识。

Claims (3)

1.一种遥感卫星多光谱数据的大气订正方法，其特征在于步骤如下：

(1)读取遥感卫星多光谱数据，对该多光谱数据进行预处理，将多光谱数据转换为表观辐亮度数据和表观反射率数据；

(2)依据表观反射率数据计算决策树法中的决策因子：比值植被指数RVI、土壤可调植被指数SAVI和归一化水体指数NDWI，通过决策树法实现DDV暗目标自动提取；

(3)根据步骤(1)中的表观反射率数据获取上述提取的DDV暗目标区域的DDV表观反射率数据；并根据该DDV表观反射率数据结合DDV地表反射率计算Path反射率；

(4)根据Path反射率和Rayleigh反射率比值与气溶胶光学厚度的关系，确定多光谱各波段的气溶胶光学厚度；并依据多光谱各波段的气溶胶光学厚度和Rayleigh散射光学厚度，确定出550nm处总的光学厚度τ；

(5)依据上述550nm处总的光学厚度，利用查找表获取反演参数ρ_path(λ)、T(λ)、S(λ)，通过辐射传输方程求解地表反射率，完成大气订正。

2.根据权利要求1所述的一种遥感卫星多光谱数据的大气订正方法，其特征在于所述步骤(3)中的DDV地表反射率的计算公式如下：

ρ_DDV(λ)＝S_vegρ_veg(λ)+(1-S_veg)ρ_soil(λ)

其中：ρ_DDV(λ)对应多光谱波段λ的DDV地表反射率；

ρ_veg(λ)、ρ_soil(λ)是实测多光谱波段λ的植被反射率和土壤反射率；

S_veg是植被调整因子，

3.根据权利要求1所述的一种遥感卫星多光谱数据的大气订正方法，其特征在于所述步骤(5)中查找表的索引τ分别与ρ_path(λ)、T(λ)、S(λ)之间的关系满足下式：

ρ_Path(λ)＝c₀(λ)+c₁(λ)τ+c₂(λ)τ²

S(λ)＝d₀(λ)+d₁(λ)τ+d₂(λ)τ²

T(λ)＝e₀(λ)+e₁(λ)τ+e₂(λ)τ²

其中，c₀(λ)、c₁(λ)、c₂(λ)、d₀(λ)、d₁(λ)、d₂(λ)、e₀(λ)、e₁(λ)、e₂(λ)对应多光谱波段λ在确定的卫星成像几何参数太阳天顶角θ_S、卫星观测天顶角θ_V、太阳与卫星的方位角

下的大气路径反射率ρ_path(λ)、大气透过率T(λ)、半球反射率S(λ)与总光学厚度τ的系数。

Images