基于辐射传输过程的水质遥感模型及其在珠江口的应用

刘 冉，徐 鹏，刘心怡，黄 瑞，朱孟兰

(中国电建集团贵阳勘测设计研究院有限公司，贵州 贵阳 550081)

1 引言

当前，大部分内陆水体受到不同程度的污染，特别是在经济较发达的地区水体污染更为严重，这不仅影响了居民的居住环境，而且制约了城市进一步发展。因此，有效的水质监测和治理成为了亟待解决的问题，其中监测是基础。传统的监测方法是定期实地取水样进行分析，但是获取的只是离散的数据，时效性也不高。遥感技术自20世纪70年代开始应用于水体，其具有监测范围广、速度快、 成本低和便于进行长期动态监测的优点，为实时高效监测水质提供了一种行之有效的方法，促进了水污染的治理和生态环境的保护，对城市、社会的可持续发展具有重要的意义。

2 研究现状

随着卫星遥感技术的不断发展，水质遥感反演由海洋向内陆水体、由定性向定量发展，可反演的水质参数主要包括叶绿素a浓度、悬浮物浓度、CODM和浊度等，反演的模型主要有经验模型、半经验模型和物理模型[1]。

2.1 经验模型

经验模型主要是基于卫星多光谱波段或波段组合和水质参数之间的相关性分析，然后选择与水质参数相关性最大的波段或波段组合建立经验模型，广泛应用于内陆水质遥感监测。如：Steven M.Kloiber等[2]根据TM影像和地面同步观测数据，利用TM1和TM3的比率建立一个三参回归模型来预测SDT(Secchi disk transparency)。余丰宁等[3]利用TM各波段数据和太湖同步观测数据进行相关性分析，选择最佳波段建立了叶绿素含量的遥感定量模型。

经验模型虽然在特定条件、区域能取得比较理想的结果，但必须以大量实验为基础，而且物理学概念、机理不清，模型适应性差。

2.2 半经验模型

半经验模型[4]是针对高光谱数据源的水质遥感监测方法，对实测水体光谱或者卫星高光谱进行相关处理，选取和水体水质参数相关性最好的特征波段或变量，然后和水质参数建立反演模型，该模型具有一定的物理意义。

Koponen等[5]通过对航空高光谱影像和MERIS影像的多个波段进行组合或比值后和叶绿素a浓度建立各种模型，根据反演结果和测点数据计算RMSE，选出叶绿素a浓度最佳的反演模型；疏小舟等[6]在太湖地区采用高分辨率地物光谱仪实地测量湖水的反射光谱曲线，并同时采样分析叶绿素等水质参数，研究发现内陆水体叶绿素a浓度>5 μg/L时，光谱比R705nm/R675nm、R700nm附近反射峰位置适合用于叶绿素a浓度的遥感监测。

与经验模型相比，半经验模型具有一定的物理意义，模型精度有所提高。

2.3 分析模型

分析模型是基于太阳光在水体中的传输过程以及水体中的叶绿素、悬浮物、污染物和水分子等组分的吸收、散射特性对传输过程的影响，构建水体反射光谱与污染浓度之间函数关系的物理模型。

Hoogenboom等[7]利用水体各组分的相互作用，构建一种快速的、可行的反演水体叶绿素浓度的矩阵模型；李四海等[8]利用SeaWIFS数据，通过现场光谱测试及准同步采样，建立了遥感参数与含沙量的定量反演模式；邓孺孺等[9]在研究污染水体反射光谱成像的物理过程的基础上, 建立了叶绿素不敏感波段遥感数据像元反射率与污染物浓度的函数关系, 用多波段数据进行求解，并采用TM数据运用该方法对珠江口海域水污染进行定量分析，取得了很好的效果。

分析模型具有清晰的物理概念、机理，不需要像经验、半经验模型依赖大量的同步实测数据，适用性强，反演精度高，是目前定量遥感常采用的方法。

3 模型原理及推理

3.1 辐射传输过程分析

太阳光与水体相互作用的辐射传输过程如图1所示。

图1 太阳光在水中辐射传输过程

