CN102508226B - 一种实现浑浊水体区域modis影像大气校正的方法 - Google Patents

Abstract

一种实现浑浊水体区域MODIS影像大气校正的方法，其步骤如下：(1)读入531、551、667和678nm波段实测的标准离水反射率之间的经验关系、Rayleigh散射查找表、模型所涉及波段的MODIS影像、同步大气压强、风速和臭氧厚度数据；(2)在步骤(1)涉及的气象数据的基础上，利用Rayleigh散射查找表算法，计算MODIS各波段Rayleigh散射贡献的反射率，并加以订正；(3)利用Rayleigh散射修正后的MODIS影像数据，构建531、551、667和678nm波段的反射率之间的经验关系；(4)在步骤(4)的基础上，利用Angstrom气溶胶指数模型，计算MODIS各波段的气溶胶散射贡献的反射率；(5)在步骤(2)、(4)的基础上计算大气的程辐射贡献的反射率，并对MODIS影像进行大气校正。

Description

一种实现浑浊水体区域MODIS影像大气校正的方法

技术领域

本发明涉及一种实现浑浊水体区域MODIS(Moderate-Resolution ImagingSpectroradiometer)影像大气校正的方法，属于水色定量遥感理论与应用技术领域，适用于水色定量遥感的理论方法和应用技术研究。

背景技术

标准MODIS大气校正算法称之为“清洁”水体大气校正算法。该算法的基本假设在于研究区域内存在“清洁”水体(叶绿素a浓度＜0.5-1.0μg/l，悬浮浓度几乎0)，则MODIS近红外两个波段的离水反射率近似为零，但是由于大气散射作用，使得传感器探测到的信号不为零，因此可以将这部分信号归根为大气程辐射的贡献，并用于大气信息的提取。标准MODIS大气校正算法就是利用“清洁”水体的这种光学性质进行大气校正的。因此，标准MODIS大气校正算法成败的关键在于是否能够在研究区域内找到一种为“清洁”水体的像元。

众所周知，中国海岸带海域，尤其是黄河口、长江口和珠江口等三角洲地带的水体，黄河、长江和珠江携带大量的泥沙入海，以致入海口很大面积的水域含沙量非常之高。研究结果表明，在近红外波段，水体的离水反射率与悬浮泥沙浓度成正比，哪怕水体中仅含微量的悬浮泥沙，都将引起近红外波段离水辐亮度成倍的增长。这就意味着标准MODIS大气校正算法在这些区域将会失败。要想在浑浊的II类水体中实现大气精校正，则需要更多的先验知识参与大气校正过程。

发明内容

本发明的目的在于提供一种实现浑浊水体区域MODIS影像大气校正的方法，以解决现有技术中标准MODIS大气校正算法在浑浊水体区域无法适用的缺陷，在浑浊的II类水体中实现大气精校正。

本发明一种实现浑浊水体区域MODIS影像大气校正的方法，具体所采用的技术方案，包括如下步骤：

步骤(1)：读入531、551、667和678nm波段实测的标准离水反射率之间的经验关系、Rayleigh散射查找表、模型所涉及波段的MODIS影像、同步大气压强、风速和臭氧厚度数据。其中标准离水反射率之间的关系如下：

[ρ_w]_N(λ₂)≈a₁[ρ_w]_N(λ₁)

[ρ_w]_N(λ₄)≈a₂[ρ_w]_N(λ₃)

式中，x_m、x_n和y_m为经验常数，可以通过回归计算得到。

步骤(2)：在步骤(1)涉及的气象数据的基础上，利用Rayleigh散射查找表算法，计算MODIS各波段Rayleigh散射贡献的反射率，并加以订正。

ρ_t(λ)＝ρ_r(λ)+ρ_a(λ)+ρ_ra(λ)+t(λ)[ρ_w]_N(λ)

ρ_at(λ)＝ρ_path-ρ_r(λ)＝ρ_a(λ)+ρ_ra(λ)+t(λ)ρ_n(λ)

