集合背景誤差方差中小波閾值去噪方法研究及試驗?

劉柏年皇群博 張衛民 任開軍 曹小群 趙軍

1)(國防科技大學海洋科學與工程研究院,長沙 410073)

2)(國防科技大學計算機學院,長沙 410073)

集合背景誤差方差中小波閾值去噪方法研究及試驗?

劉柏年1)2)?皇群博1)2)張衛民1)任開軍1)曹小群1)趙軍1)

1)(國防科技大學海洋科學與工程研究院,長沙 410073)

2)(國防科技大學計算機學院,長沙 410073)

(2016年3月26日收到;2016年10月20日收到修改稿)

背景誤差方差的集合估計值中帶有大量采樣噪聲,在應用之前需進行降噪處理.區別于一般的高斯白噪聲,采樣噪聲具有空間和尺度相關性,部分尺度上的噪聲能級遠大于平均能級.本文針對背景誤差方差中采樣噪聲的特征,引入小波閾值去噪方法,并根據截斷余項的小波系數分布特征發展了一種計算代價很小,能自動修正閾值的算法.一維理想試驗結果表明,該方法能濾除大量采樣噪聲,提高背景誤差方差估計值的精度.相對于原來的小波閾值方法,修正閾值后減少了因部分尺度上噪聲能級過大導致的殘差,去噪后的RMSE減少了13.28%.將該方法應用在實際的集合資料同化系統中,結果表明,小波閾值方法優于譜方法,閾值修正后能在不影響信號的前提下增大小波去噪強度.

小波,集合資料同化,背景誤差方差,去噪

1 引 言

數值天氣預報可用一個非線性方程組表示,初始場的精度在很大程度上決定了預報的準確性.為了得到高質量的初始場,變分資料同化需要融合各種觀測資料和背景信息[1].其中背景誤差協方差矩陣(B)在資料同化中起到了誤差結構分布、信息傳播、權重分配等作用,因此資料同化的首要任務是精確定義出B[2].在信息量不足的前提下,采用靜態的、均一的和各向同性的B模型能有效節省大量計算資源并簡化運算,但無法反映出大氣的時空變化特征[3,4].因此大量研究開始嘗試采用不同方法來獲取流依賴B[5-7].高性能計算機和大規模并行計算為集合資料同化(EDA)的發展提供了重要支撐,EDA通過擾動觀測資料、邊界條件和預報模式構造一組能夠反映B統計信息的低分辨率擾動集合,并以此估算流依賴的B.EDA被認為是解決當前確定性資料同化問題的一種有效途徑.國際上幾大數值預報業務中心也正在開發和完善EDA業務系統,如歐洲中期天氣預報中心(ECMWF)[8]、法國氣象局[9,10]等.

EDA的關鍵是選擇合適的集合成員個數.由于EDA成員是通過隨機擾動方法構造的,利用有限EDA成員統計背景誤差方差時,估計值中通常包含有隨機采樣噪聲,估計精度與集合成員數的均方根成正比[11].一方面,業務時效性要求EDA在確定性資料同化之前完成所有成員的計算和B的統計,在這一前提下,成員個數不能無限擴大;另一方面,EDA通過隨機擾動方法構造出多個能夠反映B統計信息的擾動集合,集合成員個數越多,引入的隨機采樣噪聲越小,統計值就越精確,因此成員個數不能過少.綜合以上分析,一種解決方案是擴大計算資源和并行規模并利用FPGA,GPU等加速計算技術,增加盡可能多的樣本數,但由于估計精度的低收斂性限制了這一方案的實際效率;另一種行之有效的方案是發展計算量小的去噪方法,通過濾除隨機采樣噪聲來提高估計精度.

