CN105092476A - 同时反演内陆水体混浊度、cod和叶绿素浓度的方法 - Google Patents

Abstract

本发明公开了同时反演内陆水体混浊度、COD和叶绿素浓度的方法，包括以下步骤：(S1)建立内陆地表水水质遥感物理分析模型；(S2)测量水体消光系数、水体吸收系数和水体反向散射系数，以及研究区典型地物的反射率；(S3)对卫星原始图像进行精确的大气纠正并消除水面镜面反射光的影响；(S4)简化所述内陆地表水水质遥感物理分析模型，并根据简化后的物理分析模型计算出悬浮泥沙、有机污染物和叶绿素浓度；(S5)重复步骤(S4)，得到所有像元的悬浮泥沙、有机污染物和叶绿素浓度。该方法可以有效地同时反演出内陆水体悬浮泥沙浓度、有机污染物浓度和叶绿素浓度，从而有效反映区域水环境状况的空间分布。

Description

同时反演内陆水体混浊度、COD和叶绿素浓度的方法

技术领域

本发明属于卫星遥感技术领域，特别涉及一种同时反演内陆水体混浊度、COD和叶绿素浓度的遥感方法。

背景技术

应用水质遥感技术可对大区域的水环境进行监测，不受地面条件的限制。且卫星遥感数据周期性覆盖地球表面，可以很方便地对区域水环境进行动态监测。遥感监测结果可以很好反演水环境的空间分布情况，与常规监测手段相结合，可以大大提高水环境监测的效能。

目前采用的水质遥感方法主要为经验公式法和半分析模型法两大类：

其中，经验公式法是在卫星过境时在水面上进行大量的水面采样分析，然后通过回归分析等统计方法，建立混浊度、叶绿素等与某个敏感波段反射率的经验公式，再根据该经验公式计算所有水域像元的相应水质参数。

但是，采用经验公式的方法需建立在大量的卫星同步地面采样分析数据的基础之上，需耗费大量的人力物力。且这种方法多是单参数反演，而水体反射率是多种光活性因子共同影响的结果，同时还受多种环境因素的影响。故这种方法精度有限，且适应性差，不同的地区，不同的季节，所建立的经验公式皆不适用。

半分析模型法多采用由蒙特卡罗方法得到的半分析模型：其中：R_w为水体以射率，β为水体后向散射系数，α为水体昅收体系数，f为常数。

该公式假定水深很大，可以数个参数同时反演。目前能反演的参数主要为：悬浮物、叶绿数和黄色物质。其它水质参数的遥感均是通过与该三个参数的共生关系进行估算，具有很大的不确定性，离实用化还有相当的距离。

但是，采用半分析模型的方法也存在几方面的问题。一是普遍采用的公式只适用于海域的水质遥感，应用于水深有限的内陆地表水水质遥感时常误差较大；其二，常数f是一个经验值，与成像时的环境条件有关，一般通过经验方式获取，其值与成像时间和地区气象等条件有关。在一定程度上限制了模型的适应性。

目前所能反演且能达到实用化精度要求的参数为悬浮物、叶绿数和黄色物质，而对有机污染浓度等最普遍的污染参数仍未能反演。

发明内容

本发明的目的是提供同时反演内陆水体混浊度、COD和叶绿素浓度的方法，该方法可以有效地同时反演出内陆水体悬浮泥沙浓度、有机污染物浓度和叶绿素浓度，从而有效反映区域水环境状况的空间分布，且可以反映大、中型污染源，为相关部门进行水环境污染控制，水资源保护提供极为有用的信息。

为了达到上述目的，本发明采用如下技术方案：

同时反演内陆水体混浊度、COD和叶绿素浓度的方法，包括以下步骤：

(S1)建立内陆地表水水质遥感物理分析模型；

(S2)测量水体消光系数、水体吸收系数和水体反向散射系数，以及研究区典型水下地物的反射率；

(S3)对卫星原始图像进行精确的大气纠正并消除水面镜面反射光的影响

(S4)简化所述内陆地表水水质遥感物理分析模型，并根据简化后的物理分析模型计算出悬浮泥沙、有机污染物和叶绿素的浓度；

(S5)重复步骤(S4)，得到所有像元的悬浮泥沙、有机污染物和叶绿素的浓度。

作为一种具体的实施例，所述步骤(S1)中内陆地表水水质遥感物理分析模型的函数关系为：

$R_w=\frac{\mathrm{\β}P\left(\mathrm{\Θ}\right)}{4\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})}(1-e^{-\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})h})+R_b{e}^{-\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})h};$

