基于无人机多光谱影像和OPT

2023-12-23 01:00| 来源: 网络整理| 查看: 265

水资源是生态环境的重要组成和基础保障, 定期进行水质监测十分必要.传统的水质监测方法主要为定点测量, 这种方法需要消耗大量的人力物力, 且因为监测点数量的限制, 只能获取点状分布的数据, 难以获取整个水域的水质状况.遥感技术具有高效、经济成本低以及覆盖范围广等优势, 目前已被广泛应用于湖泊、水库以及河流等水域的水质监测[1~6].例如, Torbick等[7]使用Landsat TM影像对密歇根州下半岛内陆湖泊的叶绿素a、总磷、总氮以及悬浮物浓度等多个水质参数进行了测绘; 朱利等[8]利用GF-1卫星WFV数据对太湖水域进行叶绿素a和透明度的遥感监测, 以上研究均表明卫星遥感技术在水质监测方面已得到应用.

近年来, 无人机遥感技术因其更高的分辨率以及便捷性等优点被广泛应用于各个领域中, 且得到良好的效果, 如魏鹏飞等[9]使用无人机多光谱影像对不同生育期的夏玉米叶片氮含量进行了遥感估测; 张智韬等[10]使用无人机多光谱影像对土壤含水率进行了遥感反演研究; Flynn等[11]使用无人机数据对浅水非浑浊河流水下植被进行了遥感监测等.

无人机遥感技术在水质监测方面的研究也日渐展开, 如刘彦君等[12]将无人机结合多光谱影像应用于小微水域的水质要素反演, 其中总磷、悬浮物浓度、浊度的反演模型决定系数分别达到0.7829、0.7503和0.7334, 为小微水域的污染防治提供了技术支撑; Guimaraes等[13]将神经网络应用于无人机图像预测水体悬浮物浓度的回归分析中, 回归模型的质量得到了显著的提升.尽管无人机技术在水质监测中的应用已有了一些进展, 但如何制作出稳定准确的水质参数反演模型仍是目前研究的重点.为了探索更精确的基于无人机影像的水质参数反演方法和模型, Su等[14]在处理多光谱数据时提出了MPP(matching pixel-by-pixel)算法, 所得反演模型决定系数R2接近于1, 但在线性回归分析中, 过高的R2被认为是过拟合的结果, 因此在后续的实验中Su等[15]对回归模型的预测性能进行了交叉验证, 证明了训练样本点数量的增加对于减轻过拟合现象以提高回归模型的预测性能是有效的, 但这种做法并没有解决运算量过大的问题.

本文以杭州市青山湖为研究区域, 提出一种基于无人机多光谱影像的OPT-MPP(Optimize-MPP)算法, 以克服MPP算法运算量过大和所得模型过拟合的问题, 并实现悬浮物浓度和浊度的反演.

1 材料与方法 1.1 研究区域和采样点

青山湖水域(N30°12′58″~30°15′21″, E119°44′55″~119°47′58″), 位于杭州市临安区东侧, 水域面积约为10 km2.本实验共设置4个采样区域, 分别位于苕溪入口、森林公园、琴山码头以及青山湖大坝处, 共计45个采样点, 密度约为4.5个·km-2, 并在取样时分别记录经纬度坐标.于2019年7月对4个采样区域分别进行4次无人机多光谱影像获取和同步水质采样, 其中悬浮物浓度使用杭州陆程仪器有限公司LH-SS1型悬浮物测定仪测量, 浊度使用杭州陆恒生物科技有限公司WGZ-1B浊度仪测量, 共获取45组数据, 根据散点图剔除3组异常数据后, 保留剩余42组数据, 选取均匀分布在各研究区域且具有一定梯度的14组数据用于建模, 另外28组数据用于模型精度检验(一般情况下采用2/3数据建模, 1/3数据检验, 由于本实验中参与到建模的运算数据量足够多, 且为了得到更精确的误差分析数据, 增加了检验数据).表 1为青山湖水域水质参数统计, 研究区域和采样点如图 1.