从图1中可以看出，到达水面的入射光—包括太阳直射光和天空散射光(天空光)在水面发生镜面反射，一部分反射光向上传输进入传感器；其余的入射光发生折射、透射进入到水中，其中一部分与水体中的组分如：水分子、悬浮泥沙、CODM等发生吸收散射，部分散射光经过水体二次衰减后，穿出水面进入传感器；在水深不大的前提下，另一部分入射光继续向下传输到达底部，被水底反射，反射光同样经过水体二次衰减后穿过水体，进入传感器，这部分反射光与水体的散射光共同组成水中光或叫离水反射辐射[10]。

3.2 水体辐射传输方程

基于上面的分析可以得知水域之上卫星传感器接收到的行星反射率主要由镜面反射比Rm和水体出水反射率Rwo两部分组成，即：

Rw=Rwo+Rm

(1)

水面镜面反射带有少量的水体本身的信息，它的强度主要与表面粗糙度有关，但是这部分反射光几乎不携带任何水质信息，对于水质遥感来说是干扰信息，应该予以剔除。镜面反射比Rm与波长无关，水在红外波段为强吸收，Rwo近似等于零，所以Landsat8数据的红外波段OLI6(1.560～1.660μm)的反射率约等于镜面反射比Rm，每个波段减去OLI6即可消除镜面反射比Rm即得到水体出水反射率Rwo：

Rwo=Rw-Rm=Rw-ROLI6

(2)

而水体出水反射率Rwo又由水中散射率Rws和水底对水体反射率的贡献Rg组成，即：

Rwo=Rws+Rg

(3)

又由具有朗伯反射地物的反射率与辐亮度的关系可得：

(4)

式(1)～(2)中，Lwo为水体的离水辐亮度；Eo为大气层外太阳平均辐照度；θ为太阳天顶角；

而水体离水辐亮度Lwo由两部分构成，即水中散射光向上的部分Lws和水底反射光Lg；

如果只考虑一次散射，并且假设水体中的组分只有水分子、悬浮泥沙和污染物，则水深h处的薄层水的散射对离水辐亮度Lwo的贡献为：

(5)

假设水深为H，对式(5)从0到H积分，得

(6)

将式(6)代入式(4)，得到水中的散射率Rw s：

(7)

同理，根据水体反射的物理机制，可以求出水底反射光对水体离水辐亮度的贡献为：

(8)

式(8)中Rb为水底底质的反射率；将式(8)代入(4)，可得到水底对水体反射率的贡献Rg：

Rg=Rbe-μ(α+β)H

(9)

将式(7)和式(9)代入(3)，得到体的出水反射率R_wo：

(10)

3.3 污染物和悬浮泥沙提取

由于河口水体比较浑浊，污染也比较严重，一般水域(水深>2m)到达水底的反射光极其微弱，可以忽略，其效果相当于水深极大，即H→∞，则式(10)可简化为[9]：

(11)

式(11)为本文所使用的最终水质反演模型。

4 模型在珠江口水质参数定量提取中的应用

本文所采用的遥感数据是2016年2月7日广州地区Landsat8的OLI传感器的多光谱数据，处理流程如图2所示。

图2 处理流程

4.1 辐射定标

(12)

式(12)中，Qcal为影像以16位量化的亮度值；Mρ为波段λ的反射率调整因子；Aρ为波段λ的反射率调整参数，Mρ和Aρ可从头文件中读取，其数值如表1所示。

表1 各波段定标系数

(13)

式(13)中θ为太阳天顶角，与太阳高度角互余，太阳高度角可从头文件中读取。

4.2 大气校正

大气辐射校正是进行定量遥感的先决条件, 其目的是将卫星遥感数据转换为地物真实信息的表观反射率，对后续的定量反演有着非常重要的影响。一般可通过三种途径进行大气散射校正，即辐射传递方程式计算法、野外波谱测度回归分析法及多波段图像的对比分析法[12]。

本文选择基于大气传输物理机制的暗像元法进行大气校正。该方法以山区阴影部分的植被作为黑体，先假设阴影区植被的反射率为0，估算出大气散射对程辐射的贡献，再利用迭代法对程辐射进行校正[13]。

传感器接收到的像元反射率R(λ)为：

(14)

式中，T'为入射方向大气透过率；ω为大气散射系数；Rg为地面反射率；T为垂向上大气透过率；Ρ(Θ)为散射相函数，Θ为散射角；

由式(14)可得阴影区植被行星反射率Rv为：

(15)

根据式(15)对所选取的一个暗像元进行计算，所求得的校正参数如表2。