式中，ρ_r(λ)为Rayleigh散射贡献的反射率；ρ_a(λ)为气溶胶散射贡献的反射率；ρ_at(λ)为经过Rayleigh散射校正后的表观反射率；ρ_ra(λ)为气溶胶和大气分子多次散射贡献的反射率，该项可以解释光子在大气分子和气溶胶颗粒之间多次穿梭的结果；ρ_n(λ)为待求的表观反射率；t(λ)为漫射透过系数，其表达式如下：

$t\left(\mathrm{\λ}\right)=\exp \{-[0.5{\mathrm{\τ}}_r\left(\mathrm{\λ}\right)+{\mathrm{\τ}}_{\mathrm{oz}}\left(\mathrm{\λ}\right)\left](\frac{1}{\cos {\mathrm{\θ}}_{\mathrm{sat}}}+\frac{1}{\cos {\mathrm{\θ}}_{\mathrm{sol}}})\right\}$

式中，θ_sat为卫星天顶角；θ_sol为太阳高度角；τ_r(λ)为Rayleigh散射光学厚度。Rayleigh散射光学厚度与气压之间存在如下关系：

${\mathrm{\τ}}_r\left(\mathrm{\λ}\right)=\frac{P_0}{P}{\mathrm{\τ}}_{r,0}\left(\mathrm{\λ}\right)$

式中，τ_r，0(λ)为标准大气压下的Rayleigh散射的光学厚度，可从Rayleigh散射查找表中获取；P₀为标准大气压，其值取1013.25mb；P为像元所在位置的实际大气压强。臭氧光学厚度可以用下式表示：

${\mathrm{\τ}}_{\mathrm{oz}}\left(\mathrm{\λ}\right)={\mathrm{\τ}}_{\mathrm{oz},0}\left(\mathrm{\λ}\right)\frac{\mathrm{DU}}{1000}$

式中，DU为单位臭氧浓度。

步骤(3)：利用Rayleigh散射修正后的MODIS影像数据，构建531、551、667和678nm波段的反射率之间的经验关系。

Rayleigh散射修正后的MODIS影像在531、551、667和678nm波段的反射率之间的经验关系形式如下：

ρ_at(λ₂)＝x_mρ_at(λ₁)+y_m

ρ_at(λ₄)＝x_nρ_at(λ₃)

步骤(4)：在步骤(3)的基础上，利用531、551、667和678nm波段反射率的经验关系，结合Angstrom气溶胶指数模型，计算MODIS各波段的气溶胶散射和气溶胶-Rayleigh多次散射贡献的反射率。

将步骤(1)所示的经验关系，代入ρ_at(λ)的表达式，并化简可得：

${\mathrm{\ρ}}_{\mathrm{at}}\left({\mathrm{\λ}}_2\right)=a_1\frac{t\left({\mathrm{\λ}}_2\right)}{t\left({\mathrm{\λ}}_1\right)}{\mathrm{\ρ}}_{\mathrm{at}}\left({\mathrm{\λ}}_1\right)+{\mathrm{\ρ}}_{\mathrm{as}}\left({\mathrm{\λ}}_2\right)-a_1\frac{t\left({\mathrm{\λ}}_2\right)}{t\left({\mathrm{\λ}}_1\right)}{\mathrm{\ρ}}_{\mathrm{as}}\left({\mathrm{\λ}}_1\right)$

${\mathrm{\ρ}}_{\mathrm{at}}\left({\mathrm{\λ}}_4\right)=a_2\frac{t\left({\mathrm{\λ}}_4\right)}{t\left({\mathrm{\λ}}_3\right)}{\mathrm{\ρ}}_{\mathrm{at}}\left({\mathrm{\λ}}_3\right)+{\mathrm{\ρ}}_{\mathrm{as}}\left({\mathrm{\λ}}_4\right)-a_2\frac{t\left({\mathrm{\λ}}_4\right)}{t\left({\mathrm{\λ}}_3\right)}{\mathrm{\ρ}}_{\mathrm{as}}\left({\mathrm{\λ}}_3\right)$

Angstrom气溶胶指数模型ε(λ_s，λ₁)如下：