式中：P(Θ)为水体散射相函数，Θ为水体散射角；β为水体后向散射系数，α为水体吸收系数，R_b为水底底质反射率，h为水深；μ＝secθ+secφ，θ为阳光在水中传播的天顶角，φ为水中的观测天顶角；R_w＝R_ws+R_wb，R_w为水体反射率，R_ws为水体散射光对反射率的贡献，R_wb为水底反射光对反射率的贡献。

作为一种具体的实施例，步骤(S3)具体为：

(S31)从图像上提取暗像元，迭代计算大气散射系数、大气透过率和天空光辐照度。

(S32)对图像中所有像元进行大气纠正。

(S33)消除天空光镜面反射光。

(S34)消除太阳直射光镜面反射光。

进一步地，所述步骤(S31)中迭代计算大气散射系数、大气透过率和天空光辐照度的步骤包括：

(S311)读取暗像元的反射率R_vd，根据暗像元的表观反射率R_vd的计算公式：

$R_{vd}=\frac{{\mathrm{\ωTR}}_v}{2}+\frac{\mathrm{\ω}}{4}P\left(\mathrm{\θ}\right)---\left(1\right)$

式中：ω为大气散射系数，T为大气透过率，R_v为地物反射率，P(θ)为大气散射相函数，θ为散射角，

首先忽略等式(1)右边第一项的值，计算得到大气散射系数的初值

(S312)：根据大气透过率的初值的计算公式：T′＝e^-ω′，计算出大气透过率的初值T′，将大气透过率的初值T′代入等式(1)右边第一项中的T，将大气散射系数的初值ω′代入等式(1)右边第一项中的ω，计算得到大气散射系数

(S313)：计算天空光幅照度：

将步骤(S312)中计算得到的大气散射系数ω₁代入天空光幅照度计算公式计算天空光幅照度E_sky1，式中：c为大气纠正前后的大气散射系数比，为测量值；

(S314)：将步骤(S312)中计算得到的大气散射系数ω₁代入大气透过率的计算公式T＝e^-ω，计算出大气透过率T₁；

(S315)：将步骤(S312)中计算得到的大气散射系数ω₁代入等式(1)右边第一项中的ω、将步骤(S314)中计算得到的大气透过率T₁代入等式(1)右边第一项中的T，计算出精度高于大气散射系数ω₁的大气散射系数ω₂；

根据步骤(S313)计算出精度更高的天空光幅照度E_sky2；

根据步骤(S314)计算出精度高于大气透过率T₁的大气透过率T₂；

依次办法，迭代计算出精度高于大气散射系数ω₂的大气散射系数ω_k、精度高于大气透过率T₂的大气透过率T_k和精度高于天空光幅照度E_sky2的天空光幅照度E_skyk；其中，k>2。

进一步地，所述步骤(S32)中对图像中所有像元进行大气纠正的计算公式为：

$R_g=\frac{R_1-\frac{{\mathrm{\ω}}_k}{4}P\left(\mathrm{\θ}\right)}{\[T^{\′\′}+\frac{{\mathrm{\ω}}_k}{2}\]T_k};$

式中：R_g为大气纠正后得到的像元反射率，T′′＝e^-τ/cosθ,其中τ为大气光学厚度，R₁为像元表观反射率：L为像元原始辐亮度，E₀为太阳辐照度。

进一步地，所述步骤(S33)的具体步骤为：将像元各波段的原始辐亮度减去步骤(S315)得到的天空光辐照度E_skyk与天空光水面镜面反射率R_m之积，其公式为：L₁＝L-E_skykR_m，

式中：L₁为像元消除天空光镜面反射后的辐亮度。

进一步地，所述步骤(S34)的具体步骤为：将图像转化为反射率图像，像元所有波段的反射率均减去红外波段的反射率，其中，像元所有的波段的反射率减去红外波段的反射率的公式为：R₂＝R_g-R_ir

式中：R₂为像元消除光镜面反射后的反射率，R_ir为红外波段的反射率。

作为一种具体的实施例，所述步骤(S4)中简化所述内陆地表水水质遥感物理分析模型具体为：

当水深较大时，内陆地表水水质遥感的物理分析模型简化为：

$R_w=\frac{{\mathrm{\β}}_w+D_s{\mathrm{\β}}_s+D_u{\mathrm{\β}}_u+D_c{\mathrm{\β}}_c}{4\mathrm{\μ}\[{\mathrm{\α}}_w+{\mathrm{\β}}_i+D_s({\mathrm{\α}}_s+{\mathrm{\β}}_s)+D_u({\mathrm{\α}}_u+{\mathrm{\β}}_u)+D_c({\mathrm{\α}}_c+{\mathrm{\β}}_c)\]}---\left(2\right);$

当水深较浅时，内陆地表水水质遥感的物理分析模型简化为：

