基于擴散濾波的多尺度三維變分研究

李冬，王喜冬,2，張學峰,2，吳新榮,2，李威，韓桂軍

（1.國家海洋局海洋環(huán)境信息保障技術重點實驗室 國家海洋信息中心， 天津 300171；2.中國科學院南海海洋研究所LED重點實驗室，廣東 廣州 510301）

基于擴散濾波的多尺度三維變分研究

李冬1，王喜冬1,2，張學峰1,2，吳新榮1,2，李威1，韓桂軍1

（1.國家海洋局海洋環(huán)境信息保障技術重點實驗室 國家海洋信息中心， 天津 300171；2.中國科學院南海海洋研究所LED重點實驗室，廣東 廣州 510301）

將擴散方程引入三維變分分析，揭示了傳統(tǒng) 3D-VAR不能有效提取多尺度觀測信息的根本原因，即觀測導致的目標函數(shù)梯度在空間的不連續(xù)分布。將擴散濾波融入基于梯度的最優(yōu)化算法，發(fā)展了基于擴散濾波的多尺度3D-VAR。海表面溫度數(shù)據(jù)同化試驗結果表明，新方法可從長波至短波有效地提取多尺度的觀測信息。

三維變分；數(shù)據(jù)同化；各向異性；遞歸濾波；擴散方程

Abstract: The diffusion equation is introduced into 3D-VAR analysis in this paper.As it reveals, the root cause of the inefficiency of the traditional 3D-VAR scheme in extracting multi-scale information resolved by observations lies in the discontinuous distribution of gradient of the cost function caused by the inhomogeneous distribution of observations.A multi-scale 3D-VAR scheme is developed by incorporating diffusion filtering into gradient-based minimization algorithm.The new scheme is applied to a two-dimensional idealized SST assimilation experiment with real AVHRR SST observations, which shows a superior performance and the multi-scale observational information, and from longer to shorter wavelengths, such observational informationcan be extracted successively by this scheme.

Keywords: 3D-VAR; data assimilation; anisotropic; recursive filter; diffusion equation

1 引 言

氣象、海洋三維變分數(shù)據(jù)同化（3D-VAR）通?？蓺w結為如下目標函數(shù)的最小化問題（Lorenc 1986; Couriter 1997）

在遞歸濾波3D-VAR中，遞歸濾波的實質其實是高斯濾波，用以對背景誤差協(xié)方差矩陣的作用進行模擬。由于高斯濾波最早來源于熱擴散方程的求解，本文將從擴散濾波的角度認識3D-VAR，由此分析3D-VAR不能有效提取多尺度觀測信息的根本原因，并將擴散濾波融入基于梯度的最優(yōu)化算法，發(fā)展基于擴散濾波的多尺度3D-VAR。

文章第二部分簡單介紹遞歸濾波3D-VAR；第三部分引入擴散方程，從擴散濾波的角度認識3D-VAR，并分析其不能有效提取多尺度觀測信息的原因；第四部分提出基于擴散濾波的多尺度3D-VAR；第五部分將多尺度3D-VAR用于海表面溫度（SST）同化試驗，第六部分給出本文結論。

2 遞歸濾波3D-VAR簡介

3D-VAR中的背景場誤差協(xié)方差矩陣通常采用高斯函數(shù)進行描述，例如，對二維情形的各項同性的背景誤差相關，B可表達為（Daley, R.,1991）:

為了避免B的求逆及相關運算帶來的巨大的內(nèi)存與計算開銷，并加快3D-VAR中優(yōu)化算法的收斂速度，定義新的控制變量（Lorenc 1988; Derber 和Rosati 1989）：

考慮到BT=B, (3) 可改寫為：

式中：α為濾波參數(shù)，決定了低通濾波器的帶寬。α越大，帶寬越窄，只有低頻（長波）信息可以通過；α越小，帶寬越寬，使得較高頻率（短波）的信息也可通過。理論上可以證明，當 ，u為無限長序列，時，上述遞歸濾波器可逼近于一個一維的各向同性高斯濾波器。對高維情形，由各向同性高斯濾波的可分離性，可轉化為多個一維遞歸濾波，例如，對二維情形，可以先對每一行進行濾波，再對每一列濾波。