$\mathrm{\ϵ}(\mathrm{\λ},{\mathrm{\λ}}_1)\≈\frac{{\mathrm{\ρ}}_{\mathrm{as}}\left(\mathrm{\λ}\right)}{{\mathrm{\ρ}}_{\mathrm{as}}\left({\mathrm{\λ}}_1\right)}\≈\mathrm{Exp}\left[n({\mathrm{\λ}}_1-\mathrm{\λ})\right]$

式中，n为Angstrom指数；λ₁为参考波段。将ε(λ_s，λ₁)代入本步骤中的ρ_t(λ₂)和ρ_t(λ₄)的表达式可得：

$R_{\mathrm{at}}\left({\mathrm{\λ}}_2\right)=\frac{a_{1}t\left({\mathrm{\λ}}_2\right)}{t\left({\mathrm{\λ}}_1\right)}{R}_{\mathrm{at}}\left({\mathrm{\λ}}_1\right)+[\mathrm{\ϵ}({\mathrm{\λ}}_2,{\mathrm{\λ}}_1)-\frac{a_{1}t\left({\mathrm{\λ}}_2\right)}{t\left({\mathrm{\λ}}_1\right)}\mathrm{\ϵ}({\mathrm{\λ}}_1,{\mathrm{\λ}}_1)]R_{\mathrm{as}}\left({\mathrm{\λ}}_1\right)$

$R_{\mathrm{at}}\left({\mathrm{\λ}}_4\right)=\frac{a_{2}t\left({\mathrm{\λ}}_4\right)}{t\left({\mathrm{\λ}}_3\right)}{R}_{\mathrm{at}}\left({\mathrm{\λ}}_3\right)+[\mathrm{\ϵ}({\mathrm{\λ}}_4,{\mathrm{\λ}}_1)-\frac{a_{2}t\left({\mathrm{\λ}}_4\right)}{t\left({\mathrm{\λ}}_3\right)}\mathrm{\ϵ}({\mathrm{\λ}}_3,{\mathrm{\λ}}_1)]R_{\mathrm{as}}\left({\mathrm{\λ}}_1\right)$

式中，ρ_as(λ)为ρ_ra(λ)和ρ_a(λ)的和。将步骤(3)中的经验关系代入上式可得参考波段的ρ_as(λ)和n的计算方程组：

$n=\frac{\ln x_n}{{\mathrm{\λ}}_3-{\mathrm{\λ}}_4}$

${\mathrm{\ρ}}_{\mathrm{as}}\left({\mathrm{\λ}}_1\right)=\frac{y_m}{\mathrm{\ϵ}({\mathrm{\λ}}_2,{\mathrm{\λ}}_1)-x_m}$

步骤(5)：在步骤(2)、(4)的基础上计算大气的程辐射贡献的发射率，并对MODIS影像进行大气校正，计算公式如下：

${\left[{\mathrm{\ρ}}_w\right]}_N\left(\mathrm{\λ}\right)=\frac{{\mathrm{\ρ}}_{\mathrm{at}}\left(\mathrm{\λ}\right)-\mathrm{\ϵ}(\mathrm{\λ},{\mathrm{\λ}}_1){\mathrm{\ρ}}_{\mathrm{as}}\left({\mathrm{\λ}}_1\right)}{t\left(\mathrm{\λ}\right)}.$

本发明一种实现浑浊水体区域MODI S影像大气校正的方法，其优点及功效在于：本发明以波段之间的经验关系为基础，结合查找表Rayleigh散射计算方法，研究与探讨适用于混杂II水体的MODIS大气校正算法。该算法克服了浑浊I I类水体大气校正难的问题，有助于提高浑浊水体的大气校正精度，进而有利于推动水色遥感的发展，具有重要的理论和应用意义。

附图说明

图1为本发明的气溶胶散射贡献的反射率。

图2为本发明的各波段大气校正误差。

具体实施方式

