周 煒，關(guān)洪軍，童 俊

(陸軍工程大學(xué)野戰(zhàn)工程學(xué)院，江蘇 南京 210007)

水體信息是地表覆蓋物的重要組成部分[1]，是工程保障的基礎(chǔ)信息。水體邊界信息的精確提取對(duì)于檢測(cè)江河的流量流速動(dòng)態(tài)變化、估算水庫及池塘的蓄水量具有重要意義，可以為渡河保障及給水保障提供信息支持。遙感影像的水體信息提取是遙感應(yīng)用技術(shù)方法研究的重要分支，在經(jīng)濟(jì)和軍事等領(lǐng)域具有重要的應(yīng)用價(jià)值。針對(duì)不同數(shù)據(jù)源，國(guó)內(nèi)外學(xué)者提出了不同的水體信息提取方法，如多波段譜間關(guān)系法[2]、水體指數(shù)法[3]、光譜知識(shí)法、面向?qū)ο蟮乃w信息提取方法[4]等。這些方法在進(jìn)行水體提取時(shí)，將像元作為水體提取的基本單位，對(duì)混合像元的判定比較粗糙，會(huì)導(dǎo)致水體邊界識(shí)別的精度損失。而混合像元分解技術(shù)可以對(duì)水體邊界的混合像元進(jìn)行定量分析，從而提高水體提取精度[5]。目前，國(guó)內(nèi)針對(duì)混合像元分解研究的數(shù)據(jù)源主要為多光譜影像，而利用高光譜影像進(jìn)行水體邊界混合像元分解的研究目前還未見報(bào)道，算法研究還有進(jìn)一步提升的空間。

本文針對(duì)水體邊界易與周邊地物產(chǎn)生混淆而界定困難的問題，提出一種基于高光譜混合像元分解的水體邊界信息提取方法，流程如圖1所示。

1 研究區(qū)域

研究區(qū)位于西藏阿里地區(qū)扎西崗鄉(xiāng)，北緯約32°30′，東經(jīng)約79°40′，屬于高寒干旱氣候，年平均氣溫-2～0℃，年平均降水量50～100 mm。該區(qū)域?yàn)閵A在喜馬拉雅山與岡底斯山脈間的山間谷地平原，山前冰雪融水形成大量冰水洪積扇，平原臺(tái)地多由濕地發(fā)育，植被以高寒草原為主，土壤主要為高山寒漠土。

2 遙感影像處理

2.1 影像波段篩選

本文研究使用國(guó)產(chǎn)某衛(wèi)星高光譜影像數(shù)據(jù)，攝影時(shí)段為2016年1月22日12時(shí)，影像分辨率為10 m，包含波段123個(gè)，畫幅大小700×500像元。鑒于國(guó)產(chǎn)高光譜影像各波段質(zhì)量參差不齊，根據(jù)圖像灰度范圍、圖像灰度標(biāo)準(zhǔn)差、圖像熵，對(duì)數(shù)據(jù)源進(jìn)行質(zhì)量評(píng)價(jià)和波段篩選。

灰度范圍用整個(gè)圖像灰度最大值與最小值的差值表示，記為D。圖像灰度標(biāo)準(zhǔn)差代表圖像灰度值的離散程度，記為S，計(jì)算公式如下

(1)

式中，M、N為影像的行列號(hào)；f(x,y)為目標(biāo)影像；Mean為目標(biāo)影像灰度均值。

圖像熵表示為圖像中灰度分布的聚集特征所包含的信息量，也描述了圖像的平均信息量[6]，計(jì)算公式為

(2)

式中，H表示圖像一維熵；i表示灰度值；Pi表示目標(biāo)圖像中灰度值為i的像元占總像元的比例。

篩選閾值為D≥600、S≥100、H≥5.5。根據(jù)篩選結(jié)果，將波段1～8、31～33、82～83、88～123判斷為不合格波段，使用74個(gè)質(zhì)量合格波段作為試驗(yàn)對(duì)象。

