尹德威，閆新珠

(山西省地質(zhì)測繪院，山西 運(yùn)城044000)

礦山排污作為礦產(chǎn)資源開發(fā)過程的必要環(huán)節(jié)，會造成礦山生態(tài)環(huán)境的嚴(yán)重污染。開展對礦山水體的提取工作，對礦山水環(huán)境保護(hù)和預(yù)防尾礦庫事故具有重要的意義[1]。

水體提取的方法包括人工提取和遙感方法提取，人工提取方法效率低下，而遙感手段，由于具有周期短、大范圍、波譜豐富等特點(diǎn)，成為了提取水體的主流方法[2]。目前，遙感水體提取的主要方法有：單波段閾值法、水體指數(shù)法、面向?qū)ο蠓ǖ萚2]。單波段閾值法通過單一波段的光譜差異來區(qū)分水體與其他地物，對于光譜復(fù)雜的水體提取精度不高；一些研究人員使用基于對象的方法來檢測水體，并借助紋理等附加特征，但非常耗時，且改進(jìn)效果有限；而水體指數(shù)法是通過特征波段間的比值運(yùn)算，以區(qū)分水體與其他地物。

水體指數(shù)法由于提取水體的高效性，是工業(yè)生產(chǎn)中的主流方法。典型的有McFEETERS[3]提出的歸一化差分水體指數(shù)法(NDWI)，但NDWI不能很好地抑制建筑物的影響。因此，徐涵秋[4]提出了改進(jìn)的歸一化差分水體指數(shù)法(MNDWI)，提高了對建筑的區(qū)分能力，但是由于陰影的影響，提取結(jié)果的精度受到限制?？紤]到兩者的不足，Wu等人[5]提出將城市水體指數(shù)(UWI)與城市陰影指數(shù)(USI)結(jié)合，利用兩步水體指數(shù)法(TSUWI)來提取水體，取得了較好的效果。然而，TSUWI的閾值采用的是經(jīng)驗閾值，難以適用于不同的影像數(shù)據(jù)。因此，有人提出了自適應(yīng)閾值的方法，其中以大津閾值(OTSU)方法[6]最為常見，但是OTSU閾值法只考慮到水體與地物之間光譜差異性，對于光譜變異的水體不能起到很好的區(qū)分效果。Shahriari等人[7]提出的一種基于分形技術(shù)的面積分形閾值算法，由于考慮了空間光譜特性，在影像分割中起到了很好的效果。

本文針對礦區(qū)復(fù)雜水體的提取，提出了改進(jìn)的TSUWI(MTSUWI)水體提取方法。以山西新柳礦區(qū)為研究對象，首先分別獲取UWI和USI指數(shù)；然后，創(chuàng)新性改進(jìn)了面積分形閾值算法，并分別對UWI與USI指數(shù)進(jìn)行閾值確定與水體提?。蛔詈?，將兩種水體提取結(jié)果進(jìn)行合并，通過提取交集部分的水體實現(xiàn)礦區(qū)復(fù)雜環(huán)境下的水體信息的高精度自動提取。

1 研究區(qū)概況與數(shù)據(jù)源介紹

實驗區(qū)為山西新柳礦區(qū)，位于山西省呂梁市孝義市西南陽泉曲鎮(zhèn)境內(nèi)。礦區(qū)內(nèi)地勢復(fù)雜，屬于山地地形。年降水量在450mm至550mm之間，四季降水量不均衡，夏季降水量集中，易于形成降水蓄積。

數(shù)據(jù)源為2016年8月3日獲取的GF2數(shù)據(jù)，GF2數(shù)據(jù)處于可見光和近紅波段范圍，包括四個多光譜波段和一個全色波段，多光譜波段分辨率為4m，全色波段分辨率為1m。

使用ENVI5.3對數(shù)據(jù)進(jìn)行預(yù)處理，由于礦區(qū)處于山地地形，地勢復(fù)雜，所以先結(jié)合DEM數(shù)據(jù)對多光譜與全色波段進(jìn)行正射校正。然后采用FLASSH模型對數(shù)據(jù)進(jìn)行大氣校正，消除大氣的影響。最后，采用Gram-Schmidt Pan Sharpening模型進(jìn)行多光譜數(shù)據(jù)與全色波段的融合。融合后數(shù)據(jù)的分辨率達(dá)到1m，且同時具備紅、綠、藍(lán)、近紅四個波段，顯著提高了數(shù)據(jù)的光譜空間特性。本文裁剪了2645*1897像素大小的區(qū)域作為實驗數(shù)據(jù)，實驗數(shù)據(jù)如圖1所示。