遞歸濾波法，不需顯式地構造B矩陣，而是通過選取濾波參數(shù)間接的反映B的信息，可顯著減少因B矩陣的存儲和相關運算帶來的內(nèi)存和計算開銷。但在實際應用中，遞歸濾波方法最大的困難在于濾波參數(shù)的選取。以各向同性遞歸濾波為例，若選取的濾波參數(shù)α過大，則同化結果過于平滑，無法提取觀測的短波信息；相反，若α過小，又不能準確地提取長波信息，于是短波信息也不可能很好的提取。

3 從擴散濾波的角度認識和分析3D-VAR

遞歸濾波3D-VAR中，遞歸濾波的實質是高斯濾波。由于高斯濾波最早即來源于擴散方程的求解，下面我們基于擴散方程，對傳統(tǒng)3D-VAR進行分析，以探尋傳統(tǒng)3D-VAR難以提取多尺度觀測信息的根本原因。

3.1 3D-VAR中引入擴散濾波

高斯濾波可以通過求解擴散方程實現(xiàn)?？紤]如下一維擴散方程的柯西問題：

以看出，擴散系數(shù)a和積分時間t決定了濾波器的頻窗寬度（為4σ），也即代表了尺度。a越大，擴散越快，濾波越強，結果越平滑，得到的是低頻信息；同理，積分時間t越長，結果也越平滑。顯然，在積分時間固定的情形下，濾波的強弱僅取決于擴散系數(shù)的大小。另外，需要說明的是，3D-VAR中沒有時間的概念，因此這里的并不是真正的時間，而是為用擴散方程進行濾波而引入的“虛擬”的概念。

推廣至二維有限區(qū)域，擴散方程可改寫為：

3.2 伴隨方程及目標函數(shù)的梯度

在基于梯度的優(yōu)化算法中，需要知道目標函數(shù)關于控制變量的梯度。引入了熱擴散方程，便可以利用變分伴隨方法推導出梯度的表達式。通常來講，應該從離散的擴散方程直接推導出離散的伴隨方程和目標函數(shù)梯度，但限于篇幅，本文僅給出連續(xù)形式的伴隨方程及梯度表達式，并略去推導過程。同時，為便于敘述，假設觀測點位于分析網(wǎng)格點上，O為單位矩陣，并暫時略掉背景項?？梢宰C明，基于上述假設，伴隨方程有如下形式：

3) 計算 f。

5) 采用合適的最小化算法優(yōu)化w。

6) 轉至步驟2），直至滿足收斂條件。

3.3 3D-VAR不能提取多尺度觀測信息的原因分析

以下將基于擴散濾波分析3D-VAR不能有效提取多尺度觀測信息的原因。

顯然，伴隨模型（8）也是擴散方程，因此梯度的求解過程實際上是對進行高斯濾波，濾波的強弱由擴散系數(shù)所決定。由于在有觀測的分析網(wǎng)格點上而在無觀測的網(wǎng)格點上，所以在空間的分布存在“躍變”。當濾波參數(shù) a比較大時，濾波效果比較強，因此可以濾掉在空間的“躍變”，使得梯度比較平滑并保持物理上的連續(xù)性，但濾波參數(shù)大，同時意味著將濾除許多有價值的短波信息；而當a比較小時，則無法濾除的這種“躍變”，從而導致梯度在空間上的不連續(xù)性分布。顯然，這種現(xiàn)象并不是物理上真實存在的，而是由觀測數(shù)據(jù)的空間分布造成的。如果將此梯度直接應用于優(yōu)化算法（例如最速下降法、共軛梯度法、L-BFGS方法等）進行三維變分分析，在觀測稀疏的區(qū)域將得到錯誤的結果。

以最速下降法為例，其以在當前迭代點的梯度負方向作為下一次迭代的搜索方向，即：其中搜索步長可由線性搜索算法計算得到。由于梯度存在“躍變”，因此新的迭代值也會存在“躍變”，如此迭代下去，最終的分析結果在觀測稀疏的區(qū)域會失去觀測的長波特征。其他通用的優(yōu)秀的最優(yōu)化算法，如共軛梯度法、L-BFGS方法等，其設計也主要基于收斂速度的考慮，并不考慮尺度信息，因此應用于三維變分分析時存在與最速下降法同樣的問題。

因此，由觀測引起的目標函數(shù)梯度的空間不連續(xù)分布，是導致傳統(tǒng)3D-VAR（包括相關尺度法、遞歸濾波法）無法提取多尺度觀測信息的根本原因。