为了更好的说明本发明涉及的一种实现浑浊水体区域MODIS影像大气校正的方法，利用ASD野外地物光谱仪分别测量得到黄河口、长江口和太湖的光谱数据。2003年10月27日和28日，在太湖进行表光光学量测量的同时，收集了MODIS影像数据各一景，用于算法分析计算和精度评估。本发明涉及一种实现浑浊水体区域MODIS影像大气校正的方法，具体实现步骤如下：

步骤(1)：读入531、551、667和678nm波段实测的标准离水反射率之间的经验关系、Rayleigh散射查找表、模型所涉及波段的MODIS影像、同步大气压强、风速和臭氧厚度数据。

臭氧厚度数据、风速和大气压强同步气象数据可以从NASA的官方网站ftp://oceans.gsfc.nasa.gov下载得到。531、551、667和678nm波段实测的标准离水反射率之间的经验关系如下：

[ρ_w]_N(555nm)≈0.9812[ρ_w]_N(53lnm)，R²＝0.9970

[ρ_w]_N(678nm)≈0.9812[ρ_w]_N(667nm)，R²＝0.9955

步骤(2)：在步骤(1)涉及的气象数据的基础上，利用Rayleigh散射查找表算法，计算MODIS各波段Rayleigh散射贡献的反射率，并加以订正。

步骤(3)：利用Rayleigh散射修正后的2003年10月27日和28日两景MODIS影像数据，构建531、551、667和678nm波段的反射率之间的经验关系。

2003年10月27日遥感影像的经验关系：

ρ_at(555nm)＝1.1681ρ_at(531nm)-0.0043，R²＝0.9795

ρ_at(678nm)＝0.9808ρ_at(667nm)，R²＝0.9934

2003年10月28日遥感影像的经验关系：

ρ_at(555nm)＝1.1106ρ_at(531nm)-0.0012，R²＝0.9936

ρ_at(678nm)＝0.9835ρ_at(667nm)，R²＝0.9880

步骤(4)：在步骤(3)的基础上，利用531、551、667和678nm波段反射率之间的经验关系，结合Angstrom气溶胶指数模型，计算MODIS各波段的气溶胶散射和气溶胶-Rayleigh多次散射贡献的反射率。

将步骤(3)所示的经验算法，即可求得R_as(531nm)和n的值。将上述的R_as(531nm)和n的值代入Angstrom气溶胶指数模型，便可得到412、443、488、531、551、667、678、745和870nm波段的气溶胶散射贡献的反射率，如图1所示。

步骤(5)：在步骤(2)、(4)的基础上计算大气的程辐射贡献的反射率，并对MODIS影像进行大气校正。

为了检验ACES大气校正算法的稳定性，本发明分别收集了2003年10月27日和28日的MODIS影像及其同步标准离水反射率数据。利用同步的标准离水反射率与经验算法的结果进行比对分析，并将该比对结果作为该经验算法的误差。图2给出了ACES算法的误差状况。由图2可知，ACES算法在412(11.96％)、443(13.50％)、488(14.94％)、531(3.56％)、555(4.19％)和667nm(18.99％)波段具有较好的遥感反射率反演精度，但是在678、748和870nm，该算法的误差较大(＞30.53％)。

Claims (1)

1.一种实现浑浊水体区域MODIS影像大气校正的方法，其步骤如下：

步骤（1）：读入531、551、667和678nm波段实测的标准离水反射率之间的经验关系、Rayleigh散射查找表、模型所涉及波段的MODIS影像、同步大气压强、风速和臭氧厚度数据；其中标准离水反射率之间的关系如下：

[ρ_w]_N(λ₂)≈a₁[ρ_w]_N(λ₁)

[ρ_w]_N(λ₄)≈a₂[ρ_w]_N(λ₃)

式中，[ρ_w]_N(λ₁)、[ρ_w]_N(λ₂)、[ρ_w]_N(λ₃)和[ρ_w]_N(λ₄)为标准反射率；a₁和a₂为经验系数，由实测离水反射率的回归分析得到；λ₁、λ₂、λ₃和λ₄为波长；λ₁为531nm，λ₂为551nm，λ₃为667nm，λ₄为678nm；