2 MTSUWI算法提取礦區(qū)水體

礦區(qū)復(fù)雜水體具有光譜變異、空間異質(zhì)性和易受陰影與裸土等地物干擾的特點(diǎn)，針對這些提取難點(diǎn)，本文提出MTSUWI方法，包括三個主要步驟：首先構(gòu)建UWI和USI特征圖，然后基于面積閾值法提取水體，最后將提取結(jié)果合并以求得最后的水體提取結(jié)果。

2.1 UWI與USI特征圖的構(gòu)建

UWI指數(shù)考慮了近紅、紅和綠波段的信息，USI指數(shù)則考慮了近紅、紅、綠和藍(lán)四個波段的信息，相比于傳統(tǒng)水體指數(shù)NDWI只考慮了近紅和綠波段的信息，UWI與USI更加充分地利用了光譜信息，能更穩(wěn)健地應(yīng)對光譜變異和其他地物的光譜混淆問題。參考Wu等人[5]提出的指數(shù)參數(shù)，UWI和USI構(gòu)建為式(1)和式(2)：

其中，NIR、R、G、B分別表示近紅、紅、綠、藍(lán)波段，α與β表示偏置，需要根據(jù)水體所處的不同環(huán)境來確定。式(1)和式(2)所得的特征圖UWI與USI，將水體部分進(jìn)行增強(qiáng)高亮，接下來需要確定閾值將高亮的水體部分提取出來。

2.2 基于面積分形閾值法的水體提取

面積分形閾值法最初應(yīng)用于影像分割領(lǐng)域，本文將其改進(jìn)以應(yīng)用于礦區(qū)水體閾值的自動確定。其模型可表示為：

其中，r>0表示特征值，這里為指數(shù)值，C>0表示比例常數(shù)，D>0為一般分維數(shù)。

其中，N(rj)表示相同的rj特征值對應(yīng)的個數(shù)或者頻率，rn表示最大特征值，n為特征值個數(shù)。根據(jù)公式(4)可以獲得數(shù)據(jù)集(N(r≥ri)，ri)，將數(shù)據(jù)集(N(r≥ri)，ri)代入式(3)并兩邊取對數(shù)，則問題簡化為一元線性回歸方程：

如圖2所示，ln(N(r≥ri))關(guān)于ln(ri)的散點(diǎn)(以小圓表示)會呈現(xiàn)出分段分布特點(diǎn)，通過最小二乘法可以按式(5)擬合出R1、R2、R3線段(以線段表示)，分別對應(yīng)D1、D2、D3區(qū)間。其中D3的分布區(qū)間表示高亮的區(qū)間，對應(yīng)的是純水體，D2分布區(qū)間表示非典型水體，比如污染嚴(yán)重、水較淺、水體面積較小等受噪聲污染的水體，而D1分布區(qū)間則表示非水體。因此，我們將閾值確定為T1，這樣可以很好地將純水體和光譜變異與空間異質(zhì)的水體提取出來。因此，基于UWI與USI指數(shù)，我們可以分別確定閾值Tuwi和Tusi，所以水體提取結(jié)果分別為Ruwi和Rusi：

圖2 ln(N(r≥ri))關(guān)于ln(ri)的分布圖與擬合線段

2.3 水體提取結(jié)果的合并

按照式(6)和式(7)所得結(jié)果能夠提取出水體，但是會存在噪聲的影響，因此，我們將提取的結(jié)果Ruwi與Rusi進(jìn)一步融合，以消除噪聲的影響。融合方式如式(8)：

通過式(8)求兩種提取水體的交集后，就可以獲得最終的水體提取結(jié)果。

3 實驗結(jié)果與分析

3.1 實驗結(jié)果對比