表 1 (Table 1) 表 1 水质参数统计 Table 1 Statistics of water quality parameters 水质参数最小值最大值平均值标准差悬浮物浓度/mg·L-1 15.65 87.25 28.57 15.35 浊度/NTU 10.3 44.1 20.08 7.16 表 1 水质参数统计 Table 1 Statistics of water quality parameters 图 1 Fig. 1 图 1 研究区域和采样点示意 Fig. 1 Research area and sampling point 1.2 无人机影像

本研究使用大疆悟2四旋翼无人机携带MicaSense RedEdge多光谱传感器获取遥感数据, 该传感器具有质量轻、体积小的特点, 适合搭载于中小型无人机.本实验航高设置为400 m, 航速5 m·s-1, 以更加稳定地获取研究区域的多光谱影像, 相对于卫星遥感影像数据, 由于无人机的飞行高度较低, 获取的影像区域较小, 大气折射和地球曲率的影响可以忽略[16~18], 在无人机影像数据预处理时, 使用与RedEdge多光谱传感器配套的MicaSense Calibrated Reflectance Panel配合pix4Dmapper软件进行辐射校正与辐射定标, 最终在ENVI软件中读取与采样点对应的像元区域数据, 以此作为光谱反射率数据(Rrs), 所得光谱曲线如图 2, 其呈现出一般内陆水体的光谱特征[19].表 2为传感器获取的各波段名称及其对应的中心波长和波宽, 图 3为各采样区域无人机多光谱影像图.

图 2 Fig. 2 不同颜色的曲线分别表示不同采样点的光谱曲线图 2 采样点水体光谱曲线 Fig. 2 Spectral curve of water in sample points 表 2 (Table 2) 表 2 各波段中心波长及波宽 Table 2 Center wavelength and bandwidth of each band 波段波段名中心波长/nm 波宽/nm R1 Blue 475 20 R2 Green 560 20 R3 Red 668 10 R4 NIR 840 40 R5 Red Edge 717 10 表 2 各波段中心波长及波宽 Table 2 Center wavelength and bandwidth of each band 图 3 Fig. 3 (a)1号区域苕溪入口; (b)2号区域森林公园; (c)3号区琴山码头; (d)4号区青山湖大坝图 3 各区域无人机多光谱影像 Fig. 3 Unmanned aerial vehicle(UAV) multispectral imagery of various areas 1.3 MPP算法

基于多光谱影像数据反演水质参数时, 首先需要分析采样点样本的实测值与遥感影像中该点所对应的像元值之间的相关性, 并建立回归模型.由于无人机影像的分辨率极高, 每个像元所代表的实际区域很小, 很难精确匹配采样点和像元, 因此通常使用样本点周围i×j个像元的平均值参与拟合回归模型(如图 4).这种方式虽然简洁方便, 但不能最大程度规避数据误差, 很难找到最优模型, 尤其是样本量较少时.而MPP算法的目的就是通过遍历每一个可能的像元来寻找最优的模型.

图 4 Fig. 4 图 4 3×3像元取平均值 Fig. 4 Average of 3×3 pixels

MPP算法本质上是挑选不同的像元组合进行多次回归分析, 并基于遗传算法寻找最优解.设Pij(m)是第m个采样点的第(i, j)个像元的值, 其中i为列标, j为行标, 则每个采样点有i×j种可能的像元选择, m个采样点产生(i×j)m种像元组合, 将每种组合均与实测水质参数进行一次回归分析, 总计进行(i×j)m次回归分析, 将决定系数R2最大的模型作为最优回归模型.因此, MPP算法的时间复杂度主要由i、j与m的大小决定, i与j的数值越大, 每个采样点的像元选择越多, 相应地m的数值越大, 产生的像元组合数目越多.

假设实验样本点个数为8, i与j的维度均为3, 则所有可能的候选组合个数为(32)8种, 图 5展示出了其中一种可能的解结果, 其像元组合为P12(1)、P11(2)、P21(3)、P31(4)、P33(5)、P23(6)、P13(7)和P32(8), 用于与实测水质参数数据建立反演模型.

图 5 Fig. 5 图 5 3×3像元区域及一种可能解 Fig. 5 The 3×3 pixel area and one of the possible solutions

图 6为传统平均值建模方法与MPP算法建模方法的对比, 相较于传统平均法建模方式, MPP算法需构建组合矩阵以找到与水质参数具有最佳相关性的数据.

