经验模态分解滤波法用于时变重力场去噪研究

艾尚校， 肖 云

（1.长安大学地质工程与测绘学院，西安 710054； 2.地理信息工程国家重点实验室，西安 710054；3.西安测绘研究所，西安 710054）

GRACE（Gravity Recovery and Climate Experiment）重力卫星获取了大量的全球尺度、高精度、高时空分辨率的地球重力场信息，为认识和探究地球物质质量的空间分布及变化规律提供了新方法，特别是在监测和分析陆地水储量变化的研究领域取得了大量、优异的研究成果［1-7］.

受卫星观测模式、载荷仪器误差等因素的影响，GRACE时变重力场的高阶球谐位系数误差较大，在空域中表现为严重的南北条带噪声，不利于陆地水储量监测分析研究.因此，时变重力场滤波算法研究成为一个热门研究领域，很多相关研究机构和科学家投入大量精力研究时变重力场滤波方法，产出大量研究成果，实现有效削弱条带噪声目的，获得可靠、高精度重力场时变信号，用于支撑陆地水储量变化研究分析.

目前研究较成熟的时变重力场滤波方法可以分为两类［8］：第一类方法是通过引入平滑核函数以降低高阶次项位系数权重的空间平滑滤波方法，Wahr等［9］最早利用高斯平均核函数降低高阶项位系数的权重，建立了高斯平滑滤波.随后研究发现，高阶项位系数中高次项位系数的误差大于低次项位系数的误差.于是Han 等［10］提出了与位系数阶项和次项均相关的各向异性滤波，有效地提高了空间滤波方法的滤波能力.Zhang等［11］提出了Fan滤波，即对位系数阶项和次项同时进行高斯平滑滤波处理，是一种更为简单的各向异性滤波.这类空间平滑滤波方法在去除噪声的同时会损失部分真实地球物理信号，选取的平滑半径参数越大，损失的真实信号越多，空间分辨率也越低.第二类方法是去相关滤波方法，Swenson和Wahr［12］发现球谐系数中同次的奇（偶）项阶存在相关误差，并设计滑动固定窗多项式拟合消除该误差，即去相关滤波.Chambers［13］和Chen等［14］提出了对大于m次的位系数进行n次多项式拟合的PnMm方法，Duan等［15］根据球谐系数中误差分布特点提出了滑动可变窗去相关滤波.这类去相关滤波方法在低纬度地区的去噪能力较弱.在实际时变重力场滤波应用中，通常将空间滤波方法与去相关滤波方法组合使用，可获得更好的滤除噪声效果.尽管国内外大量的相关研究取得了很好的成果，但是在时变重力场滤波方法的去噪能力、保留真实信号能力等方面还有潜力可以挖掘.基于此，本文引入了经验模态分解法，应用于时变重力场滤波处理，以期获得更好的滤波效果，提高时变重力场数据的精度.

本文首先引入经验模态分解方法，构建了适用时变重力场滤波的函数模型，分析了技术可行性，进一步与其他成熟时变重力场滤波方法进行比较，评估本文提出的方法去噪能力、保留真实信号能力、反演准确度等能力指标，最后通过与陆面水文模型GLDAS（Global Land Data assimilation System）进行对比，进一步验证该方法提取陆地水储量的有效性和可靠性.

1 原理模型

1.1 经验模态分解原理

经验模态分解（Empirical Mode Decomposition，EMD）是一种自适应信号分解方法，将输入信号分解成有限个不同特征尺度的固有模态函数（Intrinsic Mode Function，IMF）分量和一个残余分量［16］.每一个IMF分量都要满足两个条件：一是IMF中极值点的个数和过零点的个数相等或相差最大不超过1；二是由极大值和极小值构造的上包络线和下包络线的均值为0［16］.给定一个待分解的信号为x(t)，其经验模态分解过程如下.

步骤1：构造包络线.寻找待分解信号x(t)中的极大值点和极小值点，并通过极大值点利用三次样条插值函数构造上包络线，用相同方法通过极小值点构造下包络线.

步骤2：信号分解.将信号x(t)减去上下包络线的均值xˉ(t)，得到新信号y(t)，即

步骤3：条件判断.若y(t)满足IMF 的两个条件，则认定y(t)为首个IMF，并记为s1；若不满足，将y(t)作为新的待分解信号，重复步骤1～3，直至得到第一个IMF分量，记为s1.

步骤4：迭代筛选.上述过程可视为IMF分量的筛选过程，此时将待分解信号x(t)减去y(t)，得到新的待分解信号g1(t)，即