4 多尺度擴散濾波三維變分

基于擴散濾波思想，我們對3D-VAR進行了改進，將觀測的尺度信息融入最優(yōu)化算法中，設計了多尺度三維變分分析方法，以下簡稱多尺度3D-VAR。算法基本步驟如下：

3) 計算 f。

7) 減小擴散系數(shù)b。

8) 轉至步驟2），直至滿足收斂條件。

下面對此算法進行簡要的分析：

需要指出的是，雖然此算法中的每一迭代步都對梯度進行各向同性擴散濾波，但由于濾波系數(shù)隨迭代過程變化，依次提取不同尺度信息，因此整個優(yōu)化過程的累積效果，實際帶有非均勻各向異性濾波的特征，即可以在觀測稀疏的區(qū)域及其周邊區(qū)域進行較強的濾波以提取長波信息，而在觀測數(shù)據(jù)密集的區(qū)域進行較弱的濾波以便不丟失短波信息，本文后面的試驗結果也證明了這一點。

5 海表面溫度數(shù)據(jù)同化試驗

為檢驗多尺度3D-VAR提取觀測信息的能力，我們分別將遞歸濾波3D-VAR及多尺度3D-VAR應用于SST二維數(shù)據(jù)同化試驗，并對兩者進行比較。

5.1 試驗數(shù)據(jù)及參數(shù)設置

（1）“真實”場、觀測數(shù)據(jù)及相關參數(shù)

分析區(qū)域為200?E～220?E，30?N~40?N，分析網(wǎng)格的分辨率為 0.25?×0.25?。圖1a為我們假設的“真實”場，來源于2001年2月2日的AVHRR SST數(shù)據(jù)（AVHRR Pathfinder Version 5.0）。起初，在分析區(qū)域內(nèi)隨機選取500個觀測, 然后人為地去掉位于 205?E～209?E，35?N～38?N區(qū)域內(nèi)的觀測，剩余466個觀測（其分布見圖1b），下面的試驗中將用這些觀測對真實場進行重建，以此檢驗多尺度3D-VAR在觀測稀疏的區(qū)域提取觀測信息的能力。假設觀測誤差不相關，故R為對角陣，這里假定為單位陣。觀測算子H取為雙線性插值算子。

（2）多尺度3D-VAR的差分格式及相關參數(shù)

多尺度 3D-VAR中擴散方程（6）的差分格式采用ADI格式（交錯方向隱格式），可以證明，離散的伴隨方程應采用LOD格式(局部一維格式)，而且擴散方程及伴隨方程的差分格式絕對穩(wěn)定，限于篇幅，證明此處略。擴散系數(shù)取0.2，擴散系數(shù)取為下面的高斯函數(shù):

（3）優(yōu)化算法

遞歸濾波3D-VAR中的優(yōu)化算法采用L-BFGS (Liu and Nocedal, 1989)方法；多尺度3D-VAR中采用的線性搜索算法基于 Jorge J.More和 David J.Thuente的研究（Jorge J.More和David J.Thuente，1992）。

5.2 試驗結果

圖2a, 2b, 2c, 2d, 2e為遞歸濾波3D-VAR的分析結果，其濾波系數(shù)分別取α=0.2，0.3，0.4，0.5，0.6；圖2f為多尺度3D-VAR的分析結果，擴散系數(shù)取為 0.2。與“真實”場（圖1a）相比較可以看出，在遞歸濾波3D-VAR中，α的選取是非常困難的，若α過小，在觀測稀疏的區(qū)域難以提取觀測的長波信息；若α過大，會犧牲掉許多短波觀測信息。然而，在多尺度3D-VAR中，無論長波還是短波觀測都能很好地被提取出來。

