球面小波背景誤差協方差模型的設計

曹小群，張衛民，宋君強，劉柏年

國防科學技術大學計算機學院，長沙 410073

球面小波背景誤差協方差模型的設計

曹小群，張衛民，宋君強，劉柏年

國防科學技術大學計算機學院，長沙 410073

在全球變分資料同化系統中設計和實現了基于球面小波的背景誤差協方差（B）模型。引入框架理論構造了球面小波函數；設計了一個基于球面小波變換的全球B矩陣模型；分別通過理想數值試驗和在全球變分同化系統中的實現對新模型的有效性進行了驗證。試驗結果表明：基于球面小波的背景誤差協方差模擬方法能夠克服由有限背景誤差樣本引入的取樣噪聲，能估計出真實的背景誤差相關函數；在全球變分同化系統中新模型能夠模擬出在物理和動力上正確和有效的背景誤差結構函數。

背景誤差協方差；球面小波；變分資料同化；數值天氣預報

1 引言

以高性能計算為基礎的數值天氣預報（NWP）技術能有效克服一般預報方法產品不夠豐富、可用預報時效短、海量氣象觀測數據處理能力弱等缺點，已經成為氣象部門制作業務天氣預報的根本科學途徑。全球大氣資料同化是把全球范圍觀測資料和對大氣運動的動力學認識進行最優結合的一類信息融合方法，目的是產生出與所有已知信息（包括觀測資料、背景場及它們對應的誤差信息）最一致的大氣狀態場。大規模并行計算機性能的不斷發展有力地推動著全球中期數值天氣預報模式分辨率的不斷提高，當前條件下發展業務化的全球中小尺度數值天氣預報系統已成為可能，從而使全球中期數值天氣預報系統能為臺風和暴雨等災害性天氣的預報和監測提供更好的科學依據[1]。全球數值天氣預報模式分辨率的提高對提供初始場的資料同化[2-9]系統也提出了新要求：資料同化系統分辨率需要相應提高，同時背景誤差協方差模型必須能刻畫不同尺度（包括行星尺度、大尺度、中尺度和中小尺度等）大氣運動的誤差統計特征[9]。

在目前業務化的全球變分同化系統中，一般都采用譜空間對角化方法來模擬背景誤差協方差[9-10]，譜空間對角方法對所有大氣運動尺度的背景誤差不加分別地統一處理，其計算結果反映的是各種尺度的平均統計誤差特征。在這種情況下，小尺度或小振幅的大氣波動的背景誤差特征完全有可能被大尺度或大振幅的運動特征所掩蓋。譜空間對角化方法的另一個缺點是將相同的相關函數應用在所有的物理位置，不能夠描述實際垂直相關或水平相關隨空間的變化[9-10]；第二種模擬全球B矩陣的方法是物理空間格點方法：分別用經驗正交函數（EOF）和遞歸濾波模擬垂直相關和水平相關函數，然后把它們應用到全球模式格點空間上[11-12]。物理格點方法能表示B-矩陣在物理空間的信息，但會把同樣的相關應用到所有運動尺度上。綜上所述，如果只使用譜方法，就不能刻畫背景協方差在物理空間的變化特征；如果在格點空間表示B，就不能刻畫協方差隨波動尺度變化的特征。總之，單獨用譜方法或者格點空間方法都不能刻畫B的所有特征。

小波函數在譜和格點空間中同時具有局部性，能夠將相關函數同時表示為尺度和位置兩者的函數。正交小波函數具有光滑性、緊支撐、濾波和正交性等許多良好的性質，其適合解決定義域為矩形區域的問題[13-15]，曹小群等采用正交小波模擬了區域變分資料同化系統中的水平誤差相關函數[16]。但將正交小波分析擴展到球面區域是困難的：因為在球面上構建正交小波一般都存在“極點問題”[17]，所以構造球面小波需要采用引入新的數學工具。本文根據Fisher等[8-9]的思想，研究利用框架理論[17-18]構造能夠應用于業務化的全球氣象變分同化系統的球面小波函數，設計了全球變分資料同化中基于球面小波函數的B矩陣模型。