对g1(t)重复IMF分量的筛选过程，得到第二个IMF分量，记为s2.如此反复进行，直至gn(t)为单调信号时记为r(t)，经验模态分解结束.则待分解信号x(t)的函数表达式为：

式中：si表示信号x(t)分解出的第i个IMF分量，频率由高到低排列，r(t)为分解残余项.

1.2 构建适用时变重力场的滤波方法

未经滤波处理的时变重力场数据在空域中存在明显的南北条带噪声.其任意纬度带剖面数据，可以看作是随经度λ变化的纬度带信号x(λ) .对于纬度带信号，噪声主要集中在相对高的频率波段，而真实的重力场信号集中在相对低频率波段［12］.由经验模态分解原理可知，IMF分量按照频率由高到低排序，则时变重力场高频噪声集中在排序靠前的模态分量中，而排序靠后的分量是低频真实重力场信号主导的模态分量.因此，只需确定区分高频噪声和低频信号的分界模态分量sk，然后重构信号分量就可以实现纬度带信号的滤波处理.

本文使用模态相关分选准则［17］作为选取分界模态分量的依据，首先计算纬度带的原始信号x(λ)与每个模态分量之间的互相关系数，计算式如下：

式中：N为经度采样点的个数，对于1°×1°的时变重力场格网数据而言，N=360；λ表示经度；x(λ)表示以经度为参量的原始信号；xˉ为原始信号的平均值；si表示第i个模态函数分量；sˉi为第i个模态函数分量的平均值；R为互相关系数.

上述所得的互相关系数中，第一个局部极小值点对应的模态分量即为分界模态分量sk，则对于sk之后的模态分量即为纬度带中的低频信号.重构低频信号分量，并引入包含着部分细节信号信息sk和sk-1分量，就得到滤波后的纬度带信号x̂(λ)：

2 时变重力场经验模态分解试验

以2004 年11 月的时变重力场数据为例，分析时变重力场经验模态分解的可行性和有效性.在时变重力场信号中选取0°N纬度带、30°N纬度带、40°N 纬度带等几个具有代表性的时变信号x(λ)，其信号特征见图1.图1（a）～图1（c）分别表示0°N 纬度带、30°N 纬度带、40°N 纬度带的信号，其中左侧图为纬度带信号随经度的变化，右侧为相应的归一化功率谱（下同）.由图1 可见，不同纬度带的噪声呈现出不同特征，纬度越低，噪声的频率越高.

图1 不同纬度带信号及其归一化功率谱Fig.1 The signals and normalized power spectra of different latitudes

以40°N 纬度带为例分析上述方法的可行性和有效性.将信号x(λ)进行经验模态分解，得到的结果如图2所示.

图2 经验模态分解结果Fig.2 Empirical mode decomposition results

从图2 中可以看出，信号x(λ)分解后得到了9 个固有模态函数分量和1 个残余项.根据式（4）计算x(λ)与各分量之间的互相关系数，如图3所示.

图3 互相关系数曲线Fig.3 The cross correlation coefficient curve

由图3 可知，固有模态函数分量s6为分界模态，则选择s5及其之后的模态分量与残余分量进行重构，就得到滤波后信号.以本文方法对不同纬度带信号进行滤波，滤波结果如图4所示.图4（a）～图4（c）分别表示经过经验模态分解滤波后的0°N纬度带、30°N纬度带、40°N纬度带的信号.对比图1和图4可以发现，图4 左侧图呈现出低频率特性，右侧图的信号也主要集中在低频带，因此，本文提出的经验模态分解滤波对于不同纬度带信号的高频噪声均有明显的滤除效果.

图4 不同纬度带信号滤波后的特征Fig.4 The characteristics of denoised signals in different latitudes

3 应用试验

利用本文方法对实际的时变重力场数据进行滤波处理，分析其滤波效果.以2007年10月的数据为例，滤波前后结果如图5所示，其中图5（a）为未经滤波处理的时变重力场数据，图5（b）为EMD滤波后的结果.由图5可见EMD滤波有效滤除了时变重力场的噪声，滤波后可以明显辨识出陆地上亚马逊流域、奥里诺科河流域、恒河流域等地区的重力场时变信号.

图5 EMD滤波前后结果对比Fig.5 Comparison of results before and after EMD filtering

