基于多源遙感技術的水庫特征參數反演研究

李闖　李兵　王寧

摘要：以丹江口水庫為例，采用Sentinel-1A合成孔徑雷達（SAR）和Sentinel-2A/B多光譜成像儀（MSI）聯合提取水域面積，通過建立多源傳感器（SAR和MSI）所提取水域面積差異的定量函數關系，修正二者水域面積之間的差異，并將修正后的水域面積與Sentinel-3A雷達高度計（RA）所提取的水位相結合，開展水庫蓄水量反演研究，最后以“全國水雨情網”公布的實測數據等為參考，評估水庫特征參數的反演精度。結果表明：由RA數據提取的遙測水位與實測水位相比，R2達到0.99以上，RMSE達到0.093 3，遙測水位與實測水位的平均誤差僅有9.33 cm;相對于單一MSI影像或SAR影像所提取的水域面積，二者水域面積聯合參與建立的“遙測水位-多源遙感水域面積”模型的R2和RMSE均有所優化，其中，R2分別提高了0.26%和0.60%，RMSE分別降低了2.63和1.29;由遙測水位和多源遙感水域面積結合水庫蓄水量模型反演的蓄水量與標準蓄水量相比，R2達到0.99，RMSE達到1.177，相對誤差平均為1.8%。研究結果表明：利用多源遙感技術開展水庫特征參數的反演研究，能夠更加準確地反映水庫的運行狀態。

關鍵詞：水庫特征參數; 多源遙感; 雷達高度計; 多光譜成像儀; 合成孔徑雷達

中圖法分類號： TN958

文獻標志碼： A

DOI：10.16232/j.cnki.1001-4179.2021.11.015

0引 言

水庫是陸地水資源的聚集地，水庫特征參數主要包括水庫的水位、水域面積和蓄水量等。這些參數直接反映著水庫健康運行的狀態，對區域的氣候、濕度、降水等有著重要的影響。研究水庫特征參數有助于掌握水庫的運行狀態，對攔截洪水、保障供水、調配水資源等具有積極作用，因此，開展水庫特征參數的反演研究具有十分重要的現實意義[1]。

傳統方法主要依靠人力到現場布設各種監測儀器并繪制庫區地圖等來反演水庫特征參數;單一遙感方法主要利用一種遙感技術對水庫的一種參數進行提取，其他參數仍采用現場調查的方式。李書杰[2]利用人力進行外業測量獲得水庫的水位并繪制了庫區地圖，通過對內業資料的整理，利用等高線、方格網等求得水庫的水域面積，受制圖過程中人為操作產生的精度損失影響，所反演的蓄水量只是粗略的結果。蘇業助等[3]為了科學化管理水庫，減少因外業測量和內業資料整理而造成的數據偏差，根據蓄水量平衡原理修正了蓄水量反演曲線，提高了蓄水量反演精度。孫建蕓等[4]基于我國GF-1衛星數據和實測水文數據，利用歸一化水體指數（NDWI）提取了丹江口水庫的水域面積，并根據水位、水域面積和蓄水量三者之間的相關性，分析了丹江口水庫的供水情況，是遙感數據與實測數據相結合的應用。王慶健[5]利用Sentinel-1A提取的水域面積、“全國水雨情網”公布的實測水位和標準蓄水量反演了水庫參數，并通過蓄水量平衡公式計算水庫的滲透量等，最后，對影響蓄水量的各種因素進行了分析。這些研究采集信息的時間長、工作量大、耗資多，所布設的儀器也常常會因為惡劣天氣和日常消耗而損壞，從而導致數據的缺失;一段時期內，獲取的可用遙感數據少，并且在相關模型建立之后仍需輸入實測數據才能獲得參數結果，對實測數據的依賴較大，精度相對較低。

