基于改进TSUWI算法在山西矿山复杂水体提取中的应用研究

尹德威，闫新珠

(山西省地质测绘院，山西 运城044000)

矿山排污作为矿产资源开发过程的必要环节，会造成矿山生态环境的严重污染。开展对矿山水体的提取工作，对矿山水环境保护和预防尾矿库事故具有重要的意义[1]。

水体提取的方法包括人工提取和遥感方法提取，人工提取方法效率低下，而遥感手段，由于具有周期短、大范围、波谱丰富等特点，成为了提取水体的主流方法[2]。目前，遥感水体提取的主要方法有：单波段阈值法、水体指数法、面向对象法等[2]。单波段阈值法通过单一波段的光谱差异来区分水体与其他地物，对于光谱复杂的水体提取精度不高；一些研究人员使用基于对象的方法来检测水体，并借助纹理等附加特征，但非常耗时，且改进效果有限；而水体指数法是通过特征波段间的比值运算，以区分水体与其他地物。

水体指数法由于提取水体的高效性，是工业生产中的主流方法。典型的有McFEETERS[3]提出的归一化差分水体指数法(NDWI)，但NDWI不能很好地抑制建筑物的影响。因此，徐涵秋[4]提出了改进的归一化差分水体指数法(MNDWI)，提高了对建筑的区分能力，但是由于阴影的影响，提取结果的精度受到限制。考虑到两者的不足，Wu等人[5]提出将城市水体指数(UWI)与城市阴影指数(USI)结合，利用两步水体指数法(TSUWI)来提取水体，取得了较好的效果。然而，TSUWI的阈值采用的是经验阈值，难以适用于不同的影像数据。因此，有人提出了自适应阈值的方法，其中以大津阈值(OTSU)方法[6]最为常见，但是OTSU阈值法只考虑到水体与地物之间光谱差异性，对于光谱变异的水体不能起到很好的区分效果。Shahriari等人[7]提出的一种基于分形技术的面积分形阈值算法，由于考虑了空间光谱特性，在影像分割中起到了很好的效果。

本文针对矿区复杂水体的提取，提出了改进的TSUWI(MTSUWI)水体提取方法。以山西新柳矿区为研究对象，首先分别获取UWI和USI指数；然后，创新性改进了面积分形阈值算法，并分别对UWI与USI指数进行阈值确定与水体提取；最后，将两种水体提取结果进行合并，通过提取交集部分的水体实现矿区复杂环境下的水体信息的高精度自动提取。

1 研究区概况与数据源介绍

实验区为山西新柳矿区，位于山西省吕梁市孝义市西南阳泉曲镇境内。矿区内地势复杂，属于山地地形。年降水量在450mm至550mm之间，四季降水量不均衡，夏季降水量集中，易于形成降水蓄积。

数据源为2016年8月3日获取的GF2数据，GF2数据处于可见光和近红波段范围，包括四个多光谱波段和一个全色波段，多光谱波段分辨率为4m，全色波段分辨率为1m。

使用ENVI5.3对数据进行预处理，由于矿区处于山地地形，地势复杂，所以先结合DEM数据对多光谱与全色波段进行正射校正。然后采用FLASSH模型对数据进行大气校正，消除大气的影响。最后，采用Gram-Schmidt Pan Sharpening模型进行多光谱数据与全色波段的融合。融合后数据的分辨率达到1m，且同时具备红、绿、蓝、近红四个波段，显著提高了数据的光谱空间特性。本文裁剪了2645*1897像素大小的区域作为实验数据，实验数据如图1所示。

2 MTSUWI算法提取矿区水体

矿区复杂水体具有光谱变异、空间异质性和易受阴影与裸土等地物干扰的特点，针对这些提取难点，本文提出MTSUWI方法，包括三个主要步骤：首先构建UWI和USI特征图，然后基于面积阈值法提取水体，最后将提取结果合并以求得最后的水体提取结果。

2.1 UWI与USI特征图的构建

UWI指数考虑了近红、红和绿波段的信息，USI指数则考虑了近红、红、绿和蓝四个波段的信息，相比于传统水体指数NDWI只考虑了近红和绿波段的信息，UWI与USI更加充分地利用了光谱信息，能更稳健地应对光谱变异和其他地物的光谱混淆问题。参考Wu等人[5]提出的指数参数，UWI和USI构建为式(1)和式(2)：

其中，NIR、R、G、B分别表示近红、红、绿、蓝波段，α与β表示偏置，需要根据水体所处的不同环境来确定。式(1)和式(2)所得的特征图UWI与USI，将水体部分进行增强高亮，接下来需要确定阈值将高亮的水体部分提取出来。

2.2 基于面积分形阈值法的水体提取

面积分形阈值法最初应用于影像分割领域，本文将其改进以应用于矿区水体阈值的自动确定。其模型可表示为：

其中，r>0表示特征值，这里为指数值，C>0表示比例常数，D>0为一般分维数。

