CN103499815A - 一种基于氧气和水汽吸收波段的内陆水体大气校正方法 - Google Patents

Abstract

一种基于氧气和水汽吸收波段进行内陆水体大气校正的方法，涉及遥感技术领域，所述方法包括：对内陆水体的高光谱遥感影像数据进行处理，提取清洁水体像元(暗像元)；对提取的暗像元进行去瑞利散射大气校正；对校正后的像元，利用氧气和水汽吸收波段组合计算大气校正因子；外推全波段的气溶胶散射；利用计算出的气溶胶散射值，对整景影像进行大气校正。所述方法，利用水体在氧气和水汽吸收波段反射率低的特点，选取水体暗像元，并计算气溶胶散射，该方法避免了内陆水体大气校正中，因水体浑浊度高而导致的过校正现象；提高了内陆浑浊水体大气校正的精度，为提高遥感反演精度提供了技术支持。

Description

一种基于氧气和水汽吸收波段的内陆水体大气校正方法

技术领域

本发明涉及遥感技术领域，特别涉及一种针对高光谱数据的、基于氧气和水汽吸收波段进行内陆水体大气校正的方法。

背景技术

在水色遥感中，传感器在大气层顶接收到的信号中有90％以上来自大气散射，而水体信息在可见光范围内最多只占到10％，大气散射信息对水体光谱信息造成了极大的干扰，严重影响了水环境参数的遥感提取。因此，要提高水环境遥感监测的精度，必须进行精确的大气校正。大气散射由大气分子的散射(瑞利散射)以及气溶胶散射组成，由于大气分子成分及含量比较稳定，大气分子散射贡献已能通过多次散射加偏振的精确Rayleigh散射计算得到，而气溶胶含量在空间域及时间域上变化较大，要准确计算其散射比较困难，故而不同大气校正算法的区别主要体现在如何计算气溶胶散射，以及大气分子与气溶胶的共同作用上。

传统的大气校正算法中气溶胶散射部分的计算多是建立在“黑暗像元”的假设条件上，即假设水体在近红外波段(NIR＞700nm)的离水辐亮度近似为0，近红外波段的信号值都来自于瑞利散射以及气溶胶散射，且气溶胶随着波长的变化表现出一定的规律性，在此基础上，将近红外波段的气溶胶散射外推到可见光波段范围内，从而得到各波段的离水辐亮度。然而，由于内陆水体浑浊度较高，导致水体在近红外波段的反射升高，使得近红外离水辐亮度近似为0的假设不再成立，因此该方法无法直接应用于内陆及沿岸等二类水体。对于内陆水体，重新寻找适合的光谱波段，以推算其气溶胶散射，是提高大气校正的精度的有效途径。

发明内容

本发明要解决的技术问题是：针对内陆水体，确定适宜的高光谱波段，克服通常大气校正中过校正的现象，提高高光谱遥感数据的大气校正精度。

为解决上述技术问题，本发明提供一种基于氧气和水汽吸收波段的大气校正方法，其包括步骤：

A：对内陆水体的高光谱影像数据进行预处理；

B：对所述经过预处理的影像数据，提取清洁像元(暗像元)；

C：针对所述清洁像元，计算大气分子瑞利散射，进而计算得到去除瑞利散射后的辐射亮度值；

D：针对所述去除瑞利散射后的像元，利用氧气与水汽吸收波段组合计算大气校正因子；

E：根据所述计算出的大气校正因子，推算各波段气溶胶散射；

F：根据所述计算的气溶胶散射，完成高光谱影像数据的大气校正。

优选地，所述步骤A具体包括步骤：

A1：获取所需进行大气校正的高光谱遥感影像，从中提取内陆水体水域的影像数据；

A2：利用所获取数据的定标参数(通常由传感器发射方提供)，对原始影像数据进行辐射定标和几何纠正。

优选地，所述步骤B中，通过对整景影像像元值的对比，选取出像元值最小的样点(不少于20个)，并记为清洁像元(暗像元)。

优选地，所述步骤B中，考虑到气溶胶空间分布的非均匀性，暗像元应分布在影像的各个区域。

