面向Sentinel-2影像的亞像元級水體提取方法

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

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

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

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

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

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

1 研究區及數據

1.1 研究區概況

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

圖1 研究區位置

1.2 數據及處理

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

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

2 方法構建

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

2.1 提取純水像元

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

(1)

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

2.2 提取水陸混合像元

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

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

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

2.3.1 光譜混合分析模型

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

(2)

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

(3)

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

(4)

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

2.3.2 結合空間信息的多端元動態選擇

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

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

2.3.3 端元組合

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

表1 端元組合類型

3 結果與分析

3.1 目視分析

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

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

圖4 SWES和ASWM在試驗區2估算水體豐度的結果

圖5 SWES和ASWM在試驗區3估算水體豐度的結果

3.2 精度評價

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

(5)

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

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

表2 兩種方法水體提取結果RMSE

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

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

4 結 語

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

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