基于GEE 平臺的黃河流域水體指數研究

孫根云，邵寶婕，丁孫金衍，姚延娟，付 航，王霓妮

（1.中國石油大學（華東） 海洋與空間信息學院，山東 青島 266580；2.青島海洋科學與技術試點國家實驗室 海洋礦產資源評價與探測技術功能實驗室，山東 青島 266071；3.生態環境部 衛星環境應用中心 國家環境保護衛星遙感重點實驗室，北京 100094）

1 引 言

黃河流域是我國重要的生態安全屏障［1］和經濟發展的重要區域，快速、準確地提取流域水體對流域水資源調查、綜合治理等具有重要意義。當前，基于遙感技術的水體提取方法已被廣泛應用于水資源調查、海岸線動態變化監測等領域［2］，水體提取方法主要分為水體指數方法和機器學習方法［3－5］。其中機器學習方法基于人工提取特征，使用分類器來實現水體和背景地物的區分，通常需要大量樣本訓練分類器，樣本質量會直接影響最終的水體提取精度，并且對于大范圍的水體，難以獲得優質、充足的訓練樣本［6－10］。相較而言，水體指數方法利用波段間的代數運算來凸顯水體與背景地物的差異，此類方法形式簡單、處理數據快速有效、分類精度較高，在大尺度水體提取中具有很大優勢［11］。近年來，相關研究人員提出了一系列水體指數方法，如：MCFEETERS［12］利用水體在綠色波段反射率高、在近紅外波段反射率低的特點，提出了歸一化差異水體指數（NDWI），NDWI方法可以有效提取水體，但提取過程中容易將建筑物錯分為水體。針對這一問題，徐涵秋［2］采用中紅外波段代替近紅外波段，提出了改進的歸一化差異水體指數（MNDWI），但該方法容易受山體陰影（簡稱山影）的影響。Jiang 等［13］基于Sentinel－2 影像提出了陰影水體指數（SWI），其對于純水、渾濁水、鹽水和浮冰有較好的提取效果。盡管上述水體指數方法取得了不錯的提取效果，但對于大范圍復雜水體的提取效果仍有待驗證。

對于大尺度流域水體的提取，通常需要大量的遙感影像才可實現區域全覆蓋，而傳統的遙感處理平臺難以實現在某一特定時期、較大范圍內對大量衛星影像的快速調用［14］。GEE 是專門用于處理衛星圖像和其他地理數據的平臺，現階段已被廣泛用于地理數據大規模分析［11，15］。本文基于GEE 平臺，在深入分析黃河流域不同背景地物、不同狀態水體光譜特征的基礎上，提出一種黃河流域水體指數（YRBWI），并與已有的水體指數NDWI、MNDWI和SWI進行對比，驗證其去除遙感影像固有山影的有效性以及提取精度。

2 研究區概況與數據選取

2.1 研究區概況

黃河流經黃土高原水土流失區，是我國西北、華北地區的重要生態屏障，為我國北方生態文明建設提供了載體，源區被譽為“中華水塔”［16］。流域面積為79.5萬km2，從西到東橫跨青藏高原、內蒙古高原、黃土高原和黃淮海平原4 個地貌單元，是我國第二大流域［17］。黃河源頭常年積雪，上游河水清澈，大面積沼澤地和湖泊較多；中游流經黃土高原地區，水土流失嚴重；下游泥沙量大，河道泥沙淤積現象嚴重［14，16］。流域獨特的地貌特征和水文現象決定了常規的水體提取方法無法滿足全流域水體自動提取的要求。流域水體提取主要存在以下難點：流域地物復雜、地形多變，山影嚴重干擾水體提??；流域水體復雜，高含沙量渾濁水體與湖泊、水庫等清澈水體的光譜特征差別較大，傳統水體指數大多只對清澈水體有較好的提取效果，無法滿足流域內復雜水體的高精度提取。

2.2 數據選取