2 球面小波函數構造

如果給出一個被一組有限正交小波基表示的函數，對此函數在球面上作任意旋轉后，則新函數不能被同樣的小波基表示[17]。因此構造球面小波函數只能考慮非正交變換，適合的數學工具之一是框架。對于Hilbert空間Η中的一個函數族{φn;n∈Γ}，其中Γ是離散指標的可數集，如果存在常數0＜A≤B＜+∞使得下面不等式

對任意函數F∈H成立，則此函數族構成一個離散框架[17]。在式（1）中如果常數A=B，則稱此函數族為緊框架[17]。如果指標集Γ在某一個測度μ下是可測的，同時滿足下面的條件：

上面的敘述表明在框架的意義下，非正交函數族能夠精確地表示信號，能對函數進行分解和重構。利用球面上的緊框架條件可構造具體的球面小波函數。選擇一個波數節點序列{Nj;j∈Z}，且N0=0和Nj＜Nj+1。定義基函數Wj的譜系數為：

圖1 不同球面小波函數在譜空間的取值

同時，在上面描述的緊框架解釋中，可以認為fj在一個特定點(λ，φ)的值代表了小波基函數W(λ，φ，j)的系數。為了使這個解釋有意義，要求函數W(λ，φ，j)在點(λ，φ)周圍是局部化的。在圖2～3中分別顯示了第11個和第15個球面小波函數在譜空間和物理空間中值的分布圖。譜空間中的小波函數是總波數的函數，只在有限區間中非零，而在區間以外為零。物理空間中的小波函數是大圓距離（這里用角度表示）的函數，從圖中可以看出，函數呈現了很好的局部性。圖4顯示的是第10個和第11個球面小波函數在二維球面格點空間上的分布圖，從圖中可知，球面小波函數只在局部區域非零，呈現了很好的局部性。因此，每個球面小波函數的系數場只由球面上的有限格點上的值正確決定，例如，對于高斯網格，將由2Nj+1個經度格點和(2Nj+1+1)/2個緯度格點上值決定。

圖2 第11個球面小波函數在譜空間和物理空間中值的分布圖

圖3 第15個球面小波函數在譜空間和物理空間中值的分布圖

圖4 第10和11個球面小波函數在二維球面格點空間中的等值線分布圖

獲得球面小波函數后，球面上的任何一個函數或者全球大氣物理場F（例如，溫度、壓強、水平風場和濕度等）都可以通過小波函數族Wj進行分解和重構。F和每一個小波函數的卷積?將得到小波系數函數fj。由于總數波和空間尺度之間存在對應關系，因此每一個小波系數函數fj都可以被認為代表了一個特定的空間尺度。與球面小波函數卷積的結果使物理場中不同空間尺度的部分被分離，因為小波函數同時在物理空間和譜空間上是局部化的，從而在對應的系數fj所描述的特征中保留了它們的局部性質。每個函數fj和對應小波函數Wj的卷積之和能實現對函數或物理場F的重建。

3 基于球面小波的全球B矩陣模型

利用上面部分定義的球面小波函數Wj可以構造一個背景誤差協方差B的模型。因為球面小波函數同時在譜空間和物理空間中具有局部化特性，因此基于Wj的背景誤差協方差模型能提供同時在格點空間和尺度上變化的背景誤差特征量。又因為球面小波函數具有多尺度分析的能力，因此基于球面小波函數的B矩陣模型能刻畫不同尺度大氣運動的背景誤差特征。

全球變分資料同化是通過極小化目標泛函來估計大氣真實狀態的離散向量x[20-23]。

