熊龍海,何穎清,劉茉默,李 俊

(1. 水利部珠江河口治理與保護重點實驗室,廣東 廣州 510611; 2.珠江水利委員會珠江水利科學(xué)研究院,廣東 廣州 510611)

準確獲取地表水水體范圍及面積對水資源保護和合理開發(fā)利用、海岸線變化監(jiān)測、水環(huán)境監(jiān)測、水資源調(diào)查、濕地調(diào)查、洪澇災(zāi)害評估等方面具有重要意義[1]。衛(wèi)星遙感技術(shù)具有快速、大范圍、周期性成像等特點,已經(jīng)成為一種有效的地表水水體范圍提取技術(shù)。不同空間分辨率、時間分辨率和光譜分辨率的衛(wèi)星傳感器已經(jīng)被廣泛應(yīng)用于地表水水體范圍的提取[2-8]。在地表水遙感的研究中,長序列、高頻次、大范圍對水面率、蓄水量、生態(tài)流量等水資源要素進行監(jiān)測通常需要結(jié)合多源中高分辨率衛(wèi)星數(shù)據(jù)實現(xiàn)。但是在結(jié)合高分辨率遙感影像與Sentinel-2的研究中,往往因為空間分辨率差異大導(dǎo)致水體提取結(jié)果存在空間尺度效應(yīng)。因此,為解決空間尺度效應(yīng),研究面向Sentinel-2的亞像元級水體提取方法十分必要。

雖然目前有較多遙感水體提取方法,但絕大多數(shù)都是像元精度的,主要包括:①通過目視判讀進行數(shù)字化;②單波段閾值法[9];③多波段光譜水指數(shù)法[10-12];④監(jiān)督或非監(jiān)督分類算法[13];⑤圖像處理方法,如面向?qū)ο蠓治?object-oriented analysis, OOA)[14];⑥目標檢測算法,如基于約束能量最小化 (constrained energy minimization, CEM) 的目標檢測法[15];⑦上述方法的組合法[16-17]。上述方法將每個像元分為水體像元或陸地像元,即使是水陸混合像元也被分類為水體或陸地像元,無法提取水陸混合像元中水體豐度。然而,水陸混合像元普遍存在于衛(wèi)星遙感影像上,現(xiàn)有像元級水體提取方法均不能準確獲取混合像元中水體豐度,導(dǎo)致水體提取精度不足。

一些學(xué)者提出通過軟分類、機器學(xué)習(xí)和光譜混合分析的方法解決水體提取過程中水陸混合像元的問題。軟分類的方法并不計算混合像元中水體豐度,其水體提取結(jié)果還是像元級。機器學(xué)習(xí)方法需要大量的訓(xùn)練集,且依賴訓(xùn)練集的質(zhì)量,限制了該方法的適用性。光譜混合分析方法是一種有效的、不需要大量訓(xùn)練數(shù)據(jù)的計算混合像元中水體豐度的方法,文獻[18—19]基于Landsat 8 OLI數(shù)據(jù),利用光譜混合分析方法提出了亞像元級水體提取技術(shù),但該技術(shù)尚未在Sentinel-2數(shù)據(jù)上應(yīng)用,Sentinel-2數(shù)據(jù)比Landsat 8 OLI數(shù)據(jù)波段更多,直接將上述技術(shù)應(yīng)用在Sentinel-2數(shù)據(jù)上無法充分利用Sentinel-2數(shù)據(jù)的優(yōu)勢,難以取得最優(yōu)效果。目前尚缺乏有效的面向Sentinel-2數(shù)據(jù)的亞像元級水體提取方法。

因此,本文在利用植被紅邊水體指數(shù)(vegetation red edge based water index,RWI)[20]提取純水體像元的基礎(chǔ)上,結(jié)合動態(tài)選擇端元的多端元光譜混合分析方法構(gòu)建面向Sentinel-2影像的亞像元級水體提取方法(subpixel water extraction method for Sentinel-2,SWES),以提升水體提取精度,獲取更加精確的水域面積,將基于Sentinel-2數(shù)據(jù)的水體提取進行降尺度,解決多源中高分辨衛(wèi)星數(shù)據(jù)提取水體時的空間尺度效應(yīng)問題。

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