优选地，所述步骤C中，针对所述步骤B中提取的暗像元，计算大气分子瑞利散射，进而计算去瑞利散射的辐射亮度值，具体计算方法如下：

首先，计算瑞利散射L_r，用式(1)计算：

式(1)中各参数及计算方法如下：

F₀(λ)是大气层外对应不同波长λ的太阳辐照度，其大小随着日地距离而变化。

μ₀和μ分别代表太阳天顶角以及卫星天顶角的余弦，α_↓和α_↑分别表示入射光和反射光的散射相位角。入射光的散射相函数P_r(α_↓)定义为：

$P_r\left({\mathrm{\α}}_{\↓}\right)=\frac{3}{4}(1+{\cos }^2\left({\mathrm{\α}}_{\↓}\right))---\left(2\right)$

其中，

反射光的瑞利散射相函数P_r(α_↑)定义为：

$P_r\left({\mathrm{\α}}_{\↑}\right)=\frac{3}{4}[1+{({2\mathrm{\μ\μ}}_0+\cos \left({\mathrm{\α}}_{\↑}\right))}^2]---\left(4\right)$

其中，

ρ(μ₀)、ρ(μ)为菲涅尔反射率，计算公式如下：

$\begin{array}{c}\mathrm{\ρ}\left(x\right)=1-2\mathrm{xyn}[{(x+\mathrm{ny})}^{-2}+{(\mathrm{nx}+y)}^{-2}]\\ y=\frac{\sqrt{n^2+x^2-1}}{n}\end{array}---\left(6\right)$

式中，x为μ₀或μ，n为纯水的折射率，一般取值为1.34。

经过臭氧校正的大气透过率T_OZ(λ)计算公式如下：

T_OZ(λ)=exp[-τ_OZ(λ)*(μ^-1+μ₀ ^-1)]       (7)

其中，臭氧光学厚度计算式如下：

τ_OZ(λ)=a_OZ*U_OZ              (8)

式中，a_OZ为单位臭氧吸收系数，单位是cm^-1，U_OZ为大气臭氧含量，单位为cm。

瑞利光学厚度τ_r(λ)采用下式计算：

${\mathrm{\τ}}_r\left(\mathrm{\λ}\right)=0.008569{\mathrm{\λ}}^{-4}(1+0.0113{\mathrm{\λ}}^{-2}+0.00013{\mathrm{\λ}}^{-4})*\frac{P}{P_0}---\left(9\right)$

式中，P为实测大气压强，单位为百帕(hPa)，而P₀为标准大气压强，其值为1013.25hPa。

根据以上计算得到的中间过程参数，根据式(1)完成瑞利散射计算。

其次，大气顶层辐亮度减去瑞利散射辐亮度，即可得到去除瑞利散射后的辐射亮度值。

优选地，所述步骤D中，针对所述去除瑞利散射后的像元，利用氧气与水汽吸收波段组合计算大气校正因子ε，具体计算方法如下：

$\mathrm{\ϵ}({\mathrm{\λ}}_i,{\mathrm{\λ}}_j)=\frac{{\mathrm{\ρ}}_{\mathrm{as}}\left({\mathrm{\λ}}_i\right)}{{\mathrm{\ρ}}_{\mathrm{as}}\left({\mathrm{\λ}}_j\right)}=\frac{{\mathrm{\ω}}_a\left({\mathrm{\λ}}_i\right){\mathrm{\τ}}_a\left({\mathrm{\λ}}_i\right)p_a(\mathrm{\θ},{\mathrm{\θ}}_0,{\mathrm{\λ}}_i)}{{\mathrm{\ω}}_a\left({\mathrm{\λ}}_j\right){\mathrm{\τ}}_a\left({\mathrm{\λ}}_j\right)p_a(\mathrm{\θ},{\mathrm{\θ}}_0,{\mathrm{\λ}}_j)}---\left(10\right)$

其中，λ_j位于氧气吸收带，λ_j位于水汽吸收带。对于同一种气溶胶模型来说，ω_a可以认为是定值，而单次散射相函数P_a是观测几何与波长λ的函数，对大量实测数据的统计研究发现，单次散射相函数之比也可以看做常数，而在Junge谱下，气溶胶厚度是关于波长的负指数函数，即