相對于傳統方法和單一遙感方法，多源遙感方法能夠使各遙感技術優勢互補，數據源豐富，水庫特征參數反演精度高。本文利用星載的雷達高度計（Rader Altimetry，RA）提取水庫的水位，以及星載的多光譜成像儀（Multi-Spectral Instrument，MSI）和合成孔徑雷達（Synthetic Aperture Radar，SAR）聯合提取水庫的水域面積。提取水域面積利用的是衛星有效載荷生產的多源遙感數據，會存在一定的誤差。誤差來源主要有3個方面：① 不同衛星所搭載的傳感器本身會有一定誤差;② 不同衛星過境同一水庫的時間不同，會增加以小時、天等為時間單位的計算誤差;③ 不同衛星傳感器的成像機理不同，所生產的數據需要采用特定的處理方法，這些處理方法精度不一，會造成一定的誤差[6-7]。雖然可以采用軌跡間和衛星間等偏差的先驗調整，但不同衛星的修正方法使得水庫特征參數的計算變得繁瑣與復雜，因此，綜合考慮多源遙感技術在水庫特征參數上的研究并不多[8]。本文不去深究修正上述3種誤差的具體方法，而是在驗證水域面積精度的基礎上，通過建立差異函數直接對提取的水域面積結果進行修正，減少了中間過程中的復雜計算，節約了大量資源，為利用多源遙感技術開展水庫特征參數的反演研究提供了有益參考。

1研究區域與數據資料

1.1研究區域

丹江口水庫（32°36′～33°48′N、110°59′～111°49′E）位于漢江中上游，如圖1所示，是國家一級水源保護區和國家級生態文明示范區，作為我國南水北調中線工程的水源地，保障了河南、河北、北京、天津四?。ㄊ校┑纳a和生活用水[9]。

丹江口水庫從1958年9月開工，到1973年結束初期工程的建設，自建庫以來發揮了巨大的經濟、社會和生態效益。2013年8月完成對二期工程的驗收，設計壩頂高程為176.6 m，正常水位為170 m，水域面積可達1 050 km2，蓄水量可達290.5億m3，2014年12月正式向南水北調中線一期工程輸水[10-11]。

1.2數據資料

采用歐洲航天局的Sentinel系列衛星作為遙感數據源，獲取地址為：https：∥scihub.copernicus.eu/dhus/，采用“全國水雨情網”上公布的實測數據作為地面驗證的數據源，獲取地址為：http：∥xxfb.mwr.cn/，所有數據的獲取日期為2019年1月1日至2019年12月31日，其中，Sentinel-1A SAR影像29景、Sentinel-2A/B MSI影像14景和Sentinel-3A RA數據13份等。

重訪周期為12 d的Sentinel-1A衛星和重訪周期為27 d的Sentinel-3A衛星是發射電磁波的主動工作模式，所搭載的合成孔徑雷達，不受云霧和天氣的干擾，具有全天時、全天候的特點。重訪周期為5 d的Sentinel-2A/B衛星是反射光波的被動工作模式，所搭載的多光譜成像儀，可提供媲美人眼的光學數據，相對于灰度值的成像模式，其對地物目標的光譜差異性具有更加細微的捕捉能力[12-13]。

2研究方法

2.1水位反演

利用Sentinel-3A雷達高度計所生產的RA數據反演水庫水位，理想情況下，RA測量結果為衛星到水面的瞬時距離[14]，如式（1）所示。

h=cτ2（1）

式中：h為RA到水面的高度;c為光速;τ為從RA發射信號到接收信號的往返時間。

RA在測量過程中會受許多因素的影響，例如大氣、電離層、波浪等，實際計算時必須消除這些影響，主要有干對流層Rangedt、濕對流層Rangewt、電離層Rangei、固體潮Rangeset、極潮Rangept以及水面橢球高GeoEGM96的修正，如式（2）所示。受誤差影響和參考水準面的選取，該文通過Ice-1波形重跟蹤算法處理Sentinel-3A RA數據獲得以中國黃海水準面為基準的遙測水位，如式（3）所示，Range Cor為總誤差修正量，Height為水位，Altitude為衛星到地球橢球面距離，Range為衛星到水面距離。

2.2水域面積反演

2.2.1MSI影像反演水域面積

由MSI影像反演水域面積的處理流程如圖2所示，包括影像預處理、水體信息增強、水陸分割、后處理，以及水域面積計算等步驟。

其中，水體信息增強是為了加強水體與非水體之間的對比度，以便提高后續處理的精度;水陸分割用于水體提取，將影像特征相似的非水體區分開來;后處理主要用于對分割結果中可能存在的干擾信息進行去除，例如橋梁、浮島、池塘等地物的影響，并獲得完整水體輪廓。

