武浩文，趙艷玲，韓桂軍*，李威*，曹力戈，武曉博，李超亮，李云東，周功賦

(1.天津大學(xué) 海洋科學(xué)與技術(shù)學(xué)院，天津 300072；2.中國人民解放軍31010 部隊，北京 100081)

1 引言

隨著計算機(jī)技術(shù)的不斷發(fā)展，海洋數(shù)值模型越來越成為人們研究和預(yù)測海洋的重要工具。對于所有的海洋數(shù)值模式而言，參數(shù)化方案以及精確的參數(shù)值的給定對數(shù)值模擬結(jié)果有著重要的影響。在海洋數(shù)值模式調(diào)試的過程中，參數(shù)值的給定通常采用試錯的方法，以使模擬結(jié)果接近觀測。由于海洋數(shù)值模型的復(fù)雜性，這是一個極其耗費計算資源與人力資源的過程。因此，替代這種主觀調(diào)試方法的途徑之一是采用海洋數(shù)據(jù)同化技術(shù)，通過將觀測數(shù)據(jù)同化到海洋數(shù)值模型中，在模型狀態(tài)場調(diào)整的同時，對參數(shù)進(jìn)行估計，從而得到一個合理的參數(shù)值。這樣做的同時，也可以在一定程度上得以緩解試錯法在海洋數(shù)值模型調(diào)試中所帶來的繁重任務(wù)。利用基于伴隨模型的四維變分同化方法與集合卡爾曼濾波方法均可以開展海洋數(shù)值模型中參數(shù)的最優(yōu)估計。

以潮汐潮流的數(shù)值模擬為例，模型中不確定的參數(shù)包括水深、低摩擦系數(shù)和開邊界條件等。此前，研究者們大多采用伴隨方法來進(jìn)行潮汐潮流模擬中這些參數(shù)的優(yōu)化，從早期基于線性淺水方程模型的理想數(shù)據(jù)同化試驗[1-4]，到后來基于非線性二維和三維潮波數(shù)值模型的實際數(shù)據(jù)同化試驗[5-12]，均取得了很好的研究成果。相比較而言，基于集合卡爾曼濾波方法開展潮汐潮流模擬中的參數(shù)估計研究尚較少。Ngodock等[13]利用狀態(tài)增廣集合卡爾曼濾波（Augmented State Ensemble Kalman Filter，ASEnKF）對一個全球三維海洋環(huán)流模型的潮汐外強(qiáng)迫進(jìn)行修正。Wilson 和?zkan-Haller[14]、Wilson 等[15-16]、Landon 等[17]和Moghimi等[18]利用基于集合的數(shù)據(jù)同化方法對近岸、河流和潮汐汊道的水深參數(shù)進(jìn)行最優(yōu)估計。類似本文所開展的基于集合卡爾曼濾波方法對渤海和黃海這樣的陸架淺海M2分潮數(shù)值模擬中的水深參數(shù)進(jìn)行估計研究尚鮮有報道。

本文所關(guān)注的渤海和黃海，由于其獨特的地理位置和海底地形，使其潮波運動獨具特征。前人的研究表明，黃海的潮能消耗約占本海區(qū)總能量的80%以上[19-20]，潮汐、潮流成為本海區(qū)海水運動最重要的過程。因此，在上述海域建立高精度的潮汐潮流模擬數(shù)值模型有著重要意義。本研究擬采用集合卡爾曼濾波方法，以集中反映該海域潮波運動特征的優(yōu)勢分潮—M2分潮為例，對數(shù)值模型構(gòu)建過程中所涉及的不確定參數(shù)—水深進(jìn)行最優(yōu)估計，以期獲得更好的M2分潮模擬結(jié)果。

2 數(shù)據(jù)和方法

2.1 數(shù)據(jù)