ln(R_b)-ln(R_w)＝μh[α_w+β+D_s(α_s+β_s)+D_u(α_u+β_u)+D_c(α_c+β_c)]——(3)。

进一步地，步骤(S4)中根据简化后的物理分析模型计算出悬浮泥沙、有机污染物和叶绿素浓度具体为：

当水深较大时：

将三个波段的遥感数据代入(2)式，建立方程组：

$\underset{i=1}{\overset{n}{\Σ}}{R}_{wi}=\frac{{\mathrm{\β}}_{wi}+D_s{\mathrm{\β}}_{si}+D_u{\mathrm{\β}}_{ui}+D_c{\mathrm{\β}}_{ci}}{4\mathrm{\μ}\[{\mathrm{\α}}_{wi}+{\mathrm{\β}}_{wi}+D_s({\mathrm{\α}}_{si}+{\mathrm{\β}}_{si})+D_u({\mathrm{\α}}_{ui}+{\mathrm{\β}}_{ui})+D_c({\mathrm{\α}}_{ci}+{\mathrm{\β}}_{ci})\]}---\left(4\right);$

解该方程组(4)，得出悬浮泥沙浓度D_s、有机污染物浓度D_u和叶绿素浓度D_c；

当水深较浅时：

(A)：读取一像元，将四个波段的遥感数据代入(3)式，,建立方程组：

$\underset{i=1}{\overset{n}{\Σ}}ln\left(R_{bi}\right)-ln\left(R_{wi}\right)=\mathrm{\μ}h\[{\mathrm{\α}}_{wi}+{\mathrm{\β}}_{wi}+D_s({\mathrm{\α}}_{si}+{\mathrm{\β}}_{si})+D_u({\mathrm{\α}}_{ui}+{\mathrm{\β}}_{ui})+D_c({\mathrm{\α}}_{ci}+{\mathrm{\β}}_{ci})\]-$ $--\left(5\right);$

(B)：将一典型地物的水下地物反射率代入(5)式，用数值逼近法求解出悬浮泥沙浓度D_s、有机污染物浓度D_u、叶绿素浓度D_c以及水深h；

(C)：将步骤(B)中求得的悬浮泥沙浓度D_s、有机污染物浓度D_u和叶绿素浓度D_c以及水深h代入(5)式，反计算水体反射率R′_w；

(D)：对比计算水体反射率R′_w与原始图像水体反射率R_w的误差ε；

(E)：代入所有典型地物的水下地物反射率，重复计算步骤(B)-(D)，得到误差ε_i，其中i＝0、1、2···n；

误差最小的典型地物的水下地物反射率所对应的典型地物为该像元的水下地物类型，根据该典型地物的水下地物反射率计算出的悬浮泥沙浓度D_s、有机污染物浓度D_u、叶绿素浓度D_c以及水深h为该像元的悬浮泥沙浓度、有机污染物浓度、叶绿素浓度和水深。

与现有技术相比，本发明的有益效果在于：

本发明的同时反演内陆水体混浊度、COD和叶绿素浓度的方法，首先通过建立的水质遥感模型具有清晰的物理意义，且考虑了水底反射的影响，适合于内陆水体的水质遥感，实现了重要有机污染参数COD的遥感反演，大大提高的水质遥感的环境意义；得到的COD遥感结果可以直观反映区域水环境状况的空间分布，且可以反映大、中型污染源。为相关部门进行水环境污染控制，水资源保护提供极为有用的信息。

通过对珠江三角洲有机污染分布遥感反演和东江中游河段有机污染遥感反演，很好地验证了本发明方法在进行内陆水体的水质遥感反演的过程中具有较高的精度。

附图说明

图1是本发明同时反演内陆水体混浊度、COD和叶绿素浓度的方法流程图。

图2是本发明的水深与水下地物遥感模型的辐射传输原理示意图。

图3是本发明有机污染物的吸收系数的测量曲线图。

图4是是2001年11月20日珠江三角洲有机污染(COD)遥感分布图。

图5是2002年9月4日计算COD与测量值对比。

图6是东江潼湖—惠州—博罗河段COD遥感分布图。

图7是东江流域第二大水库白盆珠水库的COD遥感分布图。

具体实施方式

下面将结合附图以及具体实施方法来详细说明本发明，在本发明的示意性实施及说明用来解释本发明，但并不作为对本发明的限定。

如图1所示，本发明同时反演内陆水体混浊度、COD和叶绿素浓度的方法，包括以下步骤：

(A)：建立内陆地表水水质遥感物理分析模型。

如图2所示，为本发明的水深与水下地物遥感模型的辐射传输原理示意图。卫星传感器接收到的水面光能量由三部分组成，其一为水面镜面反射光，其二为水体散射光，其三为水底反射光。