針對以上需求,Raynaud等[12]于2008年根據噪聲和信號的特征尺度發展了一種基于空間平均的低通濾波方法,以空間加權的方式濾除具有頻率高、尺度小的隨機采樣噪聲.經過上述方法處理后的10個樣本背景誤差方差估計值在精度上與50個樣本的直接統計結果相當,但這種處理噪聲方法的效能受大氣變量、高度層等因素影響,且平滑的空間尺度需要根據天氣形勢進行人為調整和優化,不適合業務化.2009年,Raynaud等[10]改進了濾波算法,基于氣候態B矩陣得到了近似噪聲能量譜,發展了一種能自動估算截斷波數的譜濾波方法.其基本思想是根據信號和噪聲能量在譜空間的分布特征定位出信噪分離的截斷波數,由于噪聲能量主要集中在波數較大的小尺度空間上,因此構造一個低通濾波器,濾除大于截斷波數的高頻小尺度噪聲,對小于截斷波數的大尺度噪聲,則采用類似于Wiener的線性處理方法[13].該方案一方面降低了對集合樣本數的要求,減少了計算量,另一方面又顯著提高了集合方差的估算精度[8,10,12,14].文獻[15]等將該方法成功應用到了我國自主開發的集合資料同化系統中,使10個成員集合估計值的濾波結果好于30個成員集合估計值.但是,譜濾波方法也存在一定的局限性.首先,噪聲能量譜是利用氣候態B估算得到的,精度和正確完全由B和集合成員數決定.由于對B采用了靜態近似處理,由此統計的隨機采樣噪聲能量譜也是靜態的定常量,不適用于快速發展的天氣系統,如臺風和深對流天氣;其次,相同尺度上作用的濾波系數是相同的,這等同于對格點空間中相同尺度的平滑處理在全球范圍都是一致的,不能反映出信號的局地特征.

針對上述方法的局限性,本文在集合背景誤差方差中引入了具有譜和空間局地化特性的小波閾值去噪方法[16].根據小波信號特征,通過迭代算法得到閾值,避免了譜濾波方法中噪聲能量譜的近似處理和靜態假設.在此基礎上,根據集合背景誤差方差中采樣噪聲具有的空間和尺度相關特征,設計了一種能自動修正閾值的改進算法,可減少因部分尺度上噪聲能級過大導致的殘差,進而改進濾波效果.最后在一維理想模型和實際的集合資料同化系統中測試了該方法的魯棒性.

2 去噪原理

2.1 譜濾波原理

EDA通過擾動觀測資料、海表溫度和模式物理傾向項來表征觀測、邊界及模式物理過程的不確定性.多個表征誤差分布特征的擾動初值,經模式積分得到的預報集合可統計出流依賴的背景誤差[17-19].但是EDA作為一種純粹的隨機方法,有限個成員統計得到的結果中包含了隨機采樣噪聲.假定背景誤差協方差服從高維的高斯分布,則采樣噪聲的協方差E[Ge(i)Ge(j)]和集合誤差協方差矩陣近似滿足以下關系[10,20]：

其中E[·]表示數學期望,是格點i的背景誤差方差的集合估計值,Ge(i)是格點i上的采樣噪聲,N是集合樣本個數.由此計算得到噪聲與背景誤差之間的Daley[21]長度尺度關系滿足：

LGe(i),Lεb(i)分別為格點i上噪聲和背景誤差長度尺度,c為相關函數,r為格點之間的距離.上式表明,隨機采樣噪聲的相關長度尺度總小于背景誤差方差的長度尺度.誤差方差傾向于在更大長度尺度上變化.譜濾波方法根據這一特性,首先利用等式(1)近似估算出隨機采樣噪聲的能量譜Se,然后將方差的集合估計值轉換到譜空間,并以(譜濾波的原始輸入信號)表示.由于隨機采樣噪聲的相關長度尺度總小于背景誤差方差的長度尺度(見等式(2)),因此,可利用以下低通濾波器[10]來削弱輸入信號中的隨機采樣噪聲的強度：

其中ρ(n)是低通濾波器的濾波系數,np為譜空間中的波數,截斷波數Ntrunc為信號和噪聲能量譜開始分離的波數.譜濾波的處理流程如圖1所示.

圖1 譜濾波模塊處理流程Fig.1.Processes of the spectral filtering.

譜濾波方法能根據譜空間中信號、噪聲能量譜的分布情況自動計算出截斷波數并實現濾波,具有計算成本低、濾波效率高等優點[8,10,12].當信號具有空間緩變、各向同性、短距離相關特性時,濾波效率更為顯著.目前ECMWF和Mete-France的業務化集合資料同化系統都采用了這種方法[8,10,12,14].但是該方法也存在一定的局限性：1)譜濾波方法采用了確定性同化系統中的簡化、靜態B模型統計噪聲能量譜,并引入了一定的假設和近似,這種處理方式降低了噪聲能量譜的估計精度;2)噪聲能量譜的統計完全獨立于輸入信號,且作為一個時常量應用到后續的濾波過程中,不能反映背景誤差方差及隨機采樣噪聲隨天氣形勢變化的特征;3)譜濾波的實質是在相同波數上,作用一個濾波系數.在格點空間,這等同于在同一尺度空間上乘上一個全球相同的加權平均系數.這種各向同性、均一的濾波方式會忽略信號的局地變化,一些十分重要的特征信號可能被平滑或弱化.