圖3展示了Sentinel-2A/B衛星所生產的MSI影像反演水域面積過程中關鍵環節的處理結果。采用NDWI指數對圖3（a）所示的MSI影像圖以地物綠波段和近紅外波段反射率的差異性增強水體信息[15]，處理結果如圖3（b）所示。水體信息增強后的二值圖中，水體呈現較低灰度值的黑色，非水體呈現較高灰度值的灰色，具有良好的區分度。水陸分割利用決策樹將具有相似性的水體與地形陰影進行區分，結果如圖3（c）所示。圖3（c）中可以明顯看出存在非水體的干擾信息，可以通過設置連通域閾值方式對小連通域干擾物進行去除，從而獲得水體輪廓的提取結果，如圖3（d）所示。

基于提取到的水體，采用如式（4）所示的像元統計方法進行計算，可獲得水域面積Smw。

式中：M為影像中像元總數目;N為提取到的水體像元數目;Sim為影像所對應的地表區域面積，km2;Smw為最終反演的水域面積，km2。

2.2.2SAR影像反演水域面積

SAR影像反演水域面積的處理流程如圖4所示。首先要進行輻射定標、多視處理、斑點濾波、幾何校正等處理，以提高影像質量，為獲取更加精準水體信息奠定基礎[16-17]。

預處理后的Sentinel-1A SAR影像如圖5（a）所示。由細節放大圖可知，該影像具有較高分辨率。水陸分割是一個區分水體與非水體的二分類過程，模糊C均值算法（Fuzzy C-Means，FCM）是一種常用分類方法，該方法以隨機聚類中心開始，經過不斷的迭代優化完成地物類型歸屬劃分[18]。經水陸分割和小連通域去除的水體輪廓提取結果如圖5（b）所示。

采用FCM算法進行水陸分割的關鍵是確定影像中水體與非水體的臨界灰度值信息，從而獲得更加準確的聚類中心。因此，該文對SAR影像像元的灰度值進行了統計，其直方圖如圖5（c）所示。

可以明顯地看到直方圖中出現了2個波峰，第一個波峰主要為SAR影像中的水體，其像元灰度值分布在0～8內，第二個波峰以及剩余灰度值主要代表SAR影像中的非水體。灰度值8為水體與非水體的臨界值，為了避免漏提淺水水體，結合灰度統計直方圖和實踐研究，應將FCM算法的一個關鍵聚類中心設置為8。

2.3多源遙感水域面積修正研究

MSI影像的觀測周期較短，影像信息豐富，但受天氣和光照的影響較大，故MSI影像實際可用的數據較少。SAR影像可全天時、全天候工作，但其重返周期較長，加上受限于軍事保密、技術瓶頸和成本控制，獲取到合適的SAR影像數據也相對較少。針對上述情況，將MSI影像和SAR影像聯合應用于研究，可以豐富數據源、提高實驗精度。

MSI影像和SAR影像屬于異源傳感器影像，前者是被動成像模式，利用光波的反射成像，后者是主動成像模式，利用電磁波的散射成像，兩者具有不同的成像機理，出現了“同物異譜”現象。為了獲得后續更精確的蓄水量反演結果，必須構建兩者水域面積差異的定量函數關系，修正上述誤差。

在構建定量函數關系之前，需要先驗證水域面積的提取精度。水體輪廓對水域面積的提取精度有著至關重要的影響，只要驗證水體輪廓的提取精度即可說明水域面積的提取精度。利用遙感數據提取水體輪廓，手動提取法精度最高，但其效率低下。以手動提取的水體輪廓為驗證數據，通過漏警率、虛警率、準確率和Kappa系數等指標，對提取的水體輪廓進行精度驗證，其定量分析結果如表1所列?？芍?，2種影像水體輪廓的漏警率和虛警率都較低，準確率優于0.97，Kappa系數高于0.91，表明提取的水體輪廓足夠用于后續研究。

對SAR水域面積和MSI水域面積進行指數、冪次、多項式等不同函數的擬合，并選取最優函數來建立多源遙感水域面積差異的定量函數關系，所得定量分析的結果如表2所列。