本文通過對比分析的方法來驗證本文方法的有效性，采用了經(jīng)典的水體提取算法NDWI[3]、UWI[5]、USI[5]、TSUWI[5]四種對比方法與本文的MTSUWI算法進(jìn)行對比。NDWI采用OTSU閾值法確定閾值，UWI、USI與TSUWI采用文獻(xiàn)[5]中的設(shè)置，本文算法MTSUWI中涉及的參數(shù)α[公式(1)]和參數(shù)β[公式(2)]分別設(shè)置為1.3和0.9。各個算法的實驗結(jié)果如圖3所示，其中參考提取結(jié)果是通過易康軟件結(jié)合Google地球影像進(jìn)行的人工標(biāo)注結(jié)果。

結(jié)果的量化評價指標(biāo)采用精確率和召回率以及綜合指標(biāo)[8]，分別計算為：

其中，TP表示正確提取水體的數(shù)量，F(xiàn)P表示錯誤提取水體的數(shù)量，F(xiàn)N表示將水體分為其他類別的數(shù)量。按公式(9)、公式(10)和公式(11)計算各個算法的評價指標(biāo)，列在表1中。

表1 不同方法的評價指標(biāo)

3.2 實驗結(jié)果分析

根據(jù)圖3的可視化結(jié)果，可以看到，NDWI的結(jié)果很差，主要原因是不能區(qū)分地物混淆的情況，因此表現(xiàn)在定量結(jié)果上，有較高的召回率，但是精確率卻非常低。UWI的提取結(jié)果定量上表現(xiàn)出較高的精確率，但召回率卻很低，體現(xiàn)在可視化結(jié)果上就是水體的漏分情況較為嚴(yán)重，主要原因是受到了礦區(qū)水體的空間異質(zhì)性影響。而USI的提取結(jié)果則表現(xiàn)出了較高的精確率和0.808的召回率，表現(xiàn)出了相對強(qiáng)地物混淆抑制能力，但一定程度上受到水體光譜變異和空間異質(zhì)性的影響。TSUWI的結(jié)果是UWI與USI的合并求交集的結(jié)果，所以，當(dāng)兩者的漏分情況嚴(yán)重時，TSUWI是不能取得較好精度的。而本文的MTSUWI算法，在可視化結(jié)果上，對于礦區(qū)水體的空間異質(zhì)性、地物混淆和光譜變異的問題，表現(xiàn)出了較好的提取效果，且綜合指標(biāo)優(yōu)于其他方法，雖然精確率與召回率不是最高的，但結(jié)合其他方法的定量和可視化結(jié)果，可以發(fā)現(xiàn)綜合指標(biāo)更具有評價方法優(yōu)劣的能力。因此，本文方法相比于對比方法，更加適用于礦區(qū)水體的提取。

圖3 不同方法的實驗結(jié)果

3.3 MTSUWI偏置參數(shù)的影響

在公式(1)和公式(2)中，設(shè)置了兩個偏置參數(shù)α和β，主要作用是為了平衡礦區(qū)水體環(huán)境因素的影響，從而更好地基于礦區(qū)環(huán)境獲取UWI與USI特征圖。為了探究偏置參數(shù)對結(jié)果的影響，本文設(shè)置了α和β的系列值，α∈[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對應(yīng)的綜合指標(biāo)的值，所得結(jié)果顯示在圖4中，可以得出當(dāng)α=1.3，β=0.9時，綜合指標(biāo)的值達(dá)到峰值。所以，針對礦區(qū)環(huán)境的偏置參數(shù)分別為α=1.3，β=0.9。

圖4 偏置參數(shù)α與β對MTSUWI方法綜合指標(biāo)的影響

4 結(jié)束語

本文提出了MTSUWI礦山水體提取算法，通過將面積分形閾值算法創(chuàng)新性地引入到UWI與USI的閾值提取過程，改進(jìn)了TSUWI水體提取算法。通過對新柳礦區(qū)的GF2影像進(jìn)行對比分析實驗，本文方法很好地解決了礦山水體的空間異質(zhì)性、地物混淆和光譜變異的問題，雖然精確率和召回率不是最高的，但是綜合指標(biāo)是最高水平，由于綜合指標(biāo)更具有算法評價能力，從而驗證了本文方法的有效性，表明本文方法可以應(yīng)用于礦區(qū)復(fù)雜水體的提取工作中。