杭州灣三維懸浮泥沙輸運模型初始場 的伴隨法反演研究*

杜云飛 張繼才 王道勝

(1. 浙江大學海洋學院 舟山 316000； 2. 中國地質(zhì)大學(武漢)海洋學院 海洋地質(zhì)資源湖北省重點實驗室 武漢 430074； 3. 南方海洋科學與工程廣東省實驗室(廣州) 廣州 511458)

由于徑流輸入和局部泥沙再懸浮， 近岸、河口海域通常具有高懸浮泥沙濃度的特征(Yanget al，2016)。杭州灣作為典型的強潮河口灣， 具有潮流急， 潮差大， 水體含沙量高的顯著特征。大量懸浮泥沙的存在和大范圍含沙量時間、空間變化的存在， 使杭州灣懸浮泥沙分布比較復雜(中國海灣志編纂委員會， 1992)。因此， 準確了解杭州灣海域懸浮泥沙的時空分布和輸運規(guī)律， 不僅具有重要的科學價值， 還對杭州灣海域及近岸地區(qū)水質(zhì)、地貌、生態(tài)環(huán)境、海岸工程、港口建設(shè)以及灘涂開發(fā)等具有重要意義。數(shù)值模擬具有成本低、預(yù)報性強、可移植性強、定量的特點， 可有效獲取懸浮泥沙的時間和空間分布(Wanget al，2018b)， 從而彌補傳統(tǒng)現(xiàn)場觀測和衛(wèi)星遙感觀測的不足。因此， 懸沙輸運模型已經(jīng)成為研究泥沙輸運的一個有力工具(Papanicolaouet al， 2008)。

隨著對泥沙輸運機制認知的深入以及科學計算能力的不斷進步， 懸沙輸運模型從一維發(fā)展到了三維。眾多學者利用三維懸沙輸運模型開展懸浮泥沙的研究工作(Normant， 2000； Dufoiset al， 2014； Chenet al， 2015)。三維懸沙輸運模型對初始場比較敏感(Wanget al，2018a)， 即 懸 浮 泥 沙 濃 度(suspended sediment concentration， SSC)在數(shù)值模擬初始時刻的空間分布直接影響數(shù)值模擬的精度， 初始場的空間分布與實際SSC 的空間分布越接近， 相應(yīng)的數(shù)值模擬精度越高， 反之則越低。因此， 對于短期的懸浮泥沙數(shù)值模擬而言， 初始場的準確給定至關(guān)重要(Leeet al，2007； Yanget al， 2016)。到目前為止， 在懸浮泥沙輸運的數(shù)值模擬中， 初始場的確定方法主要有四種： (1)假定在整個模型計算區(qū)域內(nèi)為常數(shù)分布， 以模型運行穩(wěn)定時的結(jié)果作為最終的模擬結(jié)果(Huet al，2009； Gourgueet al， 2013)； (2)根據(jù)現(xiàn)場實測水文資料插值得到(曹振軼等，2002； Duet al， 2010)； (3)由經(jīng)驗、半經(jīng)驗的公式確定(Woolnoughet al， 1995； Schwabet al， 2000)； (4)用衛(wèi)星遙感觀測資料獲得(Chenet al， 2010； Ramakrishnanet al， 2012； Yanget al， 2016)。前三種方法與實際SSC 的空間分布偏差較大， 因而數(shù)值模擬的精度較差； 隨著衛(wèi)星遙感觀測手段的不斷發(fā)展， 初始場的確定更多依賴于衛(wèi)星遙感數(shù)據(jù)及與之相關(guān)的反演方法和經(jīng)驗公式。衛(wèi)星遙感反演方法雖然能夠獲取海洋表層大面積、連續(xù)的SSC， 但由于受到云覆蓋或惡劣天氣條件的限制， 衛(wèi)星遙感數(shù)據(jù)存在不同程度的缺測， 這在一定程度上限制了衛(wèi)星遙感方法的應(yīng)用范圍和準確度。

伴隨同化方法是建立在嚴格的數(shù)學基礎(chǔ)上的一種有效的四維變分同化技術(shù)， 它將所要解決的實際問題作為最小值問題來解決， 流體力學方程組及其初邊值條件作為約束條件， 使得根據(jù)具體問題設(shè)計的代價函數(shù)達到最小(Sasaki， 1970)。伴隨同化方法能夠充分利用已有的觀測資料， 實現(xiàn)觀測資料和數(shù)值模型的有機結(jié)合， 得到更加符合實際的模擬結(jié)果， 并能對模型重要參數(shù)如初始場、邊界條件等進行優(yōu)化反演。目前為止， 伴隨同化方法已廣泛應(yīng)用于海表溫度數(shù)值模型(馬繼瑞等， 2002)、風暴潮預(yù)報模型(Penget al， 2006)、海洋生態(tài)系統(tǒng)動力學模型(Wanget al， 2013)、海洋環(huán)流模型(Songet al， 2016)以及懸沙輸運模型(Wanget al， 2018b)等海洋數(shù)值模型中初始場的優(yōu)化反演， 并取得了良好的應(yīng)用效果。基于此背景， 本文利用Wang 等(2018a)構(gòu)建的一個三維懸沙輸運伴隨同化模型， 針對模型初始場的優(yōu)化反演進行一系列孿生實驗和實際實驗， 孿生實驗驗證伴隨同化方法優(yōu)化反演模型初始場的有效性并定量評估該模型數(shù)據(jù)同化的能力； 實際實驗則根據(jù)孿生實驗的結(jié)論同化杭州灣海域典型的小潮時期和大潮時期的GOCI衛(wèi)星遙感資料所得表層SSC 數(shù)據(jù)， 獲得模型初始場的最優(yōu)估計。