本研究的海區(qū)范圍為35°～41°N，117°～127°E，包括渤海和黃海的大部分海域（圖1）。數(shù)值模型所使用的水深數(shù)據(jù)來自美國國家地球物理數(shù)據(jù)中心（U.S.National Geophysical Data Center，NGDC）的ETOPO2（https://www.ngdc.noaa.gov/mgg/global/etopo2.html），其分辨率為（1/30）°×（1/30）°。研究海區(qū)內(nèi)M2分潮網(wǎng)格化的潮汐調(diào)和常數(shù)來自日本國家天文臺的NAO.99Jb數(shù)據(jù)[21]（https://www.miz.nao.ac.jp/staffs/nao99/index_En.html），其分辨率為（1/12）°×（1/12）°；南部開邊界處M2分潮的潮流調(diào)和常數(shù)來自美國俄勒岡州立大學(xué)的TPXO9 數(shù) 據(jù)[22]（https://tpxows.azurewebsites.net/），其 分辨率為（1/12）°×（1/12）°。由于NAO.99Jb 數(shù)據(jù)在中國沿岸海域有較高的精度[23]，因此本研究將其作為數(shù)據(jù)同化試驗中的觀測數(shù)據(jù)來使用。此外，本研究還使用了研究海域內(nèi)34 個驗潮站的M2分潮調(diào)和常數(shù)[6]，其位置和編號如圖1 所示。

圖1 研究海域范圍及各驗潮站位置和編號Fig.1 Study area and locations of tide gauges with numbers assigned

2.2 海洋數(shù)值模式

本研究用于關(guān)注海域潮汐潮流模擬的數(shù)值模型為廣義坐標(biāo)系統(tǒng)的美國普林斯頓大學(xué)海洋模式（Princeton Ocean Model with Generalized Coordinate System，POMgcs）[24]的外模式，其動力學(xué)基本方程組為

式中，t為時間；x和y為水平方向的笛卡兒坐標(biāo)；U和V為水平方向的流速；D=H+ζ，其中ζ為水位，H為水深；f為科里奧利參數(shù)；g為重力加速度；ρa和 ρw為空氣和水的密度；Cd為風(fēng)應(yīng)力拖曳系數(shù)；Wx和Wy為水平方向的風(fēng)速；Cb為 底摩擦系數(shù)；AM為水平渦動系數(shù)。

在本研究中，數(shù)值模型的空間網(wǎng)格分辨率為（1/30）°，時間步長為12 min，初始條件取為0 值。底摩擦系數(shù)在整個海區(qū)取為常數(shù)1.0×10-3。研究海域南部開邊界處的水位和流速按下式給定：

式中，Y為水位或流速；A為振幅；ω為頻率；φ為遲角；V和u為天文初相角。

2.3 集合卡爾曼濾波方法

本研究所采用的集合卡爾曼濾波方法是由Anderson[25]提出的確定性方法，即基于最小二乘框架的集合調(diào)整卡爾曼濾波（Ensemble Adjustment Kalman Filter，EAKF）。對某一觀測yo而言，該方法的實施包括如下兩個步驟：

第1 步，計算觀測增量：

式中，a為經(jīng)驗給定的局地化半徑，b為觀測點與狀態(tài)變量格點之間的距離。

在數(shù)據(jù)同化過程中，由于分析集合的不確定性（標(biāo)準(zhǔn)差或集合離散度）總是小于先驗集合的不確定性，隨著同化的持續(xù)進(jìn)行，分析集合的離散度逐漸減小。但由于不可避免的模式偏差的存在，使得所刻畫的先驗場的不確定性小于實際，高估了模式的準(zhǔn)確性，從而導(dǎo)致同化作用偏弱，甚至引起濾波發(fā)散。因此，對于狀態(tài)估計，本研究引入靜態(tài)乘法膨脹方案，即確定一個常數(shù)膨脹因子，對各個集合成員相對于集合平均的擾動進(jìn)行膨脹，用以調(diào)整集合的離散度，從而避免集合離散度降低，導(dǎo)致濾波發(fā)散。而對于參數(shù)估計，則引入條件靜態(tài)膨脹方案[27-28]，即通過判斷當(dāng)前時刻集合方差的大小，決定是否對其進(jìn)行參數(shù)膨脹。參數(shù)的條件靜態(tài)乘法膨脹公式如下：

此外，為了進(jìn)一步提高EAKF 參數(shù)估計的效果，在后續(xù)實際數(shù)據(jù)試驗中采用了強(qiáng)化的參數(shù)校正數(shù)據(jù)同化方法（Data Assimilation Scheme for Enhancive Parameter Correction，DAEPC）[29]。該方法的實現(xiàn)方式是在模型狀態(tài)場的估計達(dá)到準(zhǔn)平衡態(tài)后，再啟動對模型參數(shù)的估計。這是因為在數(shù)據(jù)同化的初期，觀測與狀態(tài)變量之間的協(xié)方差由狀態(tài)誤差主導(dǎo)；在進(jìn)行一段時間的模型狀態(tài)估計后，狀態(tài)誤差降低，參數(shù)誤差則成為模式誤差的主要成因，即觀測與狀態(tài)變量之間的協(xié)方差反映的主要是由參數(shù)主導(dǎo)的信號，因此，采用DAEPC 方法有助于強(qiáng)化參數(shù)調(diào)整過程中的信噪比，從而提高參數(shù)估計的效果。