2.2 預(yù)處理

影像預(yù)處理主要包括輻射校正、幾何校正、濾波器降噪。輻射校正方法采用平場(chǎng)法，即選一塊光譜均一的高反射區(qū)取其平均值，然后對(duì)每一個(gè)像元的光譜值除以該平均值。幾何糾正通過圖像配準(zhǔn)實(shí)現(xiàn)，以高光譜數(shù)據(jù)作為標(biāo)準(zhǔn)圖像，對(duì)高分辨率數(shù)據(jù)進(jìn)行配準(zhǔn)處理。再通過重采樣，對(duì)圖像亮度值進(jìn)行插值計(jì)算，建立新的圖像矩陣。濾波降噪通過中值濾波器實(shí)現(xiàn)，濾波窗口大小為5×5。中值濾波能夠有效濾除椒鹽噪聲，盡量保護(hù)邊緣信息的完整[7]。

2.3 水體指數(shù)法

水體對(duì)小于0.6 μm波長(zhǎng)的光波吸收性和透射性相對(duì)弱，反射率相對(duì)較高。而在近紅外、短紅外波段，水體幾乎吸收全部的光波的入射能量，反射率很低。水體指數(shù)能夠代表水體在藍(lán)光、綠光、近紅外波段的典型特征，是水體信息提取的常用方法之一。

(1) 歸一化差值水體指數(shù)[8](NDWI)是水體信息提取的常用指數(shù)，其表達(dá)式為

(3)

式中，B2表示影像的綠波段；B4表示影像的近紅外波段。

(2) 陰影水體指數(shù)[9](SWI)能有效區(qū)分水體和陰影，彌補(bǔ)歸一化差值水體指數(shù)的缺陷，其表達(dá)式為

SWI=B1+B2-B4

(4)

式中，B1、B2、B4分別表示影像的藍(lán)波段、綠波段和近紅外波段。

(3) 改進(jìn)的陰影水體指數(shù)[10](MSWI)相比于SWI而言，水體于陰影的分離程度更高，區(qū)分更為明顯，其表達(dá)式為

(5)

式中，B1、B4分別表示影像的藍(lán)波段和近紅外波段。

用0.486 μm波段代表藍(lán)波段，波長(zhǎng)0.559 μm波段代表綠光波段，波長(zhǎng)0.837 μm波段代表近紅外波段，計(jì)算水體指數(shù)并構(gòu)建決策樹。設(shè)定閾值NDWI≥0.1，SWI≥50，MSWI≥0.3，提取水體要素。

2.4 水體邊界混合像元提取

水的反射光譜特征主要由水的物質(zhì)組成決定，也受水體物理化學(xué)狀態(tài)的影響。通常懸浮泥沙含量、水深及水體葉綠素含量會(huì)對(duì)水體的光譜特征造成一定的破壞。在使用多光譜影像進(jìn)行水體提取時(shí)，這些因素會(huì)影響水體信息的判斷，導(dǎo)致精度降低。而借助高光譜信息優(yōu)勢(shì)，構(gòu)建水體邊界混合像元光譜擬合曲線，可以通過光譜曲線豐富的變化判斷水體邊界信息，排除懸浮泥沙、水深等因素的影響。對(duì)水體及水體邊界混合像元進(jìn)行采樣，擬合兩者亮度值隨波長(zhǎng)的變化趨勢(shì)，如圖2所示。

圖2中實(shí)線表示水體像元擬合曲線，虛線表示混合像元擬合曲線。由圖可知，水體邊界混合像元在藍(lán)和綠波段光譜特征與水體像元大致相同，在紅及近紅外波段輻射亮度變化趨勢(shì)與水體相近，但亮度值遠(yuǎn)遠(yuǎn)高于水體像元。采用多波段譜間關(guān)系法，憑借混合像元與水體像元相似的譜間關(guān)系提取水體及其邊界混合像元，再通過混合像元與水體像元紅、近紅外波段的亮度差異，濾除水體像元，實(shí)現(xiàn)邊界混合像元的提取。

