杜 洋，李文拓，沈維娜，李懷遠(yuǎn)

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

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

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

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

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ù)準(zhǔn)確消除大氣效應(yīng)，即大氣校正，目的是消除大氣氣溶膠對雜波反射的影響，獲得真實的表面反射率。

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

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

(1)

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

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

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

?？跒澈Ｓ蜻b感影像大氣校正后影像紅色波段數(shù)據(jù)如圖1所示。

1.3 懸沙反演方法

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

(2)

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

1.4 參數(shù)率定

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

(3)

1.5 表層與垂線平均含沙量對應(yīng)關(guān)系

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

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

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

2.1 泥沙擴(kuò)散方程

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

(4)

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

FS=αω(S-S*)

(5)

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

2.2 計算條件

2.2.1 模型區(qū)域

海口灣是一個向北敞開的半圓形海灣，其東部為南渡江三角洲的西緣，西部為玄武巖臺地構(gòu)成的澄邁角，北為瓊州海峽水域，與廣東省雷州半島隔海相望，灣口東西向?qū)挾燃s12 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ù)學(xué)模型提供，波浪作用下的輻射應(yīng)力值由第三代近岸波浪SWAN模型提供。通過TK-2D軟件中的波浪、潮流聯(lián)合作用模型，將波浪作用產(chǎn)生的輻射應(yīng)力加到潮流動量控制方程中來考慮波浪對潮流的影響。

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

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

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

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

2.2.3 懸沙初始場

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

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

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

3 模擬結(jié)果及精度分析

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

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

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

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

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

4 結(jié)語

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