3 數(shù)據(jù)同化方案設(shè)計

3.1 敏感性分析

Han 等[8]研究表明，相較于底摩擦系數(shù)和開邊界條件，水深對潮汐潮流模擬結(jié)果有著較大的影響。在本研究進(jìn)行水深參數(shù)估計前，首先開展水深參數(shù)的敏感性分析試驗。按照Han 等[8]的方法，將本文研究海域按照水深進(jìn)行劃分。經(jīng)大量測試，將區(qū)域按照水深劃分為0～20 m、20～40 m、40～60 m 和60 m 以上4 個區(qū)域（如圖1 中的黑色等深線所示）。根據(jù)水深數(shù)據(jù)的精度，本研究在上述4 個區(qū)域內(nèi)原水深數(shù)據(jù)的基礎(chǔ)上，分別疊加0.5 m、2 m、4 m 和6 m，生成相應(yīng)的有偏水深數(shù)據(jù)，并保持模式其他參數(shù)不變，進(jìn)行約10 d（以M2分潮周期12.42 h 計算，共計20 個周期）自由積分。在模式積分的最后一個M2分潮周期（后文簡稱周期）內(nèi)，將水位計算結(jié)果與原模式結(jié)果進(jìn)行對比，計算各網(wǎng)格點水位的時間平均均方根誤差（Root Mean Square Error，RMSE），結(jié)果如圖2 所示。由圖可以看出，大部分區(qū)域的水位誤差在0.1 m 以上，且靠近陸地的海域，水位誤差更大，尤其在黃海的北部和東部，最大水位誤差達(dá)到了0.6 m。這一水深參數(shù)的敏感性分析試驗結(jié)果表明，在本研究所采用的水深參數(shù)區(qū)域劃分方案下，M2分潮的模擬結(jié)果能夠產(chǎn)生明顯的差異，為開展后續(xù)的數(shù)據(jù)同化試驗奠定了基礎(chǔ)。

圖2 基于有偏水深數(shù)據(jù)獲得的水位數(shù)值模擬結(jié)果的時間平均均方根誤差Fig.2 Time mean root mean square error of water level simulated with biased water depths

3.2 數(shù)據(jù)同化試驗

本研究擬同時開展理想和實際水深參數(shù)估計試驗。其中理想數(shù)據(jù)同化試驗又稱為孿生試驗。在這兩種數(shù)據(jù)同化試驗中，均將水深增量作為待估計的參數(shù)，后文統(tǒng)一簡稱為水深參數(shù)。待估計的水深參數(shù)的設(shè)置采用與敏感性分析試驗相同的分區(qū)方案，即劃分為0～20 m、20～40 m、40～60 m和60 m 這4 個區(qū)域，從而有4 個與其相對應(yīng)的水深參數(shù)待估計。以某一個模式網(wǎng)格點 (i,j)為例，估計后的水深按下式計算：

在后續(xù)的理想和實際數(shù)據(jù)同化試驗中，集合數(shù)均設(shè)定為30，每個時間步均進(jìn)行同化，同化時長設(shè)定為3 個周期。對于狀態(tài)估計，由于POMgcs 模式采用蛙跳格式的時間差分方案，因此在同化試驗中，當(dāng)前的觀測同時用于調(diào)整當(dāng)前和上一個時間步的水位和潮流狀態(tài)變量[30]。同化后，利用更新后的水深繼續(xù)積分約3 d 以達(dá)到穩(wěn)定。利用最后一個周期的水位數(shù)據(jù)進(jìn)行調(diào)和分析。需要注意的是，同化過程中，調(diào)整后的水深有可能會使得D＜0。為此，在同化過程中為水深的調(diào)整設(shè)定一個閾值，以保證這種情況不會發(fā)生。

3.2.1 理想數(shù)據(jù)同化試驗