τ_a(λ)=σλ^-β          (11)

其中，σ为大气浑浊指数，β为埃斯特朗系数，在一定的气溶胶模型下其值为常数。

优选地，所述步骤E中，根据所述步骤D中计算出的大气校正因子，推算各波段气溶胶散射：

$\mathrm{\ϵ}({\mathrm{\λ}}_n,{\mathrm{\λ}}_{\mathrm{NIR}})=\frac{{\mathrm{\ρ}}_A\left({\mathrm{\λ}}_n\right)}{{\mathrm{\ρ}}_A\left({\mathrm{\λ}}_{\mathrm{NIR}}\right)}=\exp \left[C({\mathrm{\λ}}_{\mathrm{NIR}}-{\mathrm{\λ}}_n)\right]---\left(12\right)$

其中，λ_n为待估波长，λ_NIR为参考波长，ρ_A为气溶胶多次散射反射率，由此得到各波段大气校正参数ε(λ_n，λ_NIR)，再根据外推公式(12)得到各波段的气溶胶散射。

优选地，所述步骤F中，将所述E中计算的气溶胶散射，应用于高光谱遥感影像的逐像元，完成高光谱影像数据的大气校正。

本发明的大气校正方法，利用水体在氧气和水汽吸收波段反射率非常低的特点，计算大气气溶胶散射，该方法避免了内陆水体大气校正中，因水体浑浊度高而导致的过校正现象；提高了内陆浑浊水体大气校正的精度，为提高遥感反演精度提供了技术支持。

附图说明

图1是本发明实施例所述的基于氧气和水汽吸收波段进行高光谱影像大气校正的方法流程图；

图2是2007年11月11日、2008年11月20日、2009年4月25日3景遥感影像大气校正的效果对比图，图中利用了与3景影像过境时间相差不超过2小时的6个地面样点数据进行对比。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

本实施例以太湖作为内陆水体的示例，以MERIS数据作为高光谱遥感影像的示例，详细说明利用所述方法进行内陆水体大气校正的过程。图1是本发明实施例所述的基于氧气和水汽吸收波段进行高光谱影像大气校正的方法流程图，如图1所述，所述方法包括步骤：

A：通过网站获取了2007年11月11日、2008年11月20日、2009年4月25日三天的MERIS影像数据；

B：对所述经过预处理的影像数据，提取清洁像元(暗像元)；

C：针对所述清洁像元，计算大气分子瑞利散射，进而计算得到去除瑞利散射后的辐射亮度值；

D：针对所述去除瑞利散射后的像元，利用氧气与水汽吸收波段组合计算大气校正因子；

E：根据所述计算出的大气校正因子，推算各波段气溶胶散射；

F：根据所述计算的气溶胶散射，完成高光谱影像数据的大气校正。

2、如权利要求1所述的方法，其特征在于，所述步骤A中，需对所述高光谱影像数据进行辐射定标、几何校正的预处理。

A：对所获取的ME RIS数据进行预处理，包括：太湖水域提取、几何校正、辐射定标；

B：对A中处理得到的影像，进行逐像元数据对比，提取清洁像元(暗像元)，本示例中共选取了30个清洁像元点，分别分布于太湖北部、东部和南部。

C：针对步骤B中提取的暗像元，计算大气分子瑞利散射，进而计算去瑞利散射的辐射亮度值，具体计算方法如下：

首先，计算瑞利散射L_r，用式(1)计算：

式(1)中各参数及计算方法如下：

F₀(λ)是大气层外对应不同波长λ的太阳辐照度，其大小随着日地距离而变化。

μ₀和μ分别代表太阳天顶角以及卫星天顶角的余弦，α_↓和α_↑分别表示入射光和反射光的散射相位角。入射光的散射相函数P_r(α_↓)定义为：

$P_r\left({\mathrm{\α}}_{\↓}\right)=\frac{3}{4}(1+{\cos }^2\left({\mathrm{\α}}_{\↓}\right))---\left(2\right)$