其中，水面镜面反射光没有水深和水下地物的信息，为干扰信息，通过消除水域遥感数据的水面镜面反射影响的方法，去镜面反射光予以消除。消除水面镜面反射后，水体反射率R_w可由下式表示：

R_w＝R_ws+R_wb——(1)；

式中：R_ws为水体散射光对反射率的贡献，R_wb为水底反射光对反射率的贡献。

其中，水体散射光对反射率的贡献R_ws的计算公式为：

$R_{ws}=\frac{\mathrm{\β}P\left(\mathrm{\Θ}\right)}{4\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})}(1-e^{-\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})h})---\left(2\right);$

其中：β＝β_w+D_sβ_s+D_uβ_u+D_cβ_c，为水体后向散射系数，为水分子后向散射系数、悬浮泥沙后向散射系数、有机污染物后向散射系数、和叶绿素后向散射系数之和；β_w为水分子后向散射系数，β_s为悬浮泥沙后向散射系数，β_u为有机污染物后向散射系数，β_c为叶绿素后向散射系数，D_s为悬浮泥沙浓度，D_u为有机污染物浓度，D_c为叶绿素浓度。

α＝α_w+D_sα_s+D_uα_u+D_cα_c，为水体吸收系数，水分子吸收系数、悬浮泥沙吸收系数、有机污染物吸收系数、和叶绿素吸收系数之和；α_w为水分子吸收系数、α_s为悬浮泥沙吸收系数、α_u为有机污染物吸收系数、α_c为叶绿素吸收系数；

其中，P(Θ)为水体散射相函数，Θ为水体散射角；R_b为水底底质反射率，h为水深；μ＝secθ+secφ；θ为阳光在水中传播的天顶角，φ为水中的观测天顶角。

水底反射光对反射率的贡献R_wb为：

R_wb＝R_be^-μ(α+β)h——(3)；

将(2)式、(3)式代入(1)式中，得水面遥感反射率与水质成分浓度、水深、水下地物反射率的函数关系为：

$R_w=\frac{\mathrm{\β}P\left(\mathrm{\Θ}\right)}{4\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})}(1-e^{-\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})h})+R_b{e}^{-\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})h}---\left(4\right);$

(4)式即为内陆地表水水质遥感的物理分析模型。

(B)：精确测量水体消光系数、水体吸收系数和水体反向散射系数，以及研究区水下典型地物的反射率。所述水下典型地物包括有沙、石头、淤泥等。

其中，有机污染物吸收系数的测得结果如图3所示。

(C)：对卫星原始图像进行精确的大气纠正并消除水面镜面反射光的影响。

(D)：内陆地表水水质遥感物理分析模型的简化，并根据简化后的物理分析模型计算求解出悬浮泥沙、有机污染物和叶绿素的浓度。

深水模型

当水深较大时，e^-μ(α+β)h→0，另考虑水后向散射接近朗伯体，内陆地表水水质遥感的物理分析模型，即(4)式简化为：

$R_w=\frac{{\mathrm{\β}}_w+D_s{\mathrm{\β}}_s+D_u{\mathrm{\β}}_u+D_c{\mathrm{\β}}_c}{4\mathrm{\μ}\[{\mathrm{\α}}_w+{\mathrm{\β}}_i+D_s({\mathrm{\α}}_s+{\mathrm{\β}}_s)+D_u({\mathrm{\α}}_u+{\mathrm{\β}}_u)+D_c({\mathrm{\α}}_c+{\mathrm{\β}}_c)\]}---\left(5\right);$

将三个波段的遥感数据带入(5)式，建立三个方程组成的方程组如下：

$\underset{i=1}{\overset{n}{\Σ}}{R}_{wi}=\frac{{\mathrm{\β}}_{wi}+D_s{\mathrm{\β}}_{si}+D_u{\mathrm{\β}}_{ui}+D_c{\mathrm{\β}}_{ci}}{4\mathrm{\μ}\[{\mathrm{\α}}_{wi}+{\mathrm{\β}}_{wi}+D_s({\mathrm{\α}}_{si}+{\mathrm{\β}}_{si})+D_u({\mathrm{\α}}_{ui}+{\mathrm{\β}}_{ui})+D_c({\mathrm{\α}}_{ci}+{\mathrm{\β}}_{ci})\]}---\left(6\right);$

(6)式即为水面遥感反射率与悬浮泥沙浓度、有机污染物浓度、叶绿素浓度、水深和水下地物反射率的物理分析模型。

解该方程组，即可计算出悬浮泥沙浓度D_s、有机污染物浓度D_u和叶绿素浓度D_c。