1 模型介紹

1.1 正向模型

三維懸沙輸運模型的控制方程(Wanget al，2018a)：

其中，C是懸浮泥沙濃度，t是時間，x、y是水平坐標，σ是垂向坐標，H是包含靜水深和海面起伏的總水深，u、v和w分別是x、y和σ方向的流速，ws為懸浮泥沙的沉降速度，KH和KV分別為水平擴散系數(shù)和垂向擴散系數(shù)。

初始條件：

其中，C0是初始場。

邊界條件：

其中，是垂直于邊界的外法向量，Cobc是流入開邊界處的懸浮泥沙濃度，B1代表陸地邊界，B2代表海洋開邊界，B3代表流入開邊界，B4代表海面邊界，B5代表底邊界，E和D分別是侵蝕率和沉降率， 二者的計算公式如下(Partheniades， 1965)：

其中，M0是再懸浮率，τb是底部剪切應(yīng)力值，τce是侵蝕的臨界剪切應(yīng)力值，τcd是沉降的臨界剪切應(yīng)力值，C1是近海底處的懸浮泥沙濃度值。

正向模型的離散形式和求解方式詳見Wang 等(2018a)， 此處不再贅述。

1.2 伴隨模型

根據(jù)拉格朗日乘子法， 引入伴隨變量λ， 令拉格朗日函數(shù)關(guān)于各個變量和模型參數(shù)的梯度值為0， 可以推導出伴隨模型以及代價函數(shù)關(guān)于初始場的梯度表達式， 具體的推導過程和差分格式詳見Wang 等(2018a)， 此處不再贅述。

1.3 優(yōu)化算法

由伴隨模型求得代價函數(shù)關(guān)于初始場的梯度后， 可以利用優(yōu)化算法通過迭代獲得最優(yōu)初始場， 迭代公式可以表示為：

本文中，αn是正值， 一般取為常數(shù)， 但調(diào)整步長過大或過小都會對反演結(jié)果產(chǎn)生影響。調(diào)整步長過大時， 在極小值附近會造成約化代價函數(shù)“過調(diào)整”的問題， 產(chǎn)生“鋸齒效應(yīng)”； 調(diào)整步長過小時， 會造成收斂速度慢， 計算效率低的問題。與以往調(diào)整步長取為常數(shù)不同， 本文中調(diào)整步長依賴于初始場的空間分布， 取為空間分布的形式， 即將αn取為初始場的第n個迭代步迭代結(jié)果的3%(公式13)， 這樣能夠保證迭代收斂的穩(wěn)定性， 且不需要人為地憑經(jīng)驗給定調(diào)整步長， 可以避免由調(diào)整步長過大或者過小引起的上述問題， 同時避免需要反復調(diào)試才能獲得最佳的調(diào)整步長， 提高了計算效率。

確定dn有許多可供選擇的優(yōu)化算法， 本文選擇并比較了其中的三種， 分別是最速下降法(the gradient descent algorithm， GD)、共軛梯度法(the conjugate gradient algorithm， CG)和有限記憶 BFGS 法(the limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm， L-BFGS)。

在GD 中， 搜索方向為代價函數(shù)關(guān)于模型初始場的負梯度方向：

其中，dn為搜索方向，gn為代價函數(shù)關(guān)于模型初始場的梯度方向， 根據(jù)βn的不同求法， CG 又可以分為許多種類， 其中PRP、HS、LS 方法的數(shù)值表現(xiàn)較好， 這三種方法中βn的計算公式如下：

(1) PRP 方法(CG-PRP， Polaket al，1969； Polyak， 1969)：

在初始場優(yōu)化反演的迭代過程中， 迭代收斂標準可以有很多， 當達到某一收斂標準時迭代終止。收斂標準可以是最后兩個迭代步的代價函數(shù)值充分接近或小于某個值， 也可以是最后兩個迭代步的代價函數(shù)關(guān)于初始場的梯度值充分接近或小于某值等。本文的孿生實驗中， 迭代收斂標準是最后兩個迭代步的約化代價函數(shù)值之差的絕對值小于3.0×10-5， 且最大迭代步設(shè)置為200， 即：

2 研究區(qū)域、模型設(shè)置及數(shù)據(jù)來源

本文的研究區(qū)域為杭州灣海域(圖1)， 經(jīng)緯度范圍為29.902°—30.904°N、120.602°—121.806°E。背景流場數(shù)據(jù)由FVCOM(Chenet al，2003)計算得到， 與 Wang 等(2018b)相同。模型的水平分辨率為0.006°×0.006°； 垂向分為四層， 與輸出背景流場數(shù)據(jù)的數(shù)值模型FVCOM 一致。正向模型的時間步長為180s， 總的模擬時間是31h， 小潮時期數(shù)值模擬實驗的初始時刻是2011 年6 月26 日0：28， 結(jié)束時刻是2011 年6 月27 日7：28， 大潮時期數(shù)值模擬實驗的初始時刻是2015 年3 月24 日0：28， 結(jié)束時刻是2015年3 月25 日7：28。河流流入開邊界條件在小潮時期取為1.5×10-1kg/m3(中華人民共和國水利部， 2011)， 在大潮時期取為1.0×10-1kg/m3(中華人民共和國水利部， 2015)。根據(jù)唐建華(2007)的研究， 一個潮周期內(nèi)杭州灣海域絮凝的平均沉速的量級為10-4m/s， 因此沉降速度取為1.0×10-4m/s。參照Hu 等(2009)在模擬杭州灣懸浮泥沙輸運時再懸浮率的取值， 再懸浮率取為3.0×10-6m/s， 沉降和侵蝕的臨界切應(yīng)力均取為0.4N/m2?；赪ang 等(2018a)中的尺度分析， 水平擴散系數(shù)取為80m2/s， 垂向擴散系數(shù)在小潮時期取為3.0×10-3m2/s， 在大潮時期取為3.0×10-2m2/s。以上模型參數(shù)的具體描述和取值如表1 所示。