其中xb是對x的先驗（背景場）估計，y是觀測向量，H是把模式狀態量映射到觀測空間的算子。矩陣O和B分別是觀測和背景誤差協方差矩陣。

一般地，對方程（4）中定義的目標泛函直接進行最小化，在數值計算上是壞條件的[21-22]。因此通常的做法是定義一個控制向量v，利用它對背景誤差協方差矩陣進行對角化。通過一個控制向量轉換可以確定分析增量δx：

在控制變量v空間中，目標泛函可以寫成：

在方程（6）中B矩陣無需顯式表示，而是通過控制變換L被隱式地表示。可以很容易地證明，假如B=LLT，則極小化J(v)和J(x)將得到相同的大氣狀態估計量。因此在全球變分同化中背景誤差協方差矩陣不需要直接定義，實際應用中是將其定義為一系列控制變換，即通過選擇一個合適的控制變量v和變換矩陣L來實現的，對于變換矩陣L，無需假設為可逆的，甚至不要求是方陣[9]。

為了定義一個基于小波變換的B矩陣模型，定義控制向量vT=(，，…，，…，)，j=1，2，…，K標識不同尺度的小波空間。其中：

其中K表示不同氣象變量之間的物理變換矩陣，包含了平衡約束關系。Vj(λ，φ)表示在j尺度小波空間水平位置(λ，φ)上的垂直協方差矩陣，V(λ，φ)表示它的逆均方根矩陣。原則上，對每一個Wj?K-1δx，都必須在格點空間上指定一個垂直協方差矩陣。而在實際計算中，有必要減少垂直協方差矩陣的數目，方法是許多相鄰的格點使用同一個垂直協方差矩陣。

在目標泛函的極小化過程中不需要使用式（7），而只需要使用控制向量到分析增量的轉換關系式。從上面所定義的新控制向量和方程式（7），可以導出新的基于球面小波的控制變量轉換關系式：

通過上面的分析容易知道，基于球面小波的B矩陣模型所需要的背景誤差統計量如下：

（1）物理變換K中需要用到的統計回歸系數：速度勢和流函數之間的回歸系數；溫度和流函數之間的回歸系數；地面氣壓和流函數之間的回歸系數。

（2）每個控制向量在每個尺度和每個水平位置上的局地垂直相關矩陣Vj(λ，φ)，由于這部分統計參數需要巨大的資料存儲量，因此只能在相對變分同化系統更低分辨率的格點上進行垂直統計量的計算。在分析過程中，將低分辨率的統計量插值到同化系統的格點上。

和一般B矩陣模擬方法（例如，著名的WRF變分同化系統）相比，本文設計球面小波B矩陣模型具有以下不同。首先，普通模擬方法中控制變量轉換包含的子變換按照作用順序依次為：物理變換、垂直變換和水平變換；而基于球面小波的控制變量轉換包含的子變換按照作用順序為：物理變換、水平變換和垂直變換。后者的水平變換按照作用順序又分為：對格點場進行球諧變換；在譜空間中進行球面小波變換；小波變換后，進行逆球諧變換，將譜數據轉換成格點場。因此，球諧變換、球面小波變換和逆球諧變換三部分一起構成球面小波控制變量轉換中的水平變換。其次，需要的統計量不同，普通模擬方法中需要回歸系數、垂直特征值和特征向量以及水平功率譜等統計量。而球面小波控制變量轉換中需要統計量是：平衡回歸系數、同時依賴于小波尺度和水平位置的局地垂直相關協方差的均方根矩陣及其轉置矩陣。

4 數值試驗與分析

4.1 理想試驗

