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

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

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

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

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

利用遙感影像成像范圍廣的特點(diǎn)，通過(guò)懸沙反演手段獲取大范圍海域懸沙濃度場(chǎng)，可作為補(bǔ)充現(xiàn)場(chǎng)實(shí)測(cè)數(shù)據(jù)的有效手段。本文提出將遙感影像反演得到的懸沙場(chǎng)數(shù)據(jù)作為數(shù)值模擬的補(bǔ)充驗(yàn)證資料，來(lái)改進(jìn)懸沙數(shù)值模擬的精度。以?？跒乘?yàn)槔脤?shí)地水文測(cè)量資料，對(duì)懸沙遙感反演結(jié)果進(jìn)行率定，并對(duì)懸沙數(shù)值模擬結(jié)果進(jìn)行對(duì)比驗(yàn)證，檢驗(yàn)以懸沙遙感反演數(shù)據(jù)作為驗(yàn)證資料對(duì)數(shù)值模擬結(jié)果精度的影響。

1 遙感影像反演原理

1.1 遙感影像獲取

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

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

1.2 大氣校正

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

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

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

(1)

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

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

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

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

1.3 懸沙反演方法

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

(2)

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

1.4 參數(shù)率定

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

(3)

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

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

經(jīng)過(guò)對(duì)實(shí)測(cè)站點(diǎn)數(shù)據(jù)的分析，發(fā)現(xiàn)表層含沙量與垂線平均含沙量在一定的水深區(qū)域內(nèi)存在比例關(guān)系。隨所處水深不同，垂線平均含沙量為表層含沙量的1.25～1.32倍，水深越淺，倍數(shù)越大，反之倍數(shù)越小。因此，以海圖2 m、5 m、10 m的等深線作為邊界對(duì)影像數(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ò)散方程

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

(4)

式中：t為時(shí)間；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)驗(yàn)公式法或半理論方法確定；ω為泥沙沉降速度；α為泥沙沉降幾率。

2.2 計(jì)算條件

2.2.1 模型區(qū)域

?？跒呈且粋€(gè)向北敞開(kāi)的半圓形海灣，其東部為南渡江三角洲的西緣，西部為玄武巖臺(tái)地構(gòu)成的澄邁角，北為瓊州海峽水域，與廣東省雷州半島隔海相望，灣口東西向?qū)挾燃s12 km，海岸線長(zhǎng)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 水動(dòng)力條件

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

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

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

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

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

2.2.3 懸沙初始場(chǎng)

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

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

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

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

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

8-a 2#站含沙量對(duì)比曲線

8-b 8#站含沙量對(duì)比曲線圖8 2#、8#站含沙量計(jì)算值與實(shí)測(cè)值對(duì)比圖Fig.8 Comparison of calculated and measured sediment concentration at stations 2 and 8

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

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

4 結(jié)語(yǔ)

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