浅水模型

当水较浅时，水体反射率主要为水低反射光，水体散射光很弱，即：1-e^-μ(α+β)h→0，(3)式简化为：R_w＝R_be^-μ(α+β)h；

即：ln(R_b)-ln(R_w)＝μh[α_w+β+D_s(α_s+β_s)+D_u(α_u+β_u)+D_c(α_c+β_c)]——————(7)；

将四个波段的遥感数据带入可建立四个方程组成的方程组；

$\underset{i=1}{\overset{n}{\Σ}}ln\left(R_{bi}\right)-ln\left(R_{wi}\right)=\mathrm{\μ}h\[{\mathrm{\α}}_{wi}+{\mathrm{\β}}_{wi}+D_s({\mathrm{\α}}_{si}+{\mathrm{\β}}_{si})+D_u({\mathrm{\α}}_{ui}+{\mathrm{\β}}_{ui})+D_c({\mathrm{\α}}_{ci}+{\mathrm{\β}}_{ci})\]-$ $--\left(8\right);$

将某种类型水下地物反射率代入(8)式，用数值逼近方法可求悬浮泥沙浓度D_s、有机污染物浓度D_u和叶绿素浓度D_c，和水深h。

然后将求得的结果代入(8)式，计算水体反射率R′_w，并与原始图像水体反射率R_w对比求其误差。

将所有典型水下地物类型均代入进行相同计算，取具有最小误差的水下地物类型作为该像元的水下地物类型，其算出的悬浮泥沙浓度D_s、有机污染物浓度D_u和叶绿素浓度D_c和水深h为该点的相应值。

对所有像元均进行上述运算，得到全区水域的悬浮泥沙、有机污染物和叶绿素浓度分布图。

其中，步骤(C)对卫星原始图像进行精确的大气纠正并消除水面镜面反射光的影响，包括以下步骤：

步骤(C1)：从图像上提取暗像元，迭代计算大气散射系数、大气透过率和天空光辐照度。

首先，从原始卫星图像上提取暗像元。根据提取的像元，计算获取大气散射系数、大气透过率和天空光辐照度。具体步骤如下：

(C11)读取暗像元的反射率R_vd，计算大气散射系数初值ω′：根据暗像元的反射率R_vd的计算公式： $R_{vd}=\frac{{\mathrm{\ωTR}}_v}{2}+\frac{\mathrm{\ω}}{4}P\left(\mathrm{\θ}\right)---\left(1\right),$

首先忽略等式(1)右边第一项的值，计算得到大气散射系数的初值

式中：ω为大气散射系数，T为大气透过率，R_v为地物反射率，P(θ)为大气散射相函数，θ为散射角。

(C12)：根据大气透过率的初值的计算公式：T′＝e^-ω′，计算出大气透过率的初值T′，并将大气透过率的初值T′代入等式(1)右边第一项中的T，将大气散射系数的初值ω′代入等式(1)右边第一项中的ω，计算得到大气散射系数

(C13)：计算天空光幅照度：

将步骤(C12)中计算得到的大气散射系数ω₁代入天空光幅照度计算公式计算天空光幅照度E_sky1，式中：c为大气纠正前后的大气散射系数比，通过测量得到，为小于2的测量值。

(C14)：将步骤(C12)中计算得到的大气散射系数ω₁代入大气透过率的计算公式T＝e^-ω，计算出大气透过率T₁。

(C15)：将步骤(C12)中计算得到的大气散射系数ω₁代入等式(1)右边第一项中的ω、将步骤(C14)中计算得到的大气透过率T₁代入等式(1)右边第一项中的T，计算出精度高于大气散射系数ω₁的大气散射系数ω₂。

根据步骤(C13)计算出精度更高的天空光幅照度E_sky2。

根据步骤(C14)计算出精度高于大气透过率T₁的大气透过率T₂。

依次办法，迭代计算出精度高于大气散射系数ω₂的大气散射系数ω_k、精度高于大气透过率T₂的大气透过率T_k和精度高于天空光幅照度E_sky2的天空光幅照度E_skyk；其中，k>2。

步骤(C2)：对图像中所有像元进行大气纠正。对图像中所有像元，均进行以下计算： $R_g=\frac{R_1-\frac{{\mathrm{\ω}}_k}{4}P\left(\mathrm{\θ}\right)}{\[T^{\′\′}+\frac{{\mathrm{\ω}}_k}{2}\]T_k},$

式中：R_g为大气纠正后得到的像元反射率，T′＝e^-τ/cosθ，其中τ为大气光学厚度，R₁为像元表观反射率：L为像元原始辐亮度，E₀为太阳辐照度。