1.1 研究區(qū)概況

廣州市水系發(fā)達,包含多種水體類型,有不同的水陸混合像元,是亞像元級水體提取方法較為理想的研究區(qū)。本文在廣州市分別選擇了3個試驗區(qū),每個試驗區(qū)都包含不同的水體類型、地表覆蓋類型和不同數(shù)目的混合像元。這3個試驗區(qū)包括河流、水庫和坑塘養(yǎng)殖區(qū),試驗區(qū)內(nèi)水域的水深、渾濁度及污染程度各異。選擇這3個試驗區(qū)是因為它們分別代表了不同光譜特征的水體。圖1為試驗區(qū)的地理位置示意圖及假彩色合成圖。其中,試驗區(qū)1為水口水庫,水體類型以清潔水體為主;試驗區(qū)2為珠江干流,擁有城市背景下的復(fù)雜水生環(huán)境,主要污染來自于生活污水和工業(yè)廢水;試驗區(qū)3為坑塘養(yǎng)殖區(qū),主要為水深較淺的?；~塘。每個試驗區(qū)都存在著不同數(shù)量的水陸混合像元,其中,珠江干流混合像元最少,水口水庫有較多庫灣包含較多混合像元,坑塘養(yǎng)殖區(qū)存在大量面積不一的坑塘,水陸混合像元數(shù)量最多。

圖1 研究區(qū)位置

1.2 數(shù)據(jù)及處理

水體提取采用的Sentinel-2數(shù)據(jù)為從Google Earth Engine平臺上獲取的地表反射率產(chǎn)品。產(chǎn)品級別為L2A,該產(chǎn)品已經(jīng)由歐空局預(yù)先進行了輻射定標、大氣校正等預(yù)處理。覆蓋研究區(qū)的Sentinel-2影像成像日期為2019年12月1日,并將所有空間分辨率不為10 m的波段重采樣至10 m分辨率。

選取與Sentinel-2影像同一天成像的Google EarthTM高空間分辨率(<1 m)影像作為真實水體范圍參考數(shù)據(jù)。統(tǒng)一兩種影像的坐標系,通過選擇道路交叉口、低層建筑角點等控制點對Sentinel-2影像與Google EarthTM高空間分辨率影像進行配準。將Sentinel-2數(shù)據(jù)疊置在真實的水體范圍數(shù)據(jù)上,從而得到真實的純水像元及真實的水體豐度圖,并用于驗證本文提出的SWES。

2 方法構(gòu)建

衛(wèi)星影像上的像元可分為3類:水體像元、陸地像元、水陸混合像元。對于地表水提取只需要提取水體像元及水陸混合像元,亞像元級地表水遙感提取結(jié)果由水體像元及水陸混合像元中水體豐度組成。本文構(gòu)建的亞像元級水體遙感提取方法(SWES)共包括3個主要的步驟:①提取純水像元;②提取水陸混合像元;③結(jié)合空間信息的多端元光譜混合分析提取水體豐度。

2.1 提取純水像元

采用RWI提取純水像元,RWI為面向Sentinel-2影像的細小水體提取指數(shù),在一定程度上能夠消除水陸邊界混合像元的影響,在提取純水像元上效果較好,該指數(shù)同樣適用于大范圍水體提取。RWI使用大氣校正后的遙感反射率作為輸入,其公式為

(1)

式中,Rrs(560)、Rrs(705)、Rrs(842)、Rrs(865)、Rrs(2190)分別為Sentinel-2第3、5、8、8A 和 12 波段的遙感反射率。

2.2 提取水陸混合像元

水陸混合像元的光譜非常復(fù)雜,常常表現(xiàn)出與其他地物相似的特征,因此單純地利用光譜特征提取水陸混合像元是不可行的。而假設(shè)影像中所有像元都是水陸混合像元也是不切實際的,因此本文提出了利用空間信息提取水陸混合像元的方法。水陸混合像元往往在水體和陸地的過渡帶上,即在水域的邊界上,因此可利用這一空間特征提取水陸混合像元。具體步驟為在純水像元的基礎(chǔ)上搜尋8鄰域所有的陸地像元,如圖2所示,即在位于(i,j) 處純水像元的8鄰域內(nèi),所有的陸地像元被認為是水陸混合像元。