其中，

反射光的瑞利散射相函数P_r(α_↑)定义为：

$P_r\left({\mathrm{\α}}_{\↑}\right)=\frac{3}{4}[1+{({2\mathrm{\μ\μ}}_0+\cos \left({\mathrm{\α}}_{\↑}\right))}^2]---\left(4\right)$

其中，

ρ(μ₀)、ρ(μ)为菲涅尔反射率，计算公式如下：

$\begin{array}{c}\mathrm{\ρ}\left(x\right)=1-2\mathrm{xyn}[{(x+\mathrm{ny})}^{-2}+{(\mathrm{nx}+y)}^{-2}]\\ y=\frac{\sqrt{n^2+x^2-1}}{n}\end{array}---\left(6\right)$

式中，x为μ₀或μ，n为纯水的折射率，一般取值为1.34。

经过臭氧校正的大气透过率T_OZ(λ)计算公式如下：

T_OZ(λ)=exp[-τ_OZ(λ)*(μ^-1+μ₀ ^-1)]       (7)

其中，臭氧光学厚度计算式如下：

τ_OZ(λ)=a_OZ*U_OZ           (8)

式中，a_OZ为单位臭氧吸收系数，单位是cm^-1，U_OZ为大气臭氧含量，单位为cm。

瑞利光学厚度τ_r(λ)采用下式计算：

${\mathrm{\τ}}_r\left(\mathrm{\λ}\right)=0.008569{\mathrm{\λ}}^{-4}(1+0.0113{\mathrm{\λ}}^{-2}+0.00013{\mathrm{\λ}}^{-4})*\frac{P}{P_0}---\left(9\right)$

式中，P为实测大气压强，单位为百帕(hPa)，而P₀为标准大气压强，其值为1013.25hPa。

根据以上计算得到的中间过程参数，根据式(1)完成瑞利散射计算。

其次，大气顶层辐亮度减去瑞利散射辐亮度，即可得到去除瑞利散射后的辐射亮度值。

D：针对所述去除瑞利散射后的像元，利用MERIS的761nm和900nm波段(该两个波段分别位于氧气和水汽吸收波段)组合计算大气校正因子ε，具体计算方法如下：

$\mathrm{\ϵ}({\mathrm{\λ}}_{761},{\mathrm{\λ}}_{900})=\frac{{\mathrm{\ρ}}_{\mathrm{as}}\left({\mathrm{\λ}}_{761}\right)}{{\mathrm{\ρ}}_{\mathrm{as}}\left({\mathrm{\λ}}_{900}\right)}=\frac{{\mathrm{\ω}}_a\left({\mathrm{\λ}}_{761}\right){\mathrm{\τ}}_a\left({\mathrm{\λ}}_{761}\right)p_a(\mathrm{\θ},{\mathrm{\θ}}_0,{\mathrm{\λ}}_{761})}{{\mathrm{\ω}}_a\left({\mathrm{\λ}}_{900}\right){\mathrm{\τ}}_a\left({\mathrm{\λ}}_{900}\right)p_a(\mathrm{\θ},{\mathrm{\θ}}_0,{\mathrm{\λ}}_{900})}---\left(10\right)$

其中，λ₇₆₁位于氧气吸收带，λ₉₀₀位于水汽吸收带。对于同一种气溶胶模型来说，ω_a可以认为是定值，而单次散射相函数P_a是观测几何与波长λ的函数，对大量实测数据的统计研究发现，单次散射相函数之比也可以看做常数，而在Junge谱下，气溶胶厚度是关于波长的负指数函数，即

τ_a(λ)=σλ^-β              (11)

其中，σ为大气浑浊指数，β为埃斯特朗系数，在一定的气溶胶模型下其值为常数。

E：根据步骤D中计算出的大气校正因子，推算各波段气溶胶散射：

$\mathrm{\ϵ}({\mathrm{\λ}}_n,{\mathrm{\λ}}_{900})=\frac{{\mathrm{\ρ}}_A\left({\mathrm{\λ}}_n\right)}{{\mathrm{\ρ}}_A\left({\mathrm{\λ}}_{900}\right)}=\exp \left[C({\mathrm{\λ}}_{900}-{\mathrm{\λ}}_n)\right]---\left(12\right)$