2.5 混合像元分解

混合像元是指包含多種地物混合信息的單個(gè)像元，其記錄的地面反射光譜信號(hào)是多種地物光譜信息綜合而成的復(fù)合信號(hào)[11]?；旌舷裨喾植加诘匦纹扑橹幒投喾N類型地物的交界地區(qū)。水體信息的邊界像元幾乎都屬于混合像元?；旌舷裨纸饩褪歉鶕?jù)混合像元的光譜曲線，計(jì)算混合像元中各端元所占的百分比，即端元豐度[12]。

2.5.1 端元選取

端元是指僅包含一種地物信息的單個(gè)像元，其對(duì)應(yīng)的光譜曲線可以代表該類地物的光譜曲線[13]。在遙感影像中同一種地物的光譜特征存在差異，難以做到直接提取純地物端元。因此，對(duì)水體邊界的不同地物進(jìn)行采樣，擬合地物端元光譜曲線。采集地物樣本點(diǎn)時(shí)，空間分布要盡量均勻，且需避開多種地物交界區(qū)。在擬合端元光譜曲線時(shí)，對(duì)地物樣本進(jìn)行標(biāo)準(zhǔn)差計(jì)算，排除標(biāo)準(zhǔn)差大的波段，選擇地物輻射亮度值穩(wěn)定的波段作為混合像元分解的依據(jù)。其中，水體在不同波段的亮度均值與標(biāo)準(zhǔn)差如圖3所示。

圖3中實(shí)線表示樣本像元亮度均值；虛線表示樣本亮度值標(biāo)準(zhǔn)差。鑒于水體輻射亮度在0.891 μm以外均存在較大波動(dòng)，不利于混合像元分解模型的計(jì)算。在混合像元分解時(shí)，僅將0.891 μm以內(nèi)的擬合結(jié)果作為各地物端元的光譜曲線。

2.5.2 線性光譜混合模型

線性光譜混合模型[14]是目前應(yīng)用最廣泛的混合像元分解模型，其表達(dá)式為

(6)

式中，x∈Rl為混合像元光譜；i為波段數(shù)；A∈Rl×P為端元矩陣；P為圖像中的端元數(shù)目；S∈RP為該混合像元的豐度向量；ε∈Rl為高斯隨機(jī)噪聲或模型誤差。

2.6 水體邊界提取算法

2.6.1 混合像元分解目標(biāo)判定及分解

以水體指數(shù)法提取結(jié)果為基礎(chǔ)，建立2像元的緩沖區(qū)。將緩沖區(qū)與混合像元的提取結(jié)果疊加，作為混合像元分解的目標(biāo)區(qū)域。而線性光譜混合模型可以看做是一個(gè)不等式組的求解，不等式組如下

(7)

式中，x1、x2、x3、x4、x5、x6為各端元所占百分比；p1,p2,…,pn為植被端元各波段亮度值；w1,w2,…,wn為水體端元各波段亮度值；g1,g2,…,gn為草場(chǎng)端元各波段亮度值；f1,f2,…,fn為裸地類型一端元各波段亮度值；s1,s2,…,sn為裸地類型二端元各波段亮度值；i1,i2,…,in為灘涂端元各波段亮度值；m1,m2,…,mn為混合像元各波段亮度值。鑒于不等式求解結(jié)果不完全符合線性光譜混合模型，設(shè)置5相對(duì)亮度值為緩沖。

2.6.2 混合像元分割

將高光譜影像中的每個(gè)像元拆分成25個(gè)小像元，則混合像元分割的步驟如下：

