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

杜云飛 張繼才 王道勝

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

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

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

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

1 模型介紹

1.1 正向模型

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

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

初始條件：

其中，C0是初始場。

邊界條件：

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

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

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

1.2 伴隨模型

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

1.3 優化算法

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

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

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

在GD 中， 搜索方向為代價函數關于模型初始場的負梯度方向：

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

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

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

2 研究區域、模型設置及數據來源

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

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

表1 主要模型參數 Tab.1 Main parameters of the model

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

3 孿生實驗及結果分析

3.1 初始場的誤差對模擬結果的影響

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

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

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

圖2 初始時刻實測的GOCI 衛星遙感觀測數據經過水平空間插值后所得表層懸浮泥沙濃度的初始場 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 模擬值與觀測值之間的相對誤差分別進行水平和垂向空間平均后得到整個研究區域的SSC 的平均相對誤差(mean relative error， MRE)隨模擬時間的變化情況(圖3)。可以看出， SSC 模擬結果的MRE 均隨初始場誤差的增大而增大； 隨著模擬時間的增加， 模擬結果的MRE 沒有明顯下降。這說明初始場對于三維懸沙輸運模型是敏感參數， 初始場的給定越準確， 模擬結果的MRE 越小， 反之， 則越大。因此， 對于三維懸沙輸運模型而言， 尤其是短期的SSC 數值模擬， 初始場的給定至關重要， 將直接影響SSC 的模擬精度。

3.2 不同優化算法對初始場優化反演的影響

圖3 整個研究區域的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 中約化代價函數(log10(J/J0)) (a)、代價函數關于初始場的梯度的L1 范數(Gradient of C0) (b)、SSC 模擬結果與“被同化觀測”(AOs)之間的平均絕對誤差(MAEs at AOs) (c)、SSC 模擬結果與“檢驗觀測”(COs)之間的平均絕對誤差(MAEs at COs) (d)、初始場的反演結果與給定分布之間的平均相對誤差(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 中， 數據同化前后的約化代價函數(log10(J/J0))、代價函數關于初始場的梯度的L1 范數(Gradient of C0)、SSC 模擬結果與AOs 之間的平均絕對誤差(MAEs at AOs)、SSC 模擬結果與COs 之間的平均絕對誤差(MAEs at COs)和初始場的反演結果與給定分布之間的平均相對誤差(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

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

3.3 初始猜測值對初始場優化反演的影響

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

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

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

3.6 背景水動力場的誤差對初始場優化反演的的影響

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

圖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 中log10(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 中， 數據同化前后的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 伴隨同化方法和插值方法重構初始場的能力比較

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

4 實際實驗及結果分析

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

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

圖16 實際實驗中， 小潮(a)和大潮(b)時期初始場的反演結果 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

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

5 結論

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

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

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

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