其中，λ_n为待估波长，λ₉₀₀为900nm波长，ρ_A为气溶胶多次散射反射率，由此得到各波段大气校正参数ε(λ_n，λ₉₀₀)，再根据外推公式(12)得到各波段的气溶胶散射。

F：将步骤E中计算的气溶胶散射，应用于预处理后的MERIS影像数据的逐像元，完成整景影像的大气校正。

本发明实施例所述基于氧气与水汽吸收波段的大气校正方法，利用MERIS全分辨率的卫星影像数据，结合其位于761nm和900nm的氧气和水汽吸收波段，估算暗像元的大气气溶胶散射，进而对整景影像进行大气校正，得到水面遥感反射率数据。将该方法应用于2007年11月11日、2008年11月20日、2009年4月25日的MERIS卫星影像，与地面准同步实测样点相比，其数据拟合结果如图2所示，结果表明，在MERIS39波段的校正精度较高，其相对误差均小于20％，符合内陆水色遥感的精度要求。该发明为提高内陆水体高光谱数据的大气校正精度，提供了一种技术方法，此方法的推广，将促进遥感技术在水环境监测中的应用。

以上实施方式仅用于说明本发明，而并非对本发明的限制，有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型，因此所有等同的技术方案也属于本发明的范畴，本发明的专利保护范围应由权利要求限定。

Claims (7)

1.一种基于氧气和水汽吸收波段的内陆水体大气校正方法，其特征在于，包括步骤：

A：对内陆水体的高光谱影像数据进行预处理；

B：对所述经过预处理的影像数据，提取清洁像元即暗像元；

C：针对所述清洁像元，计算大气分子瑞利散射，进而计算得到去除瑞利散射后的辐射亮度值；

D：针对所述去除瑞利散射后的像元，利用氧气与水汽吸收波段组合计算大气校正因子；

E：根据所述计算出的大气校正因子，推算各波段气溶胶散射；

F：根据所述计算的气溶胶散射，完成高光谱影像数据的大气校正。

2.如权利要求1所述的方法，其特征在于，所述步骤A中，需对所述高光谱影像数据进行辐射定标、几何校正的预处理。

3.如权利要求1所述的方法，其特征在于，所述步骤B中，需通过对整景影像的数据对比，提取清洁像元即暗像元，考虑到气溶胶空间分布的非均匀性，暗像元应分布在影像的各个区域。

4.如权利要求1所述的方法，其特征在于，所述步骤C中，针对所述步骤B中提取的暗像元，计算大气分子瑞利散射，进而计算去瑞利散射的辐射亮度值，具体计算方法如下：

首先，计算瑞利散射L_r，用式(1)计算：

式(1)中各参数及计算方法如下：

F₀(λ)是大气层外对应不同波长λ的太阳辐照度，其大小随着日地距离而变化，

μ₀和μ分别代表太阳天顶角以及卫星天顶角的余弦，α_↓和α_↑分别表示入射光和反射光的散射相位角，入射光的散射相函数P_r(α_↓)定义为：

$P_r\left({\mathrm{\α}}_{\↓}\right)=\frac{3}{4}(1+{\cos }^2\left({\mathrm{\α}}_{\↓}\right))---\left(2\right)$

其中，

反射光的瑞利散射相函数P_r(α_↑)定义为：

$P_r\left({\mathrm{\α}}_{\↑}\right)=\frac{3}{4}[1+{({2\mathrm{\μ\μ}}_0+\cos \left({\mathrm{\α}}_{\↑}\right))}^2]---\left(4\right)$

其中，

ρ(μ₀)、ρ(μ)为菲涅尔反射率，计算公式如下：

$\begin{array}{c}\mathrm{\ρ}\left(x\right)=1-2\mathrm{xyn}[{(x+\mathrm{ny})}^{-2}+{(\mathrm{nx}+y)}^{-2}]\\ y=\frac{\sqrt{n^2+x^2-1}}{n}\end{array}---\left(6\right)$

