劉柏年皇群博 張衛(wèi)民 任開軍 曹小群 趙軍

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

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

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

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

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

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

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

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

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

1 引 言

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

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

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

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

2 去噪原理

2.1 譜濾波原理

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

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

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

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

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

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

2.2 小波閾值去噪原理

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

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

3)進行小波閾值截斷

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

圖2 噪聲方差和閾值迭代算法示意圖 粗實線圓代表T0,虛線代表T1,箭頭表示小波系數(shù)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.

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

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

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

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

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

其中截斷余項的標(biāo)準(zhǔn)差σw′計算方法為

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

3 濾波試驗

3.1 試驗設(shè)置

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

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

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

3.2 試驗結(jié)果分析

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

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

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

圖5 譜濾波(綠線)、小波濾波(藍(lán)線)、改進后的小波濾波(紅線)以及最優(yōu)閾值小波濾波(點短線)結(jié)果的比較,黑粗線為真值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隨截斷波數(shù)的變化TP=3.7為最優(yōu)閾值,濾波后的RMSE最小為0.57值;TS=2.65為修正后的閾值,對應(yīng)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.

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

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

圖8 小波系數(shù)的概率密度分布 紅、藍(lán)、綠線分別對應(yīng)噪聲,方差的集合估計值和絕對值小于TA的噪聲小波系數(shù)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對比了不同信噪比情況下,閾值修正前后的濾波效果.方差集合估計值的精度(黑線)正比于集合成員個數(shù)的開方.當(dāng)集合成員數(shù)增加,背景誤差方差估計值中的信噪比也將增加.但是這種收斂性非常緩慢,通過增加樣本個數(shù)來提高估計值精度的代價是十分昂貴的,這也凸顯了去噪方法的重要性.可以看出,在不同信噪比的條件下,閾值修正后都能顯著減少RMSE.當(dāng)集合成員數(shù)N取40和100時,改進幅度已經(jīng)接近或達到極限值(綠色星形線).由這10組試驗統(tǒng)計得出,改進濾波方法后RMSE由原來的1.11減少為0.92,相對于集合估計值X的RMSE(1.40),濾波效率提升了13.28%.

圖9 原集合方差估計值(黑色菱形線)和改進前后濾波RMSE隨集合成員的變化,藍(lán)色方框線為改進前的小波濾波結(jié)果,紅色圓圈線為改進后的小波濾波結(jié)果,綠色星形線為最優(yōu)濾波結(jié)果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 在實際系統(tǒng)中的初步應(yīng)用

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

圖10 第91個模式層上渦度標(biāo)準(zhǔn)偏差,時間對應(yīng)2013年8月2日0900 UTC (a),(b),(c),(d)分別對應(yīng)10個集合樣本的估計值、譜濾波結(jié)果、閾值未修正時小波閾值去噪結(jié)果和閾值修正后的去噪結(jié)果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 結(jié) 論

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

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

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

[2]Wang S C,LI Y,Zhang W M,Zhao J,Cao X Q 2011Acta Phys.Sin.60 099203(in Chinese)[王舒暢,李毅,張衛(wèi)民,趙軍,曹小群2011物理學(xué)報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)[劉柏年,張衛(wèi)民,曹小群,趙延來,皇群博,羅雨 2015地球物理學(xué)報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)[張衛(wèi)民,劉柏年,曹小群,趙延來,朱孟斌,趙文靜2016氣象學(xué)報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

?國家自然科學(xué)基金(批準(zhǔn)號：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