(1) 建立一個(gè)與試驗(yàn)影像等大的空矩陣A。根據(jù)混合像元分解的結(jié)果，將水體端元所占百分比在25%以上的像元換算成小像元個(gè)數(shù)N，記錄到矩陣A的相同位置。

(8)

式中，x1表示水體端元所占百分比；N表示拆分后純水體像元的個(gè)數(shù)。

(2) 若在矩陣A中，一個(gè)非零值像元四鄰域皆非零，則該像元為純水體像元，賦值25。若一非零像元四鄰域皆為零，則將該像元為非水體像元，賦值0。

(3) 判斷水體在像元中的聚集方位。像元四鄰域中，若僅有一個(gè)鄰域非零，則聚集方位為該非零鄰域方向；若有兩個(gè)鄰域非零且非零鄰域相鄰，如左鄰域和上鄰域，則聚集方位為值等于25的鄰域方向；若有兩個(gè)鄰域非零且非零鄰域相對(duì)，如左鄰域和右鄰域，則聚集方位為中央方向且連接非零鄰域；若有3個(gè)鄰域非零，則聚集方位位于零值鄰域的相反方向。

2.6.3 邊界提取

Laplacian算子對(duì)噪聲敏感，對(duì)平滑后的圖像邊緣有良好的檢測(cè)能力。以3×3窗口進(jìn)行運(yùn)算，獲得連續(xù)的水體邊界。

3 水體邊界提取效果分析

選擇影像中最復(fù)雜的河段作為對(duì)比區(qū)域，即區(qū)域中存在大量寬度小于10 m，表現(xiàn)為混合像元的河段。對(duì)比的方法分別為水體指數(shù)法、支持向量機(jī)法[15]及高光譜像元分解法，結(jié)果如圖4所示。表1為通過混淆矩陣對(duì)3種方法的精度評(píng)價(jià)。

表1 不同水體提取方法的精度驗(yàn)證

從表1看出，在水體多為混合像元的前提下，高光譜像元分解法在精度上遠(yuǎn)優(yōu)于水體指數(shù)法，略優(yōu)于SVM法。由圖4(a)可知，水體指數(shù)法對(duì)混合像元的區(qū)分度低，對(duì)于清澈淺水體像元難以確定合適的閾值進(jìn)行區(qū)分，不能有效識(shí)別。由圖4(b)、4(c)可知，SVM法與高光譜像元分解法提取效果相近，但高光譜像元分解法在支流連接及水體數(shù)學(xué)形態(tài)優(yōu)化方面更具優(yōu)勢(shì)。由圖4(d)可見，高光譜像元分解法能夠根據(jù)混合像元分解結(jié)果預(yù)測(cè)灘涂區(qū)域是否位于水位以下，從而對(duì)灘涂區(qū)域地表覆蓋屬性作出有效判斷。

4 結(jié) 論

本文以國(guó)產(chǎn)高光譜影像為數(shù)據(jù)源，提出了基于高光譜像元分解的水體邊界提取方法，實(shí)現(xiàn)了扎西崗周邊水體邊界信息的精確提取。結(jié)果表明，該方法水體提取總精度為93.86%，Kappa系數(shù)為0.87，高于水體指數(shù)法和SVM法，能為基于中分影像的水資源儲(chǔ)量評(píng)估以及灘涂水位的預(yù)測(cè)方面提供一定的科學(xué)參考。針對(duì)試驗(yàn)內(nèi)容，得出的結(jié)論如下：

(1) 憑借高光譜信息優(yōu)勢(shì)，可以通過擬合光譜曲線，排除懸浮泥沙、水深等因素對(duì)水體提取的干擾，提高水體邊界混合像元的解譯精度。

(2) 水體指數(shù)對(duì)混合像元的區(qū)分度低。高光譜像元分解法和SVM法均能有效區(qū)分混合像元的主成分類別，高光譜像元分解法在水體面積估算以及水體形態(tài)優(yōu)化方面要優(yōu)于SVM法。