同时反演内陆水体混浊度、cod和叶绿素浓度的方法
【技术领域】
[0001] 本发明属于卫星遥感技术领域，特别涉及一种同时反演内陆水体混浊度、COD和叶 绿素浓度的遥感方法。
【背景技术】
[0002] 应用水质遥感技术可对大区域的水环境进行监测，不受地面条件的限制。且卫星 遥感数据周期性覆盖地球表面，可以很方便地对区域水环境进行动态监测。遥感监测结果 可以很好反演水环境的空间分布情况，与常规监测手段相结合，可以大大提高水环境监测 的效能。
[0003] 目前采用的水质遥感方法主要为经验公式法和半分析模型法两大类：
[0004] 其中，经验公式法是在卫星过境时在水面上进行大量的水面采样分析，然后通过 回归分析等统计方法，建立混浊度、叶绿素等与某个敏感波段反射率的经验公式，再根据该 经验公式计算所有水域像元的相应水质参数。
[0005] 但是，采用经验公式的方法需建立在大量的卫星同步地面采样分析数据的基础之 上，需耗费大量的人力物力。且这种方法多是单参数反演，而水体反射率是多种光活性因子 共同影响的结果，同时还受多种环境因素的影响。故这种方法精度有限，且适应性差，不同 的地区，不同的季节，所建立的经验公式皆不适用。
[0006] 半分析模型法多采用由蒙特卡罗方法得到的半分析模型：
库中：RW 为水体以射率，β为水体后向散射系数，α为水体吸收体系数，f为常数。
[0007] 该公式假定水深很大，可以数个参数同时反演。目前能反演的参数主要为：悬浮 物、叶绿数和黄色物质。其它水质参数的遥感均是通过与该三个参数的共生关系进行估算， 具有很大的不确定性，离实用化还有相当的距离。
[0008] 但是，采用半分析模型的方法也存在几方面的问题。一是普遍采用的公式只适用 于海域的水质遥感，应用于水深有限的内陆地表水水质遥感时常误差较大；其二，常数f是 一个经验值，与成像时的环境条件有关，一般通过经验方式获取，其值与成像时间和地区气 象等条件有关。在一定程度上限制了模型的适应性。
[0009] 目前所能反演且能达到实用化精度要求的参数为悬浮物、叶绿数和黄色物质，而 对有机污染浓度等最普遍的污染参数仍未能反演。
【发明内容】
[0010] 本发明的目的是提供同时反演内陆水体混浊度、COD和叶绿素浓度的方法，该方法 可以有效地同时反演出内陆水体悬浮泥沙浓度、有机污染物浓度和叶绿素浓度，从而有效 反映区域水环境状况的空间分布，且可以反映大、中型污染源，为相关部门进行水环境污染 控制，水资源保护提供极为有用的信息。
[0011] 为了达到上述目的，本发明采用如下技术方案：
[0012] 同时反演内陆水体混浊度、COD和叶绿素浓度的方法，包括以下步骤：
[0013] (SI)建立内陆地表水水质遥感物理分析模型；
[0014] (S2)测量水体消光系数、水体吸收系数和水体反向散射系数，以及研究区典型水 下地物的反射率；
[0015] (S3)对卫星原始图像进行精确的大气纠正并消除水面镜面反射光的影响
[0016] (S4)简化所述内陆地表水水质遥感物理分析模型，并根据简化后的物理分析模型 计算出悬浮泥沙、有机污染物和叶绿素的浓度；
[0017] (S5)重复步骤（S4)，得到所有像元的悬浮泥沙、有机污染物和叶绿素的浓度。
[0018] 作为一种具体的实施例，所述步骤（SI)中内陆地表水水质遥感物理分析模型的 函数关系为：
[0019]
[0020] 式中：Ρ(Θ)为水体散射相函数，Θ为水体散射角；β为水体后向散射系数，α为 水体吸收系数，R b为水底底质反射率，h为水深；μ = sec Θ +sec φ，Θ为阳光在水中传播 的天顶角，Φ为水中的观测天顶角；RW= Rws+Rwb，lU%水体反射率，Rws为水体散射光对反射 率的贡献，R wb为水底反射光对反射率的贡献。
[0021] 作为一种具体的实施例，步骤（S3)具体为：
[0022] (S31)从图像上提取暗像元，迭代计算大气散射系数、大气透过率和天空光辐照 度。
[0023] (S32)对图像中所有像元进行大气纠正。
[0024] (S33)消除天空光镜面反射光。
[0025] (S34)消除太阳直射光镜面反射光。
[0026] 进一步地，所述步骤（S31)中迭代计算大气散射系数、大气透过率和天空光辐照 度的步骤包括：
[0027] (S311)读取暗像元的反射率Rvd，根据暗像元的表观反射率Rvd的计算公式：
[0028] 式中：ω为大气散射系数，T为大气透过率，&为地物反射率，Ρ( Θ )为大气散射 相函数，Θ为散射角， ,4R'.,
[0029] 首先忽略等式（1)右边第一项的值，计算得到大气散射系数的初值W
[0030] (S312):根据大气透过率的初值的计算公式：T' = e ，计算出大气透过率的初 值T'，将大气透过率的初值T'代入等式（1)右边第一项中的T，将大气散射系数的初值
ω'代入等式（1)右边第一项中的ω，计算得到大气散射系数
[0031] (S313):计算天空光幅照度：
[0032] 将步骤（S312)中计算得到的大气散射系数Co1代入天空光幅照度计算公式 匕.，计算天空光幅照度Eskyl，式中：c为大气纠正前后的大气散射系数比，为测量值；
[0033] (S314):将步骤（S312)中计算得到的大气散射系数Co1代入大气透过率的计算公 式T = e ω，计算出大气透过率T1;
[0034] (S315):将步骤（S312)中计算得到的大气散射系数％代入等式⑴右边第一项 中的ω、将步骤（S314)中计算得到的大气透过率1\代入等式（1)右边第一项中的Τ，计算 出精度高于大气散射系数大气散射系数ω 2;
[0035] 根据步骤（S313)计算出精度更高的天空光幅照度Esky2;
[0036] 根据步骤（S314)计算出精度高于大气透过率T1的大气透过率T 2;
[0037]依次办法，迭代计算出精度高于大气散射系数ω2的大气散射系数ω k、精度高于 大气透过率T2的大气透过率T k和精度高于天空光幅照度E _2的天空光幅照度E skyk;其中， k>2〇
[0038] 进一步地，所述步骤（S32)中对图像中所有像元进行大气纠正的计算公式为：
[0039]
[0040] 式中：Rg为大气纠正后得到的像元反射率，T' ' = e~°s0，其中τ为大气光学 厚度，R1为像元表观反射率
L为像元原始辐亮度，Ε。为太阳辐照度。
[0041] 进一步地，所述步骤（S33)的具体步骤为：将像元各波段的原始辐亮度减去步 骤（S315)得到的天空光辐照度E skyk与天空光水面镜面反射率Rni之积，其公式为：L 1 = L-EskykRm,
[0042] 式中：L1为像元消除天空光镜面反射后的辐亮度。
[0043] 进一步地，所述步骤（S34)的具体步骤为：将图像转化为反射率图像，像元所有波 段的反射率均减去红外波段的反射率，其中，像元所有的波段的反射率减去红外波段的反 射率的公式为：R 2= Rg-Ru
[0044] 式中=R2为像元消除光镜面反射后的反射率，R U为红外波段的反射率。
[0045] 作为一种具体的实施例，所述步骤（S4)中简化所述内陆地表水水质遥感物理分 析模型具体为：
[0046] 当水深较大时，内陆地表水水质