本文采用由靜止水色衛(wèi)星GOCI 的觀測數(shù)據(jù)反演所得SSC 作為同化數(shù)據(jù)， 由He 等(2013)提供。GOCI 衛(wèi)星能夠?qū)Τ毕?、赤潮、河流羽流、泥沙輸運等短期的區(qū)域性海洋現(xiàn)象實現(xiàn)高時空分辨率的監(jiān)測， 在水色遙感探測中具有很大的發(fā)展?jié)摿?Heet al，2013)。GOCI 數(shù)據(jù)的空間分辨率為500m， 時間分辨率為 1h， 每天可以獲得八幅觀測圖像(Ryuet al， 2011)。He 等(2013)構(gòu)建了一個大氣校正算法和經(jīng)驗關(guān)系模型， 可以反演得到小潮時期和大潮時期長江口和杭州灣海域表層SSC， 該數(shù)據(jù)在Wang 等(2018b)中有具體描述。

表1 主要模型參數(shù) Tab.1 Main parameters of the model

數(shù)據(jù)同化的有效性需要用未被同化的獨立觀測進行檢驗(Elbernet al，2007)。因此， 在本文的孿生實驗和實際實驗中， 為了檢驗伴隨同化方法優(yōu)化反演模型初始場的有效性， 將隨機選取10%的觀測數(shù)據(jù)不參與同化而作為獨立觀測來檢驗數(shù)據(jù)同化的結(jié)果， 將這些觀測數(shù)據(jù)稱為“檢驗觀測”(check observations， COs)； 其余90%的觀測數(shù)據(jù)參與同化， 稱其為“被同化觀測”(assimilated observations， AOs)。孿生實驗和實際實驗中觀測數(shù)據(jù)的唯一區(qū)別是， 孿生實驗中的觀測數(shù)據(jù)是給定模型參數(shù)值運行正向模型生成的表層“觀測”， 實際實驗中的觀測數(shù)據(jù)則是實際GOCI 衛(wèi)星遙感數(shù)據(jù)反演所得SSC。

3 孿生實驗及結(jié)果分析

3.1 初始場的誤差對模擬結(jié)果的影響

初始場的不確定性是模型誤差的主要來源之一?；谌S懸沙輸運伴隨同化模型， 以小潮時期為例， 利用初始時刻實測的GOCI 衛(wèi)星遙感觀測數(shù)據(jù)經(jīng)過水平空 間插值獲得表層SSC 的初始場(圖2a)， 插值方法采用MATLAB 軟件自帶的自然鄰近插值法(natural)， 其余各層初始場由表層初始場線性插值得到， 插值公式與Wang 等(2018b)中相同， 該插值公式是由SSC 原位觀測數(shù)據(jù)( 旸楊 等， 2008)線性擬合得到， 較好地體現(xiàn)了杭州灣SSC 由表及底逐漸增大的特點。運行正向模型獲得整個研究區(qū)域的SSC“觀測值”， 每小時輸出一次結(jié)果。

為了評估初始場的誤差對模擬結(jié)果的影響， 對初始場分別添加最大不超過10%、20%、30%的隨機相對誤差， 計算公式如下：

其中，C0+%為添加隨機誤差后的初始場，C0為給定初始場，p為最大隨機相對誤差百分比，r為介于-1—1之間的隨機數(shù)。

圖2 初始時刻實測的GOCI 衛(wèi)星遙感觀測數(shù)據(jù)經(jīng)過水平空間插值后所得表層懸浮泥沙濃度的初始場 Fig.2 Initial field of the surface SSC obtained by horizontal spatial interpolation of GOCI satellite remote sensing observation data measured at initial time during neap tide

將垂向四層的SSC 模擬值與觀測值之間的相對誤差分別進行水平和垂向空間平均后得到整個研究區(qū)域的SSC 的平均相對誤差(mean relative error， MRE)隨模擬時間的變化情況(圖3)?？梢钥闯觯?SSC 模擬結(jié)果的MRE 均隨初始場誤差的增大而增大； 隨著模擬時間的增加， 模擬結(jié)果的MRE 沒有明顯下降。這說明初始場對于三維懸沙輸運模型是敏感參數(shù)， 初始場的給定越準確， 模擬結(jié)果的MRE 越小， 反之， 則越大。因此， 對于三維懸沙輸運模型而言， 尤其是短期的SSC 數(shù)值模擬， 初始場的給定至關(guān)重要， 將直接影響SSC 的模擬精度。

3.2 不同優(yōu)化算法對初始場優(yōu)化反演的影響

圖3 整個研究區(qū)域的SSC 模擬值與觀測值之間的平均相對誤差隨模擬時間的變化情況 Fig.3 The change of mean relative errors of SSC between the simulation results and observations for the total study area with the model time increasing

圖4 TE11—TE15 中反演后的表層初始場(a—e)、反演后的表層初始場與給定初始場的誤差(f—j)、反演后的表層初始場與給定初始場的相對誤差(k—o)的空間分布 Fig.4 The spatial distributions of the surface initial field after inversion in TE11—TE15 (a—e)； the spatial distributions of errors between the surface initial field after inversion and the given initial field in TE11—TE15 (f—j)； and the spatial distributions of relative errors between the surface initial field after inversion and the given initial field in TE11—TE15 (k—o)