根據R2越大越好、RMSE越小越好和多項式次數越低計算越簡單的定量分析原則，四次多項式函數具有最小的RMSE、較大的R2和較低的多項式次數。因此，采用四次多項式函數表示MSI影像和SAR影像所提取水域面積之間的差異。

四次多項式函數的表達式如式（5）所示。

F（x）=ax4+bx3+cx2+dx+e （5）

式中：x為水位，m;F（x）為同一水位下MSI影像和SAR影像提取水域面積的差異值，m2;a、b、c、d、e為多項式系數，擬合后的值分別為-0.000 476，0.302 6，-72.1，7 632，-0.000 030 28。

丹江口水庫位于北半球和南半球冷暖氣流的交匯處，易受副熱帶高壓北進擺動和穩定南退的影響，天氣變化劇烈，水庫上空云霧天氣的概率較大。對于選取的部分MSI影像，不可避免地存在少量的云霧，其遮擋效應在一定程度上影響了水域面積提取的精度。相比之下，SAR影像具有全天時、全天候的特點，不存在上述誤差。另外，從水體輪廓提取結果來說，SAR影像所提取水體輪廓精度要優于MSI影像。因此，考慮到上述情況，該文以SAR影像為基準，修正水域面積差異，即將MSI影像提取的水域面積值P（x），加上面積誤差修正值F（x），得到修正后的MSI水域面積S（x），如式（6）所示。修正后的MSI水域面積與SAR水域面積聯合在一起稱為多源遙感水域面積。

S（x）=P（x）+F（x）（6）

2.4水庫蓄水量反演

水庫蓄水量是無法直接觀測的量，一般通過水位和水域面積進行間接反演。該文把水庫理想化為棱臺，獲得鄰水位或者相鄰橫斷面之間的一小部分蓄水量，并把這些小蓄水量累加起來，再加上死水位對應的蓄水量就能求得整個水庫的總蓄水量[19]，如式（7）～（9）所示。

3結果與分析

3.1水位結果及分析

根據Sentinel-3A RA數據提取遙測水位的獲取日期查找到對應的實測水位，然后將遙測水位和實測水位進行對比，結果如圖6所示，日期范圍是2019年1月1日至12月31日。

由圖6可知，2019年丹江口水庫的水位變化范圍為150～167m。隨著旱季和汛季的轉換，水位呈現規律性波動，在7～10月的汛期，水位上漲迅速，11月之后將逐漸進入枯水期，水位開始慢慢下降，并且由Sentinel-3A RA提取的遙測水位與實測水位的波動趨勢相吻合。

采用R2和RMSE對遙測水位的精度進行評價，得到R2=0.99，RMSE=0.0933，表明提取的遙測水位具有很高的精度，遙測水位與實測水位的平均誤差僅有9.33 cm。

3.2水域面積結果及分析

利用星載MSI影像和星載SAR影像分別提取的MSI水域面積和SAR水域面積，結果如圖7～8所示。

由圖7～8可知：2019年丹江口水庫的水域面積在5月初達到最小值，7月開始明顯增大，在11月中旬達到最大值，然后再次回落，最大水域面積與最小水域面積相差約250 km2。隨著時間推移，整體上MSI水域面積和SAR水域面積變化趨勢較為吻合，MSI水域面積和SAR水域面積的準確率分別達到了97.95%和98.95%。

3.3建立“水位-水域面積”模型

由2.4小節可知：只要獲得水位和水域面積就能在不依靠地面實測數據的情況下獲得蓄水量。對于水位和水域面積，已知其中一個參數，想要快速獲得另外一個參數，需要建立“水位-水域面積”關系曲線。該關系曲線建立之后，再通過多源遙感技術獲得水位或水域面積任何一個參數，就能獲得包括蓄水量在內的另外兩個參數。

以多源遙感技術獲取的遙測水位、MSI水域面積、SAR水域面積以及修正后的多源遙感水域面積分別建立“遙測水位-MSI水域面積”關系曲線、“遙測水位-SAR水域面積”

關系曲線和“遙測水位-多源遙感水域面積”關系曲線。以“遙測水位-多源遙感水域面積”關系曲線為例，建立過程如下：