在眾多光學數據中，Sentinel－2 影像數據是唯一在紅邊范圍內含有3 個波段（B5、B6 和B7）的數據。選用GEE 平臺中2018 年所有可用的黃河流域Sentinel－2 影像，首先通過QA60 波段對云陰影等觀測質量較差的像素進行掩膜，再對所有影像進行中值合成，最終將其用于2018 年黃河流域水體提取。結合Sentinel－2 影像和Google Earth 高分辨率遙感影像，選取2 種典型水體（清澈水體、渾濁水體）和6 種典型背景地物（農田、沙漠、山影、植被、建筑物、道路），對每種水體和背景地物分別選取100 個樣本點用于繪制光譜曲線。通過分層隨機抽樣，結合Sentinel－2 和Google Earth 影像進行目視解譯生成940 個水體樣本和1 073 個非水體樣本，用于水體提取結果的精度評價。Sentinel－2 影像光譜覆蓋了430～2 300 nm 范圍內的13個波段，可見光、近紅外波段的空間分辨率為10 m，紅邊（Red＿edge）、短波紅外（SWIR）波段的空間分辨率為20 m，衛星重訪周期為5 d。

3 技術與方法

3.1 YRBWI 的構建

黃河流域地表水體的反射特性經常受到周圍環境（山影、泥沙等）的影響。為研究不同水體及背景地物的光譜特征，分別繪制其光譜曲線（見圖1），參考已有研究［13，18］，選擇Sentinel－2 影像的B5 波段作為紅邊波段。可以看出，從藍色波段到紅邊波段，清澈水體和渾濁水體的反射率變化呈現相反趨勢，且二者與背景地物的反射率相近，無法同時實現兩種水體的提取。從紅邊波段到近紅外（NIR）波段，背景地物的反射率均表現為上升趨勢，而清澈水體和渾濁水體的反射率均表現為下降趨勢，水體與背景地物的反射率差異明顯，容易進行區分。從近紅外波段到SWIR2 波段，水體與部分背景地物的反射率相近且均表現為下降趨勢，無法進行區分。

圖1 黃河流域不同水體和背景地物的光譜曲線

對圖1 分析可知，從紅邊波段到近紅外波段，水體與背景地物的反射率差異明顯、變化趨勢相反，因此利用這兩個波段構造水體指數，以達到同時提取清澈水體、渾濁水體并消除背景地物影響的效果。據此構建的初步水體指數（PWI）如下：

式中：ρRed＿edge、ρNIR分別為影像在紅邊和近紅外波段的反射率。

水體與背景地物相反的反射率變化趨勢是區分水體與背景地物的關鍵。然而，根據圖1 可知，清澈水體和渾濁水體的反射率在紅邊波段和近紅外波段的下降幅度較小，這可能導致水體提取效果不明顯。為凸顯這種變化趨勢，進一步對PWI進行差異增強處理。對于式（1）中PWI的分子項，使用斜率大于1 的單調遞增函數可以增強水體反射率的下降趨勢。考慮到各種背景地物在兩個波段的反射率都在0～1 之間，選取常用的y＝sinx、y＝tanx、y＝ex、y＝lnx這4 種單調遞增函數作為增強函數。由于斜率為1 的y＝x函數能夠反映出背景地物在兩個波段原本的反射率差值，因此以y＝x作為對照函數，5 種函數在反射率為0～1 之間的函數曲線見圖2。以渾濁水體為例，5 種函數處理后渾濁水體在紅邊和近紅外波段間的反射率差值見圖3。

圖2 5 種函數在反射率為0～1 之間的函數曲線

圖3 5種函數處理后渾濁水體在紅邊和近紅外波段間的反射率差值

分析圖2 可知，反射率在0～1 之間時y＝lnx函數的斜率（增長速率）大于其他函數，且反射率小于0.3時y＝lnx函數的增長趨勢更加明顯，更加有利于加大水體在紅邊和近紅外波段之間的反射率差異。對比圖3中不同函數的反射率拉伸效果可知，y＝sinx、y＝tanx函數幾乎沒有拉伸效果，y＝ex函數的效果稍好，y＝lnx函數的拉伸效果最好。因此，選取y＝lnx作為拉伸函數，對水體在紅邊和近紅外波段的反射率差異進行增強處理。對于式（1）中PWI的分母項，考慮到PWI分母變大會導致其整體減小，從而造成反射率差異效果的抵消，因此對分母項保持原有形式，不做拉伸處理，最終得到黃河流域水體指數（YRBWI）：

