基于改進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影像進行對比分析實驗，本文方法很好地解決了礦山水體的空間異質性、地物混淆和光譜變異的問題，雖然精確率和召回率不是最高的，但是綜合指標是最高水平，由于綜合指標更具有算法評價能力，從而驗證了本文方法的有效性，表明本文方法可以應用于礦區復雜水體的提取工作中。