在理想數(shù)據(jù)同化試驗中，將ETOPO2 的水深數(shù)據(jù)視為真實水深，利用該水深自由積分得到的水位數(shù)據(jù)作為真實場。在此基礎(chǔ)上，將真實場疊加上標(biāo)準(zhǔn)差為0.1 m 的高斯白噪聲生成與模式網(wǎng)格點和時間積分步一致的水位觀測場。

水深參數(shù)集合的形成與敏感性分析試驗相同，即在0～20 m、20～40 m、40～60 m 和60 m 這4 個區(qū)域內(nèi)，分別以0.5 m、2 m、4 m 和6 m 為先驗水深參數(shù)(m=1,2,3,4)，將其疊加高斯白噪聲形成相應(yīng)的先驗水深參數(shù)集合，標(biāo)準(zhǔn)差分別取為各先驗水深參數(shù)值的5%。

數(shù)值模型基于上述有偏的水深參數(shù)集合積分約7 d以達(dá)到穩(wěn)定。需要說明的是，在測試中發(fā)現(xiàn)：同時進(jìn)行狀態(tài)和水深參數(shù)估計已經(jīng)能夠獲得很好的試驗結(jié)果，因此在理想數(shù)據(jù)同化試驗中并未采用DAEPC方法。

3.2.2 實際數(shù)據(jù)同化試驗

在實際數(shù)據(jù)同化試驗中，利用NAO.99Jb 的調(diào)和常數(shù)計算出模式積分的所有時間步、每隔3 個模式水平網(wǎng)格點的水位值，連同驗潮站處所計算得到的水位值一并作為本研究的觀測數(shù)據(jù)。令4 個先驗水深參=0(m=1,2,3,4)，然后采用疊加高數(shù)均為0，即斯白噪聲的方式各自形成先驗水深參數(shù)集合，其標(biāo)準(zhǔn)差與理想試驗相同。模型積分約7 d 達(dá)到穩(wěn)定。之后引入DAEPC 方法進(jìn)行水深參數(shù)估計。在測試中發(fā)現(xiàn)，進(jìn)行模型狀態(tài)場估計的1 個周期后，狀態(tài)場估計的誤差已達(dá)到準(zhǔn)平衡態(tài)。因此，進(jìn)行實際數(shù)據(jù)同化試驗時，在第1 個周期僅進(jìn)行狀態(tài)場估計，之后同時進(jìn)行狀態(tài)場和水深參數(shù)估計。在實際數(shù)據(jù)同化試驗中，進(jìn)行如下兩組試驗：（1）僅利用來自NAO.99Jb 的觀測數(shù)據(jù)，設(shè)置觀測誤差為0.2 m。（2）同時利用來自NAO.99Jb和驗潮站的觀測數(shù)據(jù)，設(shè)置NAO.99Jb 數(shù)據(jù)的觀測誤差為0.2 m，驗潮站數(shù)據(jù)的觀測誤差為0.1 m。這兩組試驗中所給定的觀測誤差是通過多次數(shù)據(jù)同化試驗后獲得的最優(yōu)結(jié)果。

4 結(jié)果

4.1 理想數(shù)據(jù)同化試驗

本節(jié)主要討論當(dāng)數(shù)值模式的誤差僅來自水深參數(shù)時，EAKF 方法的參數(shù)估計效果。通過多組試驗結(jié)果對比，發(fā)現(xiàn)局地化半徑為40 個網(wǎng)格點，不進(jìn)行狀態(tài)場膨脹，而參數(shù)膨脹系數(shù)為1.3 時，水深參數(shù)估計的效果最好。因此，這里基于該組試驗結(jié)果，分析EAKF方法對前述4 個參數(shù)的估計效果。

圖3a 和圖3b 分別為同化后空間平均的水位RMSE 時間序列和4 個水深參數(shù)估計值與真值偏差的時間序列。為方便進(jìn)行比較，將圖3b 中水深參數(shù)估計偏差進(jìn)行標(biāo)準(zhǔn)化處理后形成圖3c，即將4 個水深參數(shù)估計偏差的時間序列分別除以初始時刻設(shè)置的水深參數(shù)值（0.5 m、2 m、4 m 和6 m）。