2.2 小波閾值去噪原理

小波閾值去噪是一種常用的信號處理方法,它將小于閾值的小波系數置零后,再重構出信號,因此閾值的選取是一個十分關鍵的問題.Donoho和Johnstone[16]將閾值表達為集合成員個數和噪聲方差的函數.在缺乏信號的先驗信息條件下,小波閾值能減少定義在Holder和Besov空間中的誤差,MAD方法則是將小波的最小尺度系數模的中位數絕對偏差(MAD)作為噪聲的能級[20].對于高斯白噪聲,可采用以下迭代算法得到最佳閾值[22],具體步驟如下：

1)將原始信號X轉換到小波空間,并劃分為噪聲e和信號?兩部分,在循環迭代開始時置e=,即將所有信號當成噪聲;

3)進行小波閾值截斷

迭代算法如圖2所示,各個方向的白噪聲相互正交于圓內.圓對應最大的噪聲,即閾值.由于迭代開始時噪聲方差非常大,大部分小波系數都囊括在圓內.經第一次去噪后,排除了那些模大于閾值T0的系數,余下的系數則構成新的噪聲,用于計算下一個閾值T1.

圖2 噪聲方差和閾值迭代算法示意圖 粗實線圓代表T0,虛線代表T1,箭頭表示小波系數Fig.2.Illustration of the recursive algorithm for estimating the noise variance and the threshold.The estimated thresholdsTandT1are represented by the bold and dashed spheres respectively.The arrows represent a selection of wavelet coefficientsi,j.

小波閾值方法等價于用一個適合輸入信號局地變化特征的濾波器來估計真實信號.其事實依據是,函數f在尺度j和位置xj(i)的小波轉換,度量了函數f在xj(i)附近的變率,xj(i)的尺度與j成正比.快速變化的信號能使小尺度上的小波系數變大.小波閾值方法通過將小于閾值T的系數置0,構造一個依賴于小波系數的自適應濾波器.較大的系數表明函數f在小尺度范圍內的變化劇烈,這部分系數予以保留避免平滑掉.而滿足的系數則表示f變化較為平緩,通過置0后濾除.

背景誤差方差的集合估計值中,噪聲具有空間和尺度相關性,不能簡單地當作高斯白噪聲處理.譜空間中信噪比隨波數的變化可以用簡單的正弦函數表示[10],基于此構造的形如等(2)式的低通濾波器可以有效地減少采樣噪聲在各個尺度上的絕對量.在小波空間中為了濾除具有相關性的噪聲,一種處理方法是對每個尺度使用不同的閾值[23]：

其中σ(j)表示尺度上的標準偏差,nj=2j表示尺度j上小波系數個數.由于個別尺度對應的nj過小,統計得到的σ(j)具有較大的誤差.另一種處理方法是采用新的全局閾值,如對迭代閾值TA進行調整：

通過改變α值,使滿足σ(j)＞σw′上的部分過大的噪聲小波系數也能被閾值TP置0.Pannekoucke等[24]給出了初步實驗結果驗證了該方法的有效性,但由于噪聲具有非高斯性,理論建模和推導都十分復雜,目前沒有很好的方法估算α.

閾值TP在調整過程中,需要綜合考慮對立因素：既要保證部分能級較大的噪聲能被閾值置0,這要求TP要足夠大;同時TP又不能無限大,確保信號的小波系數不受TP的影響.如圖3所示,小圓表示高斯白噪聲,橢圓表示具有空間和尺度相關的非高斯噪聲,在某些尺度上噪聲能級大于平均能級,表現為橢圓的兩端溢出了小圓范圍.大圓則是以最大噪聲能級為半徑.如果簡單地將背景誤差方差中的采樣噪聲當作高斯白噪聲處理,采用σ=σw計算得到的閾值將明顯偏弱(小圓只包含了橢圓的一部分),大量噪聲仍未被濾除;但如果采用σ=max(σ(j)),則濾波偏強,部分信號被濾除,導致失真.一種可行的處理方式如下：

