尹德威,閆新珠
(山西省地質(zhì)測(cè)繪院,山西 運(yùn)城044000)
礦山排污作為礦產(chǎn)資源開發(fā)過(guò)程的必要環(huán)節(jié),會(huì)造成礦山生態(tài)環(huán)境的嚴(yán)重污染。開展對(duì)礦山水體的提取工作,對(duì)礦山水環(huán)境保護(hù)和預(yù)防尾礦庫(kù)事故具有重要的意義[1]。
水體提取的方法包括人工提取和遙感方法提取,人工提取方法效率低下,而遙感手段,由于具有周期短、大范圍、波譜豐富等特點(diǎn),成為了提取水體的主流方法[2]。目前,遙感水體提取的主要方法有:?jiǎn)尾ǘ伍撝捣ā⑺w指數(shù)法、面向?qū)ο蠓ǖ萚2]。單波段閾值法通過(guò)單一波段的光譜差異來(lái)區(qū)分水體與其他地物,對(duì)于光譜復(fù)雜的水體提取精度不高;一些研究人員使用基于對(duì)象的方法來(lái)檢測(cè)水體,并借助紋理等附加特征,但非常耗時(shí),且改進(jìn)效果有限;而水體指數(shù)法是通過(guò)特征波段間的比值運(yùn)算,以區(qū)分水體與其他地物。
水體指數(shù)法由于提取水體的高效性,是工業(yè)生產(chǎn)中的主流方法。典型的有McFEETERS[3]提出的歸一化差分水體指數(shù)法(NDWI),但NDWI不能很好地抑制建筑物的影響。因此,徐涵秋[4]提出了改進(jìn)的歸一化差分水體指數(shù)法(MNDWI),提高了對(duì)建筑的區(qū)分能力,但是由于陰影的影響,提取結(jié)果的精度受到限制??紤]到兩者的不足,Wu等人[5]提出將城市水體指數(shù)(UWI)與城市陰影指數(shù)(USI)結(jié)合,利用兩步水體指數(shù)法(TSUWI)來(lái)提取水體,取得了較好的效果。然而,TSUWI的閾值采用的是經(jīng)驗(yàn)閾值,難以適用于不同的影像數(shù)據(jù)。因此,有人提出了自適應(yīng)閾值的方法,其中以大津閾值(OTSU)方法[6]最為常見,但是OTSU閾值法只考慮到水體與地物之間光譜差異性,對(duì)于光譜變異的水體不能起到很好的區(qū)分效果。Shahriari等人[7]提出的一種基于分形技術(shù)的面積分形閾值算法,由于考慮了空間光譜特性,在影像分割中起到了很好的效果。
本文針對(duì)礦區(qū)復(fù)雜水體的提取,提出了改進(jìn)的TSUWI(MTSUWI)水體提取方法。以山西新柳礦區(qū)為研究對(duì)象,首先分別獲取UWI和USI指數(shù);然后,創(chuàng)新性改進(jìn)了面積分形閾值算法,并分別對(duì)UWI與USI指數(shù)進(jìn)行閾值確定與水體提??;最后,將兩種水體提取結(jié)果進(jìn)行合并,通過(guò)提取交集部分的水體實(shí)現(xiàn)礦區(qū)復(fù)雜環(huán)境下的水體信息的高精度自動(dòng)提取。
實(shí)驗(yàn)區(qū)為山西新柳礦區(qū),位于山西省呂梁市孝義市西南陽(yáng)泉曲鎮(zhèn)境內(nèi)。礦區(qū)內(nèi)地勢(shì)復(fù)雜,屬于山地地形。年降水量在450mm至550mm之間,四季降水量不均衡,夏季降水量集中,易于形成降水蓄積。
數(shù)據(jù)源為2016年8月3日獲取的GF2數(shù)據(jù),GF2數(shù)據(jù)處于可見光和近紅波段范圍,包括四個(gè)多光譜波段和一個(gè)全色波段,多光譜波段分辨率為4m,全色波段分辨率為1m。
使用ENVI5.3對(duì)數(shù)據(jù)進(jìn)行預(yù)處理,由于礦區(qū)處于山地地形,地勢(shì)復(fù)雜,所以先結(jié)合DEM數(shù)據(jù)對(duì)多光譜與全色波段進(jìn)行正射校正。然后采用FLASSH模型對(duì)數(shù)據(jù)進(jìn)行大氣校正,消除大氣的影響。最后,采用Gram-Schmidt Pan Sharpening模型進(jìn)行多光譜數(shù)據(jù)與全色波段的融合。融合后數(shù)據(jù)的分辨率達(dá)到1m,且同時(shí)具備紅、綠、藍(lán)、近紅四個(gè)波段,顯著提高了數(shù)據(jù)的光譜空間特性。本文裁剪了2645*1897像素大小的區(qū)域作為實(shí)驗(yàn)數(shù)據(jù),實(shí)驗(yàn)數(shù)據(jù)如圖1所示。
礦區(qū)復(fù)雜水體具有光譜變異、空間異質(zhì)性和易受陰影與裸土等地物干擾的特點(diǎn),針對(duì)這些提取難點(diǎn),本文提出MTSUWI方法,包括三個(gè)主要步驟:首先構(gòu)建UWI和USI特征圖,然后基于面積閾值法提取水體,最后將提取結(jié)果合并以求得最后的水體提取結(jié)果。
UWI指數(shù)考慮了近紅、紅和綠波段的信息,USI指數(shù)則考慮了近紅、紅、綠和藍(lán)四個(gè)波段的信息,相比于傳統(tǒng)水體指數(shù)NDWI只考慮了近紅和綠波段的信息,UWI與USI更加充分地利用了光譜信息,能更穩(wěn)健地應(yīng)對(duì)光譜變異和其他地物的光譜混淆問題。參考Wu等人[5]提出的指數(shù)參數(shù),UWI和USI構(gòu)建為式(1)和式(2):
其中,NIR、R、G、B分別表示近紅、紅、綠、藍(lán)波段,α與β表示偏置,需要根據(jù)水體所處的不同環(huán)境來(lái)確定。式(1)和式(2)所得的特征圖UWI與USI,將水體部分進(jìn)行增強(qiáng)高亮,接下來(lái)需要確定閾值將高亮的水體部分提取出來(lái)。
面積分形閾值法最初應(yīng)用于影像分割領(lǐng)域,本文將其改進(jìn)以應(yīng)用于礦區(qū)水體閾值的自動(dòng)確定。其模型可表示為:
其中,r>0表示特征值,這里為指數(shù)值,C>0表示比例常數(shù),D>0為一般分維數(shù)。
其中,N(rj)表示相同的rj特征值對(duì)應(yīng)的個(gè)數(shù)或者頻率,rn表示最大特征值,n為特征值個(gè)數(shù)。根據(jù)公式(4)可以獲得數(shù)據(jù)集(N(r≥ri),ri),將數(shù)據(jù)集(N(r≥ri),ri)代入式(3)并兩邊取對(duì)數(shù),則問題簡(jiǎn)化為一元線性回歸方程:
如圖2所示,ln(N(r≥ri))關(guān)于ln(ri)的散點(diǎn)(以小圓表示)會(huì)呈現(xiàn)出分段分布特點(diǎn),通過(guò)最小二乘法可以按式(5)擬合出R1、R2、R3線段(以線段表示),分別對(duì)應(yīng)D1、D2、D3區(qū)間。其中D3的分布區(qū)間表示高亮的區(qū)間,對(duì)應(yīng)的是純水體,D2分布區(qū)間表示非典型水體,比如污染嚴(yán)重、水較淺、水體面積較小等受噪聲污染的水體,而D1分布區(qū)間則表示非水體。因此,我們將閾值確定為T1,這樣可以很好地將純水體和光譜變異與空間異質(zhì)的水體提取出來(lái)。因此,基于UWI與USI指數(shù),我們可以分別確定閾值Tuwi和Tusi,所以水體提取結(jié)果分別為Ruwi和Rusi:
圖2 ln(N(r≥ri))關(guān)于ln(ri)的分布圖與擬合線段
按照式(6)和式(7)所得結(jié)果能夠提取出水體,但是會(huì)存在噪聲的影響,因此,我們將提取的結(jié)果Ruwi與Rusi進(jìn)一步融合,以消除噪聲的影響。融合方式如式(8):
通過(guò)式(8)求兩種提取水體的交集后,就可以獲得最終的水體提取結(jié)果。
本文通過(guò)對(duì)比分析的方法來(lái)驗(yàn)證本文方法的有效性,采用了經(jīng)典的水體提取算法NDWI[3]、UWI[5]、USI[5]、TSUWI[5]四種對(duì)比方法與本文的MTSUWI算法進(jìn)行對(duì)比。NDWI采用OTSU閾值法確定閾值,UWI、USI與TSUWI采用文獻(xiàn)[5]中的設(shè)置,本文算法MTSUWI中涉及的參數(shù)α[公式(1)]和參數(shù)β[公式(2)]分別設(shè)置為1.3和0.9。各個(gè)算法的實(shí)驗(yàn)結(jié)果如圖3所示,其中參考提取結(jié)果是通過(guò)易康軟件結(jié)合Google地球影像進(jìn)行的人工標(biāo)注結(jié)果。
結(jié)果的量化評(píng)價(jià)指標(biāo)采用精確率和召回率以及綜合指標(biāo)[8],分別計(jì)算為:
其中,TP表示正確提取水體的數(shù)量,F(xiàn)P表示錯(cuò)誤提取水體的數(shù)量,F(xiàn)N表示將水體分為其他類別的數(shù)量。按公式(9)、公式(10)和公式(11)計(jì)算各個(gè)算法的評(píng)價(jià)指標(biāo),列在表1中。
表1 不同方法的評(píng)價(jià)指標(biāo)
根據(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)出了相對(duì)強(qiáng)地物混淆抑制能力,但一定程度上受到水體光譜變異和空間異質(zhì)性的影響。TSUWI的結(jié)果是UWI與USI的合并求交集的結(jié)果,所以,當(dāng)兩者的漏分情況嚴(yán)重時(shí),TSUWI是不能取得較好精度的。而本文的MTSUWI算法,在可視化結(jié)果上,對(duì)于礦區(qū)水體的空間異質(zhì)性、地物混淆和光譜變異的問題,表現(xiàn)出了較好的提取效果,且綜合指標(biāo)優(yōu)于其他方法,雖然精確率與召回率不是最高的,但結(jié)合其他方法的定量和可視化結(jié)果,可以發(fā)現(xiàn)綜合指標(biāo)更具有評(píng)價(jià)方法優(yōu)劣的能力。因此,本文方法相比于對(duì)比方法,更加適用于礦區(qū)水體的提取。
圖3 不同方法的實(shí)驗(yàn)結(jié)果
在公式(1)和公式(2)中,設(shè)置了兩個(gè)偏置參數(shù)α和β,主要作用是為了平衡礦區(qū)水體環(huán)境因素的影響,從而更好地基于礦區(qū)環(huán)境獲取UWI與USI特征圖。為了探究偏置參數(shù)對(duì)結(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],并分別計(jì)算α和β不同組合下,MTSUWI對(duì)應(yīng)的綜合指標(biāo)的值,所得結(jié)果顯示在圖4中,可以得出當(dāng)α=1.3,β=0.9時(shí),綜合指標(biāo)的值達(dá)到峰值。所以,針對(duì)礦區(qū)環(huán)境的偏置參數(shù)分別為α=1.3,β=0.9。
圖4 偏置參數(shù)α與β對(duì)MTSUWI方法綜合指標(biāo)的影響
本文提出了MTSUWI礦山水體提取算法,通過(guò)將面積分形閾值算法創(chuàng)新性地引入到UWI與USI的閾值提取過(guò)程,改進(jìn)了TSUWI水體提取算法。通過(guò)對(duì)新柳礦區(qū)的GF2影像進(jìn)行對(duì)比分析實(shí)驗(yàn),本文方法很好地解決了礦山水體的空間異質(zhì)性、地物混淆和光譜變異的問題,雖然精確率和召回率不是最高的,但是綜合指標(biāo)是最高水平,由于綜合指標(biāo)更具有算法評(píng)價(jià)能力,從而驗(yàn)證了本文方法的有效性,表明本文方法可以應(yīng)用于礦區(qū)復(fù)雜水體的提取工作中。