首先構建“遙測水位-多源遙感水域面積”關系曲線的指數函數、冪函數和多項式函數的數學模型，由于四次及四次以上多項式函數出現過擬合現象，選擇指數函數、冪函數和四次以下多項式函數對擬合效果作進一步驗證，如表3所列。

由表3可知：指數函數的R2最小，RMSE最大;冪函數與指數函數的R2和RMSE最接近;三次多項式函數的R2最大，但RMSE不是最小，其R2只比二次多項式函數的R2多0.0003，二次多項式函數的RMSE最小。根據R2越大越好，RMSE越小越好的定量分析原則，二次多項式函數最適合描述遙測水位與多源遙感水域面積之間的關系。

得到“遙測水位-多源遙感水域面積”模型如式（10）所示，關系曲線如圖9（a）所示;同理可得“遙測水位-MSI水域面積”模型如式（11）所示，關系曲線如圖9（b）所示;“遙測水位-SAR水域面積”模型如式（12）所示，關系曲線如圖9（c）所示。

“水位-水域面積”關系曲線建立之后，各模型最優函數的R2和RMSE如表4所列。

由表4可知：各模型本身也已經具有較高的R2和較低的RMSE，相比單獨使用MSI影像或SAR影像，聯合使用兩種遙感影像后，“遙測水位-多源遙感水域面積”模型的R2分別提高了0.26%和0.60%，RMSE分別降低了2.63和1.29，即R2和RMSE均有所優化。

3.4蓄水量結果及分析

利用多源遙感技術反演得到水庫的遙測水位和多源遙感水域面積，進而得到水庫的蓄水量，將反演的蓄水量與“全國水雨情網”上公布的標準蓄水量進行對比，對比結果如圖10所示，日期范圍是2019年1月1日至12月31日。

由圖10可知，反演蓄水量和標準蓄水量的波動趨勢具有較高吻合度。在7～10月的汛期，相關地區降雨較多，丹江口水庫的蓄水量增長迅速，在11月17日左右達到最高峰，然后開始慢慢減少，逐漸進入枯水期，在5月初達到最低值。整體上，隨著旱季和汛季的轉換，其蓄水量也呈現規律性的變化。

采用R2和RMSE進行定量分析。分析結果表明，R2=0.99，說明反演蓄水量與標準蓄水量具有極高的吻合度，能夠真實地反映水庫蓄水量的變化情況;另外，RMSE=1.177，表明反演的水庫蓄水量標準差在1.177億m3以內，相對誤差平均為1.8%。

4結 語

本文以丹江口水庫為例，利用多源遙感技術反演了水庫參數，其中，利用RA數據提取了水庫的水位，利用多源傳感器生產的MSI影像和SAR影像提取水庫的水域面積。通過構建多源傳感器水域面積之間的定量函數關系，修正了二者水域面積之間的誤差，最后，基于修正后多源遙感水域面積，結合蓄水量模型反演水庫的蓄水量。

實驗結果表明，利用多源遙感技術反演的水庫特征參數能夠真實地反映水庫健康運行的狀態。在水庫特征參數反演過程中，除了利用地面實測數據對反演結果進行驗證之外，不再依賴實測數據的支持就能達到對水庫特征參數的反演，降低了成本，節約了資源，有效補充了實測數據缺失情況下的研究工作。同時，利用多源遙感技術反演水庫特征參數能夠達到對水庫異常預警和日常管理的最大效益化，特別是對堰塞湖或出現大面積違法的圍湖造田等的預警。該文關于水庫特征參數的反演方法不僅能夠克服當前傳統方法和單一遙感方法的不足，也可以進一步推廣到其它大型水庫特征參數的反演工作當中。未來還可以考慮采用更多類型的遙感數據和更加優良的水庫蓄水量模型等參與反演。

參考文獻：

[1]周建中，頓曉晗，張勇傳.基于庫容風險頻率曲線的水庫群聯合防洪調度研究[J].水利學報，2019，50（11）：1318-1325.

[2]李書杰.水庫庫容的量算精度[J].東北水利水電，1997（7）：34-36.

[3]蘇業助，洪為善.根據水量平衡原理修正水庫庫容曲線方法[J].人民長江，1999，30（12）：39-41.

