李健鋒， 葉虎平， 張宗科， 魏顯虎

(1.陜西地建土地工程技術(shù)研究院有限責(zé)任公司，西安 710021；2.中國科學(xué)院遙感與數(shù)字地球研究所，北京 100094；3.中國科學(xué)院中國-斯里蘭卡水技術(shù)研究與示范聯(lián)合中心，北京 100085；4.陜西省土地工程建設(shè)集團(tuán)有限責(zé)任公司，西安 710075；5.中國科學(xué)院地理科學(xué)與資源研究所 資源與環(huán)境信息系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室，北京 100101)

0 引言

隨著航空航天技術(shù)的快速發(fā)展，遙感技術(shù)為資源調(diào)查、環(huán)境監(jiān)測、區(qū)域分析和全球宏觀研究等諸多領(lǐng)域提供了先進(jìn)手段[1-3]?？焖?、準(zhǔn)確地從遙感衛(wèi)星影像中提取水體信息對水資源調(diào)查、流域規(guī)劃、水資源的監(jiān)測與保護(hù)具有重要的意義[4-5]。

近年來，針對不同衛(wèi)星傳感器的影像提取水體，國內(nèi)外學(xué)者開展了廣泛的研究。Shih[6]基于Landsat MSS影像的近紅外波段，采用密度分割法，有效提取了平原地區(qū)的水體信息。周成虎等[7]研究發(fā)現(xiàn)，TM影像中水體具有獨(dú)特的譜間關(guān)系特征，即(TM2+TM3)>(TM4+TM5)，該方法特別適合山區(qū)水體的提取。Mcfeeters[8]根據(jù)TM2與TM4波段間的關(guān)系，提出了歸一化水體指數(shù)(normalized difference water index,NDWI)提取模型，該模型可以有效抑制土壤和植被信息，但難以抑制建筑物信息；徐涵秋[9]在NDWI基礎(chǔ)上，提出改進(jìn)的歸一化水體指數(shù)(modified NDWI,MNDWI)，尤其在城市水體提取方面，MNDWI表現(xiàn)出了較高的精度。Zhang等[10]針對Landsat-8 OLI影像，提出了基于LBV變換的地表水體提取方法，并與NDWI、MNDWI和AWEI提取結(jié)果相比，具有更好的精度。張景奇等[11]利用TM影像經(jīng)K-T轉(zhuǎn)換后水體的綠度(greenness)分量和第四分量小于濕度(wetness)分量的特征，有效地對水體信息進(jìn)行了提取。Sentinel-2衛(wèi)星是歐盟委員會(huì)(European Commission,EC)和歐洲航天局(European Sapce Agency,ESA)共同倡議的全球環(huán)境與安全監(jiān)測系統(tǒng)“哥白尼計(jì)劃”中的第二顆衛(wèi)星。Sentinel-2衛(wèi)星的主要任務(wù)是對全球陸地表面進(jìn)行高分辨率多光譜成像，分為2A和2B 2顆衛(wèi)星。相比于Landsat-8衛(wèi)星影像，Sentinel-2影像擁有更高的空間分辨率和時(shí)間分辨率[12]。雖然國內(nèi)外學(xué)者針對不同傳感器提出了許多水體提取模型，但這些模型在復(fù)雜的水體環(huán)境下，提取精度普遍不高，無法有效地消除云陰影、山體陰影及水田泥地對水體提取精度的影響，存在一定的誤分現(xiàn)象?！?1世紀(jì)海上絲綢之路”沿線國家斯里蘭卡位于南亞次大陸南端，受印度洋季風(fēng)系統(tǒng)的影響，導(dǎo)致強(qiáng)降雨一年中在該地區(qū)系統(tǒng)性地遷移，存在明顯季節(jié)性缺水，水資源時(shí)空分布十分不均衡[13]。地表水體面積精確提取是評估水資源量的一個(gè)重要手段。斯里蘭卡屬于熱帶地區(qū)，近中午前后常有各種云覆蓋，致使對地觀測衛(wèi)星過境影像常年受云霧覆蓋及云陰影影響，并且斯里蘭卡有很多種植水稻的水田，收割完水稻的泥地與水體的光譜特征較為接近，也增加了遙感圖像水體提取的難度。

本文選取斯里蘭卡中東部為研究區(qū)，推導(dǎo)了用于Sentinel-2影像的LBV變換方程，在分析了水體與植被、建筑、陰影等典型地物經(jīng)LBV和K-T變換后特征的基礎(chǔ)上，提出了基于LBV和K-T變換的水體提取模型，實(shí)現(xiàn)了復(fù)雜水環(huán)境下水體的精確提取，并與水體指數(shù)法(NDWI、MNDWI)的提取結(jié)果進(jìn)行了對比分析。

1 數(shù)據(jù)與方法