图 6 Fig. 6 图 6 建模方法对比 Fig. 6 Comparison diagram of the modeling method

虽然MPP算法相较于传统平均值法已有了较大改进, 但本文认为MPP算法仍存在如下缺陷：①组合数量过多, 计算量大, 导致运算速度过慢, 仅能适用于较少的采样点个数; ②将R2作为模型选择的唯一标准欠妥, 容易产生模型过拟合现象.

1.4 OPT-MPP算法

针对MPP算法计算量大运算速度缓慢的缺陷, 本研究将3×3像元区域内的数据按列从小到大排序, 选取大小排在2、5和8位的数据参与后续运算, 相比MPP算法而言将每个像元区域内所需处理的数据量从9个减少到了3个, 将9m种组合运算减少到3m种, 其中m表示采样点的个数, 通过减少单个采样点参与运算的数据量, 以使更多的其他样本数据参与运算, 处理过程如图 7.

图 7 Fig. 7 图 7 OPT-MPP算法处理示意 Fig. 7 Processing diagram of the OPT-MPP algorithm

针对MPP算法可能产生过拟合的现象, OPT-MPP算法在筛选模型时设置双重条件, 即决定系数ROPT2大于RAVE2(传统平均值建模方法所得模型决定系数), 同时平均相对误差EOPT小于EAVE(传统平均值建模方法所得模型相对误差), 以此产生模型集合S：

式中y表示拟合模型, 最终从集合S中筛选出误差最小模型M：

式中, yE表示拟合模型y的平均相对误差.这种双重条件的筛选模式, 去除了决定系数较高但精度较低的拟合模型, 这样既能避免模型产生过拟合现象, 又能保证反演模型的精度, 建模过程如图 8.

图 8 Fig. 8 图 8 OPT-MPP算法建模方法 Fig. 8 Modeling method for the OPT-MPP algorithm 2 结果与讨论 2.1 相关性分析

根据悬浮物浓度和浊度的光谱响应特征, 参考已有的特殊波段组合方式, 构建了8种满足本实验需求的波段组合, 包括R4/R3、R4/R2、R4/R1、R4/(R1+R2+R3)、(R2+R3)/R4、(R2+R3)×R4、(R3+R4)/R2和R3/R4[20~22].悬浮物浓度和浊度进行相关性及显著性分析后, 选用与悬浮物浓度相关性最高的R4/R1波段组合, 与浊度相关性最高的R4/(R1+R2+R3)波段组合进行实验, 其显著性均满足小于0.05的统计学要求[23].表 3为相关性及显著性分析结果.

表 3 (Table 3) 表 3 相关性及显著性分析结果1) Table 3 Analysis result for correlation and significance 编号波段或波段组合悬浮物浓度浊度相关性显著性相关性显著性 B1 R1 -0.530 0* -0.376 0.011* B2 R2 -0.353 0.017* -0.265 0.079 B3 R3 -0.343 0.021* -0.321 0.031* B4 R4 0.360 0.015* 0.478 0.001* B5 R5 0.057 0.711 0.089 0.560 B6 R4/R3 0.575 0* 0.639 0* B7 R4/R2 0.649 0* 0.660 0* B8 R4/R1 0.701 0* 0.663 0* B9 R4/(R1+R2+R3) 0.685 0* 0.692 0* B10 (R2+R3)/R4 -0.536 0* -0.566 0* B11 (R2+R3)×R4 -0.010 0.947 0.110 0.470 B12 (R3+R4)/R2 0.606 0* 0.566 0* B13 R3/R4 -0.501 0* -0.559 0* 1)*表示通过0.05水平的显著性检验, 黑体字表示相关性最大值表 3 相关性及显著性分析结果1) Table 3 Analysis result for correlation and significance 2.2 模型结果分析

本实验使用决定系数R2、平均相对误差ARE、均方根误差RMSE、相对均方根误差rRMSE以及综合误差CE这5种指标来评价模型精度[24, 25], 模型决定系数R2越大, ARE、RMSE、rRMSE和CE越小, 则模型的估算效果越好.相关公式如下：

(1) (2) (3) (4) (5)