表2 第一次计算后的校正参数

表3 第二次计算后的校正参数

由于第一次计算时，我们先假设在阴影区植被的反射率为0，把在阴影区植被处所读得的数值直接作为程辐射反射率处理，没有考虑植被对下行的天空光的反射，会存在一定误差，因此需要进行一步订正[13]。利用第一次计算的参数进行第二次迭代，计算得到的校正参数如表3所示。

为了检验大气校正的效果，从影像上选取了植被这一典型地物样本，读取植被某点的反射率值，并绘制大气校正前后反射率曲线的对比图，如图3所示。

波段

从图中可以看出，校正前植被反射率曲线与标准植被反射率曲线有很大偏差，特别是在反射峰附近偏差更大，而大气校正后植被反射率有了明显的改正，反射率曲线则与标准植被反射率曲线相似，这也说明了暗像元法是一种比较有效的大气校正方法。

4.3 污染物和悬浮泥沙浓度的提取

可见光波段能够很好地反映水体信息，对水质信息比较敏感。OLI2波段为蓝光波段，对水的穿透力最大，包含更多的水下细节和水深信息，不利于水质信息的提取；而OLI3波段为绿光波段，对水体具有一定的穿透能力，可反映水体浑浊度，对水体污染的研究效果好，OLI4波段为红光波段，是信息量最丰富的一个波段，对水中的悬浮泥沙反应敏感，因此本文选用OLI3和OLI4这两个波段进行计算。

本文研究的对象是水域的水质信息，所以对遥感图像进行水陆分离，剔除各波段的陆地信息部分，提取水体部分，可加快运算速度。由于水体在近红外波段强烈吸收，为低反射，可通过观察图像在近红外波段水陆交界处的实际反射率值，选取一个合适的值作为水陆分界的阈值，再根据植被在近红外波段的反射率高于红光波段这一特征将植被阴影区域分离出来[14]。

所使用的水体各组分的散射和吸收参数如表4所示。

表4 水质参数

将表4水体各组分的参数代入式(11)，建立一个二元一次方程组，可求解出污染物浓度和悬浮泥沙浓度的解析解，再利用PCI Geomatica V9.0图像处理软件建模对影像进行处理，得到珠江口水域的污染物和悬浮泥沙浓度分布图，最后使用ArcMap进行分级渲染，结果如图4和图6所示。

图4 珠江口污染物浓度

从图4中可以直观地看出珠江口西部海域污染物浓度总体比东侧的高，对比图5的实际影像中，也可以看出西岸水体呈现灰黄色，确实比东岸的浑浊。

但是其中东侧水域中，在三角电厂、东宝河入海口、深圳机场福永码头和深圳湾附近污染物浓度明显比周围海域的都高，从影像中也可以观察到该区域的水体呈现黑黄色，污染物浓度高。

图5 珠江口部分水域影像(真彩色合成)

图6 珠江口悬浮泥沙浓度分布

从图6中，可以看出珠江口西部海域悬浮泥沙浓度明显比东侧的高，特别是在南沙十九涌以南的水域，从影像中也可以清晰地看出该区域水呈泥黄色，比较浑浊；一些河流的入海口处水体也比较浑浊，悬浮泥沙浓度较高。

5 结论

(1)从最终的反演结果来看，珠江口污染物和悬浮泥沙浓度的分布与实际情况基本相符，西部水域这两种污染物的浓度大体比东岸的高，但是东岸存在着典型的污染源，这与附近的地区的生产活动相关，这就为水污染治理提供了有效的信息，有助于治理方案的制定和实施，这也是水质定量遥感所带来的社会效益和经济效益。

(2)本文使用暗像元法进行大气校正，假设大气是均质的，而实际上大气的结构和成分是很复杂的，在时空分布上也很不均匀，而且只选取单一暗像元对整幅影像进行校正，会产生一定的误差。另外最终的结果也缺乏卫星过境时地面同步实测数据来进行精度检验，明显还存在不足之处。

(3)总的来说，该模型物理机理清晰，是在分析入射光在水体中的传输过程的基础上建立的，模型中的参数也比较少，也不像经验半经验模型需要做大量的实验，而且有些参数，如清洁水体的吸收、散射系数等适用性比较高，不需要花费大量人力物力进行重复测量，可以应用于其他数据或区域的反演，模型有较好的适用性，可行性也比较高。