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

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

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

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

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

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

1 遙感影像反演原理

1.1 遙感影像獲取

表1 實驗數(shù)據(jù)參數(shù)介紹Tab.1 Introduction of experimental data parameters

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

1.2 大氣校正

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

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

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

(1)

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

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

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

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

1.3 懸沙反演方法

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

(2)

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

1.4 參數(shù)率定

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

(3)

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

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

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

2 懸移質(zhì)泥沙數(shù)值模型

2.1 泥沙擴散方程

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

(4)

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

FS=αω(S-S*)

(5)

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

2.2 計算條件

2.2.1 模型區(qū)域

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

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

2.2.2 水動力條件

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

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

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

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

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

2.2.3 懸沙初始場

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

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

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

3 模擬結果及精度分析

基于上述懸移質(zhì)泥沙數(shù)值模型，在水動力條件不變的情況下，采用兩組工況條件對懸沙運動過程進行了模擬。通過對遙感影像成像時刻后24 h內(nèi)的懸沙運移變化過程的計算，得到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。可見通過融合遙感影像反演數(shù)據(jù)作為數(shù)值模擬的初始含沙量場，可以提高數(shù)值模擬結果的精度。

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

4 結語

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