圖5 實驗TE11—TE15 中約化代價函數(shù)(log10(J/J0)) (a)、代價函數(shù)關(guān)于初始場的梯度的L1 范數(shù)(Gradient of C0) (b)、SSC 模擬結(jié)果與“被同化觀測”(AOs)之間的平均絕對誤差(MAEs at AOs) (c)、SSC 模擬結(jié)果與“檢驗觀測”(COs)之間的平均絕對誤差(MAEs at COs) (d)、初始場的反演結(jié)果與給定分布之間的平均相對誤差(MREs of C0) (e)隨迭代步的變化 Fig.5 The values of log(J/J0) versus the iteration steps (a)； the variations of Gradient of C0 (b)； the variations of MAEs at AOs (c)； the variations of MAEs at COs (d)； and the variations of MREs of C0 (e) in TE11—TE15

表2 實驗TE11—TE15 中， 數(shù)據(jù)同化前后的約化代價函數(shù)(log10(J/J0))、代價函數(shù)關(guān)于初始場的梯度的L1 范數(shù)(Gradient of C0)、SSC 模擬結(jié)果與AOs 之間的平均絕對誤差(MAEs at AOs)、SSC 模擬結(jié)果與COs 之間的平均絕對誤差(MAEs at COs)和初始場的反演結(jié)果與給定分布之間的平均相對誤差(MREs of C0) Tab.2 log10(J/J0)， gradient of C0， MAEs at AOs， MAEs at Cos， and MREs of C0 before and after data assimilation in TE11—TE15

本節(jié)探討 GD、CG-PRP、CG-HS、CG-LS 和L-BFGS 五種算法在三維懸沙輸運伴隨同化模型中優(yōu)化反演初始場的表現(xiàn)。設(shè)計孿生實驗TE11—E15， 以圖2a 所示的初始場作為給定初始場， 初始場的初始猜測值取為給定初始場的水平空間平均值的一半， 即常數(shù)0.03325kg/m3， 其余模型參數(shù)均為默認值。從反演結(jié)果和誤差分布(圖4)來看， 五種算法的反演結(jié)果與給定分布較為吻合， 利用GD、CG-PRP 和CG-LS算法反演所得初始場的誤差主要集中分布在杭州灣東側(cè)灣口的海洋開邊界處， 而在其余的大部分區(qū)域誤差較小， 均取得了良好的反演效果； 然而利用CG-HS 和L-BFGS 算法反演所得初始場的誤差在灣頂、灣口和南部的大部分區(qū)域都很大， 反演效果較差， 且CG-HS 算法要比L-BFGS 算法更差。從統(tǒng)計結(jié)果(圖5、表2)來看， 五組實驗中約化代價函數(shù)、SSC 模擬結(jié)果與“觀測數(shù)據(jù)”之間的平均絕對誤差(mean absolute error， MAE)均有不同程度的下降， 可以定量說明數(shù)據(jù)同化的效果(Jinet al，2015； Liuet al， 2017； Zonget al， 2018)。以約化代價函數(shù)為例， 實驗TE15的約化代價函數(shù)下降速度最快， 但只下降了2 個數(shù)量級； 實驗TE11、TE12和TE14的約化代價函數(shù)下降速度相當， 均下降了3 個數(shù)量級； 實驗TE13的約化代價函數(shù)下降速度最慢， 在第3 步已達到收斂準則， 迭代停止。此外， 代價函數(shù)關(guān)于初始場的梯度的L1范數(shù)、反演結(jié)果與給定分布之間的MRE 兩個統(tǒng)計量可以直觀地說明初始場的反演效果。實驗TE11、TE12和TE14中兩個統(tǒng)計量分別經(jīng)過171、181、181 步均大幅度下降， 表明初始場的反演效果較好； 而實驗TE13和TE15中反演結(jié)果的誤差較大， 反演效果較差。綜合以上分析， GD、CG-PRP 和CG-LS 三種算法反演初始場的效果相當， 但GD 算法具有易于實現(xiàn)、迭代平穩(wěn)易于控制、良好的可操作性等優(yōu)點， 在本模型中具有很大的優(yōu)勢； L-BFGS 算法的反演速率最快， 但最終的反演結(jié)果較差； CG-HS 算法無論是反演速率還是反演結(jié)果都是最差的。以上數(shù)值實驗的結(jié)論與Chen 等(2013)和Zhang 等(2015)的研究結(jié)論相似。因此， 本節(jié)之后的孿生實驗和實際實驗選用GD 算法用于模型初始場的優(yōu)化反演。

3.3 初始猜測值對初始場優(yōu)化反演的影響

本節(jié)設(shè)計孿生實驗TE21—TE23研究初始場的初始猜測值對反演結(jié)果的影響。初始場的初始猜測值分別取為給定初始場的水平空間平均值的1/2、給定初始場的水平空間平均值和給定初始場的水平空間平均值的3/2， 即常數(shù)0.03325、0.0665 和0.09975kg/m3。從反演結(jié)果和誤差分布(圖6)來看， 三組實驗的反演結(jié)果相差不大， 反演結(jié)果與給定分布吻合良好， 相對誤差在大部分區(qū)域均小于20%。從統(tǒng)計結(jié)果(圖7、表3)來看， 三組實驗中相關(guān)統(tǒng)計量均大幅度下降。其中， 約化代價函數(shù)均下降了至少3 個數(shù)量級， 代價函數(shù)關(guān)于初始場的梯度的1L范數(shù)均至少下降了2 個數(shù)量級， SSC 模擬結(jié)果與觀測之間的MAE 均至少下降了1 個數(shù)量級， 初始場的反演結(jié)果與給定分布之間的MRE 均小于16%， 表明初始場的反演結(jié)果充分收斂到了“真實值”附近。本節(jié)數(shù)值實驗結(jié)果表明， 初始場的反演效果受 初始猜測值的影響較小， 初始場的“真實值”均能得到精確反演。因此， 本節(jié)之后的敏感性理想實驗中初始猜測值均取為0.03325kg/m3。