由圖3a 和圖3b 可見，在數(shù)據(jù)同化試驗的初期，模式水位誤差由淺水區(qū)域主導(dǎo)，因此很快達(dá)到穩(wěn)定。與此同時，EAKF 的狀態(tài)估計很快將水位誤差調(diào)整到觀測誤差（0.1 m）以下。一段時間后，狀態(tài)誤差已經(jīng)很低，此時模式誤差主要來自。在同化3 個周期后，兩個參數(shù)在真實值附近穩(wěn)定，水位誤差也進(jìn)一步降低。而雖然很快達(dá)到穩(wěn)定，但穩(wěn)定后呈周期性震蕩，與真實值始終維持一定的距離。同化結(jié)束后，水位誤差保持在0.01 m 以下，說明模式狀態(tài)場誤差已經(jīng)大幅降低，且同化后模式的自由積分（3個周期后）結(jié)果與真實場也十分接近（圖形未在此給出）。

圖3 后驗水位均方根誤差（a）、后驗水深參數(shù)估計偏差（b）時間序列和標(biāo)準(zhǔn)化的后驗水深參數(shù)估計偏差時間序列（c）Fig.3 Time series of the root mean square error for the posterior water level (a),bias (b) and standardized bias (c) for the posterior water depth parameters

由標(biāo)準(zhǔn)化后的水深參數(shù)偏差（圖3c）可見，在參數(shù)估計的起始階段，很快回到了真實值附近，說明有很高的敏感性。在參數(shù)估計的第1 和第2 個周期，產(chǎn)生了較明顯的震蕩，這是由于和仍未穩(wěn)定，在參數(shù)估計的過程中不斷被調(diào)整。而到了第3 個周期，基本穩(wěn)定后，也保持在真值附近。在整個參數(shù)估計過程中趨勢比較平穩(wěn)，偏差減小的速度較慢，說明這兩處位置的敏感性相對較低，需要較長時間的參數(shù)估計才能達(dá)到穩(wěn)定。而在穩(wěn)定后一直與真值保持著一定的偏差，說明并不敏感，即使沒有回到真值，也不能產(chǎn)生很大的狀態(tài)誤差。

圖4 為基于真實場和理想數(shù)據(jù)同化試驗前后數(shù)據(jù)繪制的M2分潮同潮圖。其中，遲角以格林威治時間（對應(yīng)遲角為0°）為基準(zhǔn)繪制。由圖4b 可見，同化前，海區(qū)邊緣處的振幅和無潮點附近的遲角均出現(xiàn)一定誤差，且渤海北部的無潮點位置相較于真實場（圖4a）出現(xiàn)偏移。同化后，振幅和遲角（圖4c）相比于真實場均已無明顯差異，渤海北部無潮點的位置也得到改善。

圖4 利用真實場數(shù)據(jù)（a）、理想試驗中的先驗結(jié)果（b）和后驗結(jié)果（c）繪制的M2 分潮同潮圖Fig.4 Cotidal chart of M2 constituent from the truth (a),prior (b) and posterior (c)

表1 給出了先驗和后驗振幅與遲角的空間平均誤差。同化前，振幅的誤差約為7.6 cm，遲角誤差約為12°54′；同化后，振幅誤差僅為0.2 cm，遲角誤差僅為18′。表明同化后，模式誤差整體上得以大幅降低。

表1 理想數(shù)據(jù)同化試驗中M2 分潮的先驗和后驗振幅與遲角空間平均誤差Table 1 Spatial averaged errors of amplitude and phase lag of M2 constituent from the prior and posterior in twin experiment

綜上所述，在理想數(shù)據(jù)同化試驗中，EAKF 方法通過對水深參數(shù)進(jìn)行最優(yōu)估計，有效地降低了模式的偏差。這一結(jié)果也表明EAKF 方法可以很好地應(yīng)用于潮汐潮流數(shù)值模擬的水深參數(shù)估計中。

4.2 實際數(shù)據(jù)同化試驗

4.2.1 基于NAO.99Jb 數(shù)據(jù)的水深估計

在本試驗中，將NAO.99Jb 數(shù)據(jù)作為觀測進(jìn)行水深估計。經(jīng)過多組試驗，選取局地化半徑為15 個網(wǎng)格點、狀態(tài)膨脹系數(shù)為1.03、參數(shù)膨脹系數(shù)為1.10 的數(shù)據(jù)同化試驗結(jié)果在此討論。