其中截斷余項的標準差σw′計算方法為

圖3 高斯白噪聲和具有空間、尺度相關噪聲的示意圖Fig.3.Graphical representation of white and correlated noises.

3 濾波試驗

3.1 試驗設置

將赤道緯圈展開為一維等距的n=401個格點,構造一個簡單的B矩陣,其方差V?隨空間緩慢變化,并采用均一、各向同性的高斯函數作為相關函數：

其中i為緯圈上的格點序數,r是格點間距,相關長度尺度Lεb設為300 km.集合樣本數N取50個.為能真實模擬采樣噪聲的特征,這里采用與實際系統相同的隨機方法估計背景誤差標準偏差.首先隨機生成50個服從N(0,I)的隨機矢量αk,k=1,···,N,其次將B1/2與每個矢量αk相乘得到樣本此時滿足N(0,B1/2)分布,最后計算出的方差V,即為背景誤差方差的集合估計值,方差計算方法如下：

因Coif5具有良好的正交和雙正交特性,且在時域和頻域都具有良好的緊支撐和消失矩,試驗中選用Matlab自帶的Coif5正交小波模塊,實現信號分解和重構.

3.2 試驗結果分析

圖4給出了預定義的背景誤差方差真值,即B的對角元素,以及50個樣本的方差估計值.在大部分區域,預定義的方差真值變化較為平緩,而在第150,250個格點附近出現了陡峭變化.這種變化可表征風暴、深對流和臺風等劇烈天氣事件的背景誤差.從圖4可以看出,盡管集合平均能清晰反映出真實方差的大尺度特征,但是存在以信號為中心上下波動的小尺度采樣噪聲.在方差偏大的格點附近,采樣噪聲也明顯偏大.

采用上述方法對集合估計值進行去噪,得到圖5中的結果.對比可以看出,譜濾波(綠線)結果最為平滑,但兩個特征信號同時也被大幅度削弱.其原因在于譜濾波本質上屬于一種空間加權平均濾波,平滑范圍由噪聲和信號的特征尺度決定.權重系數是一個全局量,不能隨空間位置變化,因此譜濾波無法刻畫出信號的局地特征.當特征信號的尺度與噪聲相當時,信號被當作噪聲處理,導致信號失真.未經修正的小波去噪結果,如圖5藍線所示,與真值之間的RMSE(0.77)小于譜方法(1.44),但是去噪后信號波動最為劇烈.從某種意義上說,這是去噪強度偏弱的一種表象,即迭代方法得到的閾值TA=1.42偏小.主要原因在于尺度j=2,3上的噪聲標準偏差σj=2=6.33,σj=3=2.95均大于σw=0.41,部分噪聲的小波系數并沒有被置0.采用等式(6)對閾值進行修正得到TS=2.65,圖5中紅線即為改進后的去噪結果.可以看出,去噪后信號的平滑性有了較大的提高,RMSE也由原來的0.77減至0.65,接近最優濾波的0.57.最優濾波為圖5中點劃線,是通過手動調整閾值使RMSE達到最小得到的.最優濾波代表了這種去噪方法的最佳效果.圖6給出了不同閾值對應的RMSE值.

圖4 背景誤差方差真值(粗線),集合估計值(細線)以及噪聲(點線)隨格點的變化Fig.4.True(bold solid)and EDA-based estimated(thin solid)variances and corresponding sampling noise(dashed).

圖5 譜濾波(綠線)、小波濾波(藍線)、改進后的小波濾波(紅線)以及最優閾值小波濾波(點短線)結果的比較,黑粗線為真值Fig.5.Comparison of filtered variances through spectral filtered(green),primal wavelet filtered(blue),modified wavelet filtered(red)and optimal filtered(dotted line),the true variance is denoted by bold solid line.

圖6 RMSE隨截斷波數的變化TP=3.7為最優閾值,濾波后的RMSE最小為0.57值;TS=2.65為修正后的閾值,對應RMSE為0.65.Fig.6.RMSEs corresponding to different thresholds,TP=3.7 is the optimal threshold leading to the minimal RMSE value 0.57;TS=2.65 is obtained by equation(6),its filtered RMSE=0.65.