圖6 TE21—TE23 中反演后的表層初始場(a—c)、反演后的表層初始場與給定初始場的誤差(d—f)、反演后的表層初始場與給定初始場的相對誤差(g—i)的空間分布 Fig.6 The spatial distributions of the surface initial field after inversion in TE21—TE23 (a—c)； the spatial distributions of errors between the surface initial field after inversion and the given initial field in TE21—TE23 (d—f)； and the spatial distributions of relative errors between the surface initial field after inversion and the given initial field in TE21—TE23 (g—i)

圖7 實驗TE21—TE23 中l(wèi)og10(J/J0) (a)、Gradient of C0 (b)、MAEs at AOs(c)、MAEs at COs (d)、MREs of C0(e)隨迭代步的變化 Fig.7 The values of log(J/J0) versus the iteration steps (a)； the variations of Gradient of C0 (b)； the variations of MAEs at AOs (c)； the variations of MAEs at COs (d)； and the variations of MREs of C0 (e) in TE21—TE23

表3 實驗TE21—TE23 中， 數(shù)據(jù)同化前后的log10(J/J0)、Gradient of C0、MAEs at AOs、MAEs at COs 和MREs of C0 Tab.3 log10(J/J0)， gradient of C0， MAEs at AOs， MAEs at Cos， and MREs of C0 before and after data assimilation in TE21—TE23

3.4 衛(wèi)星遙感數(shù)據(jù)空間分布數(shù)量對初始場優(yōu)化反演的影響

在實際的衛(wèi)星遙感觀測中， 由于云覆蓋或惡劣天氣條件等因素的影響， 衛(wèi)星遙感數(shù)據(jù)往往存在一定程度的缺失， 因此有必要討論衛(wèi)星遙感數(shù)據(jù)數(shù)量對初始場優(yōu)化反演的影響。本節(jié)設(shè)計孿生實驗TE31—TE34， 四組實驗中“被同化觀測”的數(shù)量依次隨機減少1/10、1/4、1/2、3/4。從反演結(jié)果和誤差分布(圖8)可以看出， 隨著衛(wèi)星數(shù)據(jù)數(shù)量的減少， 反演精度相差不大， 反演結(jié)果在大部分區(qū)域與給定分布吻合良好， 誤差較大的區(qū)域主要集中在灣口的海洋開邊界處。從統(tǒng)計結(jié)果(圖9、表4)可以看出， 反演結(jié)果與給定分布之間的MRE 分別為10.94%、10.82%、12.71%和 13.62%， 依然取得了較高的反演精度， 實驗TE34相比于實驗TE31初始場的反演結(jié)果與給定分布之間的MRE 只增大了2.68%， 相差不大。據(jù)此我們認為， 在目前的實驗設(shè)置下， 衛(wèi)星遙感數(shù)據(jù)數(shù)量對初始場的優(yōu)化反演影響較小， 即便是分布稀疏或云覆蓋現(xiàn)象較為嚴重， 依然能夠取得良好的反演效果， 表明伴隨同化方法是還原SSC 空間分布的有效手段。

圖8 TE31—TE34 中反演后的表層初始場(a—d)、反演后的表層初始場與給定初始場的誤差(e—h)、反演后的表層初始場與給定初始場的相對誤差(i—l)的空間分布 Fig.8 The spatial distributions of the surface initial field after inversion in TE31—TE34 (a—d)； the spatial distributions of errors between the surface initial field after inversion and the given initial field in TE31—TE34 (e—h)； and the spatial distributions of relative errors between the surface initial field after inversion and the given initial field in TE31—TE34 (i—l)

圖9 實驗TE31—TE34 中l(wèi)og10(J/J0) (a)、Gradient of C0 (b)、MAEs at AOs (c)、MAEs at COs (d)、MREs of C0(e)隨迭代步的變化 Fig.9 The values of log(J/J0) versus the iteration steps (a)； the variations of Gradient of C0 (b)； the variations of MAEs at AOs (c)； the variations of MAEs at COs (d)； and the variations of MREs of C0 (e) in TE31—TE34

表4 實驗TE31—TE34 中， 數(shù)據(jù)同化前后的log10(J/J0)、Gradient of C0、MAEs at AOs、MAEs at COs 和MREs of C0 Tab.4 log10(J/J0)， gradient of C0， MAEs at AOs， MAEs at COs， and MREs of C0 before and after data assimilation in TE31—TE34

3.5 衛(wèi)星遙感數(shù)據(jù)同化時間窗口的寬度對初始場優(yōu)化反演的影響