圖5 為后驗水位均方根誤差（將數(shù)據(jù)同化試驗得到的水位與NAO.99Jb 數(shù)據(jù)對比）的時間序列分布圖。由圖可見，在同化開始后，后驗水位均方根誤差在周期性震蕩的過程中下降。這是由于模式計算結(jié)果與NAO.99Jb 數(shù)據(jù)之間存在較大的相位誤差，導(dǎo)致兩者水位變化方向有時不同，從而使得在這段時間內(nèi)，模式模擬的水位均方根誤差不斷增大。但總體上來看，水位均方根誤差保持在觀測誤差（0.2 m）附近波動。而模式狀態(tài)的調(diào)整無法改正這一相位誤差。第1 個周期后，開始同時進(jìn)行狀態(tài)和參數(shù)估計，這時由于水深參數(shù)的調(diào)整，使得水位模擬均方根誤差不斷減小，最后維持在0.1 m 附近，且小于所給定的NAO.99Jb 數(shù)據(jù)的觀測誤差（0.2 m），說明數(shù)據(jù)同化試驗是成功的。

圖5 后驗水位均方根誤差的時間序列Fig.5 Time series of the root mean square error for the posterior water level

圖6 為基于NAO.99Jb 和實際數(shù)據(jù)同化試驗前后數(shù)據(jù)繪制的M2分潮同潮圖。由圖6c 可見，整體上，同化后模式的調(diào)和分析結(jié)果更接近NAO.99Jb 數(shù)據(jù)，尤其是在黃海的東北和東南部，同化后振幅明顯減小，遲角的分布也更趨近NAO.99Jb。同時，各無潮點的位置也有一定改善：渤海北部的無潮點向左偏移，渤海南部、黃海中部的無潮點則向下移動。對比圖6a、圖6b 和圖6c 可見，江華灣部分誤差在同化后誤差較大。導(dǎo)致這一現(xiàn)象的原因可能是由于本研究中估計的水深參數(shù)個數(shù)較少，使得在當(dāng)前這樣的參數(shù)估計方案下，估計后的水深參數(shù)雖然使得總體的誤差在減小，但這種減小的方式有可能是以某些區(qū)域的模擬誤差減小、某些區(qū)域的模擬誤差增大的方式而獲得，而不是以整個研究區(qū)域的模擬誤差同時減小的方式獲得。因此，設(shè)計更好的水深參數(shù)估計方案是后續(xù)需要進(jìn)一步開展的研究工作。

圖6 利用NAO.99Jb 數(shù)據(jù)（a）和NAO.99Jb 實際數(shù)據(jù)同化試驗中的先驗（b）和后驗（c）結(jié)果繪制的M2 分潮同潮圖Fig.6 Cotidal chart of M2 constituent from the NAO.99Jb (a),prior (b) and posterior (c)

表2 給出了實際數(shù)據(jù)同化試驗中M2分潮的調(diào)和常數(shù)在同化前后與NAO.99Jb 相比的空間平均誤差。從整體上來看，無論是振幅還是遲角，后驗結(jié)果都優(yōu)于先驗結(jié)果。

表2 NAO.99Jb 實際數(shù)據(jù)同化試驗中M2 分潮的先驗和后驗振幅與遲角空間平均誤差Table 2 Spatial averaged errors of amplitude and phase lag of M2 constituent for the prior and posterior in the NAO.99Jb data assimilation experiment

4.2.2 基于NAO.99Jb 和驗潮站數(shù)據(jù)的水深估計

在4.2.1 節(jié)的實際數(shù)據(jù)同化試驗基礎(chǔ)上，加入驗潮站數(shù)據(jù)進(jìn)行水深估計。在大量同化試驗方案測試的基礎(chǔ)上得到如下兩組較優(yōu)的數(shù)據(jù)同化試驗：試驗1 為使用全部NAO.99Jb 和驗潮站數(shù)據(jù)，開展水深估計試驗；試驗2 為模式網(wǎng)格點先驗水深小于20 m 且非渤海區(qū)域內(nèi)僅使用驗潮站數(shù)據(jù)，其余部分則使用全部數(shù)據(jù)，開展水深參數(shù)估計試驗。

在上述同化試驗中，將驗潮站數(shù)據(jù)的觀測誤差設(shè)置為0.1 m，NAO.99Jb 數(shù)據(jù)的觀測誤差設(shè)置為0.2 m。其余試驗設(shè)置與4.2.1 節(jié)相同。基于試驗1 和試驗2 的結(jié)果，將4.2.1 節(jié)的NAO.99Jb 數(shù)據(jù)同化試驗作為對比試驗，進(jìn)行結(jié)果討論分析。