经过大气纠正后，大气散射光的影响被很好消除，影像清晰度得到明显提高。但水域镜面反射光的影响严重，水质及水下地物信息受到严重干扰。

本发明根据水面镜面反射的特征，消除每个像元的水面镜反射，其中，水面镜反射由水面对太阳直射光的镜面反射和对天空光的镜面反射光组成。由于太阳直射光和天空光的光谱组成不同，故对二者的镜面反射光分别消除。

步骤(C3)消除天空光镜面反射光。

通过测量或暗像元提取等方法获取天空光强度，然后像元各波段辐亮度减去天空光辐照度与镜面反射率之积即可消除天空光镜面反射的影响。

将像元各波段的原始辐亮度减去步骤(A5)得到的天空光辐照度E_skyk与天空光水面镜面反射率R_m之积，其公式为：L₁＝L-E_skykR_m。

式中：L₁为像元消除天空光镜面反射后的辐亮度。

(C4)消除太阳直射光镜面反射光。

根据大量的测量实验证明，各光学波段的水面镜面反射光的反射率近相等。对同一像元，各波段的镜面反射率近相等；长波近红外波段，由于水强吸收，从水中出来的光强为0，其像元反射率即为镜面反射光的反射率。像元所有波段的反射率均减去红外波段的反射率，即可消除太阳直射光镜面反射的影响。

将图像转化为反射率图像，像元所有波段的反射率均减去红外波段的反射率，其中，波段的反射率减去红外波段的反射率的公式为：R₂＝R-R_ir。

式中：R像元原始反射率，R₂像元消除光镜面反射后的反射率，R_ir红外波段原始反射率。

根据大气纠正后得到的像元发射率R_g，则得到R₂＝R_g-R_ir。以此消除太阳直射光镜面反射的影响。

消除镜面反射后，图像水域部分的清晰度得到显著改善，镜面反射的影响得到有效消除，水下及水质信息得到清楚显示。效果良好。

采用原始数据中阴影和清深水体等暗像元计算天空光强度的方法，无需实地测量，可通过计算机自动完成水面镜面反射光的消除，更有实用价值。用叠代方法计算可得到较精确的天空光辐照度，从而更精确地消除天空光镜面反射的影响。

通过对该方法进行大量的应用试验，包括对珠江三角洲有机污染分布遥感反演，东江中游河段有机污染遥感，结果表明具有较高的精度。以下是利用该方法的具体应用实例及部分结果。

应用实例1：珠江三角洲有机污染分布遥感

采用的遥感数据为2001年11月20日美国Landsat卫星的TM数据，空间分辨率为30米。经辐射定标和大气纠正后得到反射率图像。进行水陆分离运算后应用本发明提取出COD、悬浮泥沙和叶绿素浓度图，然后采用卫星同步水面监测数据对遥感结果进行精度验证。

图4为珠江三角洲河网区COD遥感结果图。由该图可见，污染最严重的是广州市珠江河段；西江水上游水质较好，进入珠江三角洲河网后污染程度迅速加重；东江水上游也基本没有污染，但进入东莞市区东江三角洲河网后污染程度也快增加，与实际情况完全相符。全区水体污染的空间分布情况在该图上一目了然。

采用采用卫星同步水面监测数据对遥感结果进行精度验证，结果如图5。遥感结果与实际测量结果的相关系数为0.77。经定标后，遥感结果精度可达76％。精度达到实用化水平。

应用实例2：东江中游河段有机污染遥感

采用的遥感数据为2009年11月24日的国产环境卫星数据，空间分辨率也为30米。经辐射定标和大气纠正后得到反射率图像。进行水陆分离运算后应用本发明提取出COD、悬浮泥沙和叶绿素浓度图，然后采用卫星同步水面监测数据对遥感结果进行精度验证。

图6为东江潼湖—惠州—博罗河段COD遥感分布图，图上可清楚看到，东江经惠州后，江水受到明显污染。且系列污染源均得到清楚反映。包括污染源位置，污染源排污染量及其影响范围均得到较显楚显示。

图7为东江流域第二大水库白盆珠水库的COD遥感分布图。可清楚见出，水库东部库湾A处明显存在一污染源。根据实际考查，A点为二十余个大型养鸦场。二十余万只鸦组成的鸦群的排泄物成为污染源。

以上对本发明实施例所提供的技术方案进行了详细介绍，本文中应用了具体个例对本发明实施例的原理以及实施方式进行了阐述，以上实施例的说明只适用于帮助理解本发明实施例的原理；同时，对于本领域的一般技术人员，依据本发明实施例，在具体实施方式以及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

Claims (9)

1.同时反演内陆水体混浊度、COD和叶绿素浓度的方法，其特征在于，包括以下步骤：