本節(jié)設(shè)計孿生實驗TE41—TE44研究衛(wèi)星遙感數(shù)據(jù)同化時間窗口的寬度對初始場優(yōu)化反演的影響。SSC 數(shù)值模擬的起始時刻是GOCI 衛(wèi)星遙感觀測數(shù)據(jù)同化時間窗口的起始時刻， 同化窗口依次增大， 分別為1、2、3、4h， 即分別覆蓋2 次、3 次、4 次、5 次的GOCI 衛(wèi)星遙感觀測。從反演結(jié)果和誤差分布(圖10)可以看出， 隨著同化窗口的寬度的增大， 初始場的反演效果越來越好， 當同化窗口的寬度增加到3h， 在杭州灣中部和東南部兩個誤差相對較大的區(qū)域得到明顯改善。從統(tǒng)計結(jié)果(圖11、表5)可以看出， 四組實驗的約化代價函數(shù)、代價函數(shù)關(guān)于初始場的梯度的1L范數(shù)以及模擬結(jié)果和觀測之間的MAE 分別經(jīng)過200、200、194、191 步均顯著下降， 四組實驗初始場的反演結(jié)果與給定分布之間的MRE 分別為13.53%、10.03%、8.15%和8.73%， 即增大同化窗口的寬度可以有效改善反演結(jié)果。由此可見， 初始場的優(yōu)化反演效果對同化窗口的寬度較為敏感， 當SSC 數(shù)值模擬的起始時刻是GOCI 衛(wèi)星遙感觀測數(shù)據(jù)同化窗口的起始時刻， 只需要至少3h 的同化窗口即可得到改善效果明顯的反演結(jié)果。因此， 對于SSC 的數(shù)值模擬而言， 合理設(shè)置衛(wèi)星遙感數(shù)據(jù)同化時間窗口的寬度可以有效提高初始場的優(yōu)化反演效果。

3.6 背景水動力場的誤差對初始場優(yōu)化反演的的影響

圖10 TE41—TE44 中反演后的表層初始場(a—d)、反演后的表層初始場與給定初始場的誤差(e—h)、反演后的表層初始場與給定初始場的相對誤差(i—l)的空間分布 Fig.10 The spatial distributions of the surface initial field after inversion in TE41—TE44 (a—d)； the spatial distributions of errors between the surface initial field after inversion and the given initial field in TE41—TE44 (e—h)； and the spatial distributions of relative errors between the surface initial field after inversion and the given initial field in TE41—TE44 (i—l)

圖11 實驗TE41—TE44 中l(wèi)og10(J/J0) (a)、Gradient of C0 (b)、MAEs at AOs (c)、MAEs at COs (d)、MREs of C0(e)隨迭代步的變化 Fig.11 The values of log(J/J0) versus the iteration steps (a)； the variations of Gradient of C0 (b)； the variations of MAEs at AOs (c)； the variations of MAEs at COs (d)； and the variations of MREs of C0 (e) in TE41—TE44

表5 實驗TE41—TE44 中， 數(shù)據(jù)同化前后的log10(J/J0)、Gradient of C0、MAEs at AOs、MAEs at COs 和MREs of C0 Tab.5 log10(J/J0)， gradient of C0， MAEs at AOs， MAEs at COs， and MREs of C0 before and after data assimilation in TE41—TE44

本文中背景水動力場是由FVCOM 模擬生成， 用于驅(qū)動三維懸沙輸運模型。由于背景流場的模擬結(jié)果與實際流場不可避免地存在誤差， 因而三維懸沙輸運模型的模擬精度會很大程度上依賴于背景水動力場的模擬精度(Mehtaet al， 1989， Xieet al， 2009， Wanget al，2018b)， 因此有必要討論背景水動力場的誤差對初始場優(yōu)化反演的影響。李丹(2016)利用FVCOM 對長 江口杭州灣海域的水動力場進行模擬， 將模擬結(jié)果與實測資料進行對比， 潮流流速誤差的最大值在15%以內(nèi)。Wang 等(2018b)利用FVCOM 對杭州灣海域水動力場的模擬結(jié)果顯示， 模擬的潮流流速的相對誤差大多不超過20%。因此， 本節(jié)設(shè)計孿生實驗TE51—TE53， 三組實驗中分別對背景水動力場中流速的三個分量u、v、w均添加最大不超過10%、20%和30%的隨機誤差， 誤差添加方式與3.1 節(jié)中相同。從反演結(jié)果和誤差分布(圖12)來看， 三組實驗中初始場的反演結(jié)果相差不大， 反演結(jié)果均與給定分布吻合良好。從統(tǒng)計結(jié)果(圖13、表6)來看， 三組實驗相關(guān)統(tǒng)計量分別經(jīng)過179、179、180 步都下降到了相同的數(shù)量級， 三組實驗中反演結(jié)果與給定分布之間的MRE 均小于12%。因此，背景水動力場的誤差對初始場優(yōu)化反演的影響較小， 本文所采用的FVCOM 模擬結(jié)果足以滿足研究要求。

圖12 TE51—TE53 中反演后的表層初始場(a—c)、反演后的表層初始場與給定初始場的誤差(d—f)、反演后的表層初始場與給定初始場的相對誤差(g—i)的空間分布 Fig.12 The spatial distributions of the surface initial field after inversion in TE51—TE53 (a—c)； the spatial distributions of errors between the surface initial field after inversion and the given initial field in TE51—TE53 (d—f)； and the spatial distributions of relative errors between the surface initial field after inversion and the given initial field in TE51—TE53 (g—i)

圖13 實驗TE51—TE53 中l(wèi)og10(J/J0) (a)、Gradient of C0(b)、MAEs at AOs (c)、MAEs at COs (d)、MREs of C0(e)隨迭代步的變化 Fig.13 The values of log(J/J0) versus the iteration steps (a)； the variations of Gradient of C0 (b)； the variations of MAEs at AOs (c)； the variations of MAEs at COs (d)； and the variations of MREs of C0 (e) in TE51—TE53

