基于Sentinel-2遙感影像的水庫庫容曲線快速構建
2025-09-02 00:00:00韓廣
人民珠江 2025年7期
中圖分類號：TV123；TP79 文獻標識碼：A 文章編號：1001-9235（2025）07-0086-15
Abstract：Asasetofautomated process hasnotyet been formed forconstructing reservoircapacitycurvesusingremotesensing technology，itisproposedtodeterineteptialtresholdwiththeinflectiopontmethodoftheumulativefrequncycurebsd onthewaterbodyindexmethod.Theedgeofthe waterbodyisoptimizedbytheregionalseedgrowthalgorithm，soastoconstructa completesetofautomatedprocess toquicklyconstruct thereservoircapacitycurves．Withentinel-2remotesensingimagesasthe datasourceadCheguReservoirinWuanCityHandanCityHebeiProvinceasteresearchbect，tispaperusesielandadater image metrics，amely，NDWMDWI，RWI，BWndAWEInshtoextracttheareaofthewaterbodyostructteeor capacitycurve，andcompareitwiththeresultsofthemeasurementsin2O17.Theresultsareasfollows.Theextractionaccuracyof NDWI and RWI waterbody index is high.The average error of Chegu water body area extraction is 2. 5% . The relative error is -7 3% } （20 5.8% and -6 2% ～4. 7% ， and the R2 is O.993 36 and O.990 49，respectively. The range of the relative error of the reservoir capacity is -9 4% to -1 8% and -9 3% to -1. 7% ，respectively. The study shows that the cumulative frequency curve inflection point method ishighlyaplicableandstableindeterminingtheoptimalthresholdof waterbodyindex，withasimpleandfeasibleprinciple.The regionalsedgrowthalgorithmcanefectivelyeliminatethediscontinuityattheedgeofthewaterbody，reducesystematiceo，nd furtherimprovetheaccuracy.Thecumulativefrequencycurveinflectionpointmethodprovidesanewideaforremotesensingwater bodyextraction，servingasareferenceforthebatchconstructionofnationalsmallandmedium-sizedreservoircapacitycurve extraction.
Keywords：sentinel-2;water Index; inflectionpointmethodofthecumulativefrequencycurve;growthalgorithm;storageapacity curve; Chegu Reservoir
水庫是流域防洪工程體系的重要組成部分，是國家水網重要結點，是保障國家水安全的重要基礎設施。水庫庫容是水庫有效發揮防洪、供水、生態、發電等功能的重要保障，事關工程安全、防洪安全、供水安全和生態安全。截至2023年，全國已建成各類水庫94877座，其中大型水庫836座，中型水庫4230座，小型水庫89811座[1]。中國各類水庫的80% 修建于20世紀50—70年代，受當時基本圖件、測量儀器及方法限制，庫容曲線精度偏低，加上水庫運行期間受到庫區微地貌、淤積、土地利用變化等因素，水庫庫容曲線難免發生變化，直接影響到水庫防洪計算及發電調度運行[2]。為強化庫區、壩體、下游河道監管與問題整治，維護水庫庫容安全、工程安全、河道安全，科學推進水庫減淤清淤，水利部辦公廳印發《關于加快開展承擔防洪任務的大中型水庫庫容曲線復核工作的通知》。
目前，中國水庫庫容曲線復核主要采用陸域地形測量和水域地形測量相結合的方法，其中陸域地形測量主要有航空攝影測量、機載激光掃描和全野外數字測量等方法；水域地形測量主要為GNSSRTK與水深測量相結合的方法。這項工作不僅涉及對水庫的實地測量，還需要借助專業的測繪設備和軟件，通常需要投入大量的人力、物力和財力。中國水庫眾多，分布范圍廣，各地經濟狀況和水庫管理水平各異，尋找一種即實用又能滿足測量精度的方法復核水庫庫容曲線，顯得尤為重要。
為貫徹落實《水利部印發關于加強水文情報預報工作指導意見的通知》（水文[2019]203號）《國務院辦公廳關于切實加強水庫除險加固和運行管護工作的通知》（國辦發[2021]8號）要求，各省均開展了中小型水庫雨水情測報工作，實現對庫水位實時監測，積累了大量的水位觀測資料，為使用高分遙感影像提取水庫庫容曲線提供了前提條件。
隨著遙感技術的快速發展和普及，高分遙感影像逐步被應用于水體遙感解譯中，在水利方面的應用越來越成熟廣泛[3]。陳曦等4]使用ERDASIMAGINE遙感圖像處理系統，利用CBERS-02B、資源一號02C以及資源三號衛星高分辨率遙感影像，推求桐仁橋小水庫庫容曲線。張莉芳等5采用Landsat衛星遙感影像，提取了柘林水庫庫容曲線。王金鑫等提出了多源時相遙感影像數據耦合的庫容曲線重構方法，并重構了深圳市徑心水庫庫容曲線。
雖然目前基于遙感技術對水庫庫容曲線提取已得到應用，但仍存在批量提取水體面積效率低等問題，難以形成一套自動化流程。因此，本文基于水體指數法，提出利用累積頻率曲線拐點法自動確定最優閾值，應用區域種子生長算法優化水體邊緣，以實現高效、批量地提取水體面積，并采用JAVA語言，使用GDAL庫編制相關計算程序，構建了一套完整地自動化流程來快速構建水庫庫容曲線。
1 研究方法
1.1遙感影像水體提取
目前使用遙感影像提取水體的算法很多，一般而言，這些方法大致可分為單波段法、譜間關系法、影像分類法和水體指數法等幾個類別，其中水體指數法被研究人員廣泛采用[。水體指數是一種簡單且高效地提取水體信息的方法。但由于受區域自然條件、成像時大氣條件、傳感器差異等因素影響，水體指數法的閾值隨不同地區和季節影像的不同而有差異。閾值的不穩定性和設定的繁瑣性常常是水體遙感指數被詬病之處[8]。
為實現大批量、自動化、高效地提取水體信息，本文采用水體指數法。同時為實現閾值的自動確定，從統計學角度，按聚類分類的思想，繪制水體指數累積頻率曲線，經數學分析提出最優閾值位于曲線上由凹變凸的拐點處，進而提出一種高效、易編程的自動化閾值確定方法，即累積頻率曲線拐點法。
1.1.1水體指數法
水體指數的構建基本可分為2種類型：比值型和差值型。比值型指數通過比值運算擴大水體和非水體信息的區別，分離度較強，并能抵消地形的影像。由于其變量常不帶系數，因此具有較強的普適性。差值型指數則多帶有經驗系數，在不同衛星、不同地區和不同環境中表現不太穩定[8。本文選取國內外比較有影響力的水體指數，其中比值型選取歸一化差異水體指數（NormalizedDifferenceWaterIndex，NDWI）[9]和改進的歸一化差異水體指數（Modified Normalized Diference Water Index，MNDWI）[0]；差值型指數選取多波段水體指數（Multi-bandWaterIndex，MBWI）和自動水體提取指數（Automated Water ExtractionIndex，AWEInsh）[12];以及專門用于Sentinel-2提取細小水體的植被紅邊水體指數算法（RedsideWaterIndex，RWI）[7]。基于Sentinel-2遙感影像的各水體指數計算見式（1）—（5）：
1. 1. 2 累積頻率曲線拐點法
1.1.2.1 理論分析
利用水體指數法提取水體信息時，自動閾值的確定非常關鍵。常用的方法有最大類間方差法（Otsu法）迭代法、直方圖雙峰法等，以及相關的改進算法[13-14]。每種方法的核心思想均是利用一定的原則，進行分類，且每種方法都有各自的優缺點和適用范圍。如：Otsu方法已被許多研究證明其確定的最優閾值并不可靠[15-17]；直方圖雙峰法常常會面臨單峰或多峰的情況，導致獲取閾值結果無意義。
水體指數法原理是基于水體與其他地表覆蓋類型在特定波段反射的差異，通過波段運算，經歸一化處理后，將影像二值化，以區分水體和其他地表覆蓋物。該方法本身具有較強的物理基礎和理論依據，且該方法的本質是通過指數值的差異來分類。為更好有效地利用水體指數法，建議應先對水體指數圖進行聚類再分類。統計分析中累積頻率曲線常用于分析數據的分布特征、中心趨勢和離散程度等。由于累積頻率曲線生成過程中需要先排序，在一定程度上完成了聚類。本文通過繪制水體指數圖的累積頻率曲線，提出最優閾值位于曲線上由凹變凸的拐點處。
以為橫坐標，以水體指數 I 為縱坐標，繪制累積頻率曲線圖。其中P（I?i） 表示水體指數圖中像元數值大于等于 i 的百分比； 
為指示函數，當滿足 I?i 時取1，否則取 0;N 為水體指數圖的像元總數。
完全二值化影像圖的累積頻率曲線為不連續的分段函數（圖1），曲線中 iw 代表水體像元值， io 代表非水體像元值。
圖1 累積頻率曲線
Fig.1Cumulativefrequencycurves
實際中，對于歸一化的水體指數圖，其水體指數為在某范圍內變動的離散函數，但整體呈二值化趨勢。水體指數圖的累積頻率曲線為非嚴格單調遞減函數，即 I′（p）?0 。實際應用中，經統計分析得到的累積頻率曲線圖是離散的點集。為方便分析，用光滑的曲線將點集連接，并假定曲線是連續的，二階可導的。
曲線上 （pmw，iw） 代表水體像元密度最大的點。在該點指數值 i 的小幅減少，將引起水體累積頻率 p 的大幅增加，即 
為附近鄰域[pmw-ε，pmw+ε] 內的最大值。由于 I′（p） 連續可導，因此 
為極大值點，即 
。由于
因此區間 [pmw-ε，pmw] 內 
，區間 [pmw，pmw+ε] 內
。由函數凹凸性和拐點定義可知，點 pmw 為曲線上由凸變凹的拐點。曲線上 （pmo，io） 代表非水體像元密度最大的點，在該點指數值 i 的小幅減少，將引起非水體累積頻率 p 的大幅增加。同理點pmo 也為曲線由凸變凹的拐點。
由于 
和 
均為附近鄰域內最大值（極值點）。依據連續函數的有界性定理，設 I′（p） 在區間 
內的最小值為 I′（pw） ，則存在區間[pw-ε，pw+ε] ，使得 I′（pw-ε）gt;I′（pw） ，且 
ε）gt;I′（pw） 由于 I′（p） 在閉區間 [pw-ε，pw+ε] 上是連續可導的，因此 pw 為 I′（p） 的極值點，即
。由于區間 [pw-ε，pw] 上 
，區間 
，因此 pw 為由凹變凸的拐點，即在 pw 點附近隨累積頻率的增加，指數值的變化速度先減少后增加。
對于無其他地物干擾，應用效果較好的水體指數，水體和非水信息分離度較大且不重合，即累積頻率曲線上 （pmw，iw） 和 （pm0，io） 是唯一的，區間[pmw，pm0] 內必然存在一個拐點 pw 。累積頻率曲線圖在區間 
內的變化趨勢是水體指數二值化效果的直接體現。在 pw 點附近存在某個頻率范圍，表達了水體向非水體的過渡過程。該頻率范圍越小，凹凸變化越明顯，說明在較小的累積頻率區間水體指數的變化范圍越大，二值化效果越好。
綜上所述，提出點 pw 為水體和非水體的分界點 .I（pw） 即為最優閾值。因此，可以通過尋找累積頻率曲線上由凹變凸的拐點來確定最優閾值。此外，水體指數圖的累積頻率曲線突出了水體向非水體的轉變過程（區間） 
，弱化了其他影響，即使影像中水體占比較小，也可以較好地區分水體和非水體。但同時也要求有足夠的非水體像元，以保證累積頻率曲線的完整性，才能充分體現水體向非水體的轉變過程。
1. 1. 2. 2 算法實現
水體指數最優閾值的確定可以轉換為尋找累積頻率曲線上由凹變凸的拐點 pw 。實際應用中，由于受區域自然條件、大氣條件和季節等因素影響，相同區域不同時期采用相同的水體指數法對水體提取的效果可能不同。體現在累積曲線上：一是可能不存在由凹變凸的拐點，二是可能存在多個由凹變凸的拐點。因此，本文給出以下建議： ① 若累積頻率曲線上不存在由凹變凸的拐點，說明水體和非水體的區分界面不明顯，這種情況下最優閾值的選取較為敏感，即閾值的較小變動，將引起較大的誤差。此時水體指數法已無法有效區分水體和非水體； ② 若累積頻率曲線上存在多個由凹變凸的拐點，說明水體指數圖中存在較多類水體的干擾信息，如云層、冰雪等，基于水體指數二值化的原理，累積頻率曲線上除兩端某區間范圍外，最優閾值處I′（pw） 最小，由于 I′（p）?0 ，即 |I′（pw）| 最大。此時，建議先利用簡化線算法進行處理，在保持曲線形狀的同時，減少干擾信息，之后再尋求由凹變凸的拐點。若仍存在多個凹變凸的拐點，則判斷各拐點處的 |I′（pw）| 值，選取 |I′（pw）| 值最大處的閾值作為最優閾值。
綜上所述，水體指數最優閾值確定的關鍵步驟如下。
a）采用道格拉斯-普克抽稀算法（DP算法）[18]進行簡化線處理。選擇累積頻率曲線首尾兩點連接，遍歷兩點之間的所有點，找到距離首尾連線垂線距離最大的點，若最大垂距小于設定閾值，則刪除中間所有點；若最大垂距大于設定閾值，則保留此點，并將其與首尾相連，得到2段新的線段。分別對新線段進行如上操作，直至判斷完整個曲線上的點。由于構建累積頻率曲線時，兩軸的意義并不相同，因此建議采用垂距在1軸上投影長度和設定閾值進行比較。此外，不同水體指數得到的指數范圍是不同的。為具有普適性，本文將閾值的確定轉化為簡化后至少保留的點數，即由大到小試算閾值，當保留點數滿足要求時便停止。
b）利用離散型函數拐點算法分析由凹變凸的拐點區間。設累積頻率曲線共有 n 個散點，分別為 p0，p1，p2，…，pk…，pn-1， 對應的縱坐標為 I（p0） ，I（p1），I（p2），…，I（pk）…，I（pn-1） 。如果拐點處于點pk-1 與 pk（k=3，4，5…n-3） 之間，則有：
若在 pk-1～pk 區間，累積頻率曲線由凹變凸，則有：
若存在多個由凹變凸的區間，則判斷
值，選取該值最大的區間。
Πc） 判斷 pk-pk-1 是否滿足精度要求。若滿足，則最優閾值取 I（pw）=（I（pk）+I（pk-1））/2 。若不滿足，則按精度要求，對 pk～pk-1 區間進行細化，選取
最大點對應的水體指數為最優閾值。
1.1.3區域種子生長算法
研究中發現：基于水體指數采用累積頻率曲線拐點法確定最優閥值，提取的水體面積整體偏小，存在明顯系統誤差。經初步分析主要原因為：受影像空間分辨和水體指數法自身限制（水體邊緣處指數存在突變），導致水體邊緣提取不完整。為進一步提高水體提取精度，本文采用區域種子生長算法[20]對水體邊緣進行優化。
應用區域種子生長算法是對水體指數法的補充，為實現方法間的互補，生長準則的制定不易再基于水體指數。參考人工目視的原理，考慮到水對紅外波段的吸收系數值非常大[21]，將R（紅）G（綠）B（藍）色彩模式中的R（紅）替換成NIR（近紅外），形成偽色域（NIRGB），并基于標準色差公式2來制定生長準則。
通過歸一化處理，得到每個像元的NIR、G和B的亮度值，見式（8）一（10）。區域生長時，以種子像元為基準，向外逐層生成緩沖區，搜尋連續的目標像元，并判斷目標像元和種子像元的色差，見式（11）。當色差小于給定閾值時，將目標像元歸為水體，實現增長。
式中： x 為影像柵格的行號； y 為影像柵格的列號；LNIR?Lg 和 Lb 分別為 （x，y） 處近紅外、綠和藍光的亮度值； R8，R3 和 R2 分為Sentinel-2的第8、3和2波段的遙感反射率； R8，max?R3，max 和 R2，max 分別為影像柵格中相應波段的最大遙感反射率值; R8，min?R3，min 和R2，min 分別為影像柵格中相應波段的最小遙感反射率值； R8，（x，y）?R3，（x，y） 和 R2，（x，y） 分別為 （x，y） 處相應波段的遙感反射率值； ΔLc.ΔLNR.ΔLg 和 ΔLb 分別為目標像元與種子像元的色差、近紅外、綠和藍光亮度差值。
實際應用時可能會遇到生長不可控的情況，即將非水體誤判為水體。其原因一是由于閾值設置過于寬泛；二是由于水體指法誤將邊緣處非水體確定為水體。因此區域種子生長算法運行前，需要進一步判斷水體邊緣處和孤立的水體像元是否為真實的水體。
1.2庫容曲線快速構建
1. 2. 1 庫容計算
本文庫容計算采用臺柱體積法[5]。首先把水體按不同水位分成 n 層，依據臺柱體積公式計算每層的體積，經累加得到庫容。
臺柱體積公式：
累計庫容計算公式：
式中： Vi 為相鄰水位之間的庫容差值，萬 m3;Δh 為相鄰兩水位的水位差， m;Si-1 和 Si 分別為兩相鄰水位對應的水面面積， 104m2;i 為序數； n 為累計個數； V0 為初始庫容，萬 m3 ；V為累計庫容，萬 m3 。
1.2.2庫容曲線構建
本文研究目的是利用多時序的遙感影像快速自動構建水庫庫容曲線，更注重水庫庫容曲線構建的整體精度，并不需要每幅影像都能高精度地成功提取水體面積。通常受各種條件影響，水體指數法偶爾會出現不適用的情況。因此，當時序影像足夠多時，可適當舍棄誤差較大的，以保證整體的精度。庫容曲線快速構建流程如下： ① 依據各景影像生成的時間，從水庫觀測資料中提取相應的水位； ② 使用水體指數法對預處理后的各景影像進行歸一化處理，得到水體指數圖； ③ 對水體指數圖進行統計分析，得到累積頻率曲線； ④ 對累積頻率曲線進行簡化線處理，以減少干擾； ⑤ 在處理后的累積頻率曲線上查找由凹變凸的拐點區間，并確定最優閾值，查找拐點前，可先依據水位變動情況，估算或通過人工目視解譯選取影像中水體的最大和最小面積，并作為先驗條件輸入，以減少搜尋范圍，提高效率，若不存在由凹變凸的拐點區間，則舍棄該景影像； ⑥ 依據最優閾值提取水體，利用區域種子生長算法對水體邊緣進行優化，計算水體面積，建議輸入庫區最大矢量范圍作為限制條件，以減少庫外區域噪聲的影響； ⑦ 生成水位-水面面積散點圖，并采用多項式6對水位-面積曲線進行擬合； ⑧ 判斷水位-水面面積散點和擬合曲線的誤差關系，舍棄誤差較大的點，并重新擬合，反復進行該步驟直至各散點誤差均滿足精度要求為止； ⑨ 依據擬合的水位-面積曲線得到等米間距的水面面積； ⑩ 采用臺柱體積公式計算庫容差，并經累積得到水位-庫容曲線。
2 案例分析
2.1 研究區概況
以河北省邯鄲市武安市洺河流域的車谷水庫為分析對象。該水庫是一座綜合利用的中型水利樞紐工程，總庫容3315萬 m3 ，壩址以上集水面積為124km2 ，集水區范圍為36.816～36.995N、113.724～113.857E，上游主要河流為管陶川，見圖2。
圖2研究區概況
Fig.2Overviewofresearcharea
2.2 數據來源
2017年武安市水利局對車谷水庫庫容曲線進行復核，得到水庫水位、水面面積和蓄水量等特征數據，見表1。
表1車谷水庫特征數據
Tab.1CharacteristicdataofChegu Reservoir
本文選取Sentinel-2衛星多光譜遙感數據，共有13個波段，空間分辨率為 10，20，60m， A和2B兩顆衛星同時運行，時間分辨率5d，時間和空間分辨率較高。選取2016年6月至2020年12月期間云量小于10% 的71景影像，以及云量大于 15% 但庫區無云的2019年6月2日和7月2日2景影像，共計73景，見圖3。所選影像均為經大氣校正后的L2A級產品。分析前，使用歐空局推薦的哨兵數據應用平臺（SentinelApplicationPlatform，簡稱SNAP，版本10.0.O），對各景影像進行重采樣和裁剪，得到水庫長時序的 10m 分辨率影像。分析范圍共21. 62km2 。
圖3選取衛星遙感影像的時間矩陣
Fig.3Temporalmatrixofselected satelliteremote sensing images
2.3水體面積提取評估
以武安市車谷水庫為研究對象，采用JAVA語言，使用GDAL庫編制相關計算程序，分別計算73景影像的NDWI、MNDWI、RWI、MBWI和AWEInsh指數，先利用累積頻率曲線拐點法提取水體，再采用區域種子生長算法對水體邊緣進行優化，最后計算水體面積，并和2017年實測水位-面積進行比較，分析提取精度，見表2和圖4。依據水庫蓄水情況，選取2016-06-22（水位 681.01m ）、2017-05-18（水位 704.51m ）、2018-04-08（水位 698.21m 和2020-09-04（水位 688.99m 作為典型時期進行分析，見表3、圖5—7。
表2水體面積提取成果
Tab.2Resultsofwaterareaextraction
圖4遙感解譯水體面積和實測面積對比
Fig.4Comparison of remote sensing interpreted water bodyarea and actual measuredarea
表3典型時期水體面積提取成果
Tab.3Resultsofwaterareaextractionfortypical periods
圖5典型時期累積頻率曲線
Fig.5Cumulative frequencycurves for typical periods
圖62016年6月22號水體提取對比
Fig.6Comparison of waterbody extraction on June22，2016
圖72020年9月4號水體提取對比
Fig.7Comparison of water body extractionon September 4，2020
程序運行過程中有3處主要參數需要設置：一是將累積頻率曲線等分成多少份，均采用500；二是由DP算法處理后，累積頻率曲線至少保留多少個點，以水位-面積擬合較好，且影像利用率盡量高為條件進行試算，本文NDWI、MNDWI、RWI、MBWI和AWEInsh指數分別采用200、250、150、250和200；三是生長準則的閾值，車谷水庫水質較好，庫區顏色比較一致，均采用15。
經分析，采用頻率曲線拐點法能夠較好地確定最優閾值，各水體指數均能較好地利用遙感影像，NDWI指數利用率最高為 100% ，MNDWI利用率最少為 90.4% 。NDWI指數和RWI指數提取水體比較完整，穩定性好，且RWI指數對細小水體的提取有更好的效果；MNDWI指數法雖然在描述水體邊緣具有明顯優勢，但水體邊緣指數的明顯變化，往往導致水體提取不完整，水體具有明顯的不連續現象，導致誤差較大。同樣，MBWI和AWEInsh指數也存在不同程度的邊緣不連續現象，見圖6、7。
利用區域種子生長算法對提取水體邊緣進行優化后，可明顯提高水體提取質量，各指數法水體面積提取的平均誤差由 4.6%～9.4% 降低至 2.5% ～3.1% 。優化后NDWI指數提取水體效果最好，有效利用73景，水體面積相對誤差 -7.3%～5.8% ，平均誤差 2.5% R2 為0.99336；其次是RWI指數，有效利用72景，水體面積相對誤差 -6.2～4.7 ，平均誤差2.5% R2 為0.99049；MNDWI指數僅利用69景，水體面積相對誤差 -7.3～4.2，R2 為 0.992 19 。
通過對典型時期水庫水體提取分析可知，2016年6月22日應用區域種子生長算法后水體面積變化不大，但與實測值比較相對誤差在 -3.7%～-6.9% 。經人工目視檢查RWI指數提取成果比實際水面略大，其他指數和實際水面基本吻合。因此，考慮該景影像提取的誤差主要是由于水位觀測導致的，其原因為：本次采用水位為日均水位，而影像生成時間為一天中某個時段。車谷水庫承擔下游供水任務，當水位較低蓄水量較少時，供水將導致一天內水位出現明顯的波動。
綜上所述，頻率曲線拐點法對各類型指數最優閾值的確定均有較好的適用性，結合區域種子生長算法可以有效消除水體邊緣不連續的現象，減少系統誤差，進一步提高水體提取精度。
2.4水庫水位-面積曲線擬合
水庫水位-面積曲線常采用多項式進行擬合。利用不同時期的水位和提取的水體面積繪制散點圖，選取效果較好的二次多項式進行擬合，見圖8。使用不同水體指數所擬合的水位-面積曲線的相關系數均大于0.97，說明得到的關系式能較好反映水庫水位和水面面積之間的關系。
在選取73景影像中，車谷水庫水位變動范圍為681.01～706.47m? 。因此，在 681～707m 取等米水位，從水位-面積擬合曲線上提取面積，并和2017年測量面積進行對比可知：隨著水位的抬高，擬合曲線提取的面積相對誤差絕對值逐漸減少，特別是應用區域種子生長算法對水體邊緣進行優化后，水位690m 以上提取水體面積和實測面積基本一致。各水體指數水位面積擬合誤差見表4。
2.5水庫庫容曲線構建
取 681～707m 等米水位，在考慮水體邊緣優化得到的水位-面積擬合曲線上，逐米查取相應的水面面積，采用臺柱體積計算區間容積，經累加得到車谷水庫 681m 水位以上的水位-庫容關系曲線（圖9），并和2017年實測成果進行對比，分析生成庫容曲線的精度，見表5。
表4水位面積擬合效果
Tab.4Effect of water level-area fitting
圖9車谷水庫庫容曲線提取成果
Fig.9ResultsofCheguReservoircapacitycurves
選取73景影像對應的最低水位為 681.01m ，且681～686m 無對應的影像數據，因此該水位下遙感影像提取的水體面積對低水位的庫容影響較大。雖然人工目視檢查各水體指數均將水體完全提取，但由于水位選取誤差，導致提取面積仍較實測小。因此構建的水庫庫容隨水位抬高，誤差逐漸減少。
水位 707m 時各水體指數法構建的庫容相對誤差為-1.7%～-2.5% ，和實測庫容非常接近。
3 驗證分析
武安市車谷水庫位于山區，周邊建筑物較少，應用本文提出的方法提取水庫水體面積取得了比較滿意的結果。為進一步驗證本文提出方法的普適性，對武安市大洺遠水庫進行驗證分析。大洺遠水庫位于武安市城區南側，周邊建筑物較多，屬湖泊型水庫。
武安市大洺遠水庫和車谷水庫位于同一景遙感影像中，因此同樣選用73景進行分析。基于NDWI和RWI指數法提取水體。分析時將累積頻率曲線等分為500份；DP算法處理時，至少保留200個點；考慮大洺遠水庫自庫尾至壩前有明顯的色差，生長準則的閾值相比車谷水庫較大，取20。經分析，NDWI和RWI指數法分別有效提取71和70景。
表5水位庫容提取效果
Tab.5Effectof extracting storagecapacity fromwater level
重點對典型時期（2016-06-22、2018-12-14、2019-08-16和2020-04-27）進行分析，結果見圖10、11。通過人工目視檢查，2種方法均將庫區水體柵格提取出來。2020年4月27日2種指數對庫區外左上角的水體提取存在明顯差異。主要原因是RWI將該區域提取為較小的孤立水體，并且該區域偽色域空間和庫區均值的色差大于給定的閾值，被判斷為非水體，從而未實現邊緣生長。
整體來看，本文提出采用累積頻率曲線拐點法和區域種子生長算法相結合的方法提取水體對位于城區的水庫同樣適用，具有較好的普適性。
圖10NDWI水體指數提取成果Fig.10ExtractionresultsofNDWI
圖11RWI水體指數提取成果Fig.11ExtractionresultsofRWI
4結論
中國中小型水庫眾多，傳統的庫容曲線復核工作通常需要投入大量的人力、物力和財力。隨著中小型水庫水雨情測報系統完善，積累了大量的水位觀測資料，為使用高分遙感影像提取水庫庫容曲線提供了前提條件。本文基于Sentinel-2遙感影像，針對水體指數法，提出利用累積頻率曲線拐點法自動識別最優閾值，采用區域種子生長算法對水體邊緣進行優化，進而構建了一套快速構建水庫庫容曲線的完整自動化流程。以武安市車谷水庫為研究對象，并在大洺遠水庫進行驗證。
a）水體指數法的原理是通過波段運算，經歸一化處理后，將影像二值化，以區分水體和其他地表覆蓋物。從統計學角度，基于聚類分類的思路，通過分析歸一化處理后水體指數圖的累積頻率曲線，提出最優閾值位于該曲線上由凹變凸的拐點處，進而提出采用累積頻率曲線拐點法確定最優閾值。該方法的提出為遙感水體的提取提供了新思路。
b）針對基于水體指數法，采用累積頻率曲線拐點法提取水體邊緣不連續，存在明顯系統誤差的問題，提出基于偽色域（NIRGB）的生長準則，并利用區域種子生長算法優化水體邊緣，有效提高了提取精度。
c）本文構建了一套完整地自動化流程來快速構建水庫庫容曲線，并對細節問題進行深入分析。水庫庫容曲線提取的核心為快速構建水位-面積關系曲線，重點在于保證曲線整體的精度，當影像景數足夠多時，不必苛求每景影像都具有較高的提取精度。通常由于受區域自然條件和成像時大氣條件等環境條件影響，以及水體指數本身的局限性，有時并不存在最優閾值。本文在對車谷水庫水體面積進行提取時，共選取2016年6月至2020年12月期間共73景影像。NDWI水體指數的影像利用率達到 100% ，但MNDWI水體指數利用率為 94.5% 。
d）選取比值型水體指數NDWI、MNDWI和RWI，差值型水體指數MBWI和AWEInsh，共5種指數進行分析。結果表明：頻率曲線拐點法對各類型指數最優閾值的確定均有較好的適用性和穩定性，且算法簡單，可以充分利用先驗條件，以提高搜尋效率。比值類指數水面提取效果整體優于差值型指數。在車谷水庫水體面積提取中，NDWI和RWI水體指數水體面積提取效果較好，平均誤差均為2.5% ，相對誤差分別為 -7.3%～5.8% 和 -6.2%. 4.7% R2 分別為0.99336和0.99049。誤差主要受水位觀測、水體面積提取效果和水位面積關系擬合精度等因素影響。
e）在對車谷水庫庫容曲線提取中，采用NDWI和RWI水體指數提取精度較高，庫容相對誤差范圍分別在 -9.4%～-1.8% 和 -9.3%～-1.7% ，且隨水位抬高誤差越小，高水位提取的庫容值和實測值基本一致。說明基于Sentinel-2遙感影像提取的水庫庫容曲線具有一定的精度，對大批量構建全國中小型水庫庫容曲線具有一定的參考意義。
T選取緊鄰武安市城區的大洺遠水庫進行了驗證分析，表明本文基于水體指數法提出的累積頻率曲線拐點法和區域種子生長算法相結合的方法提取水體具有較好的普適性。
參考文獻：
[1］2023年全國水利發展統計公報[M].北京：中國水利水電出版社，2024.
[2］王國興，李士鴻.應用遙感資料獲取庫區水下地形的方法研究[J].河海大學學報（自然科學版），1998，26（6）：91-94.
[3]PHILIPAT，STEPHENJW.Modeling floodplainin undationusinganintegrated GIS with radar and optical remote sensing[J].Geomorphology，1998，21（3）：295-312.
[4]陳曦，裴毅，姚幫松，等.水位庫容曲線的衛星影像測定方法研究[J].人民長江，2013，44（20）：25-28.
[5］張莉芳，潘華海，單定軍，等.基于遙感技術的柘林水庫庫容曲線復核[J].水利水電技術，2017，48（6）：1-6
[6]王金鑫，桑學鋒，常家軒，等.多源時相遙感影像數據耦合的水庫庫容曲線重構[J].測繪通報，2021（11）：42-47.
[7]吳慶雙，汪明秀，申茜，等.Sentinel-2遙感圖像的細小水體提取[J].遙感學報，2022，26（4）：781-794.
[8］徐涵秋.水體遙感指數研究進展[J].福州大學學報（自然科學版），2021，49（5）：613-625.
[9]MCFEETERSSK.The use of the normalized difference waterindex（NDWI）in the delineation ofopen water features[J].International Journal ofRemote Sensing，1996，17（7）：1425-1432.
[10］徐涵秋.利用改進的歸一化差異水體指數（MNDWI）提取水體信息的研究[J].遙感學報，2005，9（5）：589-595.
[11]WANGXB，XIE SP，ZHANGXL，et al.A robust multi-bandwaterindex（MBWI） forautomated extractionofsurfacewaterfromLandsat8 OLI imagery[J]. International Journal of Applied EarthObservation and Geoinformation，2018，68：73-91.
[12]FEYISAGL，MEILBYH，FENSHOLTR，etal.Automatedwaterextraction index：a new technique for surface water mapping usingLandsat imagery[J].Remote Sensing of Environment，2O14，140：23-35.
[13］王茹月，周航宇.直方圖雙峰法在遙感影像水體提取中的研究與改進[J].測繪通報，2023（9）：77-81.
[14]張珂，吳星宇，吳南，等.基于高分一號遙感影像的水體提取方法對比分析與改進[J].水資源保護，2024，40（4）：9-16.
[15]SEKERTEKIN A.A survey on global thresholding methods formapping open water body using Sentinel-2 Satelite imagery andnormalized difference water index[J].Archives of ComputationalMethods in Engineering，2021，28：1335-1347.
[16]MUKHERJEE A，KUMAR A A，RAMACHANDRAN P.Development of new index-based methodology for extraction ofbuilt-up area from Landsat7 imagery：comparison of performancewithSVM，ANN，and existing indices[J].IEEE TransactionsonGeoscience and Remote Sensing，2021，59（2）：1592-1603.
[17]LICM，SHAO ZF，ZHANGL，etal.A comparative analysis ofindex-based methods for impervious surface mapping using multi-Seasonal Sentinel-2 satellitedata[J].IEEEJournal ofSelectedTopics in Applied Earth Observations and Remote Sensing，2021，14：3682-3694.
[18]DOUGLASDH，PEUCKER TK.Algorithmsfor thereduction ofthenumber ofpointsrequired torepresentadigitizedlineoritscaricature［J].The Canadian Cartographer，1973，10（2）：112-122.
[19］李濤，田曉君.離散型函數拐點算法及應用[J].微計算機信息，2007，23（6）：248-249.
[20]ADAMS R，BISCHOF L.Seeded region growing[J].IEEETransactionson Pattern Analysisand Machine Intelligence，1994，16（6）：641-647.
[21]鄧孺孺，何穎清，秦雁，等.近紅外波段（900—2500nm）水吸收系數測量[J].遙感學報，2012，16（1）：192-206.
[22]李貴俊，劉正熙，游志勝，等.一種基于色差和彩色歸一化的車身顏色識別算法[J].計算機應用，2004（9）：47-49.
（責任編輯：向飛）