為了驗證球面小波全球B矩陣模型的有效性，首先將進行理想數值試驗：使用上面構造的球面小波函數模擬大圓上的協方差函數，從而說明球面小波的濾波性質。研究的地理區域是半徑為a的地球大圓。坐標x表示距離，變化范圍是[0，2πa]。試驗中只考慮大圓上一個物理場，例如一維溫度場。假設圓周上有均勻分布的121個格點xj，j=0，1，…，120，同時使用具有變化特征長度尺度的高斯函數來指定格點上真實的背景誤差相關函數Ct。特征長度尺度l(x)在個格點1和121上具有最大值：6.5倍格點距；在格點61上具有最小值：1.5倍格點距，其他格點上的特征長度尺度在最大值和最小值之間是線性變化的。Ct表示如下：

在已知Ct的情況下，可以通過隨機模擬方法產生背景誤差樣本，公式如下：

（1）對隨機模擬的誤差樣本{εn，n=1，2，…，50}進行統計，得到格點標準偏差D，然后對誤差樣本進行規則化=Dεn，n=1，2，…，50。

獲得統計量后可以計算出隱含的協方差矩陣，用公式表示為：

圖5 第30個格點上的相關函數比較

圖6 第90個格點上的相關函數比較

總之，通過上面的數值試驗說明了基于球面小波的背景誤差協方差模擬方法能夠克服由有限樣本數引入的取樣噪聲，表明了球面小波形式的背景誤差協方差在模擬真實誤差相關函數的有效性。雖然這里只對地球大圓上的一維問題進行了試驗，但結論同樣適應于球面上的二維和三維情形。

4.2 新模型在全球變分同化系統中的實現

為了檢驗球面小波背景誤差協方差矩陣模型的正確性，利用已經實現的全球四維變分資料同化系統（YH4DVAR）[22]進行了單個觀測資料同化試驗。首先將新的B矩陣模型實現到YH4DVAR中，同時采用NMC方法[23-25]估計了全球背景誤差統計量，在此基礎上利用最優化算法對目標泛函進行極小化求解。

在資料同化系統中，由單點觀測生成的分析增量隱含地表現了背景場誤差協方差矩陣B的作用，因此Thepaut等在1991年[20]引入了單點試驗方法來反映變分同化中B矩陣的結構函數。在變分資料同化系統中，如果假設觀測為單個模式變量，就有：

其中的yi是第i個網格點上的單個觀測，xi是從背景場計算得到的相關等價量，σo和σb分別是觀測和背景場誤差，Bi是背景場誤差協方差的第i列。從上式可以明顯看出，分析增量xa-xb是正比于背景場誤差協方差。單點觀測的信息傳播依賴于背景場誤差協方差矩陣，B矩陣在變分同化系統的信息傳播中起著重要作用。

圖7顯示的是YH4DVAR系統在采用新的背景誤差協方差模型后對北半球中緯度單個溫度觀測（經緯度為（122.78°W，53.90°N），高度為800 Hpa）進行同化后，產生的地面氣壓（圖7（a）），模式第7層上的溫度（圖7（b））、緯向風（圖7（c））、經向風（圖7（d））和風矢量（圖7（e））的水平分析增量。單點觀測試驗通過全球變分同化系統YH4DVAR可以產生全球范圍的增量，但在距離觀測較遠格點上的增量為零，即單個觀測資料對分析場的影響在空間上局部化的，因此上面所有的圖只顯示靠近觀測點區域上的增量。各個分析增量的分布在物理是一致的，具體解釋如下：在觀測位置上溫度的升高，引起氣壓降低，從而形成輻合風場，但由于受地轉科氏力的作用，形成在觀測位置以左為北風，觀測位置以右為南風的增量分布；低層的輻合將產生上升氣流，在高層引起氣體質量堆積，從而形成輻散風場，但受地轉科氏力的作用，形成觀測位置以左為南風，觀測位置以右為北風的增量分布，子圖圖7（e）顯示在模式第7層位于觀測氣壓層797.0 Hpa以上，從而形成了一個反氣旋式的風場。分析結果表明新的背景誤差協方差模型在物理和動力上是正確和有效的。