3.2 對比水體指數

為驗證YRBWI水體指數的有效性，選取NDWI、MNDWI和SWI水體指數進行對比，這3 種水體指數公式分別為

式中：ρGreen、ρSWIR1分別為影像在綠色、SWIR1 波段的反射率。

3.3 精度評價指標的確定

采用總體精度（OA）和Kappa系數評定各水體指數提取水體的精度。總體精度指被正確分類的類別像元數（樣本數）與總類別像元數的比值；Kappa系數代表正確分類與完全隨機分類產生錯誤減少的比例。OA和Kappa的計算公式如下：

式中：pkk為被正確分為k類的樣本數；p為樣本總數；pii為被正確分類為i類的樣本數；pi＋為分類所得到的第i類的樣本總數；p＋i為目視解譯得到的第i類的樣本總數。

4 結果分析

4.1 精度評價

各水體指數的最佳閾值以人工方式確定，將0 作為基本閾值，在此基礎上多次調整，選取水體提取精度最高時對應的閾值作為該指數的最佳閾值（T）。經計算，YRBWI、NDWI、MNDWI、SWI的最佳閾值分別為－0.05、－0.10、0.25、0.15，4 種水體指數在最佳閾值下的OA值分別為96.87%、93.64%、95.78%、96.08%，Kappa系數分別為0.937 0、0.872 4、0.914 9、0.921 1，因此YRBWI的OA和Kappa值最大，精度最高。

4 種水體指數在最佳閾值下對渾濁水體、清澈水體和山影的正確分類率見表1，可以看出，與NDWI、MNDWI和SWI相比，YRBWI對渾濁水體和山影有最高的正確分類率，盡管對于清澈水體來說，YRBWI的正確分類率略低于NDWI，但仍取得了99%以上的效果。說明YRBWI在準確提取黃河流域高含沙量水體的同時，還能有效去除流域內大量的山影噪聲，YRBWI能夠很好地適應該流域的廣闊范圍和特殊的地理環境。

表1 4 種水體指數對渾濁水體、清澈水體和山影的正確分類率

4.2 水體指數增強性能

為比較不同水體指數對水體與背景地物之間差異的增強性能，不同水體和背景地物的各水體指數均值見圖4（a），圖4（b）是對圖4（a）中NDWI、MNDWI和SWI均值的放大圖。與其他3 種指數相比，不同水體和背景地物的YRBWI均值的差異最大，說明YRBWI有最好的增強性能。另外，清澈水體、渾濁水體的YR?BWI均值均大于0，背景地物均值均小于0，說明YRBWI的分類效果較好。分析圖4（b）可知，NDWI、MNDWI和SWI對于水體與背景地物之間的差異有一定的增強效果，但渾濁水體的NDWI均值小于0，說明NDWI存在把渾濁水體錯分為背景地物的情況；山影的MNDWI和SWI均值大于0，說明MNDWI和SWI存在把山影錯分為水體的情況。由此可知，與其他3 種指數相比，YRBWI對水體與背景地物之間差異的增強效果明顯，具有較好的分類效果。

圖4 不同地物類型的各水體指數平均值

4.3 不同水體指數的提取效果對比