式中, yi表示第i个点的实测值,ŷi表示第i个点的预测值, y表示水质参数实测值的平均值, n表示验证样本数.

将检测得到的悬浮物浓度、浊度数据以及经OPT-MPP算法处理后的水体多光谱数据结合, 依次构建线性(U)、多项式(PL)、幂函数(P)和指数函数(E)这4项基本模型, 并计算决定系数R2及各项模型精度指标.悬浮物浓度和浊度建模结果与精度评价指标分别如表 4和表 5, 其中表 4的x为R4/R1, y表示悬浮物浓度, 表 5的x为R4/(R1+R2+R3), y表示浊度.

表 4 (Table 4) 表 4 悬浮物浓度回归模型1) Table 4 Regression model of the suspended sediment concentration 方法编号拟合模型 R2 ARE/% RMSE rRMSE CE 传统平均值法 ASS(U) y=30.28x+0.49 0.701 7 17.87 5.101 2 0.178 5 0.178 6 ASS(PL) y=35.06x2-35.43x+28.31 0.775 1 13.11 3.820 9 0.133 6 0.132 4 ASS(P) y=30.88x1.032 0.702 2 18.35 5.115 0.179 0.181 3 ASS(E) y=9.26e1.12x 0.757 3 15.46 4.321 1 0.151 2 0.152 9 未设置双重筛选条件 RSS(U) y=35.72x-3.65 0.838 9 20.35 5.525 8 0.193 4 0.198 5 RSS(PL) y=38.16x2-36.35x+27.09 0.910 2 14.16 4.030 2 0.141 1 0.141 4 RSS(P) y=31.55x1.219 0.849 8 20.16 5.333 7 0.186 7 0.194 2 RSS(E) y=8.33e1.28x 0.9 16.68 4.482 4 0.156 9 0.161 9 OPT-MPP算法 ISS(U) y=27.51x+2.35 0.702 1 16.53 4.935 2 0.172 7 0.169 ISS(PL) y=37.12x2-41.03x+30.44 0.787 0 12.35 3.944 3 0.138 0 0.130 8 ISS(P) y=30.15x0.95 0.702 6 17.15 5.004 1 0.175 2 0.173 4 ISS(E) y=10.27e1.03x 0.761 6 14.52 4.300 1 0.150 5 0.147 9 1)黑体字表示综合误差最小值表 4 悬浮物浓度回归模型1) Table 4 Regression model of the suspended sediment concentration 表 5 (Table 5) 表 5 浊度回归模型1) Table 5 Regression model of turbidity 方法编号拟合模型 R2 ARE/% RMSE rRMSE CE 传统平均值法 ATU(U) y=106.3x-9.34 0.680 4 23.02 5.264 8 0.262 2 0.246 2 ATU(PL) y=625.4x2-275.8x+45.57 0.803 2 16.87 4.409 4 0.219 6 0.194 2 ATU(P) y=146.2x1.58 0.714 6 22.48 5.388 4 0.268 3 0.246 6 ATU(E) y=4.55e5.09x 0.761 8 20.54 5.002 0.249 1 0.227 3 未设置双重筛选条件 RTU(U) y=116.42 x-12.13 0.820 0 24.75 5.745 6 0.286 1 0.266 8 RTU(PL) y=524.62x2-211.36x+35.71 0.891 6 17.82 4.532 4 0.225 7 0.202 RTU(P) y=161.07x1.66 0.848 8 23.32 5.643 3 0.281 0.257 1 RTU(E) y=4.6e5.1x 0.874 5 20.78 5.065 9 0.252 3 0.230 1 OPT-MPP算法 ITU(U) y=97.76x-8.07 0.681 7 19.45 4.611 9 0.229 7 0.212 1 ITU(PL) y=572.5x2-265.9x+46.19 0.8043 13.29 3.365 4 0.167 6 0.150 3 ITU(P) y=126.38x1.51 0.715 5 19.25 4.668 6 0.232 5 0.212 5 ITU(E) y=5.09e4.58x 0.762 5 17.13 4.144 9 0.206 4 0.188 9 1)黑体字表示综合误差最小值表 5 浊度回归模型1) Table 5 Regression model of turbidity