表3 給出了同化前、試驗1、試驗2 和對比試驗與驗潮站振幅，遲角的空間平均誤差。由表3 可以看出，各試驗的整體誤差相較于同化前均有所降低。其中，試驗2 的誤差最小，相比于同化前，振幅誤差減小了40.27%，遲角誤差降低了49.19%；而試驗1 的振幅和遲角誤差只分別降低了28.52%和29.39%。

圖7 給出了試驗1、試驗2 和對比試驗所得到的M2分潮同潮圖以及在驗潮站位置處同化前后振幅和遲角誤差的變化情況。圖8 為先驗、試驗1、試驗2 和對比試驗在各驗潮站位置處的調(diào)和常數(shù)誤差統(tǒng)計直方圖。由圖7 和圖8 可見，相較于同化前，試驗1中驗潮站位置處的振幅和遲角誤差減小的站點數(shù)分別為21 個和26 個，試驗2 中分別為25 個和30 個。這表明雖然試驗2 相較于試驗1 所使用的觀測信息更少，但驗潮站位置處振幅和遲角的整體誤差以及誤差分布都優(yōu)于試驗1。這可能是由于在同化試驗中的水深參數(shù)估計中出現(xiàn)了“少數(shù)服從多數(shù)”的現(xiàn)象，即EAKF 方法根據(jù)大部分觀測信息來調(diào)整水深，使觀測信息較多的地方水位模擬結(jié)果得到改進(jìn)，而觀測信息較少的地方水位模擬結(jié)果反而變差。具體來說，在試驗1 中，黃海北部的模擬結(jié)果更好，渤海的模擬誤差則較大，其余驗潮站位置處接近對比試驗。而試驗2 中江華灣區(qū)域改進(jìn)明顯，渤海和黃海北部的誤差則接近對比試驗。這表明在先驗水深小于20 m 且非渤海區(qū)域的NAO.99Jb 數(shù)據(jù)使EAKF 更傾向于優(yōu)化黃海北部地區(qū)。這也說明去除這些觀測信息后，其余觀測信息的比重被調(diào)整，反而使更多區(qū)域的水深參數(shù)被優(yōu)化。

在圖8a 和圖8b 中，試驗2 相較試驗1 得到優(yōu)化的站點共20 個。其中，渤海和江華灣部分的振幅得到改進(jìn)，尤其是江華灣區(qū)域的改進(jìn)尤為明顯；而黃海北部和南部反而變差。在圖8c 和圖8d 中，試驗2 的遲角優(yōu)于試驗1的驗潮站同樣是20 個。其中，黃海北部、渤海和黃海東南部的提升較為明顯，江華灣、黃海西部和南部則有所下降。因此，試驗2 雖然在整體上振幅和遲角的誤差更小，但誤差的空間分布沒有明顯改進(jìn)，說明在后續(xù)的研究中仍需要進(jìn)一步探索如何改進(jìn)參數(shù)估計的同化試驗方案。

此外，由圖8c 和圖8d 可見，試驗1 中葫蘆島站（圖1）中4 號的遲角存在較大誤差，達(dá)到了175°10′，幾乎是反位相。對照其地理位置可以看出，該處距離渤海北部的無潮點較近（圖7b）。由于無潮點模擬位置的偏移，在該點處產(chǎn)生了較大的誤差，從而也導(dǎo)致了試驗1 模擬的遲角整體誤差大幅增大（表3）。如果誤差統(tǒng)計中不考慮該驗潮站的遲角誤差，則試驗1 模擬的遲角誤差為16°36′。相比較而言，對比試驗的遲角誤差則為16°56′。因此，從這個角度來說，試驗1 相較于對比試驗，遲角的模擬仍然有所改進(jìn)。這從對比圖7b和圖7f 中驗潮站誤差變小的個數(shù)上也可以得到佐證。

圖7 試驗1（a，b）、試驗2（c，d）和對比試驗（e，f）的振幅（a，c，e；單位：cm）和遲角（b，d，f；單位：（°））及其在各個驗潮站位置的數(shù)值相較于同化前的變化情況Fig.7 Amplitude (a,c,e;unit:cm) and phase lag (b,d,f;unit:(°)) from experiment 1 (a,b),experiment 2 (c,d) and NAO.99Jb data assimilation experiment (e,f),and change of errors at each tide gauge