式中，x为μ₀或μ，n为纯水的折射率，一般取值为1.34，

经过臭氧校正的大气透过率T_OZ(λ)计算公式如下：

T_OZ(λ)=exp[-τ_OZ(λ)*(μ^-1+μ₀ ^-1)]       (7)

其中，臭氧光学厚度计算式如下：

τ_OZ(λ)=a_OZ*U_OZ                (8)

式中，a_OZ为单位臭氧吸收系数，单位是cm^-1，U_OZ为大气臭氧含量，单位为cm，

瑞利光学厚度τ_r(λ)采用下式计算：

${\mathrm{\τ}}_r\left(\mathrm{\λ}\right)=0.008569{\mathrm{\λ}}^{-4}(1+0.0113{\mathrm{\λ}}^{-2}+0.00013{\mathrm{\λ}}^{-4})*\frac{P}{P_0}---\left(9\right)$

式中，P为实测大气压强，单位为百帕(hPa)，而P₀为标准大气压强，其值为1013.25hPa。

根据以上计算得到的中间过程参数，根据式(1)完成瑞利散射计算，

其次，大气顶层辐亮度减去瑞利散射辐亮度，即可得到去除瑞利散射后的辐射亮度值。

5.如权利要求1所述的方法，其特征在于，所述步骤D中，针对所述去除瑞利散射后的像元，利用氧气与水汽吸收波段组合计算大气校正因子ε，具体计算方法如下：

$\mathrm{\ϵ}({\mathrm{\λ}}_i,{\mathrm{\λ}}_j)=\frac{{\mathrm{\ρ}}_{\mathrm{as}}\left({\mathrm{\λ}}_i\right)}{{\mathrm{\ρ}}_{\mathrm{as}}\left({\mathrm{\λ}}_j\right)}=\frac{{\mathrm{\ω}}_a\left({\mathrm{\λ}}_i\right){\mathrm{\τ}}_a\left({\mathrm{\λ}}_i\right)p_a(\mathrm{\θ},{\mathrm{\θ}}_0,{\mathrm{\λ}}_i)}{{\mathrm{\ω}}_a\left({\mathrm{\λ}}_j\right){\mathrm{\τ}}_a\left({\mathrm{\λ}}_j\right)p_a(\mathrm{\θ},{\mathrm{\θ}}_0,{\mathrm{\λ}}_j)}---\left(10\right)$

其中，λ_j位于氧气吸收带，λ_j位于水汽吸收带。对于同一种气溶胶模型来说，ω_a可以认为是定值，而单次散射相函数P_a是观测几何与波长λ的函数，对大量实测数据的统计研究发现，单次散射相函数之比也可以看做常数，而在Junge谱下，气溶胶厚度是关于波长的负指数函数，即

τ_a(λ)=σλ^-β                (11)

其中，σ为大气浑浊指数，β为埃斯特朗系数，在一定的气溶胶模型下其值为常数。

6.如权利要求1所述的方法，其特征在于，所述步骤E中，根据所述D中计算出的大气校正因子，推算各波段气溶胶散射：

$\mathrm{\ϵ}({\mathrm{\λ}}_n,{\mathrm{\λ}}_{\mathrm{NIR}})=\frac{{\mathrm{\ρ}}_A\left({\mathrm{\λ}}_n\right)}{{\mathrm{\ρ}}_A\left({\mathrm{\λ}}_{\mathrm{NIR}}\right)}=\exp \left[C({\mathrm{\λ}}_{\mathrm{NIR}}-{\mathrm{\λ}}_n)\right]---\left(12\right)$

其中，λ_n为待估波长，λ_NIR为参考波长，ρ_A为气溶胶多次散射反射率，由此得到各波段大气校正参数ε(λ_n，λ_NIR)，再根据外推公式(12)得到各波段的气溶胶散射。

7.如权利要求1所述的方法，其特征在于，所述步骤F中，将所述E中计算的气溶胶散射，应用于高光谱遥感影像的逐像元，完成高光谱影像数据的大气校正。

Images