由表 4可以看出, 在悬浮物浓度各拟合模型中, ISS(PL)拥有OPT-MPP算法建模方式所得模型中最高决定系数R2和最低综合误差CE, 相较于传统平均值法最优模型ASS(PL)(最高决定系数R2及最低综合误差CE), 决定系数R2从0.775 1提升到0.787 0, 说明在建模过程中拟合效果更好.此外平均相对误差ARE从13.11%降低到12.35%, 均方根误差RMSE与相对均方根误差rRMSE略有上升, 但总体综合误差CE从0.132 4下降到0.130 8, 说明OPT-MPP算法得到的拟合模型检验精度更高, 更适用于悬浮物浓度反演.未设置双重筛选条件所得4种模型虽然拥有几种建模方式下较高的决定系数R2, 但其综合误差均相对较大, 说明这4种模型产生了过拟合现象, 不适用于悬浮物浓度反演.

浊度各拟合模型如表 5, 在浊度各拟合模型中, ITU(PL)拥有OPT-MPP算法建模方式所得模型中最高决定系数R2和最低综合误差CE, 相较于传统平均法最优模型ATU(PL)(最高决定系数R2及最低综合误差CE), 决定系数R2从0.803 2提升到0.804 3, 说明在建模过程中拟合效果更好.此外平均相对误差ARE从16.87%降低到13.29%, 均方根误差RMSE从4.409 4下降到3.365 4, 相对均方根误差rRMSE从0.219 6下降到0.167 6, 综合误差CE从0.194 2下降到0.150 3, 说明OPT-MPP算法得到的拟合模型检验精度更高, 更适用于浊度反演.未设置双重筛选条件所得4种模型虽然拥有几种建模方式下较高的决定系数R2, 但其综合误差均相对较大, 说明这4种模型产生了过拟合现象, 不适用于浊度反演.

表 6展示了ISS(PL)模型与ITU(PL)模型所对应的波段组合值, 即构建拟合模型所采样用的14个采样点对应波段组合值.图 9为ISS(PL)模型与ITU(PL)模型的拟合图像, 从中可以看出用于构建拟合模型的14个采样点均匀分布在拟合曲线两侧, 证明悬浮物浓度和浊度在这两种模型的拟合程度较好.

图 9 Fig. 9 图 9 ISS(PL)与ITU(PL)模型拟合 Fig. 9 Model fitting of ISS(PL)and ITU(PL) 表 6 (Table 6) 表 6 OPT-MPP算法所得最优解1) Table 6 Optimal solution by the OPT-MPP algorithm 回归模型波段组合值 y=37.12x2-41.03x+30.44y:悬浮物浓度x: R4/R1 m=1 m=2 m=3 m=4 m=5 m=6 m=7 i=3 i=3 i=3 i=3 i=3 i=3 i=3 1.075 1 0.850 1 0.792 5 0.667 3 0.712 0 0.691 2 0.861 2 m=8 m=9 m=10 m=11 m=12 m=13 m=14 i=3 i=3 i=1 i=2 i=3 i=3 i=2 0.670 1 0.555 6 0.440 0 0.972 5 0.915 2 1.406 5 1.362 1 y=572.5x2-265.9x+46.19y:浊度x: R4/(R1+R2+R3) m=1 m=2 m=3 m=4 m=5 m=6 m=7 i=3 i=3 i=3 i=3 i=3 i=3 i=1 0.350 3 0.251 5 0.295 8 0.310 3 0.202 4 0.188 0 0.233 8 m=8 m=9 m=10 m=11 m=12 m=13 m=14 i=3 i=2 i=1 i=1 i=3 i=3 i=3 0.269 6 0.278 7 0.217 8 0.250 9 0.436 5 0.421 8 0.414 4 1)m表示采样点编号, i表示光谱反射率编号表 6 OPT-MPP算法所得最优解1) Table 6 Optimal solution by the OPT-MPP algorithm

根据构建的拟合模型, 结合28个被用于模型检验的采样点水质参数实测数据, 分别计算ISS(PL)与ITU(PL)模型的反演估测值, 并制作实测值与估测值分布曲线图, 由图 10(a)可以看出ISS(PL)模型估测的悬浮物浓度值与对应的实测值在绝大部分下比较接近.由图 10(b)可以看到ITU(PL)模型估测的浊度值除15号点(推测为异常点)以外的其他点均与对应实测值较为接近, 说明这两个模型适用于悬浮物浓度与浊度反演.