(S1)建立内陆地表水水质遥感物理分析模型；

(S2)测量水体消光系数、水体吸收系数和水体反向散射系数，以及研究区典型水下地物的反射率；

(S3)对卫星原始图像进行精确的大气纠正并消除水面镜面反射光的影响；

(S4)简化所述内陆地表水水质遥感物理分析模型，并根据简化后的物理分析模型计算出悬浮泥沙、有机污染物和叶绿素的浓度；

(S5)重复步骤(S4)，得到所有像元的悬浮泥沙、有机污染物和叶绿素的浓度。

2.根据权利要求1所述的同时反演内陆水体混浊度、COD和叶绿素浓度的方法，其特征在于，所述步骤(S1)中内陆地表水水质遥感物理分析模型的函数关系为：

$R_w=\frac{\mathrm{\β}P\left(\mathrm{\Θ}\right)}{4\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})}(1-e^{-\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})h})+R_b{e}^{-\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})h};$

式中：P(Θ)为水体散射相函数，Θ为水体散射角_；β为水体后向散射系数，α为水体吸收系数，R_b为水底底质反射率，h为水深；μ＝secθ+secφ，θ为阳光在水中传播的天顶角，φ为水中的观测天顶角；R_w＝R_ws+R_wb，R_w为水体反射率，R_ws为水体散射光对反射率的贡献，R_wb为水底反射光对反射率的贡献。

3.根据权利要求1所述的同时反演内陆水体混浊度、COD和叶绿素浓度的方法，其特征在于，步骤(S3)具体为：

(S31)从图像上提取暗像元，迭代计算大气散射系数、大气透过率和天空光辐照度；

(S32)对图像中所有像元进行大气纠正；

(S33)消除天空光镜面反射光；

(S34)消除太阳直射光镜面反射光。

4.根据权利要求3所述的同时反演内陆水体混浊度、COD和叶绿素浓度的方法，其特征在于，所述步骤(S31)中迭代计算大气散射系数、大气透过率和天空光辐照度的步骤包括：

(S311)读取暗像元的反射率R_vd，根据暗像元的表观反射率R_vd的计算公式：

$R_{vd}=\frac{{\mathrm{\ωTR}}_v}{2}+\frac{\mathrm{\ω}}{4}P\left(\mathrm{\θ}\right)---\left(1\right)$

式中：ω为大气散射系数，T为大气透过率，R_v为地物反射率，P(θ)为大气散射相函数，θ为散射角，

首先忽略等式(1)右边第一项的值，计算得到大气散射系数的初值

(S312)：根据大气透过率的初值的计算公式：T′＝e^-ω′，计算出大气透过率的初值T′，将大气透过率的初值T′代入等式(1)右边第一项中的T，将大气散射系数的初值ω′代入等式(1)右边第一项中的ω，计算得到大气散射系数

(S313)：计算天空光幅照度：

将步骤(S312)中计算得到的大气散射系数ω₁代入天空光幅照度计算公式计算天空光幅照度E_sky1，式中：c为大气纠正前后的大气散射系数比，为测量值；

(S314)：将步骤(S312)中计算得到的大气散射系数ω₁代入大气透过率的计算公式T＝e^-ω，计算出大气透过率T₁；

(S315)：将步骤(S312)中计算得到的大气散射系数ω₁代入等式(1)右边第一项中的ω、将步骤(S314)中计算得到的大气透过率T₁代入等式(1)右边第一项中的T，计算出精度高于大气散射系数ω₁的大气散射系数ω₂；

根据步骤(S313)计算出精度更高的天空光幅照度E_sky2；

根据步骤(S314)计算出精度高于大气透过率T₁的大气透过率T₂；

依次办法，迭代计算出精度高于大气散射系数ω₂的大气散射系数ω_k、精度高于大气透过率T₂的大气透过率T_k和精度高于天空光幅照度E_sky2的天空光幅照度E_skyk；其中，k>2。

5.根据权利要求4所述的同时反演内陆水体混浊度、COD和叶绿素浓度的方法，其特征在于，所述步骤(S32)中对图像中所有像元进行大气纠正的计算公式为：

$R_g=\frac{R_1-\frac{{\mathrm{\ω}}_k}{4}P\left(\mathrm{\θ}\right)}{\[T^{\′\′}+\frac{{\mathrm{\ω}}_k}{2}\]T_k}$

式中：R_g为大气纠正后得到的像元反射率，T″＝e^-τ/cosθ,其中τ为大气光学厚度，R₁为像元表观反射率：L为像元原始辐亮度，E₀为太阳辐照度。