為驗證YRBWI在全流域內不同區域的提取效果，選取3 個代表性區域（見圖5），不同區域的各水體指數提取結果見圖6～圖8。圖6 中細長的線狀部分是水體，YRBWI提取的水體清晰明亮且與周圍背景地物差異明顯；MNDWI和SWI提取的水體較為細小，存在部分水體的缺失；NDWI提取的水體雖然較為完整，但是周圍山影噪聲較多，水體邊界不清晰。圖6 中紅框部分是山影，YRBWI提取的水體幾乎不存在山影噪聲，而NDWI、MNDWI和SWI提取的水體均存在大量山影噪聲，錯分情況嚴重。圖7 中綠色橢圓部分是清澈的湖泊，藍色橢圓部分是渾濁的黃河水體，YRBWI、NDWI和SWI提取的清澈水體、渾濁水體都較為完整，相比之下，MNDWI提取的渾濁水體略有缺失；圖7 中的紅框部分是城市區域，YRBWI、MNDWI和SWI能夠在錯分極少的情況下很好地提取城市中較為細碎的水體，而NDWI存在嚴重的背景地物錯分現象。圖8 中綠色橢圓部分是黃河支流，YRBWI和NDWI提取的黃河支流水體清晰明亮，MNDWI和SWI能夠識別黃河支流但是水體缺失嚴重；圖8 中的紅框部分是大量積雪和部分山影，YRBWI對于積雪和山影噪聲的去除效果很好，NDWI提取的水體存在少量積雪，而MNDWI和SWI對于積雪的去除效果極差，積雪的明亮程度甚至超過綠色橢圓部分的水體，說明積雪錯分情況嚴重。綜上來看，與其他3 種水體指數相比，YRBWI對黃河流域的水體提取效果最好。

圖5 黃河流域內選取的3 個代表性區域

圖6 A 區域的各水體指數提取結果

圖7 B 區域的各水體指數提取結果

圖8 C 區域的各水體指數提取結果

4.4 閾值敏感性分析

為比較各水體指數對于閾值的敏感性，計算4 種水體指數在不同閾值下的精度值。根據上文確定的YRBWI、NDWI、MNDWI、SWI的最佳閾值分別為－0.05、－0.10、0.25、0.15，以0.05 為間隔選取最佳閾值附近的8 個閾值。其中：YRBWI的閾值范圍是－0.10～0.25，NDWI的閾值范圍是－0.20～0.15，MNDWI的閾值范圍是－0.05～0.30，SWI的閾值范圍是－0.10～0.25。各水體指數在不同閾值下的OA和Kappa值的變化見圖9?？梢钥闯觯S閾值的增大，4 種水體指數的OA和Kappa值都呈現先增大后減小的趨勢，但精度值的變化范圍有所不同。以OA為例，YRBWI的OA值在最佳閾值后雖有所下降，但總體上變化較為平緩，在－0.10～0.25閾值范圍內其OA的最大值和最小值相差僅5%左右。而NDWI在－0.20～0.15 閾值范圍內其OA的最大值和最小值相差約20%，MNDWI和SWI在各閾值范圍內其OA的最大值和最小值均相差15%以上。由此可見，YRBWI的水體提取精度隨閾值變化的波動較小，對閾值變化不敏感，抗干擾性強。而其他3種水體指數尤其是NDWI的水體提取精度隨閾值變化波動較大，對閾值變化較敏感，抗干擾性弱。

圖9 4 種水體指數在不同閾值下的精度變化

此外，在－0.10～0.10 閾值范圍內，YRBWI的精度最高，而閾值分別大于0.10 和0.15 時，其精度分別低于SWI和MNDWI的精度。以SWI為例，結合圖4（b）可知，山影的SWI均值大于0，說明大部分山影的SWI值大于0，因此在閾值大于0 的基礎上適度增大閾值可以去除SWI提取結果中山影噪聲的影響，從而提高SWI的提取精度。而對于YRBWI，結合圖4（a）可知，水體的YRBWI均值大于0，背景地物的YRBWI均值小于0，因此在大于0 的基礎上增大閾值會使YRBWI存在水體漏分現象，從而導致精度降低。

5 結 論

針對黃河流域的復雜水體信息，提出一種新的水體指數YRBWI用于黃河流域的大尺度水體提取，結論如下。

（1）使用紅邊和近紅外波段構造PWI，在對比多個增長函數后選擇y＝lnx函數對PWI進行增強處理得到YRBWI。與NDWI、MNDWI和SWI這3 種水體指數相比，YRBWI在黃河流域水體提取中表現優異，既能有效增強渾濁水體和清澈水體與背景地物的光譜差異，準確提取高含沙量渾濁水體，也能解決大部分指數存在的山影問題，可實現快速、自動、高精度的水體提取。

（2）相比于其他指數，YRBWI對閾值變化不敏感，具有良好的抗干擾性。在未來研究中，可考慮結合水體的時間序列信息來進行黃河流域水體提取和變化研究。