表6 實驗TE51—TE53 中， 數(shù)據(jù)同化前后的log10(J/J0)、Gradient of C0、MAEs of AOs、MAEs of COs 和MREs of C0 Tab.6 log10(J/J0)， gradient of C0， MAEs at AOs， MAEs at COs， and MREs of C0 before and after data assimilation in TE51—TE53

3.7 伴隨同化方法和插值方法重構(gòu)初始場的能力比較

為了比較伴隨同化方法和插值方法重構(gòu)模型初始場的能力， 我們構(gòu)建一個新的初始場(圖14a)， 該初始場由圖 2a 所示初始場運行正向模型所得與實際GOCI 衛(wèi)星遙感觀測時刻對應(yīng)的“觀測數(shù)據(jù)”也即3.1節(jié)中的“觀測數(shù)據(jù)”經(jīng)過時間平均后獲得。利用新的初始場運行正向模型， 獲得與實際觀測數(shù)據(jù)一一對應(yīng)的受云覆蓋影響的“被同化觀測”和“檢驗觀測”用于本節(jié)的數(shù)據(jù)同化和初始場重構(gòu)。在伴隨同化方法中， 初始場的初始猜測值取為新的初始場的水平空間平均值的一半， 即常數(shù)0.0535kg/m3， 衛(wèi)星遙感數(shù)據(jù)同化時間窗口的寬度設(shè)置為3h， 其余模型設(shè)置不變。在插值方法中， 分別利用雙三次插值法、自然鄰近插值法和最近鄰點插值法插值數(shù)值模擬初始時刻的“觀測數(shù)據(jù)”， 獲得重構(gòu)的初始場。從初始場的重構(gòu)結(jié)果(圖14b—e)和誤差分布(圖14f—m)可以看出， 伴隨同化方法重構(gòu)的初始場與給定分布最接近， 在絕大部分區(qū)域二者吻合良好， 很好地重現(xiàn)了SSC 變化劇烈的區(qū)域， 且保持了很好的平滑性， 重構(gòu)結(jié)果與給定分布之間的MRE 只有11.53%； 雙三次插值法和自然鄰近插值法重構(gòu)初始場的效果較差， 雖然總體保持了與給定分布一致的空間變化規(guī)律， 但在杭州灣中西部、南部陸地邊界處和中東部存在三個相對誤差局部高值區(qū)， 插值結(jié)果對三處SSC 的劇烈變化刻畫不夠細致； 最近鄰點插值法重構(gòu)的初始場與給定分布偏差最大， 在上述三個局部區(qū)域SSC 均遠高于給定初始場， 重構(gòu)結(jié)果與給定分布之間的MRE 高達46.95%， 且插值結(jié)果的平滑性較差， 與實際偏差較大。因此， 與傳統(tǒng)的插值方法相比， 伴隨同化方法是重構(gòu)模型初始場更有效的手段。

4 實際實驗及結(jié)果分析

圖14 給定的(a)和重構(gòu)的(b—e)表層初始場的空間分布； 重構(gòu)的表層初始場與給定初始場的誤差(f—i)以及相對誤差(j—m)的空間分布 Fig.14 The spatial distribution of the given surface initial field (a)； the spatial distributions of initial field of the surface layer reconstructed by adjoint assimilation method， bicubic interpolation method， natural neighbor interpolation method， and nearest neighbor interpolation method (b—e)； the spatial distributions of errors between the surface initial field of the reconstructed surface initial field and the given initial field (f—i)； and the spatial distributions of relative errors of the reconstructed surface initial field and the given initial field (j—m)

實際實驗中， 分別同化典型的小潮時期和大潮時期GOCI 衛(wèi)星遙感資料所得杭州灣海域表層SSC 數(shù)據(jù)， 利用伴隨同化方法對模型初始場進行優(yōu)化反演。根據(jù)孿生實驗的結(jié)論， 優(yōu)化算法選用GD 算法， 初始場的初始猜測值取為初始時刻實測的GOCI 衛(wèi)星遙感觀測數(shù)據(jù)經(jīng)過空間插值獲得的SSC(圖2)， 衛(wèi)星遙感數(shù)據(jù)同化時間窗口的寬度設(shè)置為3h， 其余模型設(shè)置不變。從圖15 可以看出， 實際初始場缺測嚴重， 在缺測區(qū)域優(yōu)化反演后的初始場(圖16)與反演前的插值初始場(圖2)有很大不同。在小潮時期， 反演前的初始場的SSC 高值區(qū)出現(xiàn)在杭州灣的西部、中部和東南部， 而反演后的初始場的SSC 高值區(qū)出現(xiàn)在杭州灣的東南部及南部沿岸的大部分地區(qū)。在大潮時期， 反演前的初始場的SSC 高值區(qū)出現(xiàn)在杭州灣的中南部、東部和東南部的大部分區(qū)域， 而反演后的初始場的SSC高值區(qū)出現(xiàn)在杭州灣的東南部及南部沿岸的大部分地區(qū)。無論是小潮時期還是大潮時期， 反演后的初始場的SSC 均呈現(xiàn)出由南到北逐漸遞減的趨勢， 這與基于一個潮汐周期得到的杭州灣GOCI 衛(wèi)星遙感觀測中SSC“北低南高”的分布特征相似(Heet al， 2013； 江彬彬等， 2015； Liuet al， 2018)， 顯然優(yōu)化反演后的初始場更符合實際， 證明了伴隨同化方法優(yōu)化模型初始場的有效性。從圖17 可以看出， 優(yōu)化模型初始場可以很大程度上提高SSC 的模擬精度。小潮(大潮)時期， 經(jīng)過260(102)步迭代， 滿足收斂標準， 約化代價函數(shù)、代價函數(shù)關(guān)于初始場的梯度的L1范數(shù)均大幅下降并趨于平穩(wěn)， “被同化觀測”和模擬結(jié)果之間的MAE 從8.72×10-2kg/m3(3.66kg/m3)下降到了2.72×10-2kg/m3(7.25×10-1kg/m3)， 下降了68.81%(80.19%)； “檢驗觀測”和模擬結(jié)果之間的 MAE 從 9.02×10-2kg/m3(3.63kg/m3)下降到了2.88×10-2kg/m3(7.39×10-1kg/m3)， 下降了68.07%(79.64%)?？梢姡?反演后SSC的模擬結(jié)果和觀測之間的誤差大幅度降低， 模擬結(jié)果與觀測更接近。以上結(jié)果表明， 初始場對SSC 的模擬結(jié)果有重要的影響， 利用伴隨同化方法可以獲得模型的最優(yōu)初始場， 從而很大程度上改善SSC 的模擬結(jié)果。