譜濾波方法也可以通過調整截斷波數達到最優化,結果如圖7所示.可以看出,調整截斷波數能夠減少譜濾波的誤差,但是RMSE值均在1.25以上,總大于小波去噪結果.這也進一步說明在背景誤差方差濾波中,小波去噪方法要整體優于譜方法.需要指出的是,背景誤差方差的集合估計值在進入到同化系統之前還需要進行降分辨率處理,如ECMWF的業務化集合資料同化系統,需要將T399分辨率的方差估計值進行譜截斷,得到T65分辨率的方差.這種譜截斷處理,會消除小波濾波結果中存在的高頻小尺度振蕩信息,而保留大尺度特性,使濾波結果更加平滑.

圖7 譜濾波RMSE隨截斷波數的變化Fig.7.The variation of RMSE with truncation wave number.

圖8 小波系數的概率密度分布 紅、藍、綠線分別對應噪聲,方差的集合估計值和絕對值小于TA的噪聲小波系數Fig.8.Probability density distribution of the wavelet coefficients,the red,blue and green line corresponding to noise,EDA-based estimate variance and filtered noise respectively whose coefficients magnitude are small than thresholdTA.

為進一步測試該方法的魯棒性,圖9對比了不同信噪比情況下,閾值修正前后的濾波效果.方差集合估計值的精度(黑線)正比于集合成員個數的開方.當集合成員數增加,背景誤差方差估計值中的信噪比也將增加.但是這種收斂性非常緩慢,通過增加樣本個數來提高估計值精度的代價是十分昂貴的,這也凸顯了去噪方法的重要性.可以看出,在不同信噪比的條件下,閾值修正后都能顯著減少RMSE.當集合成員數N取40和100時,改進幅度已經接近或達到極限值(綠色星形線).由這10組試驗統計得出,改進濾波方法后RMSE由原來的1.11減少為0.92,相對于集合估計值X的RMSE(1.40),濾波效率提升了13.28%.

圖9 原集合方差估計值(黑色菱形線)和改進前后濾波RMSE隨集合成員的變化,藍色方框線為改進前的小波濾波結果,紅色圓圈線為改進后的小波濾波結果,綠色星形線為最優濾波結果Fig.9.EDA-based variance estimates(black with diamond),filtered results processed by primal wavelet filter(blue with square)and modified wavelet filter(red with cycle)and the optimal filter(green with star).

3.3 在實際系統中的初步應用

YH4DVAR是我國自主研發的全球四維變分資料同化系統,能夠為全球中期和區域短期數值天氣預報模式提供高質量的初值場[1,2].文獻[15]已經基于YH4DVAR初步構建了集合資料同化系統,提高了預報效果[25].以下采用與文獻[15]中相同的試驗平臺和設置,選擇2013年8月2日0900 UTC(國際時)第91個模式層上的渦度的10個樣本估計值作為去噪對象.此時2013年第九號臺風“飛燕”位于海南省文昌市的東南部,臺風渦旋導致渦度背景誤差方差出現局部最大值.由于實際系統缺乏真值,圖10僅給出了原集合估計值、譜濾波和改進前后小波濾波的結果.可以看出譜濾波和小波閾值方法都能有效濾除集合估計值中的小尺度噪聲,但是與30個樣本估計值(這里未給出結果)相比,譜濾波嚴重弱化了臺風中心渦度方差的極值(由原來的8.31×10-5變為6.20×10-5),且中心位置出現小幅度的漂移,效果低于小波閾值方法,與一維理想試驗的結論相同.小波閾值1.87×10-5經修正后變為2.62×10-5,去噪后臺風中心極值分別為7.65×10-5和7.30×10-5,兩者整體結構相似.閾值修正后濾除了部分小尺度信息,形狀更加平滑,去噪強度有所增強.由于缺乏真值和可信的參考值,該方法的有效性及對系統的影響還有待進一步分析,如修正閾值后對預報和分析場的貢獻.

圖10 第91個模式層上渦度標準偏差,時間對應2013年8月2日0900 UTC (a),(b),(c),(d)分別對應10個集合樣本的估計值、譜濾波結果、閾值未修正時小波閾值去噪結果和閾值修正后的去噪結果Fig.10.Standard deviations of vorticity at ML=91,corresponding to 2 August 2013 at 0900 UTC(a).Unit：10-5s-1.Filtered result based on 10 members raw with：(b)spectral filter,(c)raw wavelet filter and(d)modified wavelet filter.