為了更清楚地顯示多尺度3D-VAR提取多尺度觀測的過程，圖4給出了迭代過程中分析場及下降方向的變化，迭代步分別為第10，20，40，60，80，160步。為便于比較，圖3還給出了遞歸濾波3D-VAR的相應結果，其優(yōu)化算法采用最速下降法，迭代步分別為第10，20，40，60步，濾波系數(shù)α取為0.2?？梢钥闯觯谶f歸濾波 3D-VAR 中，當α比較小時，在無觀測（或觀測稀疏）的區(qū)域，由于下降方向（采用最速下降法時其下降方向即為負梯度方向）缺乏長波結構（圖3b，3d，3f，3h），從而導致分析場（圖3a，3c，3e，3g）同樣在此區(qū)域不能反映觀測場的長波。事實上，若遞歸濾波3D-VAR中采用其它的優(yōu)化算法，例如 L-BFGS、共軛梯度法等，依然會存在同樣的問題。而圖4a，4c，4e，4g，4k，4m，表明，多尺度3D-VAR方法可以從長波至短波依次提取觀測信息，這是由于算法中的下降方向（圖4b，4d，4f，4h，4l，4n）是對梯度進行高斯濾波得到的，彌補了梯度在觀測稀疏的區(qū)域存在的缺陷。

圖1 真實海溫場（a）及觀測的分布（b）Fig.1 (a): 圖 True temperature field (b): Distribution of observations

圖2 遞歸濾波3D-VAR及多尺度3D-VAR的分析結果。(a), (b), (c), (d), (e)分別為濾波系數(shù)α=0.2, 0.3, 0.4, 0.5, 0.6時遞歸濾波3D-VAR的分析結果；(f)為多尺度3D-VAR當擴散系數(shù)a= 0.2時的分析結果Fig.2 Analysis results of the recursive filter 3D-VAR and the multi-scale 3D-VAR.(a), (b), (c), (d) and (e) are for the recursive filter 3D-VAR with α=0.2, 0.3, 0.4, 0.5, and 0.6 respectively; (f) is for the multi-scale 3D-VAR with a=0.2.

圖3 采用最速下降法的遞歸濾波3D-VAR的分析場（左列）及目標函數(shù)下降方向（右列）。其中，濾波系數(shù)α=0.2，迭代步數(shù)分別為10, 20, 40, 60Fig.3 Analysis results (left) and the descent direction (-g) (right) at iteration 10, 20, 40 and 60 using the steepest descent algorithm for the recursive filter 3D-VAR

圖4 多尺度3D-VAR的分析場（左列）及目標函數(shù)下降方向（右列）。其中，擴散系數(shù)a=0.2，迭代步數(shù)分別為10, 20, 40, 60, 80, 160Fig.4 Analysis results (left) and the descent direction (right) at iteration 10, 20, 40, 60, 80 and 160 using the multi-scale 3D-VAR ( a=0.2)

6 結 語

本文將擴散濾波引入三維變分分析，并與最優(yōu)化算法相結合，提出了基于擴散濾波的多尺度三維變分方法，克服了傳統(tǒng)三維變分難以提取多尺度觀測信息的缺陷。應用于二維海溫SST同化試驗，效果良好。本文主要結論概括如下：

(1)傳統(tǒng)三維變分難以很難有效提取多尺度的觀測信息。由背景誤差的非均勻和各向異性特征，理論上講，3D-VAR中的相關尺度或濾波系數(shù)也應為非均勻及各向異性，且應根據(jù)誤差統(tǒng)計確定。但在實際應用中，由于我們無從確切地知曉物理量的真實狀態(tài)，因而據(jù)此進行的誤差統(tǒng)計也很難保證其準確性和合理性。而若簡單地以各向同性濾波器（如各向同性的遞歸濾波器）來模擬背景誤差協(xié)方差的作用，則不能根據(jù)觀測自動調整濾波行為，無法提取多尺度觀測信息。

(2)3D-VAR中求解梯度的過程實際是濾波過程。當濾波系數(shù)很小時，在優(yōu)化算法迭代的每一步，梯度在物理意義上代表當前解的觀測余量。因而，縱觀優(yōu)化算法的整個迭代過程，所有的觀測信息實際上都包含在這些梯度之中，從而使得我們可將濾波與優(yōu)化算法相結合，對梯度進行濾波以提取需要的觀測信息，這是本文算法的基礎。

(3)由觀測數(shù)據(jù)引起的目標函數(shù)梯度在空間的不連續(xù)分布，是導致傳統(tǒng)三維變分方法無法提取多尺度觀測信息的根本原因。如果將此梯度直接應用于傳統(tǒng)的最優(yōu)化算法（例如最速下降法、共軛梯度法、L-BFGS方法等）進行三維變分分析，在觀測稀疏的區(qū)域有可能導致錯誤的結果，其原因是這些優(yōu)秀算法的設計主要基于收斂速度的考慮，并不考慮尺度信息。