圖15 小潮(a)和大潮(b)時期實際初始場的空間分布 Fig.15 The spatial distribution of actual initial field during neap tide (a) and spring tide(b)

圖16 實際實驗中， 小潮(a)和大潮(b)時期初始場的反演結(jié)果 Fig.16 The inversion result of initial field during neap tide (a) and spring tide (b) in practical experiments

圖17 實際實驗小潮時期的log10(J/J0)(a)、Gradient of C0(b)、MAEs at AOs(c)、MAEs at COs (d) 隨迭代步的變化； 實際實驗大潮時期的log10(J/J0) (e)、Gradient of C0 (f)、MAEs at AOs (g)、MAEs at COs (h) 隨迭代步的變化 Fig.17 The values of log(J/J0) (a， e)， the variations of Gradient of C0 (b， f)， the variations of MAEs at AOs (c， g)， and the variations of MAEs at COs (d， h) during neap tide and spring tide in practical experiments

根據(jù)以上分析， 實際初始場缺測嚴重時， 利用插值方法獲得的初始場(圖2)不能體現(xiàn)杭州灣南部沿岸高SSC 的特征， 且插值結(jié)果的平滑性較差。而利用伴隨同化方法反演所得初始場在有觀測區(qū)域與實際初始場較為一致， 且在初始場的觀測數(shù)據(jù)缺測區(qū)域， 伴隨同化方法的反演效果要更好， 很好地反演出了SSC“北低南高”的分布特征， 且總體上保持了很好的 平滑性， 更符合實際。這是因為插值方法僅考慮了SSC 的空間分布， 未考慮時間變化， 不能充分地利用觀測資料， 因而不能較好地重現(xiàn)SSC 的時空變化特征， 特別是在云覆蓋較為嚴重的區(qū)域， 傳統(tǒng)的插值方法往往是失效的， 模擬所得SSC 與實際偏差較大； 而伴隨同化方法既考慮SSC 的空間分布， 又考慮其時間變化， 可以獲得與實際觀測更加接近的初始場， 因而模擬結(jié)果更接近觀測值， 且更符合實際。

5 結(jié)論

本文基于一個三維懸沙輸運伴隨同化模型， 通過一系列孿生實驗和實際實驗， 對模型初始場進行優(yōu)化反演研究。

孿生實驗中， 首先驗證了三維懸沙輸運模型對于初始場的相對敏感性， 引出了利用伴隨同化方法優(yōu)化初始場的必要性。其次， 探討了初始場對于不同影響因素的相對敏感性， 包括優(yōu)化算法、初始猜測值、衛(wèi)星遙感數(shù)據(jù)數(shù)量、衛(wèi)星遙感觀測數(shù)據(jù)同化時間窗口的寬度和背景流場誤差。結(jié)果表明， 在五種優(yōu)化算法中， GD 算法是最優(yōu)的； 對本模型設(shè)置而言， 初始猜測值對反演效果影響較小， 初始場的“真實值”均能得到精確反演； 初始場對衛(wèi)星遙感數(shù)據(jù)數(shù)量較不敏感， 即便是分布稀疏或云覆蓋現(xiàn)象較為嚴重， 依然能夠取得良好的反演效果； 初始場對同化窗口的寬度較為敏感， 合理設(shè)置同化窗口的寬度可以有效改善初始場的優(yōu)化反演效果； 初始場對背景流場的誤差較不敏感， 即便是30%的隨機誤差， 反演結(jié)果與給定分布之間的MRE 依然小于12%。最后， 比較了伴隨同化方法和插值方法重構(gòu)模型初始場的能力。結(jié)果表明， 與傳統(tǒng)的插值方法相比， 伴隨同化方法是重構(gòu)模型初始場更有效的手段。

實際實驗中， 利用孿生實驗的相關(guān)結(jié)論， 在杭州灣海域同化典型的小潮時期和大潮時期的GOCI 衛(wèi)星遙感資料所得表層SSC 數(shù)據(jù)， 優(yōu)化反演了模型初始場。結(jié)果表明， 無論是小潮時期還是大潮時期， 利用伴隨同化方法均得到了更加符合實際的最優(yōu)初始場， 使得SSC 的數(shù)值模擬結(jié)果和實測數(shù)據(jù)更加接近。

該研究進一步證實了伴隨同化方法優(yōu)化模型初始場的有效性， 并為改進懸沙輸運模型及其他海洋數(shù)值模型的初始化方案提供了一種有效的方法。未來我們將嘗試同時同化更多、更長時間尺度的觀測資料， 為懸沙輸運的數(shù)值模擬提供更為精確的初始場。