4 結 論

集合資料同化中方差濾波能減少因采樣導致的隨機噪聲,提高背景誤差方差集合估計值的精度.本文針對譜濾波方法不能表征信號在格點空間局地變化的局限,引入了同時具有譜和空間局地化能力的小波閾值方法,顯著改善了方差的濾波效果.尤其當方差信號局地變化較大時,小波濾波方法的優勢更加明顯.當隨機噪聲為均勻的高斯白噪聲時,通過迭代方法計算得到的小波閾值能夠消除大部分噪聲的小波系數.但是,集合資料同化中隨機噪聲具有尺度相關、非均勻、非高斯等特性,個別尺度上的噪聲能級過高,導致迭代得到的小波閾值偏小.對此,本文根據截斷余項的分布特征,對閾值的計算方法進行了修正和改進.一維理想試驗結果表明,改進后的小波閾值能進一步減少噪聲比重,提升背景誤差方差估計值的精度.10組試驗的統計結果表明,采用改進的小波閾值方法,能使背景誤差方差的RMSE值減少13.28%,同時該方法具有很好的魯棒性.

將該方法應用在了實際的業務系統中,以臺風中心的渦度背景誤差方差的集合估計值作為試驗對象,對比了三種去噪方法的效果.初步結果表明,小波閾值方法的去噪效果優于譜方法.雖然閾值修正前后的去噪結果相似,但修正后能減少更多的小尺度信息.后續將進一步評估研究該方法對系統整體性能的影響,以及研究基于球面小波框架的去噪方法.

[1]Zhang W M,Cao X Q,Song J Q 2012Acta Phys.Sin.61 249202(in Chinese)[張衛民,曹小群,宋君強 2012物理學報61 249202]

[2]Wang S C,LI Y,Zhang W M,Zhao J,Cao X Q 2011Acta Phys.Sin.60 099203(in Chinese)[王舒暢,李毅,張衛民,趙軍,曹小群2011物理學報60 099203]

[3]Laroche S,Gauthier G 1998Tellus Ser.A50 557

[4]Derber J,Bouttier F 1999Tellus Ser.A51 195

[5]Buizza R,Houtekamer P L,Pellerin G,Toth Z,Zhu Y J,Wei M Z 2005Mon.Weather Rev.133 1076

[6]Houtekamer P L,Mitchell H L 1998Mon.Weather Rev.126 3

[7]Evensen G 1994J.Geophys.Res.99 10143

[8]Bonavita M,Raynaud L,Isaksen L 2011Q.J.R.Meteorol.Soc.137 423

[9]Berre L,Varella H,Desroziers G 2015Q.J.R.Meteorol.Soc.141 2803

[10]Raynaud L,Berre L,Desroziers G 2009Q.J.R.Meteorol.Soc.135 1177

[11]Pereira M B,Berre L 2006Mon.Weather Rev.134 2466

[12]Raynaud L,Berre L,Desroziers G 2008Q.J.R.Meteorol.Soc.134 1003

[13]WienerN 1949Extrapolation,Interpolation,and Smoothing of Stationary Time Series(Cambridge：Massachusetts Institute of Technology)pp86-90

[14]Bonavita M,Isaksen L,Hólm E 2012Q.J.R.Meteorol.Soc.138 1540

[15]Liu B N,Zhang W M,Cao X Q,Zhao Y L,Huang Q B 2015China J.Geophys.58 1526(in Chinese)[劉柏年,張衛民,曹小群,趙延來,皇群博,羅雨 2015地球物理學報58 1526]

[16]Donoho D L,Johnstone J M 1994Biometrical81 425

[17]Parrish D F,Derber J C 1992Mon.Weather Rev.120 1747

[18]Fisher M 2003Proceedings ECMWF Seminar on“Recent Developments in Data Assimilation for Atmosphere and Ocean”Reading,September 8-12,2003 p45

[19]Isaksen L,Fisher M,Berner J 2006ECMWF Tech.Memo.492

[20]Moore S,Wood S,Davies P 1998Annals of Statistics26 1

[21]Daley R 1993Atmospheric Data Analysis1993(Cambridge：Cambridge University Press)pp46-50

