利用遙感反演數據提高懸沙數值模擬精度的方法研究

杜 洋，李文拓，沈維娜，李懷遠

(1.中海油海南能源有限公司，海口 570105；2.交通運輸部天津水運工程科學研究所 工程泥沙交通行業 重點實驗室，天津 300456)

基于水動力理論的泥沙數值模擬技術經過多年的發展，已形成模式化的計算理論和方法[1-2]，然而模擬精度不理想，計算結果一般僅能提供工程設計參考。一方面由于水動力及泥沙運動過程本身極為復雜，理論與實際情況存在差異，另一方面驗證資料的豐富性和準確性對模擬結果精度有較大影響[3]。然而，獲取海上水沙實測數據的成本較高，通過開展大范圍水文測量獲取驗證數據一般不具備現實可行性。現行做法大多通過特征點位的短期測量獲取驗證資料，或利用以往類似工況下的測量數據作為模型驗證資料，或基于專家經驗在模型中融入部分可調參數[4]。實測資料的豐富性已成為限制泥沙數值模擬精度進一步提升的主要因素。

隨著遙感技術的快速發展，遙感影像在海洋監測領域的應用日益廣泛[5-6]。國內外有關水體懸浮泥沙質量濃度遙感定量反演的研究已有大量成果，例如溫令平[7]采用TM影像的多波段組合構建了懸浮泥沙遙感定量模式，黎夏等[8]建立了懸浮泥沙遙感反演統一模式。總結前人對懸沙反演的研究成果，大多采用不同的遙感數據源，開展地面同步測量試驗，并基于理論分析法或統計經驗分析法的遙感定量反演模式來計算水體中懸浮泥沙濃度[9-11]。一般而言，由于近岸水體的光學特性遠比大洋水體復雜，不同水色成分的吸收和散射光譜有重疊區域，組合起來的光譜特性是非線性的、復雜多變的，因此懸浮泥沙反演的通用模式和高精度模式尚未真正建立起來，實際應用中大多根據率定資料選擇最佳擬合的模型和參數[12-13]。

利用遙感影像成像范圍廣的特點，通過懸沙反演手段獲取大范圍海域懸沙濃度場，可作為補充現場實測數據的有效手段。本文提出將遙感影像反演得到的懸沙場數據作為數值模擬的補充驗證資料，來改進懸沙數值模擬的精度。以海口灣水域為例，利用實地水文測量資料，對懸沙遙感反演結果進行率定，并對懸沙數值模擬結果進行對比驗證，檢驗以懸沙遙感反演數據作為驗證資料對數值模擬結果精度的影響。

1 遙感影像反演原理

1.1 遙感影像獲取

表1 實驗數據參數介紹Tab.1 Introduction of experimental data parameters

本文使用的遙感衛星影像數據是Landsat 8 OLI 1級矯正產品。Landsat 8衛星由美國國家航空航天局(NASA)于2013年2月11號發射成功，該衛星攜帶有兩個主要載荷：OLI陸地成像儀和TIRS熱紅外傳感器。OLI傳感器主要包括可見光和近紅外波段，其空間分辨率為30 m，其波段參數如表1所示。

1.2 大氣校正

為了從衛星影像中獲取真實的地表反射率信息，需要對影像數據準確消除大氣效應，即大氣校正，目的是消除大氣氣溶膠對雜波反射的影響，獲得真實的表面反射率。

利用基于6S模型的大氣校正查找表進行離線計算[14]，可有效提高大氣校正的效率。該查找表包含不同的大氣氣溶膠光學厚度、太陽天頂角、觀測天頂角和相對方位角，可提供大氣校正系數，實現Landsat 8數據的大氣校正。在去除圖像中大氣效應的同時使用絕對大氣校正方法，將DN值轉換為表面反射率。

在 6S 模型中，大氣校正基于以下方程

(1)

式中：xa，xb，xc分別為3個大氣校正系數；Li為傳感器第i波段數據的輻射值。

1-a 左側為校正前 1-b 右側為校正后圖1 大氣校正前后對比圖Fig.1 Comparison map before and after atmospheric correction

輸入參數時，在6S模型中選擇了預定義的標準模式作為輸入。大氣模型選擇低緯度熱帶模式，氣溶膠模式選擇近海型，并假定地表具有均勻的蘭伯特反射特征。