图 10 Fig. 10 图 10 实测值和预测值对比 Fig. 10 Comparison diagram of estimated values and measured values

表 7展示了本实验中各种建模方法所需的时间, 在用于模型制作的样本点个数均为14个的条件下, 传统平均值法运用MATLAB软件中的curve fitting拟合工具箱, 拥有较快的运算速度; MPP算法因运算量过大导致运算速度缓慢, 需要约一周甚至更长的运算时间; 而本实验所使用的OPT-MPP算法在运算时间上做出了较大的优化, 在保障模型精度的同时可以于2 h左右完成回归模型的制作.

表 7 (Table 7) 表 7 建模方法所需时间 Table 7 Time required for the modeling method 建模方法样本点个数/个所需时间/h 传统平均值法 14 约0.1 MPP算法 14 约168 OPT-MPP算法 14 约2 表 7 建模方法所需时间 Table 7 Time required for the modeling method 2.3 水质参数空间分布

将最优模型ISS(PL)与ITU(PL)分别代入Band Math波段运算工具, 制作实验各区域的悬浮物与浊度空间分布[26, 27], 悬浮物空间分布如图 11, 浊度空间分布如图 12.从中可以看到, 区域1和区域4中, 悬浮物和浊度分布近岸边略高, 推测是因为近岸边设有绿道, 与人类活动有关; 区域3的近岸区域两种参数在进行波段运算后的结果多呈红偏绿色, 说明浓度值较高, 与近期正在进行的青山湖绿道建设存在一定关系; 而区域2为青山湖国家森林公园附近, 近岸边多为树林, 具有较好的生态环境.

图 11 Fig. 11 (a)1号区域苕溪入口; (b)2号区域森林公园; (c)3号区琴山码头; (d)4号区青山湖大坝; 图例自下向上表示浓度值从低到高, 下同图 11 各区域悬浮物浓度空间分布 Fig. 11 Spatial distribution of the suspended sediment concentration in various regions 图 12 Fig. 12 (a)1号区域苕溪入口; (b)2号区域森林公园; (c)3号区琴山码头; (d)4号区青山湖大坝图 12 各区域浊度空间分布 Fig. 12 Spatial distribution of the turbidity in various regions 3 结论

(1) 本实验过程中, OPT-MPP算法经过前期数据处理, 将用于计算的数据量减少为MPP算法的(1/3)m, 其中m为采样点个数, 使得计算量大大降低.

(2) OPT-MPP算法计算所得最佳悬浮物反演模型ISS(PL)相较于传统平均值法所得最优模型ASS(PL), 决定系数R2从0.775 1提升到0.787 0, 综合误差CE从0.132 4下降到0.130 8, 表明模型ISS(PL)更适用于悬浮物反演.最佳浊度模型ITU(PL)相较于模型ATU(PL), 决定系数R2从0.803 2提升到0.804 3, 综合误差CE从0.194 2下降到0.150 3, 表明模型ITU(PL)更适用于浊度反演.

(3) 本方法不仅仅适用于悬浮物浓度与浊度两项水质参数指标, 同样适用于其他参数, 也可推广到其他相关的无人机影像领域, 特别是样本量较小的研究.

(4) 利用无人机遥感影像结合建模算法进行水质参数浓度反演, 可以有效提高反演模型的精度与反演质量, 得到更加准确的水质参数空间分布情况, 对于更准确地评估水域水质参数的时空变化和更好地进行水域治理等具有重要的应用价值和实际意义.但无人机遥感反演仍处于探索阶段, 在实际研究中也存在着一些问题和困难, 由于实验环境、仪器和时间等因素的限制, 本实验存在个别异常点且未能获取到更长时间跨度的数据, 在今后研究中将着重考虑时间因素, 构建多时段反演模型, 对青山湖水域进行长时序的水质参数变化状况分析.

致谢: 临安区五水共治工作领导小组办公室对本项目提供了各方面的协助与支持, 特此致谢！

【本文地址】

公司简介

联系我们