[22]Azzalini A,Farge M,Schneider K 2005Appl.Comput.Harmon.Anal.18 177

[23]Yen R N V,Farge M,Schneider K 2012Physica D Nonlinear Phenomena241 186

[24]Pannekoucke O,Raynaud L,Farge M 2014Q.J.R.Meteorol.Soc.140 316

[25]Zhang W M,Liu B N,Cao X Q,Zhao Y L,Zhu M B,Zhao W J 2016Acta Meteorol Sin.74 410(in Chinese)[張衛民,劉柏年,曹小群,趙延來,朱孟斌,趙文靜2016氣象學報74 410]

PACS:05.40.-a,05.40.Ca,02.60.-x,43.60.Hj DOI:10.7498/aps.66.020505

Invesitgation and experiments of wavelet thresholding in ensemble-based background error variance?

Liu Bai-Nian1)2)?Huang Qun-Bo1)2)Zhang Wei-Min1)Ren Kai-Jun1)Cao Xiao-Qun1)Zhao Jun1)

1)(Academy of Ocean Science and Engineering,National University of Defense Technology,Changsha 410073,China)
2)(College of Computer,National University of Defense Technology,Changsha 410073,China)

26 March 2016;revised manuscript

20 October 2016)

A large amount of sampling noise which exists in the ensemble-based background error variance need be reduced effectively before being applied to operational data assimilation system.Unlike the typical Gaussian white noise,the sampling noise is scaled and space-dependent,thus making its energy level on some scales much larger than the average.Although previous denoising methods such as spectral filtering or wavelet thresholding have been successfully used for denoising Gaussian white noise,they are no longer applicable for dealing with this kind of sampling noise.One can use a different threshold for each scale,but it will bring a big error especially on larger scales.Another modified method is to use a global multiplicative factor,α,to adjust the filtering strength based on the optimization of trade-o ffbetween removal of the noise and averaging of the useful signal.However,tuningαis not so easy,especially in real operational numerical weather prediction context.It motivates us to develop a new nearly cost-free filter whose threshold can be automatically calculated.

According to the characteristics of sampling noise in background error variance,a heterogeneous filtering method similar to wavelet threshold technology is employed.The threshold,TA,determined by iterative algorithm is used to estimate the truncated remainder whose norm is smaller thanTA.The standard deviation of truncated remainder term is regard as first guess of sampling noise.Non-Guassian term of sampling noise,whose coefficient modulus is aboveTA,is regarded as a small probability event.In order to incorporate such a coefficient into the domain of[-T,T],a semi-empirical formula is used to calculate and approach the ideal threshold.

Investigations are first conducted in a one-dimensional(1D)framework:several methods such as spectral filter,primal wavelet filter,optimal wavelet filter,and proposed filter are used to denoise the scale and space-dependent sampling noise in variance estimations.Their validity and performance are compared and examined with different ensemble sizes.Results show that the proposed method can efficiently filter out a large amount of sampling noise efficiently and improve the estimation accuracy of background error variance.Compared with previous filters,the modified threshold can help to reduce some residual noise arising from the scales where the noise level is above the average level,and the filtering performance of the new method is improved by almost 13.28%.Application to the real ensemble four-dimensional variational data assimilation system is then considered.Results show that the wavelet threshold method outperforms the spectral method.Modified threshold can enhance denosing without deforming the signal of interest.

A new nearly cost-free filter is proposed to reduce the scale and space-dependent sampling noise in ensemble-based background error variance.It is able to remove most of the sampling noises,while extracting the signal of interest.Compared with those of primal wavelet filter and spectral filter,the performance and efficiency of proposed method are improved in 1D framework and real data assimilation system experiments.Further work should focus on the sphere wavelets,which is appropriate for analysing and reconstructing the signals on the sphere in global spectral models.

wavelet,ensemble data assimilation,background error variance,filter

:05.40.-a,05.40.Ca,02.60.-x,43.60.Hj

10.7498/aps.66.020505

?國家自然科學基金(批準號：41375113,41475094,41305101,41605070)資助的課題.

?通信作者.E-mail：bnliu@nudt.edu.cn

*Project supported by the National Natural Science Foundation of China(Grant Nos.41375113,41475094,41305101,41605070).

?Corresponding author.E-mail：bnliu@nudt.edu.cn