圖2 低分辨率影像上水陸邊界

2.3 結(jié)合空間信息的多端元光譜混合分析

2.3.1 光譜混合分析模型

光譜混合分析能夠用來獲取水陸混合像元中水體的豐度。常用的線性光譜混合模型假設(shè)混合像元的光譜是由構(gòu)成該像元的地物端元的反射率以其所占像元面積比例為權(quán)重系數(shù)的線性組合,即

(2)

式中,Riλ為第i個混合像元的反射率;N為端元的數(shù)量;fji為第i個混合像元中第j個端元的豐度;Rjλ為第j個端元的反射率;eiλ為殘差。在全約束條件的混合像元分解模型中,端元的豐度通常是非負而且滿足和為1的約束,即

(3)

對于第i個混合像元,地物的豐度可通過最小二乘法得到,也即是求取使各個波段上殘差的均方根誤差(root mean squared error,RMSE)最小的解?？杀硎緸?/p>

(4)

式中,T為波段數(shù);eik為第i像元上第k個波段的殘差。

2.3.2 結(jié)合空間信息的多端元動態(tài)選擇

由于地表水體的光譜隨水質(zhì)變化而變化,淺水情況下還受水深及底質(zhì)影響。一般情況下,鄰近的水體像元的水質(zhì)、水深、底質(zhì)條件最為接近,因此,本文使用結(jié)合空間信息的水體端元動態(tài)提取方案。對于一個水陸混合像元,所有3×3的窗口8鄰域內(nèi)的純水像元都被認為是這個混合像元的候選水體端元。如圖2(d)所示,位于(k,l)處的水陸混合像元,其8鄰域內(nèi)所有純水像元均被認為是該混合像元的候選水體端元。

3類陸地端元(植被、土壤、不透水面)通過利用純凈像元指數(shù)(pixel purity index,PPI)建立端元光譜庫,并使用Count-based Endmember Selection(CoB),Endmember Average RMSE(EAR)和Minimum Average Spectral Angle(MASA)優(yōu)化建立的端元光譜庫。

2.3.3 端元組合

解混過程中考慮了多種陸地端元與水體端元的組合,主要分為2端元、3端元、4端元和5端元組合,具體的組合方式見表1。選擇陰影端元的目的是消除光照的變化。根據(jù)前人的研究[21],分解結(jié)果同時滿足下列條件的混合像元才被分解:①所有端元組分的豐度在-0.05到1.05之間;②陰影的豐度低于0.8;③各個波段上殘差的均方根誤差RMSE小于0.025。如果分解的結(jié)果無法同時滿足上述條件則該混合像元不分解,并將水體豐度設(shè)置為0。對于一個水陸混合像元,可能有多個端元組合分解的結(jié)果滿足上述條件,選擇其中RMSE最小的組合作為最佳的端元組合,其分解結(jié)果作為最優(yōu)結(jié)果。

表1 端元組合類型

3 結(jié)果與分析

3.1 目視分析

由于缺少專門面向Sentinel-2的亞像元級水體提取方法,本文選擇可移植到Sentinel-2的自動亞像元水體制圖方法(automatic subpixel water mapping method,ASWM)作為對比。在ASWM中,首先利用NDWI提取純水像元和純陸地像元,然后根據(jù)閾值提取水陸混合像元,最后利用固定端元組合的混合像元分解算法求解水陸混合像元中水體的豐度,純水像元和水陸混合像元中水體豐度共同構(gòu)成了最終的水體提取結(jié)果。在這一過程中純水像元的豐度被設(shè)置為1,純陸地像元的豐度被設(shè)置為0。