6.根据权利要求5所述的同时反演内陆水体混浊度、COD和叶绿素浓度的方法，其特征在于，所述步骤(S33)的具体步骤为：将像元各波段的原始辐亮度减去步骤(S315)得到的天空光辐照度E_skyk与天空光水面镜面反射率R_m之积，其公式为：L₁＝L-E_skykR_m，

式中：L₁为像元消除天空光镜面反射后的辐亮度。

7.根据权利要求6所述的同时反演内陆水体混浊度、COD和叶绿素浓度的方法，其特征在于，所述步骤(S34)的具体步骤为：将图像转化为反射率图像，像元所有波段的反射率均减去红外波段的反射率，其中，像元所有的波段的反射率减去红外波段的反射率的公式为：R₂＝R_g-R_ir

式中：R₂为像元消除光镜面反射后的反射率，R_ir为红外波段的反射率。

8.根据权利要求1所述的同时反演内陆水体混浊度、COD和叶绿素浓度的方法，其特征在于，所述步骤(S4)中简化所述内陆地表水水质遥感物理分析模型具体为：

当水深较大时，内陆地表水水质遥感的物理分析模型简化为：

$R_w=\frac{{\mathrm{\β}}_w+D_s{\mathrm{\β}}_s+D_u{\mathrm{\β}}_u+D_c{\mathrm{\β}}_c}{4\mathrm{\μ}\[{\mathrm{\α}}_w+{\mathrm{\β}}_i+D_s({\mathrm{\α}}_s+{\mathrm{\β}}_s)+D_u({\mathrm{\α}}_u+{\mathrm{\β}}_u)+D_c({\mathrm{\α}}_c+{\mathrm{\β}}_c)\]}---\left(2\right);$

当水深较浅时，内陆地表水水质遥感的物理分析模型简化为：

ln(R_b)-ln(R_w)＝μh[α_w+β+D_s(α_s+β_s)+D_u(α_u+β_u)+D_c(α_c+β_c)]——(3)。

9.根据权利要求8所述的同时反演内陆水体混浊度、COD和叶绿素浓度的方法，其特征在于，步骤(S4)中根据简化后的物理分析模型计算出悬浮泥沙、有机污染物和叶绿素浓度具体为：

当水深较大时：

将三个波段的遥感数据代入(2)式，建立方程组：

$\underset{i=1}{\overset{n}{\Σ}}{R}_{wi}=\frac{{\mathrm{\β}}_{wi}+D_s{\mathrm{\β}}_{si}+D_u{\mathrm{\β}}_{ui}+D_c{\mathrm{\β}}_{ci}}{4\mathrm{\μ}\[{\mathrm{\α}}_{wi}+{\mathrm{\β}}_{wi}+D_s({\mathrm{\α}}_{si}+{\mathrm{\β}}_{si})+D_u({\mathrm{\α}}_{ui}+{\mathrm{\β}}_{ui})+D_c({\mathrm{\α}}_{ci}+{\mathrm{\β}}_{ci})\]}---\left(4\right);$

解该方程组(4)，得出悬浮泥沙浓度D_s、有机污染物浓度D_u和叶绿素浓度D_c；

当水深较浅时：

(A)：读取一像元，将四个波段的遥感数据代入(3)式，,建立方程组：

$\underset{i=1}{\overset{n}{\Σ}}ln\left(R_{bi}\right)-ln\left(R_{wi}\right)=\mathrm{\μ}h\[{\mathrm{\α}}_{wi}+{\mathrm{\β}}_{wi}+D_s({\mathrm{\α}}_{si}+{\mathrm{\β}}_{si})+D_u({\mathrm{\α}}_{ui}+{\mathrm{\β}}_{ui})+D_c({\mathrm{\α}}_{ci}+{\mathrm{\β}}_{ci})\]---\left(5\right);$

(B)：将一典型地物的水下地物反射率代入(5)式，用数值逼近法求解出悬浮泥沙浓度D_s、有机污染物浓度D_u、叶绿素浓度D_c以及水深h；

(C)：将步骤(B)中求得的悬浮泥沙浓度D_s、有机污染物浓度D_u和叶绿素浓度D_c以及水深h代入(5)式，反计算水体反射率R′_w；

(D)：对比计算水体反射率R′_w与原始图像水体反射率R_w的误差ε_；

(E)：代入所有典型地物的水下地物反射率，重复计算步骤(B)-(D)，得到误差ε_i，其中i＝0、1、2···n；

误差最小的典型地物的水下地物反射率所对应的典型地物为该像元的水下地物类型，根据该典型地物的水下地物反射率计算出的悬浮泥沙浓度D_s、有机污染物浓度D_u、叶绿素浓度D_c以及水深h为该像元的悬浮泥沙浓度、有机污染物浓度、叶绿素浓度和水深。