集合背景误差方差中小波阈值去噪方法研究及试验∗

刘柏年皇群博 张卫民 任开军 曹小群 赵军

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