其中，N(rj)表示相同的rj特征值对应的个数或者频率，rn表示最大特征值，n为特征值个数。根据公式(4)可以获得数据集(N(r≥ri)，ri)，将数据集(N(r≥ri)，ri)代入式(3)并两边取对数，则问题简化为一元线性回归方程：

如图2所示，ln(N(r≥ri))关于ln(ri)的散点(以小圆表示)会呈现出分段分布特点，通过最小二乘法可以按式(5)拟合出R1、R2、R3线段(以线段表示)，分别对应D1、D2、D3区间。其中D3的分布区间表示高亮的区间，对应的是纯水体，D2分布区间表示非典型水体，比如污染严重、水较浅、水体面积较小等受噪声污染的水体，而D1分布区间则表示非水体。因此，我们将阈值确定为T1，这样可以很好地将纯水体和光谱变异与空间异质的水体提取出来。因此，基于UWI与USI指数，我们可以分别确定阈值Tuwi和Tusi，所以水体提取结果分别为Ruwi和Rusi：

图2 ln(N(r≥ri))关于ln(ri)的分布图与拟合线段

2.3 水体提取结果的合并

按照式(6)和式(7)所得结果能够提取出水体，但是会存在噪声的影响，因此，我们将提取的结果Ruwi与Rusi进一步融合，以消除噪声的影响。融合方式如式(8)：

通过式(8)求两种提取水体的交集后，就可以获得最终的水体提取结果。

3 实验结果与分析

3.1 实验结果对比

本文通过对比分析的方法来验证本文方法的有效性，采用了经典的水体提取算法NDWI[3]、UWI[5]、USI[5]、TSUWI[5]四种对比方法与本文的MTSUWI算法进行对比。NDWI采用OTSU阈值法确定阈值，UWI、USI与TSUWI采用文献[5]中的设置，本文算法MTSUWI中涉及的参数α[公式(1)]和参数β[公式(2)]分别设置为1.3和0.9。各个算法的实验结果如图3所示，其中参考提取结果是通过易康软件结合Google地球影像进行的人工标注结果。

结果的量化评价指标采用精确率和召回率以及综合指标[8]，分别计算为：

其中，TP表示正确提取水体的数量，FP表示错误提取水体的数量，FN表示将水体分为其他类别的数量。按公式(9)、公式(10)和公式(11)计算各个算法的评价指标，列在表1中。

表1 不同方法的评价指标

3.2 实验结果分析

根据图3的可视化结果，可以看到，NDWI的结果很差，主要原因是不能区分地物混淆的情况，因此表现在定量结果上，有较高的召回率，但是精确率却非常低。UWI的提取结果定量上表现出较高的精确率，但召回率却很低，体现在可视化结果上就是水体的漏分情况较为严重，主要原因是受到了矿区水体的空间异质性影响。而USI的提取结果则表现出了较高的精确率和0.808的召回率，表现出了相对强地物混淆抑制能力，但一定程度上受到水体光谱变异和空间异质性的影响。TSUWI的结果是UWI与USI的合并求交集的结果，所以，当两者的漏分情况严重时，TSUWI是不能取得较好精度的。而本文的MTSUWI算法，在可视化结果上，对于矿区水体的空间异质性、地物混淆和光谱变异的问题，表现出了较好的提取效果，且综合指标优于其他方法，虽然精确率与召回率不是最高的，但结合其他方法的定量和可视化结果，可以发现综合指标更具有评价方法优劣的能力。因此，本文方法相比于对比方法，更加适用于矿区水体的提取。

图3 不同方法的实验结果

3.3 MTSUWI偏置参数的影响

在公式(1)和公式(2)中，设置了两个偏置参数α和β，主要作用是为了平衡矿区水体环境因素的影响，从而更好地基于矿区环境获取UWI与USI特征图。为了探究偏置参数对结果的影响，本文设置了α和β的系列值，α∈[0.7，0.8，0.9，1.0，1.1，1.2，1.3，1.4，1.5，1.6]，β∈[0.3，0.4，0.5，0.6，0.7，0.8，0.9，1.0，1.1，1.2]，并分别计算α和β不同组合下，MTSUWI对应的综合指标的值，所得结果显示在图4中，可以得出当α=1.3，β=0.9时，综合指标的值达到峰值。所以，针对矿区环境的偏置参数分别为α=1.3，β=0.9。

图4 偏置参数α与β对MTSUWI方法综合指标的影响

4 结束语

本文提出了MTSUWI矿山水体提取算法，通过将面积分形阈值算法创新性地引入到UWI与USI的阈值提取过程，改进了TSUWI水体提取算法。通过对新柳矿区的GF2影像进行对比分析实验，本文方法很好地解决了矿山水体的空间异质性、地物混淆和光谱变异的问题，虽然精确率和召回率不是最高的，但是综合指标是最高水平，由于综合指标更具有算法评价能力，从而验证了本文方法的有效性，表明本文方法可以应用于矿区复杂水体的提取工作中。