圖7中的另外一個顯著特點是緯向風和經向風的增量形狀和地轉協方差中的增量分布比較相似，但明顯具有各向異性特征，即相關結構函數的值隨兩點的剛性平移或旋轉而變化。這主要是因為基于球面小波的背景誤差協方差矩陣能夠模擬局地的、多尺度大氣運動的相關特征。從圖7中可以明顯看到，球面小波背景誤差協方差具有模擬細微相關特征的能力，而基于譜方法的背景誤差協方差模型只能均勻和各向同性地向觀測位置周圍傳播信息。從而表明，新模型能夠正確地刻畫YH4DVAR中實際背景誤差協方差矩陣的物理結構函數。

5 結束語

超高分辨率的全球數值天氣預報模式將具有多尺度模擬和預報能力，預報產生中將包含行星尺度、大尺度和中小尺度等天氣系統的運動信息。為了提供適合該類數值天氣預報模式的初始場，有必要在變分資料同化系統中構造能對不同尺度大氣波動的背景誤差特征分別進行刻畫的B矩陣模型。本文利用球面小波方法設計了一個新的全球B矩陣模型，首先引入框架理論構造了球面小波函數；其次構造了一個基于球面小波變換的全球B-矩陣模型；最后分別通過理想數值試驗和在全球變分同化系統YH4DVAR中的實現對新B矩陣模型的有效性進行了驗證。主要結論是：基于球面小波的背景誤差協方差模擬方法能夠克服由有限背景誤差樣本引入的取樣噪聲，能估計出真實的背景誤差相關函數；在全球變分同化系統中新B矩陣模型能夠模擬出在物理和動力上是正確和有效的背景誤差結構函數。

圖7 采用新模型后，對于單個溫度觀測資料，YH 4DVAR所產生的分析增量分布圖

[1]Bader D A.Petascale computing：algorithms and applications[M].New York：Chapman & Hall，2007.

[2]Zou X，Navon I M，Le Dimet F X.An optimal nudging data assimilation scheme using parameter estimation[J].Q J R Meteor Soc，1992，118：1163-1193.

[3]Navon I M.Practical and theoretical aspects for adjoint parameter estimation and identifiability in meteorology and oceanography[J].Dyn Atmos Oceans，1997，27：55-79.

[4]黃思訓，韓威，伍榮生.結合反問題技巧對一維海溫模式變分資料同化的理論分析及數值試驗[J].中國科學，2003，47：903-911.

[5]Talagrand O，Courtier P.Variational assimilation of meteorological observations with the adjoint and vorticity equation.Part I：theory[J].Q J R Meteorol Soc，1987，113：1311-1328.

[6]黃思訓，蔡其發，項杰，等.臺風風場分解[J].物理學報，2007，56（5）：3022-3027.

[7]鐘劍，費建芳，黃思訓，等.模式誤差弱約束四維變分同化研究[J].物理學報，2012，61（14）.

[8]Fisher M，Andersson E.Developments in 4D-Var and Kalman filtering[J].ECMWF Tech Memo，2001，347：38-65.

[9]Fisher M.Background error covariance modelling[C]//Proc of ECMWF Seminar on Recent Developments in Data Assimilation for Atmosphere and Ocean，2004：45-64.

[10]Derber J，Bouttier F.A reformulation of the background error covariance in the ECMWF global data assimilation system[J].Tellus，1999，51A：195-221.

[11]Purser R，Wu W S，Parrish D，et al.Numerical aspects of the application of recursive filters to variational statistical analysis.Part I：spatially homogeneous and isotropic Gaussian covariances[J].M on Wea Rev，2003，131：1524-1535.

[12]Purser R，Wu W S，Parrish D，et al.Numerical aspects of the application of recursive filters to variational statistical analysis.Part II：spatially inhomogeneous and anisotropic general covariances[J].Mon Wea Rev，2003，131：1536-1548.

[13]Jawerth B，Sweldens W.An overview of wavelet based multiresolution analyses[J].SIAM Rev，1994，36（3）：377-412.