步骤（2）：在步骤（1）涉及的气象数据的基础上，利用Rayleigh散射查找表算法，计算MODIS各波段Rayleigh散射贡献的反射率，并加以订正；

ρ_t(λ)=ρ_r(λ)+ρ_a(λ)+ρ_ra(λ)+t(λ)[ρ_w]_N(λ)

ρ_at(λ)=ρ_path-ρ_r(λ)=ρ_a(λ)+ρ_ra(λ)+t(λ)ρ_n(λ)

式中，ρ_r(λ)为Rayleigh散射贡献的反射率；ρ_a(λ)为气溶胶散射贡献的反射率；ρ_at(λ)为经过Rayleigh散射校正后的表观反射率；ρ_ra(λ)为气溶胶和大气分子多次散射贡献的反射率，该项解释为光子在大气分子和气溶胶颗粒之间多次穿梭的结果；ρ_n(λ)为待求的表观反射率；t(λ)为漫射透过系数，其表达式如下：

$t\left(\mathrm{\λ}\right)=\exp \{-[0.5{\mathrm{\τ}}_r\left(\mathrm{\λ}\right)+{\mathrm{\τ}}_{\mathrm{oz}}\left(\mathrm{\λ}\right)\left](\frac{1}{\cos {\mathrm{\θ}}_{\mathrm{sat}}}+\frac{1}{\cos {\mathrm{\θ}}_{\mathrm{sol}}})\right\}$

式中，θ_sat为卫星天顶角；θ_sol为太阳高度角；τ_r(λ)为Rayleigh散射光学厚度；Rayleigh散射光学厚度与气压之间存在如下关系：

${\mathrm{\τ}}_r\left(\mathrm{\λ}\right)=\frac{P_0}{P}{\mathrm{\τ}}_{r,0}\left(\mathrm{\λ}\right)$

式中，τ_r，0(λ)为标准大气压下的Rayleigh散射的光学厚度，从Rayleigh散射查找表中获取；P₀为标准大气压，其值取1013.25mb；P为像元所在位置的实际大气压强；臭氧光学厚度用下式表示：

${\mathrm{\τ}}_{\mathrm{oz}}\left(\mathrm{\λ}\right)={\mathrm{\τ}}_{\mathrm{oz},0}\left(\mathrm{\λ}\right)\frac{\mathrm{DU}}{100}$

式中，DU为单位臭氧浓度；

步骤（3）：利用Rayleigh散射修正后的MODIS影像数据，构建531、551、667和678nm波段的反射率之间的经验关系；

Rayleigh散射修正后的MODIS影像在531、551、667和678nm波段的反射率之间的经验关系形式如下：

ρ_at(λ₂)=x_mρ_at(λ₁)+y_m

ρ_at(λ₄)=x_nρ_at(λ₃)

式中，x_m、x_n和y_m为经验常数，通过回归计算得到；

步骤（4）：在步骤(3)的基础上，利用531、551、667和678nm波段反射率的经验关系，结合Angstrom气溶胶指数模型，计算MODIS各波段的气溶胶散射和气溶胶Rayleigh多次散射贡献的反射率；

将步骤(1)所示的经验关系，代入ρ_at(λ)的表达式，并化简可得：

${\mathrm{\ρ}}_{\mathrm{at}}\left({\mathrm{\λ}}_2\right)=a_1\frac{t\left({\mathrm{\λ}}_2\right)}{t\left({\mathrm{\λ}}_1\right)}{\mathrm{\ρ}}_{\mathrm{at}}\left({\mathrm{\λ}}_1\right)+{\mathrm{\ρ}}_{\mathrm{as}}\left({\mathrm{\λ}}_2\right)-a_1\frac{t\left({\mathrm{\λ}}_2\right)}{t\left({\mathrm{\λ}}_1\right)}{\mathrm{\ρ}}_{\mathrm{as}}\left({\mathrm{\λ}}_1\right)$