海口灣海域遙感影像大氣校正后影像紅色波段數據如圖1所示。

1.3 懸沙反演方法

基于多光譜遙感影像的表層懸浮泥沙反演經過多年發展已有大量研究成果，一般認為紅、綠波段對二類水體懸浮顆粒濃度變化較為敏感，反演模型構建多選擇紅綠波段組合的指數模型、多項式模型，并依據實測數據對模型系數進行率定。本研究即采用紅色和綠色兩波段數據計算作為反演模型構建參數，并進行懸沙反演分析，反演模型可表示為

(2)

式中：S為含沙量，kg/m3；Red為紅色波段表面反射率；Green為綠色波段表面反射率；a、b為反演模型待定參數，通過同步光譜測量數據對待定參數進行率定。

1.4 參數率定

2018年9月22日，在海口灣海域現場采用ASD光譜測量儀和濁度計，分別對光譜和含沙量進行了測量。共設置了8個同步測量點，點位布置如圖2所示。其中ASD光譜儀的波長范圍為350～2 500 nm，參考反射率板為50%。為了消除隨機噪聲的影響，在Landsat 8衛星過境前后1 h內完成光譜測量得到各測點的平均表面光譜數據，同時采集太陽天頂角。含沙量觀測為每隔10 min記錄一次，共連續觀測了25 h。

(3)

1.5 表層與垂線平均含沙量對應關系

從遙感影像上直接反演得到的結果為表層含沙量場，為泥沙數學模型中模擬的垂線平均含沙量的運動過程，為了使反演數據符合數學模型的使用要求，從實測站點數據中找出表層含沙量與垂線平均含沙量的關系，從而將表層含沙量反演結果轉換為垂線平均含沙量反演結果，以便用于模型驗證。

經過對實測站點數據的分析，發現表層含沙量與垂線平均含沙量在一定的水深區域內存在比例關系。隨所處水深不同，垂線平均含沙量為表層含沙量的1.25～1.32倍，水深越淺，倍數越大，反之倍數越小。因此，以海圖2 m、5 m、10 m的等深線作為邊界對影像數據的水域分化，換算比例分別取：1.32倍(0～2 m)、1.3倍(2～5 m)、1.28倍(5～10 m)以及1.25倍(>10 m)，從而將表層含沙量換算為垂線平均含沙量。

2 懸移質泥沙數值模型

2.1 泥沙擴散方程

泥沙運動根據運動規律的不同可分為懸移質和推移質，本文主要模擬垂線平均含沙量的運動過程，即以懸移質運動過程為主要模擬對象。懸移質的二維泥沙擴散方程如下

(4)

式中：t為時間；x、y為與靜止海面重合的直角坐標系坐標；u、v分別為x、y方向的垂線平均流速；h為基準面到床面的距離；ζ為基準面到自由水面的距離；S為水體含沙量；Dx、Dy分別為x、y方向的泥沙擴散系數；FS為泥沙源匯函數或床面沖淤函數，按下面方法確定

FS=αω(S-S*)

(5)

式中：S*為水體的挾沙力，一般采用經驗公式法或半理論方法確定；ω為泥沙沉降速度；α為泥沙沉降幾率。

2.2 計算條件

2.2.1 模型區域

海口灣是一個向北敞開的半圓形海灣，其東部為南渡江三角洲的西緣，西部為玄武巖臺地構成的澄邁角，北為瓊州海峽水域，與廣東省雷州半島隔海相望，灣口東西向寬度約12 km，海岸線長20.5 km，沿岸分布著2～22 m高的沙壩，屬于沙質海岸。灣內水深為2～6 m，海灣南北向縱深自灣頂到灣口5 m等深線的距離為4.5 km，至10 m等深線的距離約6 km。海灣水域的面積約55 km2。

本次模擬的區域范圍為20.095 8°N～20.333 3°N、110.2°E～110.333 3°E，區域位置如圖4所示。

2.2.2 水動力條件

本文采用TK-2D軟件模擬研究區域的水動力變化情況，其中潮汐潮流邊界由全球潮汐模型及海南島南部海域潮汐數學模型提供，波浪作用下的輻射應力值由第三代近岸波浪SWAN模型提供。通過TK-2D軟件中的波浪、潮流聯合作用模型，將波浪作用產生的輻射應力加到潮流動量控制方程中來考慮波浪對潮流的影響。