1.1 研究區(qū)數(shù)據(jù)

本文以斯里蘭卡中東部為研究區(qū)，研究區(qū)東部沿海，西南部為高山地區(qū)，區(qū)內(nèi)環(huán)境復(fù)雜且影響水體提取精度的因素較多，包含湖泊、水庫、坑塘、瀉湖、河流等多種類型的水體。研究區(qū)Sentinel-2假彩色(8、4、3波段)合成圖像如圖1所示。實(shí)驗(yàn)選取覆蓋研究區(qū)2018年9月的1景Sentinel-2B影像，含云量為33.92%，主要集中在圖像的東南部分。利用SNAP軟件和Sen2Cor插件對影像進(jìn)行輻射定標(biāo)、大氣校正以及裁剪處理。

圖1 研究區(qū)Sentinel-2假彩色合成圖像(8、4、3波段)

1.2 研究方法

1) Sentinel-2影像的LBV變換。LBV變換是由Zeng[14]提出的一種通用的光學(xué)遙感數(shù)據(jù)變換方法。LBV從地物的物理意義出發(fā)，分析不同典型地物在多光譜影像上輻射值的變化，總結(jié)出地物的3個(gè)主要輻射特性：地物總輻射水平L(general radiance level)、地物可見光-紅外光輻射平衡B(visible-infrared radiation balance)和地物輻射隨波段變化的向量-方向和速度V(band radiance variation vector-direction and speed)。LBV變換的輸入波段組合和方程隨衛(wèi)星的不同而不同。

根據(jù)Sentinel-2A影像的波段特征并參考文獻(xiàn)[10]，選取3(green)、4(red)、8(NIR)和11(SWIR)波段作為輸入波段組合，其LBV變換方程可通過線性回歸和二次回歸得到，其過程如式(1)～式(3)所示。

(1)

(2)

(3)

地物可見光-紅外光輻射平衡B的變換方程通過線性回歸推導(dǎo)，采用二次回歸進(jìn)行推導(dǎo)L和V的變換方程。根據(jù)公式(2)和λi可得到式(4)～式(7)。

R3=a+0.559 8b+0.559 82c

(4)

R4=a+0.664 6b+0.664 62c

(5)

R8=a+0.832 8b+0.832 82c

(6)

R11=a+1.613 7b+1.613 72c

(7)

確定3個(gè)方程求解二次回歸中3個(gè)未知的系數(shù),如式(8)～式(10)所示。

(8)

(9)

(10)

結(jié)合式(4)～式(10)得到二次回歸方程系數(shù),如式(11)～式(13)所示。

a=4.408R3+0.394 9R4-4.367 6R8+0.564 7R11

(11)

b=-8.437 1R3+0.049R4+10.057R8-1.669 2R11

(12)

c=3.538 3R3-0.187 4R4-4.552 2R8+1.201 4R11

(13)

(14)

線性回歸的推導(dǎo)類似與二次回歸，最終得到線性回歸方程系數(shù)如式(15)～式(16)所示。

m=0.730 3R3+0.589 7R4+0.364R8-0.684R11

(15)

n=-0.523 4R3-0.370 2R4-0.124 2R8+1.669 2R11

(16)

地物可見光-紅外光輻射平衡B定義為線性回歸的斜率負(fù)值[13]，B的表達(dá)式如式(17)所示。

B=-n=0.523 4R3+0.370 2R4+
0.124 2R8-1.669 2R11

(17)

地物輻射隨波段變化的向量-方向和速度V定義為二次回歸的組合殘差[13]，V的計(jì)算過程如式(18)～式(19)所示。

(18)

V=-v3+v4-v8+v11=0.741 1R3-
1.338 7R4+0.626 7R8-0.029 1R11

(19)

式(14)、式(17)和式(19)為Sentinel-2A影像的LBV變換方程，同理可推導(dǎo)出Sentinel-2B影像的LBV變換方程,如式(20)所示。

(20)

2) Sentinel-2影像的K-T變換。K-T變換(又稱穗帽變換)是一種特殊的主成分分析，由Kauth和Thomas于1976年提出，通過旋轉(zhuǎn)坐標(biāo)空間使坐標(biāo)軸指向與特征密切相關(guān)的方向，特別與植物生長過程和土壤密切相關(guān)[15]。K-T變換有助于對作物性狀進(jìn)行解釋和分析，具有重要的現(xiàn)實(shí)意義。K-T變換的第一分量亮度(brightness)表示地物總的輻射能量水平，第二分量綠度(greenness)反映植被的生長狀況，第三分量濕度(wetness)反映了地物的濕度信息，對土壤濕度信息最為敏感，是較好的水體信息識別的特征波段。Nedkov[15]在2017年推導(dǎo)了Sentinel-2影像K-T變換的亮度、綠度和濕潤度指標(biāo)的變換系數(shù)。