圖8 模式先驗、對比試驗、試驗1 和試驗2 的振幅（a,b）和遲角（c,d）在驗潮站位置處的均方根誤差Fig.8 Amplitude (a,b) and phase lag (c,d) root mean square errors of the prior,NAO.99Jb data assimilation experiment,experiment 1,and experiment 2 with respect to the tide gauges

表3 同化前、試驗1、試驗2 和對比試驗與驗潮站振幅、遲角的平均空間誤差Table 3 Spatial averaged errors of amplitude and phase lag from the model,experiment 1,experiment 2 and NAO.99Jb data assimilation experiment with respect to those from tide gauges

在上述實際數(shù)據(jù)同化試驗中，4 個水深參數(shù)的調(diào)整量分別約為6 m、-4 m、15 m 和10 m。需要說明的是，不同于理想數(shù)據(jù)同化試驗有“真實”的水深參數(shù)可供比較，對于實際數(shù)據(jù)同化試驗來說，并沒有這樣的“真實”水深存在。而對于數(shù)值模式來說，其模擬結(jié)果的誤差主要來自模式的動力框架本身以及不確定的開邊界條件、底摩擦系數(shù)和水深。因此，在本研究中，由于只進(jìn)行了水深的參數(shù)估計，上述模式的動力框架本身以及不確定的開邊界條件和底摩擦系數(shù)所導(dǎo)致的模擬誤差都涵蓋在水深參數(shù)的調(diào)整中。因此，調(diào)整后的水深其實并不代表實際海域所觀測的水深。

5 結(jié)論

本研究基于集合調(diào)整卡爾曼濾波（EAKF）方法，采用廣義坐標(biāo)系統(tǒng)的美國普林斯頓大學(xué)海洋模式（POMgcs）的外模式開展了渤海和部分黃海海域M2分潮模擬中的水深估計研究。通過理想數(shù)據(jù)同化試驗（又稱孿生試驗）和基于NAO.99Jb 與驗潮站數(shù)據(jù)，并采用強(qiáng)化的參數(shù)校正數(shù)據(jù)同化方法（DAEPC）的實際數(shù)據(jù)同化試驗結(jié)果的討論分析，得出主要結(jié)論如下。

（1）EAKF 方法在理想數(shù)據(jù)同化試驗中可以實現(xiàn)對潮汐潮流模擬中的水深參數(shù)的最優(yōu)估計，還原水深參數(shù)的“真實場”。這表明，EAKF 方法可以很好地應(yīng)用于潮汐潮流模擬中的狀態(tài)場和參數(shù)估計。

（2）在實際數(shù)據(jù)同化試驗中，利用DAEPC 方法開展了狀態(tài)和水深參數(shù)估計。當(dāng)僅使用NAO.99Jb 數(shù)據(jù)作為觀測時，EAKF 方法可以通過實際的水深估計，提高數(shù)值模型對M2分潮振幅和遲角的模擬精度。相較于NAO.99Jb 數(shù)據(jù)，水深參數(shù)估計后的振幅誤差減小了30.65%，遲角誤差減小了33.86%。

（3）在NAO.99Jb 數(shù)據(jù)基礎(chǔ)上加入驗潮站數(shù)據(jù)的水深估計試驗中，通過研究改進(jìn)觀測系統(tǒng)，進(jìn)一步提高了水深估計的效果。對比驗潮站數(shù)據(jù)，模擬的M2分潮振幅與遲角誤差相較于同化前分別減小了40.27%和49.19%。在先驗水深小于20 m 且非渤海的區(qū)域僅使用驗潮站的觀測信息后，水深參數(shù)估計可以在整體上進(jìn)一步提高M(jìn)2分潮的模擬效果，但誤差的空間分布沒有明顯改進(jìn)。

值得注意的是，在實際數(shù)據(jù)同化試驗中，江華灣區(qū)域的誤差最大。這是由于該區(qū)域海底地形復(fù)雜，且水深較淺，數(shù)值模式本身模擬結(jié)果的準(zhǔn)確度不高。即使進(jìn)行了數(shù)據(jù)同化，對狀態(tài)場和水深參數(shù)估計后，該區(qū)域的誤差仍然較大。因此，需要在本研究的基礎(chǔ)上，進(jìn)一步開展數(shù)值模式調(diào)試的同時，設(shè)計更為適合的數(shù)據(jù)同化方案，以提高狀態(tài)場和參數(shù)估計的效果。