采用2018年7月～8月12個測站連續25 h的潮流、潮位測量數據作為驗證資料，測點位置如圖5所示。

以1#測點為例，觀測值與模擬值的潮流、潮位驗證結果如圖6所示。

圖7 表層含沙量反演結果圖(2018-09-22)Fig.7 Inversion results of surface sediment concentration(2018-09-22)

由潮位驗證曲線可見，潮位計算值和實測值在量值、位相上吻合良好。本次驗證高低潮時間的潮位相位偏差都在0.5 h以內，高、低潮位值偏差亦基本在10 cm以內，說明數學模型模擬的研究區域潮波運動與天然潮波運動基本相似，能夠反映本海域內潮波傳遞和潮波變形。由潮流驗證曲線可見，各驗證點計算流速和實測資料基本吻合，大小潮的流速及流向的周期性變化刻畫較為準確，漲、落潮段平均流速偏差小于10%，計算流速過程線的形態與實測值情況基本一致，能夠客觀反映項目附近海域的潮流運動規律。

2.2.3 懸沙初始場

為加入懸沙反演數據作為驗證資料對數值模擬結果的影響，考慮兩種工況。

工況一：直接利用遙感影像成像時刻的8個實測站實測數據作為驗證擬合資料。

工況二：依據率定后的遙感影像反演公式，對2018年9月22日Landsat 8 OLI數據進行表層含沙量分布場反演，結果如圖7所示。進一步根據海域表層含沙量與垂線平均含沙量的換算關系得到垂線平均含沙量反演數據，作為數值模擬的驗證擬合資料。

3 模擬結果及精度分析

基于上述懸移質泥沙數值模型，在水動力條件不變的情況下，采用兩組工況條件對懸沙運動過程進行了模擬。通過對遙感影像成像時刻后24 h內的懸沙運移變化過程的計算，得到8個實測點處含沙量模型計算值，以外側2#站和近岸8#站為例，同時刻含沙量模型計算值與實測值對比曲線如圖8所示。可見各驗證點含沙量計算值和實測資料吻合較好，計算的含沙量變化趨勢與實測值一致，該海域整體含沙量較小，隨漲、落潮略有變化。

8-a 2#站含沙量對比曲線

8-b 8#站含沙量對比曲線圖8 2#、8#站含沙量計算值與實測值對比圖Fig.8 Comparison of calculated and measured sediment concentration at stations 2 and 8

不同觀測時刻工況一、工況二計算值與實測值的對比結果見表2。由于各測站每個時刻的實測值均不相同，為便于比較計算值精度，以同位置、同時刻的3個值(實測值、工況1計算值、工況2計算值)為一組，進行歸一化處理，然后通過計算標準差對工況1和工況2提取的含沙量模擬結果精度進行比較，其中工況1的標準差為0.318，工況2標準差為0.231。可見通過融合遙感影像反演數據作為數值模擬的初始含沙量場，可以提高數值模擬結果的精度。

表2 各時刻各點位含沙量模型計算值與實測值對比(以成像時刻后6時、12時、18時、24時為例)Tab.2 Comparison of calculated and measured sediment concentration at each point at each time kg·m-3

4 結語

懸沙數值模擬結果的精度受驗證資料的影響較大，由于一般的常規現場觀測只能獲得離散的少量觀測點數據，很難反映工程大范圍實際含沙量分布和變化過程，限制了數值模擬精度的提升。衛星遙感影像具有周期性強、圖像直觀的特點，采用遙感定量反演技術能迅速地獲得大面積海域懸浮泥沙含量，是補充和豐富數值模擬驗證資料的有力手段。本文通過將海口灣水域的懸沙遙感反演結果作為驗證資料應用到泥沙數值模擬過程中，提升了數值模擬結果精度，驗證了本文所提出方法的可行性，對于泥沙數值模擬的工程應用具有一定參考價值。然而由于每個海區的水色、泥沙特征不同，在實際工程應用中遙感反演數據并不能完全取代現場測量，而應針對具體應用場景進行補充觀測，以便更好地匹配率定參數。本文提出的方法在不同海區的應用模式值得進一步深入研究。