2 結(jié)果與討論

2.1 典型地物光譜特征分析

在研究區(qū)內(nèi)針對水體、植被、建筑、陰影和水田泥地5種地物分別選取50個(gè)樣本，統(tǒng)計(jì)其在Sentinel-2B影像中2(blue)、3(green)、4(red)、8(NIR)和11(SWIR)波段的反射率均值，得到典型地物的光譜特征曲線如圖2所示。從圖2可以看出，在可見光波段，水體的反射率在藍(lán)色和綠色波段較高，紅色波段較低。在近紅外和短波紅外波段，水體反射率較低，幾乎所有入射能量都被吸收。水體反射率呈遞減的趨勢，即：藍(lán)波段>綠波段>紅波段>近紅外波段>短波紅外。影像中陰影的光譜特征與水體極為相似，二者的反射率在紅波段相交，在其他波段較為接近，區(qū)分度較低。研究區(qū)內(nèi)主要存在山體陰影和云陰影，影響水體提取的精度。水田泥地含水量高，整體的反射率較低，在遙感影像中呈暗色調(diào)，會(huì)對水體的提取結(jié)果產(chǎn)生一定的影響。建筑和植被的光譜響應(yīng)曲線與水體差別較大，易于區(qū)分。

圖2 典型地物光譜響應(yīng)曲線

2.2 基于LBV和K-T變換的水體提取模型

針對處理后的實(shí)驗(yàn)區(qū)Sentinel-2B影像，根據(jù)推導(dǎo)出的LBV變換方程和K-T變換系數(shù)[15]，基于IDL 8.5平臺編譯Sentinel-2B影像的LBV和K-T變換算法，得到研究區(qū)的LBV和K-T變換影像，分別統(tǒng)計(jì)典型地物經(jīng)LBV和K-T變換后的波譜采樣均值,如圖3所示。

圖3 典型地物L(fēng)BV和K-T變換后波譜均值直方圖

由圖3(a)可以看出，水體經(jīng)LBV變換后L波段波譜采樣均值與陰影、水田泥地的波譜采樣均值較為接近，植被的波譜采樣均值最低；水體B波段波譜采樣均值最高，為13.09，其他4種典型地物B波段的波譜采樣均值均為負(fù)值；5種典型地物V波段波譜采樣均值整體較為接近，其中植被的波譜采樣均值最高。分析5種典型地物經(jīng)LBV變換后各波段之間的關(guān)系發(fā)現(xiàn)，只有水體經(jīng)LBV變換后的波譜采樣均值呈現(xiàn)出B>V的特征，植被、建筑、陰影和水田泥地等地物均呈現(xiàn)出B波段V的灰度值可進(jìn)行水體的粗提取，提取結(jié)果如圖4(b)所示。從圖中可以看出，影像無云區(qū)域水體提取結(jié)果較好，能夠準(zhǔn)確地識別水體的邊界，影像有云區(qū)域水體提取結(jié)果包含了部分云和云陰影，并且還存在水體和水田泥地的混分現(xiàn)象。

由圖3(b)可以看出，水體經(jīng)K-T變換后brightness(亮度)波段波譜采樣均值最小，為39.29，與其他典型地物相差較大；陰影、水田泥地和水體brightness(亮度)波段波譜采樣均值較為接近，均趨于0；5種典型地物wetness(濕度)波段的波譜采樣均值差異明顯，水體的最大，為19.54。因此，根據(jù)水體在brightness和wetness波段與其他典型地物差異明顯的特點(diǎn)，通過對K-T變換后的圖像進(jìn)行多尺度分割確定brightness波段的閾值為75，wetness波段的閾值為11，進(jìn)行水體信息的提取，提取結(jié)果如圖4(c)所示。相比于利用LBV變換后波段B的灰度值>波段V的灰度值的水體提取結(jié)果，brightness和wetness波段閾值法水體提取的結(jié)果有效地消除了云、陰影以及水田泥地的影響，但存在水體與部分植被的混分現(xiàn)象，特別是水體附近的植被。

綜上所述，K-T變換brightness和wetness波段閾值法的水體提取結(jié)果能有效地消除云、陰影以及水田泥地的影響，利用LBV變換波段B的灰度值>波段V的灰度值能有效地區(qū)分水體與植被。因此，結(jié)合2種方法，建立如式(21)所示的水體提取模型。

(21)

式中：G(x，y)表示圖像水體提取后的二值化影像；1表示為水體；0表示為非水體；B和V分別表示圖像經(jīng)LBV變換后的B波段和V波段；Brig.和Wet.分別表示圖像經(jīng)K-T變換后的brightness和wetness波段；m和n分別為brightness和wetness波段的閾值。

2.3 研究區(qū)水體提取模型比較