圖3—圖5分別為利用SWES和ASWM在3個研究區(qū)進行亞像元級水體提取的結(jié)果。從結(jié)果中可以看出:①整體而言,SWES在3個研究區(qū)的水體提取結(jié)果較好,明顯優(yōu)于ASWM。②對于水陸混合像元,SWES利用空間關(guān)系能夠準確提取出來,并且獲得相應(yīng)的水體豐度,在結(jié)果圖上可以明顯地看出由水體到陸地的過渡,通過與真實水體豐度對比,可以看出SWES獲取的水陸混合像元水體豐度也能與真實值保持較好的一致性。而ASWM根據(jù)光譜信息提取水陸混合像元帶來較大誤差,將較多陸地像元誤提為水陸混合像元,而且難以準確提取3個試驗區(qū)純水體邊界處的水陸混合像元。③對于純水體像元,SWES能夠較準確提取出來,漏提和誤提較少,水體提取效果優(yōu)于ASWM,尤其在水陸混合較嚴重的區(qū)域,如試驗區(qū)3,SWES比ASWM取得更好的效果,在這些區(qū)域由于存在較多水陸混雜,水質(zhì)、水深條件不一,水體光譜變化強烈,ASWM采用NDWI提取純水體像元,很難取得較好效果。④通過圖4可以看出,建筑及其陰影對ASWM的影響比SWES的大,SWES幾乎沒有將建筑及其陰影誤提為水體,而ASWM將很多建筑陰影識別為水陸混合像元,也將圖5中的道路識別為水陸混合像元。⑤兩種方法在提取寬度小于3個像元的水體方面都存在一些誤差。

圖4 SWES和ASWM在試驗區(qū)2估算水體豐度的結(jié)果

圖5 SWES和ASWM在試驗區(qū)3估算水體豐度的結(jié)果

3.2 精度評價

RMSE是亞像元分類中廣泛使用的誤差指標,量化了地物豐度分解的相對誤差。本文選擇RMSE對建立的亞像元級水體提取模型SWES進行精度評價。計算公式為

(5)

式中,fi(ref)和fi(esti)分別表示真實的和計算的水體豐度圖中第i個像元的水體豐度;N為計算的水體豐度圖上所有像元的數(shù)目。

表2為SWES和ASWM在3個試驗區(qū)提取水體時的RMSE。整體上看,SWES在3個試驗區(qū)的平均RMSE為0.147,而ASWM的平均RMSE為0.175,SWES在所有試驗區(qū)都比ASWM的RMSE更低,特別是在有很多水陸混合像元的坑塘養(yǎng)殖區(qū),相比于ASWM,本文提出的SWES的RMSE降低最多。

表2 兩種方法水體提取結(jié)果RMSE

表3為SWES和ASWM在試驗區(qū)提取的水域面積誤差。通過表3可以看出,SWES在3個試驗區(qū)提取的水體面積精度均較高,平均相對誤差為8.03%,遠低于ASWM的平均相對誤差。兩種方法提取的水體面積均小于真實面積,這是由于試驗區(qū)有些區(qū)域水體寬度小于3個像元,兩種方法皆難以提取。

表3 兩種方法提取水體面積誤差表

4 結(jié) 語

本文首先利用RWI提取純水體像元;然后根據(jù)水陸混合與純水像元的空間關(guān)系,提出用膨脹算法提取水陸邊界混合像元;最后在對水陸邊界混合像元進行分解時,為了解決地物的類內(nèi)光譜變化問題,采取了多端元光譜混合分析算法,在提取水體端元時,基于地理相似定律提出了動態(tài)選擇水體端元的方法,進而構(gòu)建了面向Sentinel-2數(shù)據(jù)的亞像元級水體提取方法SWES。

SWES提取水體的精度達到亞像元級別,充分考慮了水體光譜的變化,結(jié)合空間信息提取水體端元,減少了解混計算成本。在3個不同水體類型的研究區(qū)開展了水體提取試驗,結(jié)果均表明SWES取得較好效果,平均RMSE為0.147,在3個試驗區(qū)上的水體提取效果均優(yōu)于ASWM,尤其在水陸混合像元較多的坑塘養(yǎng)殖區(qū)。SWES在試驗區(qū)獲取的水體面積也有較高精度。由此可見,SWES既能有效提取純水體像元又能有效提取水陸混合像元及水體豐度,準確獲取水體面積,將基于Sentinel-2數(shù)據(jù)的水體提取進行降尺度,解決多源中高分辨衛(wèi)星數(shù)據(jù)提取水體時的空間尺度效應(yīng)問題。