(4)本文提出的多尺度3D-VAR方法，從另一個角度看實際也是一種適用于多尺度變分分析的優(yōu)化算法。在迭代的每一步，對梯度進行高斯濾波，作為新的線性搜索方向，而高斯濾波算子的正定性保證了此方向必為下降方向。濾波系數(shù)隨迭代次數(shù)的增加而減小，從而可從長波至短波依次提取多尺度的觀測信息。雖然，與現(xiàn)在通用的高效優(yōu)化算法（如L-BFGS）相比，收斂速度尚有欠缺，但至少啟發(fā)我們：適用于變分分析的優(yōu)化算法既要考慮算法的收斂速度，也應考慮尺度信息，需要兩者之間的平衡。

(5)多尺度 3D-VAR雖然在每一迭代步采用的是均勻的各向同性濾波，但從整個優(yōu)化過程的累積效果看，實際帶有非均勻各向異性的特征，因此即使在觀測稀疏的區(qū)域也能捕捉到觀測的長波信息，同時又不會損失觀測密集區(qū)域的短波信息。

[1]Behringer D W，Leetmaa M J.A.An improved coupled model for ENSO prediction and implications for ocean initialization.Part I: the ocean data assimilation system[J].Mon.Wea.Rev，1998，126：1013-1021.

[2]Courtier P.Variational methods[J].Meteor.Soc.Japan，1997，75：211–218.

[3]Derber J，Rosati A.A global oceanic data assimilation system[J].Phys.Oceanogr，1989，19：1333–1347.

[4]Daley R.Atmospheric Data Assimilation[M].Cambridge University Press，1991：457.

[5]More J J，Thuente D J.On line search algorithms with guaranteed sufficient Decrease[R].Mathematics and Computer Science Division Preprint MCS-P330-1092. Argonne National Laboratory(Argonne,IL)，1992.

[6]Hayden C M，Purser R J.Recursive filter for objective analysis of meteorological fields: Applications to NESDIS operational processing[J].Appl.Meteor，1995，34：3–15.

[7]Huang X Y.Variational analysis using spatial filters[J].Mon.Wea.Rev，2000，128：2588–2600.

[8]Liu D C，Nocedal J.On the limited memory BFGS method for large scale optimization methods[J].Mathematical Programming，1989，45：503-528.

[9]Lorenc A C.Analysis methods for numerical weather prediction[J].Quart.Roy.Meteor.Soc，1986，112：1177–1194.

[10]Lorenc A C.Optimal nonlinear objective analysis[J].Quart.J.Roy.Meteor.Soc，1988， 114：205-240.

[11]Lorenc A C.Interactive analysis using covariance functions and filters[J].Quart.J.Roy.Meteor.Soc，1992，118：569–591.

[12]Purser R J，Wu W S，Parrish D.Numerical aspects of the application of recursive filters to variational statistical analysis.Part I: Spatially homogeneous and isotropic Gaussian covariances[J].Mon.Wea.Rev，2003a，131：1524–1535.

[13]Purser R J，Wu W S，Parrish D.Numerical aspects of the application of recursive filters to variational statistical analysis.Part II: Spatially inhomogeneous and anisotropic general covariances[J].Mon.Wea.Rev，2003b，131：1536–1548.

Multi-scale 3D-VAR based on diffusion filter

Li Dong1, Wang Xi-dong1,2, Zhang Xue-feng1,2, Wu Xin-rong1,2, Li Wei1, Han Gui-jun1

（1.Key Laboratory of Marine Environmental Information Technology, SOA, National Marine Data and Information Service, Tianjin 300171, China; 2.Key Laboratory of Tropical Marine Environmental Dynamics(LED),South China Sea Institute of Oceanology, Chinese Academy of Sciences, Guangzhou 510301,China)

P731.22

A

1001-6932(2011)02-0165-08

2009-11-23；收修改稿日期：2010-12-07

國家重點基礎研究發(fā)展計劃課題(2007CB816001)和國家自然科學基金項目（40776016、41030 854和40906016）聯(lián)合資助。

李冬(1974-)，男，研究員，現(xiàn)主要從事海洋數(shù)據(jù)同化、并行計算及科學數(shù)據(jù)可視化研究。電子郵箱：lidong2003@gmail.com。