[4]孫建蕓，袁琳，王新生，等.基于GF_1衛星的丹江口水庫水面面積-蓄水量-水位相關性研究[J].南水北調與水利科技，2017，15（5）：89-96.

[5]王慶健.基于SAR監測的水庫主要參數計算系統的設計與實現[D].開封：河南大學，2018.

[6]BONNEFOND P，EXERTIER P，LAURAIN O，et al.Absolute calibration of Jason-1 and Jason-2 altimeters in Corsica during the for- mation flight phase[J].Marine Geodesy，2010，33（S1）：80-90.

[7]MERTIKAS S P，IOANNIDES R T，TZIAVOS I N，et al.Statistical models and latest results in the determination of the absolute bias for the radar altimeters of Jason satellites using the gavdos facility[J].Marine Geodesy，2010，33（S1） ：114-149.

[8]王文種，黃對，劉九夫，等.基于Landsat與Sentinel-3A衛星數據的當惹雍錯1988～2018年湖泊水位-水量變化及歸因[J].湖泊科學，2020，32（5）：1552-1563.

[9]王元超，王旭，雷曉輝，等.丹江口水庫入庫徑流特征及其演變規律[J].南水北調與水利科技，2015，13（1）：15-19.

[10]吳川，張玉龍，許秀貞，等.基于Landset TM/ETM和HJ-1A/B影像的丹江口水庫水域變化監測研究[J].長江流域資源與環境，2013，22（9）：108-114.

[11]尹志杰，王容，趙蘭蘭，等.丹江口水庫以上流域夏秋汛期來水及相關性分析[J].水文，2019，39（4）：74-77.

[12]楊魁，楊建兵，江冰茹.Sentinel-1衛星綜述[J].城市勘測，2015（2）：26-29.

[13]田穎，陳卓奇，惠鳳鳴，等.歐空局哨兵衛星Sentinel-2A/B數據特征及應用前景分析[J].北京師范大學學報（自然科學版），2019，55（1）：57-65.

[14]郭華東.雷達對地觀測理論與應用[M].北京：科學出版社，2000.

[15]GAO B C.NDWI—A normalized difference water index for remote sensing of vegetation liquid water from space[J].Remote Sensing of Environment，1995，58（3）：257-266.

[16]潘習哲.星載SAR圖像處理[M].北京：科學出版社，1996.

[17]王超.SAR圖像處理新進展[J].中國圖象圖形學報，2009，14（1）：1-2.

[18]BEZDEK J C，EHRLICH R，FULL W.FCM：The fuzzy c-means clustering algorithm[J].Computers & Geoences，1984，10（2-3）：191-203.

[19]張國慶.青藏高原湖泊變化遙感監測及其對氣候變化的響應研究進展[J].地理科學進展，2018，（2）：214-223.

（編輯：黃文晉）

Abstract：Take Danjiangkou Reservoir as an example，the water area was extracted by using the Sentinel-1A Synthetic Aperture Radar （SAR） and Sentinel-2A/B Multi-Spectral Imager （MSI），the difference of water areas by SAR and by MSI was corrected by establishing a quantitative function relation of these areas，and the revised water area was combined with the inverted water level of Sentinel-3A Radar Altimeter （RA） to calculate the reservoir volume.Finally，the retrieval precision of the reservoir parameters was evaluated according to the measured data published by the National Water and rain regime network.The results showed that R2 was above 0.99 and RMSE was 0.0933，which meant that the average error between inverted water level and measured water level was only 9.33 cm.Compared with the water area extracted from single MSI image or SAR image，the R2 and RMSE of the "remote sensing water level-Multi-Source Remote Sensing water area" model were optimized，especially R2 increased by 0.26% and 0.60% ，while RMSE decreased by 2.63 and 1.29，respectively.Compared with the standard water storage，for the novel model，R2 reached 0.99 and RMSE reached 1.177，respectively，and the relative error was 1.8% on average.The result shows that the reservoir characteristic parameters inversion using multi-source remote sensing technology can more accurately reflect the operation of the reservoir.

Key words：reservoir characteristic parameters;multi-source remote sensing;Rader Altimetry;Multi-Spectral Instrument;Synthetic Aperture Radar