${\mathrm{\ρ}}_{\mathrm{at}}\left({\mathrm{\λ}}_4\right)=a_2\frac{t\left({\mathrm{\λ}}_4\right)}{t\left({\mathrm{\λ}}_3\right)}{\mathrm{\ρ}}_{\mathrm{at}}\left({\mathrm{\λ}}_3\right)+{\mathrm{\ρ}}_{\mathrm{as}}\left({\mathrm{\λ}}_4\right)-a_2\frac{t\left({\mathrm{\λ}}_4\right)}{t\left({\mathrm{\λ}}_3\right)}{\mathrm{\ρ}}_{\mathrm{as}}\left({\mathrm{\λ}}_3\right)$

Angstrom气溶胶指数模型ε(λ,λ_l)如下：

$\mathrm{\ϵ}(\mathrm{\λ},{\mathrm{\λ}}_1)\≈\frac{{\mathrm{\ρ}}_{\mathrm{as}}\left(\mathrm{\λ}\right)}{{\mathrm{\ρ}}_{\mathrm{as}}\left({\mathrm{\λ}}_1\right)}\≈\mathrm{Exp}\left[n({\mathrm{\λ}}_1-\mathrm{\λ})\right]$

式中，n为Angstrom指数；将ε(λ,λ_l)代入本步骤中的ρ_t(λ₂)和ρ_t(λ₄)的表达式可得：

${\mathrm{\ρ}}_{\mathrm{at}}\left({\mathrm{\λ}}_2\right)=\frac{a_{1}t\left({\mathrm{\λ}}_2\right)}{t\left({\mathrm{\λ}}_1\right)}{\mathrm{\ρ}}_{\mathrm{at}}\left({\mathrm{\λ}}_1\right)+[\mathrm{\ϵ}({\mathrm{\λ}}_2,{\mathrm{\λ}}_1)-\frac{a_{1}t\left({\mathrm{\λ}}_2\right)}{t\left({\mathrm{\λ}}_1\right)}\mathrm{\ϵ}({\mathrm{\λ}}_1,{\mathrm{\λ}}_1)]{\mathrm{\ρ}}_{\mathrm{as}}\left({\mathrm{\λ}}_1\right)$

${\mathrm{\ρ}}_{\mathrm{at}}\left({\mathrm{\λ}}_4\right)=\frac{a_{2}t\left({\mathrm{\λ}}_4\right)}{t\left({\mathrm{\λ}}_3\right)}{\mathrm{\ρ}}_{\mathrm{at}}\left({\mathrm{\λ}}_3\right)+[\mathrm{\ϵ}({\mathrm{\λ}}_4,{\mathrm{\λ}}_1)-\frac{a_{2}t\left({\mathrm{\λ}}_4\right)}{t\left({\mathrm{\λ}}_3\right)}\mathrm{\ϵ}({\mathrm{\λ}}_3,{\mathrm{\λ}}_1)]{\mathrm{\ρ}}_{\mathrm{as}}\left({\mathrm{\λ}}_1\right)$

式中，ρ_as(λ)为ρ_ra(λ)和ρ_a(λ)的和；将步骤(3)中的经验关系代入上式可得参考波段的ρ_as(λ)和n的计算方程组：

$n=\frac{\ln x_n}{{\mathrm{\λ}}_3-{\mathrm{\λ}}_4}$

${\mathrm{\ρ}}_{\mathrm{as}}\left({\mathrm{\λ}}_1\right)=\frac{y_m}{\mathrm{\ϵ}({\mathrm{\λ}}_2,{\mathrm{\λ}}_1)-x_m}$

步骤(5)：在步骤(2)、(4)的基础上计算大气的程辐射贡献的发射率，并对MODIS影像进行大气校正，计算公式如下：

${\mathrm{\ρ}}_n\left(\mathrm{\λ}\right)=\frac{{\mathrm{\ρ}}_{\mathrm{at}}\left(\mathrm{\λ}\right)-\mathrm{\ϵ}(\mathrm{\λ},{\mathrm{\λ}}_1){\mathrm{\ρ}}_{\mathrm{as}}\left({\mathrm{\λ}}_1\right)}{t\left(\mathrm{\λ}\right)}.$

Images