[14]Fournier A.Introduction to orthonormal wavelet analysis with shift invariance：application to observed atmospheric blocking spatial structure[J].J Atmos Sci，2000，57：3856-3880.

[15]Deckmyn A，Berre L.Wavelet approach to representing backgrond error covariances in a limited area model[J].Mon Wea Rev，2004，133：1279-1294.

[16]曹小群，黃思訓，杜華棟.變分同化中水平誤差函數的正交小波模擬新方法[J].物理學報，2008，57（3）：1984-1989.

[17]Kaiser G A.Friendly guide to wavelets[M].Boston：Birkhauser，1994.

[18]Daubechies I.Ten lectures on wavelets，CBMS-NSF series on applied mathematics[M].Philadelphia，PA：SIAM，1992.

[19]Bannister R N.Can wavelets improve the representation of forecast error covariances in variational data assimilation?[J] Mon Wea Rev，2007，135：387-408.

[20]Thepaut J N，Courtier P.Four-dimensional data assimilation using the adjoint of a multi-level primitive equation model[J].Quart J Roy Meteor Soc，1991，117：1225-1254.

[21]Courtier P，Andersson E，Heckley W，et al.The ECMWF implementation of three dimensional variational assimilation（3D-Var）.Part I：formulation[J].Quart J Roy Meteor Soc，1998，124：1783-1808.

[22]張衛民，曹小群，宋君強.以全球譜模式為約束的四維變分資料同化系統YH4DVAR的設計和實現[J].物理學報，2012，61（24）.

[23]Parrish D F，Derber J C.The national meteorological center’s spectral statistical-interpolation analysis system[J].Mon Wea Rev，1992，120：1747-1763.

[24]Berre L.Estimation of synoptic and mesoscale forecast error covariances in a limited area model[J].Mon Wea Rev，2000，128：644-667.

[25]Fisher M，Courtier P.Estimating the covariance matrices of analysis and forecast error in variational data assimilation[J].ECMWF Tech Memo，1995，220：1-28.

CAO Xiaoqun, ZHANG Weimin, SONG Junqiang, LIU Bainian

School of Computer Science, National University of Defense Technology, Changsha 410073, China

Based on spherical wavelet, a new model of background error covariance（B）in global variational data assimilation system is designed and implemented. The frame theory is introduced, which is used to construct the basis functions of spherical wavelet. A model of B-matrix is designed based on spherical wavelet transform. The validity of the new method to model global B-matrix with spherical wavelet is demonstrated by numerical and simple experiments. The results of experiments show that the method has the ability to filter the sample noise from limited number of background error member effectively and to estimate the background error correlation functions correctly.

background error covariance; spherical wavelet; variational data assimilation; numerical weather prediction

CAO Xiaoqun, ZHANG Weimin, SONG Junqiang, et al. Design of background error covariance’s model with spherical wavelet. Computer Engineering and Applications, 2014, 50（17）：49-55.

A

TN911.72

10.3778/j.issn.1002-8331.1308-0296

國家自然科學基金（No.41105063）；高分青年創新基金（No.GFZX 04060103-5-19）。

曹小群（1980—），男，博士，副研究員，研究領域為計算機應用技術、數值天氣預報技術、反問題；張衛民（1966—），男，博士，研究員，研究領域為計算機應用技術、數值天氣預報技術；宋君強（1962—），男，院士，研究員，研究領域為數值天氣預報技術、高性能計算；劉柏年（1985—），男，助理研究員，研究領域為計算機應用技術、數值天氣預報技術。E-mail：caoxiaoqun@nudt.edu.cn

2013-08-22

2013-11-15

1002-8331（2014）17-0049-07

CNKI網絡優先出版：2014-01-26,http://www.cnki.net/kcms/doi/10.3778/j.issn.1002-8331.1308-0296.htm l