基于IDL 8.5平臺，對預(yù)處理后的研究區(qū)Sentinel-2B影像分別進(jìn)行LBV轉(zhuǎn)換和K-T轉(zhuǎn)換，并合成為6波段圖像(包含L、B、V波段和brightness、greenness、wetness波段)，利用本文提出的基于LBV和K-T變換的水體提取模型進(jìn)行水體的提取。通過對K-T變換后的圖像進(jìn)行多尺度分割確定brightness波段閾值設(shè)為75，wetness波段閾值設(shè)為11，進(jìn)行水體信息的提取，提取結(jié)果如圖4(d)所示。為了進(jìn)一步對比分析基于LBV和K-T變換的水體提取模型提取精度，采用水體指數(shù)法(NDWI、MNDWI)進(jìn)行水體提取，通過大津法(Otsu)[16]確定NDWI閾值為-0.002 7，MNDWI閾值為0.008 4，提取結(jié)果如圖4(e)和圖4(f)所示。通過提取結(jié)果與影像疊加對比分析發(fā)現(xiàn)，在研究區(qū)的東南有云部分，NDWI和MNDWI水體提取模型提取精度較差，存在嚴(yán)重的誤分現(xiàn)象。NDWI模型的提取結(jié)果主要存在云陰影和水田泥地誤分為水體的現(xiàn)象；MNDWI模型的提取結(jié)果主要存在云和水田泥地誤分為水體的現(xiàn)象。在研究區(qū)的西南高山區(qū)，NDWI模型將一部分山體陰影誤分為水體。在研究區(qū)的無云平原區(qū)，NDWI和MNDWI模型提取精度較高，能比較準(zhǔn)確地提取出水體的邊界。基于LBV和K-T變換的水體提取模型在整個(gè)研究區(qū)內(nèi)都表現(xiàn)出了較高的精度，有效地消除了云、云陰影、山體陰影和水田泥地的影響。

圖4 不同模型水體提取結(jié)果圖

為了進(jìn)一步對比分析3種模型的提取結(jié)果，以研究區(qū)水體的目視解譯結(jié)果為基準(zhǔn)，在研究區(qū)內(nèi)選取200個(gè)樣本，包括100個(gè)水體樣本和100個(gè)非水體樣本，通過計(jì)算混淆矩陣[17]，分別從漏提率、誤提率和總體精度3個(gè)方面反映水體的提取精度，其統(tǒng)計(jì)結(jié)果如表1所示。從表1中可以看出，NDWI和MNDWI水體提取模型的誤提率較高，存在明顯的非水體提取為水體的現(xiàn)象，總體精度較低，均低于90%，同時(shí)也存在一定的漏提現(xiàn)象?；贚BV和K-T變換的水體提取模型總體精度最高，達(dá)到了98.13%，存在少量的漏提和誤提現(xiàn)象。

表1 3種水體提取方法的精度對比

3 結(jié)束語

本文推導(dǎo)了用于Sentinel-2影像的LBV變換方程，以水環(huán)境復(fù)雜的斯里蘭卡中東部為研究區(qū)，在分析了水體與植被、建筑、陰影等典型地物相比經(jīng)LBV和K-T變換后特征的基礎(chǔ)上，提出了基于LBV和K-T變換的水體提取模型，從目視判讀與定量分析2個(gè)角度評價(jià)了本文模型與NDWI、MNDWI模型提取水體的精度，發(fā)現(xiàn)基于LBV和K-T變換的水體提取模型總體精度最高，達(dá)到了98.13%，NDWI和MNDWI 模型的總體精度較低，低于90%，存在較為明顯的誤提現(xiàn)象。NDWI模型將一部分云陰影、山體陰影和水田泥地誤分為水體，MNDWI模型將一部分云和水田泥地誤分為水體，同時(shí)2個(gè)模型還存在一定的漏提現(xiàn)象。總地來說，相比于NDWI、MNDWI模型，基于LBV和K-T變換的水體提取模型實(shí)現(xiàn)了復(fù)雜水環(huán)境下水體的精確提取，有效地消除了云、云陰影、山體陰影和水田泥地的影響。

本文提出的基于LBV和K-T變換的水體提取模型，實(shí)現(xiàn)了斯里蘭卡中東部地區(qū)水體的精確提取，為在復(fù)雜水環(huán)境下水體的精確提取提供了有效的解決途徑，對水資源調(diào)查、監(jiān)測與保護(hù)有重要的現(xiàn)實(shí)意義。但模型閾值的選取需要對圖像進(jìn)行多尺度分割后人工確定與干預(yù)，尚未實(shí)現(xiàn)自動(dòng)化水體的精確提取，在后續(xù)的研究中應(yīng)聚焦于閾值的自動(dòng)確定，實(shí)現(xiàn)復(fù)雜水環(huán)境下水體的自動(dòng)化精確提取。