为进一步分析本文提出的滤波方法在时变重力场滤波的有效性，选取P3M6多项式拟合（P3M6）和滑动可变窗多项式拟合（Duan）两种广泛应用的去相关滤波作为对比方法，在相同的高斯平滑半径下对比分析各方法的滤波效果.本文采用两个评价指标比较各滤波方法效果，一是选取滤波后海洋上残余信号的均方根（Root Mean Square，RMS）作为评价滤波方法去噪能力的指标［1］，二是选择滤波后陆地信号和海洋残余信号的RMS比值（信噪比）作为评价滤波方法保留真实信号能力的指标［18］.计算三种方法滤波结果的均方根指标和信噪比指标，结果如表1所示.

由表1可知，在相同的平滑半径下，本文方法滤波后的海洋残余信号均方根值最小，且信噪比值最高，说明该方法的去噪能力和保留真实信号的能力均优于其他两种方法.另外，由表可见平滑半径是影响滤波效果的重要因素，当平滑半径为300 km时，三种方法均取到了最高信噪比值，分别为2.90（P3M6）、3.01（Duan）和3.14（EMD）.选择最优信噪比条件，即平滑半径设置为300 km，利用三种滤波方法分别对图5（a）数据进行处理，结果如图6所示.

图6（a）～图6（c）分别为当平滑半径为300 km时P3M6、Duan和EMD滤波的结果.对比图5（a）和图6，三种滤波方法均较好地滤除了南北条带噪声，可以清晰地辨识时变重力场信息.其中本文滤波方法结果中残余噪声低于其他两种方法，说明本文滤波方法的去噪能力更强，与表1结论相吻合.为评价滤波结果之间的一致性，计算三种方法滤波结果的两两差值，并统计差值的均方根，结果分别为1.2 cm（P3M6和Duan作差）、1.3 cm（EMD和P3M6作差）和0.92 cm（EMD和Duan作差），证明本文滤波方法具有较高的可靠性.此外，相比于其他两种方法，本文滤波方法在格陵兰岛西部和南美洲中部的信号泄漏误差最小，反演结果的准确度更高.

为了进一步验证本文滤波方法的可靠性，选取亚马逊流域、密西西比河流域和长江流域作为研究区域，比较三个区域重力卫星反演的时变重力场滤波结果与GLDAS 水文模型估计的水资源变化结果一致性.利用本文提出滤波方法及上述两种滤波方法反演选定区域陆地水储量变化，时间跨度为2003年1月至2014 年12 月，同时利用GLDAS 水文模型反演相同时段内研究区域水储量的变化，为保证空间分辨率一致性，对GLDAS 数据进行300 km 高斯滤波.各流域水储量变化的时间序列如图7 所示.通过最小二乘法拟合估算各滤波方法反演结果的周年振幅和周年相位，并计算各滤波方法与水文模型反演结果之间的相关系数，结果如表2所示.

由图7 和表2 可知，本文提出的滤波方法与其他两个传统成熟算法相比较，结果之间具有较高的一致性，时间序列曲线几乎重合.在周年振幅和相位方面，本文滤波方法与Duan滤波结果基本一致，充分说明本文方法的可靠性.进一步与水文数据GLDAS反演的水储量信号比较，发现本文方法结果与GLDAS 反演的水储量相关性最高，其相关系数均大于0.86；但是重力卫星反演的水储量信号振幅明显大于水文模型反演的水储量信号振幅，究其原因是水文模型反演的水储量仅包含浅层地表水的变化，而重力卫星反演的水储量囊括了土壤水、地下水、湖泊、湿地等更综合的水文变化信号.

图7 研究区域水储量变化Fig.7 Water storage changes in study areas

表2 各滤波方法结果的周年振幅、相位和与GLDAS的相关系数Tab.2 Anniversary amplitudes，phases and correlation coefficients with GLDAS for different filter results

4 结语

本文引入了一种新的时变重力场空域滤波方法——经验模态分解滤波法，介绍了经验模态分解基本原理，构建了匹配时变重力场滤波的滤波函数，分析了经验模态分解滤除时变重力场条带噪声的有效性，得到主要结论如下：

1）本文提出时变重力场经验模态分解滤波法是一种有效滤波方法，与现有成熟的滤波方法比较，有更优的去噪能力.试验分析表明，经验模态分解滤波方法的信噪比达到3.14，而P3M6滤波和Duan滤波的信噪比分别为2.90和3.01，因此本文滤波方法具有较好的保留真实信号能力.

2）本文提出方法与其他两种成熟方法在反演亚马逊流域、密西西比河流域和长江流域等地区的水储量方面具有很好一致性，与Duan滤波结果的周年振幅和相位基本相同.进一步与水文数据GLDAS反演的水储量信号比较，反演水储量相关性最高，其相关系数均大于0.86.