基于混合滤波架构的地震信号非线性估计与去噪方法研究

Research on Nonlinear Estimation and Denoising Methods for Seismic Signals Based on a Hybrid Filtering Architecture

摘要: 地震勘探采集的数据常受到多种噪声干扰，这些噪声会掩盖有效的地震反射信号，影响后续地质构造的识别与解释。为提高信号质量，本文提出一种基于期望最大化(EM)算法、粒子滤波(PF)与卡尔曼滤波(KF)的混合智能去噪模型(EM-KF-PF)。该模型通过EM算法增强粒子多样性，利用PF处理非线性与非高斯噪声，并引入KF以提升状态估计的收敛效率，从而在信噪比(SNR)与估计精度之间实现有效平衡。数值实验表明，相较于传统方法，本文算法在提升SNR、降低均方根误差(MSE)及状态轨迹拟合方面均表现更优，展现出良好的非线性适应性与多尺度噪声抑制能力。该模型适用于地震信号处理场景，但在参数自适应性、高维扩展性及非高斯噪声建模等方面仍有改进空间。未来研究方向包括GPU并行加速、在线学习机制及物理约束嵌入等，以提升算法的工程适用性。

Abstract: Seismic exploration data are often contaminated by various types of noise, which can obscure effective seismic reflection signals and impact the identification and interpretation of subsequent geological structures. To improve signal quality, this paper proposes a hybrid intelligent denoising model (EM-KF-PF) based on the Expectation-Maximization (EM) algorithm, Particle Filtering (PF), and Kalman Filtering (KF). This model enhances particle diversity through the EM algorithm, utilizes PF to handle nonlinear and non-Gaussian noise, and incorporates KF to improve the convergence efficiency of state estimation, thereby achieving an effective balance between signal-to-noise ratio (SNR) and estimation accuracy. Numerical experiments demonstrate that, compared to traditional methods, the proposed algorithm performs better in terms of improving SNR, reducing mean square error (MSE), and fitting state trajectories, exhibiting strong nonlinear adaptability and multi-scale noise suppression capabilities. This model is suitable for seismic signal processing scenarios, but there is still room for improvement in areas such as parameter self-adaptability, high-dimensional scalability, and non-Gaussian noise modeling. Future research directions include GPU parallel acceleration, online learning mechanisms, and the incorporation of physical constraints to enhance the engineering applicability of the algorithm.

1. 引言

地震数据作为地球物理勘探的核心数据源，其处理质量直接影响对地下介质结构的推断与解释，在能源资源勘探、地质灾害预警及地球深部研究中具有关键作用[1]。然而，实际采集的地震信号往往受到多种噪声干扰，导致有效波信号被污染，信噪比降低，进而影响成像与反演结果的准确性。因此，发展高效、自适应的地震信号去噪方法具有重要的理论意义与应用价值。

随着高密度采集与微震监测技术的进步，传统基于频域滤波或时频分析的方法在处理非平稳、频谱重叠的噪声时面临局限。近年来，稀疏表示理论与数据驱动的深度学习方法被引入地震数据处理中[2]-[4]，通过构建信号与噪声在特征空间中的差异性模型，实现噪声的自适应压制，为高精度勘探与地球物理建模提供了数据基础。尽管如此，噪声与有效信号的非线性耦合机制、算法的计算效率与泛化能力仍是当前研究的核心挑战。地震信号去噪方法的研究历程反映了信号处理技术在地球物理领域的持续演进。20世纪80至90年代，以Donoho的小波阈值去噪[5] [6]和Yilmaz的FK滤波[7]为代表的传统线性方法奠定了理论基础。进入21世纪后，非线性方法与稀疏表示技术开始兴起，如2001年Sacchi提出的Cadzow滤波[8]-[10]和2005年Herrmann引入的压缩感知理论[11]。2010年代初期，自适应与混合方法得到发展，包括EMD结合小波阈值[12]-[18]、字典学习[19]-[21]等创新。2016年起，深度学习彻底改变了该领域，CNN、WaveNet和GAN等模型相继应用于地震去噪[22]-[24]。这一演进过程呈现出从线性到非线性、从模型驱动到数据驱动、从单模态到多模态融合的明显趋势，未来将朝着实时轻量化、跨尺度建模和可解释AI等方向发展，持续推动地震信号处理技术的革新。当前地震信号去噪领域面临的核心挑战在于传统方法依赖强假设且参数调优困难，机器学习受限于特征工程和小样本问题，深度学习则面临数据需求大、实时性差和可解释性不足的瓶颈。跨领域难点包括噪声与信号的非线性耦合、实时处理与精度的权衡，以及多通道协同处理的复杂性。未来突破需聚焦物理模型与深度学习的融合、无监督学习技术、边缘计算优化和标准化噪声库建设，通过交叉创新解决算法适应性、计算效率和工程落地问题，推动地震预警和勘探技术的发展。

本文提出一种基于期望最大化(EM)算法、粒子滤波(PF)与卡尔曼滤波(KF)的三级混合滤波架构。该框架以EM算法为核心，通过极大似然估计迭代优化系统参数，有效缓解模型设定误差对状态估计的影响；粒子滤波模块通过重要性采样与重采样策略处理非高斯噪声，并设计峰值保护机制以抑制粒子退化；卡尔曼滤波则在线性高斯假设下提供最优状态估计，构成基础滤波层。通过建立基于后验概率的动态权重融合准则，该架构在非线性、非平稳系统中实现计算效率与估计精度的平衡。数值实验表明，所提方法在非高斯噪声环境下仍保持稳定的实时性能。然而，当前模型在状态转移机制的自适应性方面存在局限，未来可通过嵌入物理先验、构建多尺度状态空间模型及设计在线参数学习策略，进一步提升算法在复杂系统下的泛化能力。本研究为非线性滤波理论与智能优化方法的融合提供了有效途径，对推动现代信号处理算法的理论发展与工程应用具有参考价值。

2. 基本原理

2.1. 理论基础

2.1.1. 空间模型

设动态系统的状态空间模型为：

$\{\begin{array}{l}{x}_t=f_t\left(x_{t-1},{\omega }_t\right)\left(��方程\right)\\ y_t=h_t\left(x_t,v_t\right)\left(��方程\right)\end{array}$

其中，x_t：t时刻的状态向量；y_t：t时刻的观测向量；ω：过程噪声(协方差Q)；v_t：观测噪声(协方差R)。

2.1.2. 递推

初始状态分布：

$p\left(x_0|y_0\right)\equiv p\left(x_0\right)$

预测方程：

$p\left(x_t|y_{1:t-1}\right)=\underset{���移模型}{\underbrace{{\displaystyle \int p\left(x_t|x_{t-1}\right)}}}\cdot \underset{后�}{\underbrace{p\left(x_{t-1}|y_{1:t-1}\right)}}\text{d}{x}_{t-1}$

更新步方程：

$p\left(x_t|y_{1:t}\right)=\frac{\stackrel{似然函�}{\overbrace{p\left(y_t|x_t\right)}}\cdot \stackrel{先�}{\overbrace{p\left(x_t|y_{1:t-1}\right)}}}{\underset{�据}{\underbrace{{\displaystyle \int p\left(y_t|x_t\right)p\left(x_t|y_{1:t-1}\right)\text{d}{x}_t}}}}$

2.2. 卡尔曼滤波更新步

卡尔曼滤波更新步的推导建立在贝叶斯定理和高斯分布性质的基础上，以下是完整的数学推导过程：

系统模型(线性高斯)：

$x_t=F{x}_{t-1}+{\omega }_t,{\omega }_t~\mathcal{N}\left(0,\mathcal{Q}\right)$

$y_t=H{x}_t+v_t,v_t~\mathcal{N}\left(0,ℛ\right)$

其中，F是状态转移矩阵，H是观测矩阵，$\mathcal{Q}$ 是过程噪声协方差，$ℛ$ 是观测噪声协方差。

预测步结果：

${\widehat{x}}_t{}^{-}=F{\widehat{x}}_{t-1}$

$P_t^{-}=F{P}_{t-1}{F}^{\text{T}}+\mathcal{Q}$

其中，P是误差协方差矩阵。

更新步：

$p\left(y_t|x_t\right)=\mathcal{N}\left(H{x}_t,ℛ\right)$

定义联合随机变量：

$z=\left[\begin{array}{c}{x}_t\\ y_t\end{array}\right]$

其均值和协方差为：

${\mu }_z=\left[\begin{array}{c}{\widehat{x}}_t^{-}\\ H{\widehat{x}}_t^{-}\end{array}\right],{\sum }_z=\left[\begin{array}{cc}{P}_t^{-}& P_t^{-}{H}^{\text{T}}\\ H{P}_t^{-}& H{P}_t^{-}{H}^{\text{T}}+ℛ\end{array}\right]$

对于联合高斯分布：

$p\left(x|y\right)=\mathcal{N}\left(\left[\begin{array}{c}{\mu }_x\\ {\mu }_y\end{array}\right],\left[\begin{array}{cc}{\sum }_{xx}& {\sum }_{xy}\\ {\sum }_{yx}& {\sum }_{yy}\end{array}\right]\right)$

条件分布为：

$p\left(x|y\right)=\mathcal{N}\left({\mu }_{x|y},{\sum }_{x|y}\right)$

其中：

${\mu }_{x|y}={\mu }_x+{\sum }_{xy}{\sum }_{yy}^{-1}\left(y-{\mu }_y\right)$

${\sum }_{x|y}={\sum }_{xx}-{\sum }_{xy}{\sum }_{yy}^{-1}{\sum }_{yx}$

应用到我们的问题中：

${\sum }_{yx}=H{P}_t^{-},{\sum }_{yy}=H{P}_t^{-}{H}^{\text{T}}+ℛ$

定义卡尔曼增益：

$K_t=P_t^{-}{H}^{\text{T}}{\left(H{P}_t^{-}{H}^{\text{T}}+ℛ\right)}^{-1}$

所以状态更新方程：

${\widehat{x}}_t={\widehat{x}}_t^{-}+K_t\left(y_t-H{\widehat{x}}_t^{-}\right)$

协方差更新方程：

$P_t=\left(I-K_{t}H\right) P_t^{-}$

2.3. 粒子滤波原理的数学推导

粒子滤波(Particle Filter)是基于蒙特卡洛方法的非线性非高斯系统状态估计技术，其核心思想是用一组带权重的粒子来近似后验概率分布。以下是数学推导过程：

根据贝叶斯定理：

$p\left(x_{0:t}|y_{1:t}\right)=\frac{p\left(y_{1:t}|x_{0:t}\right)p\left(x_{0:t}\right)}{p\left(y_{1:t}\right)}$

由于直接计算该分布非常困难(涉及高维积分)，采用蒙特卡洛近似。

对于任意函数f(x)的期望：

$\mathbb{E}\left[f\left(x\right)\right]={\displaystyle \int f\left(x\right)p\left(x\right)\text{d}x}$

用样本近似：

$\mathbb{E}\left[f\left(x\right)\right]\approx \frac{1}{N}{\displaystyle \sum }_{i=1}^{N}f\left(x^{\left(i\right)}\right),x^{\left(i\right)}~p\left(x\right)$

当无法直接从p(x)采样时，引入建议分布q(x)：

$\mathbb{E}\left[f\left(x\right)\right]={\displaystyle \int f\left(x\right)\frac{p\left(x\right)}{q\left(x\right)}q\left(x\right)\text{d}x}\approx \frac{1}{N}{\displaystyle \sum _{i=1}^{N}f\left(x^{\left(i\right)}\right)\frac{p\left(x^{\left(i\right)}\right)}{q\left(x^{\left(i\right)}\right)}}$

定义重要性权重：

${\omega }^{\left(i\right)}=\frac{p\left(x^{\left(i\right)}\right)}{q\left(x^{\left(i\right)}\right)}$

对于时序问题，建议分布分解为：

$q\left(x_{0:t}|y_{1:t}\right)=q\left(x_t|x_{0:t-1},y_{1:t}\right)q\left(x_{0:t-1}|y_{1:t-1}\right)$

权重递推关系：

${\omega }^{\left(i\right)}\propto \frac{p\left(x_{0:t}^{\left(i\right)}|y_{1:t}\right)}{q\left(x_{0:t}^{\left(i\right)}|y_{1:t}\right)}=\frac{p\left(y_t|x_t^{\left(i\right)}\right)p\left(x_t^{\left(i\right)}|x_{t-1}^{\left(i\right)}\right)p\left(x_{0:t-1}^{\left(i\right)}|y_{1:t-1}\right)}{q\left(x_t^{\left(i\right)}|x_{0:t-1}^{\left(i\right)},y_{1:t-1}\right)q\left(x_{0:t-1}^{\left(i\right)}|y_{1:t-1}\right)} ={\omega }_{t-1}{}^{\left(i\right)}\cdot \frac{p\left(y_t|x_t^{\left(i\right)}\right)p\left(x_t^{\left(i\right)}|x_{t-1}^{\left(i\right)}\right)}{q\left(x_t^{\left(i\right)}|x_{0:t-1}^{\left(i\right)},y_{1:t-1}\right)}$

理论上的最优建议分布：

$q\left(x_t|x_{0:t-1},y_{1:t-1}\right)=p\left(x_t|x_{t-1},y_t\right)$

但通常难以采样，因此常用先验分布：

$q\left(x_t|x_{0:t-1},y_{1:t-1}\right)=p\left(x_t|x_{t-1}\right)$

此时权重更新简化为：

${\omega }_t^{\left(i\right)}={\omega }_{t-1}^{\left(i\right)}\cdot p\left(y_t|x_t^{\left(i\right)}\right)$

为避免数值问题，进行权重归一化：

${\tilde{\omega }}_t^{\left(i\right)}=\frac{{\omega }_t^{\left(i\right)}}{{\displaystyle {\sum }_{j=1}^N{\omega }_t^{\left(j\right)}}}$

解决粒子退化问题，重采样概率：

$P\left(��粒子i\right)={\tilde{\omega }}_t^{\left(i\right)}$

重采样后权重重置为：

${\omega }_t^{\left(i\right)}=\frac{1}{N}$

波峰区域权重调整：

${\omega }_t^{\left(i\right)}=\{\begin{array}{l}\exp \left(-\frac{{\left(y_t-H{x}_t^{\left(i\right)}\right)}^2}{2ℛ\cdot 0.5}\right)波峰�域\\ \exp \left(-\frac{{\left(y_t-H{x}_t^{\left(i\right)}\right)}^2}{2ℛ}\right)其他�域\end{array}$

相当于在波峰区域将观测噪声方差减半($R\to 0.5R$ )，使得：观测更“可信”、对匹配度要求更高、增强波峰特征的跟踪精度。

2.4. 期望最大化算法的推导

EM算法通过迭代优化解决含隐变量的最大似然估计问题。在卡尔曼滤波语境下，系统状态序列被视为隐变量，观测序列为已知数据，待估参数集完整描述了系统动力学特性。EM算法通过交替执行期望步骤(E-step)和最大化步骤(M-step)逐步提升参数估计质量，其数学本质是在无法直接观测完整数据的情况下，通过迭代逼近最大似然估计。

在E步骤中，算法基于当前参数估计计算完整数据对数似然函数的条件期望：

$Q\left(\theta |{\theta }^{\left(k\right)}\right)={\rm E}\left[\log p\left(X,Y|\theta \right)|Y,{\theta }^{\left(k\right)}\right]$

这一步骤的核心的卡尔曼平滑技术的应用，通过平滑算法，我们获得状态的后验分布：

平滑状态估计，即给定全部观测数据条件下状态变量的条件期望：$x_t^s=\text{Ε}\left[x_t|Y,{\theta }^{\left(k\right)}\right]$ ；

平滑误差协方差，表征了平滑状态估计的不确定性：$P_t^s=\text{Cov}\left[x_t|Y,{\theta }^{\left(k\right)}\right]$ ；

相邻状态互协方差，描述了不同时刻状态变量之间的相关性：$P_{t,t-1}^s=\text{Cov}\left[x_t,x_{t-1}|Y,{\theta }^{\left(k\right)}\right]$ 。

这三个量共同构成了后续M步骤中更新模型参数所需的充分统计量(Sufficient Statistics)的基础，是连接E步骤与M步骤的桥梁。

基于这些平滑结果，我们计算五个关键的充分统计量：

$S_{11}={\displaystyle {\sum }_{t=1}^T\text{Ε}\left[x_t{x}_t^{\text{T}}|Y\right]}={\displaystyle {\sum }_{t=1}^T}\left(P_t^s+x_t^s{\left(x_t^s\right)}^{\text{T}}\right)$

$S_{10}={\displaystyle {\sum }_{t=2}^T\text{Ε}\left[x_t{x}_{t-1}^{\text{T}}|Y\right]}={\displaystyle {\sum }_{t=2}^T\left(P_{t,t-1}^s+x_t^s{\left(x_{t-1}^s\right)}^{\text{T}}\right)}$

$S_{00}={\displaystyle {\sum }_{t=1}^{T-1}\text{Ε}\left[x_t{x}_t^{\text{T}}|Y\right]}={\displaystyle {\sum }_{t=1}^{T-1}\left(P_t^s+x_t^s{\left(x_t^s\right)}^{\text{T}}\right)}$

$S_{yz}={\displaystyle {\sum }_{t=1}^T{y}_t\text{Ε}\left[x_t^{\text{T}}|Y\right]}={\displaystyle {\sum }_{t=1}^T{y}_t{\left(x_t^s\right)}^{\text{T}}}$

$S_{yy}={\displaystyle {\sum }_{t=1}^T{y}_t{y}_t^{\text{T}}}$

这些统计量完整封装了在当前参数设定下，系统状态与观测数据之间的统计关联特性。

M步骤：参数更新优化。

在M步骤中，算法通过最大化Q函数来更新参数估计：

${\theta }^{\left(k+1\right)}=\arg {\max }_{\theta }Q\left(\theta |{\theta }^{\left(k\right)}\right)$

通过求解相应的极值问题，我们获得参数的解析更新公式：

状态转移矩阵更新：$A^{\left(k+1\right)}=S_{10}{S}_{00}^{-1}$ 。

这一更新确保了新的状态转移矩阵最优地拟合状态序列的动态演化规律。

过程噪声协方差更新：$Q^{\left(k+1\right)}=\frac{1}{T-1}\left(S_{11}-A^{\left(k+1\right)}{S}_{10}^{\text{T}}\right)$ 。

该公式本质上是状态预测误差的样本协方差，反映了系统模型的不确定性。

观测矩阵更新：$C^{\left(k+1\right)}=S_{yx}{S}_{11}^{-1}$ 。

观测矩阵的更新建立了状态与观测之间的最优线性映射关系。

观测噪声协方差更新：$R^{\left(k+1\right)}=\frac{1}{T}\left(S_{11}-C^{\left(K+1\right)}{S}_{yt}^{\text{T}}-S_{yx}{\left(C^{\left(k+1\right)}\right)}^{\text{T}}+C^{\left(K+1\right)}{S}_{11}{\left(C^{\left(K+1\right)}\right)}^{\text{T}}\right)$ 。

这一更新量化了观测过程中的不确定性。

通过EM算法的迭代优化，系统能够自动从观测数据中学习这些关键参数，实现了完全数据驱动的自适应滤波，无需任何关于信号和噪声统计特性的先验知识。这一特性使该方法在非平稳信号处理和时变系统建模中具有显著优势。

2.5. 基于EM算法的自适应融合框架

本节提出一种基于期望最大化算法的自适应融合滤波方法，通过建立完整的参数学习与融合决策机制，实现卡尔曼滤波与粒子滤波的智能融合。该方法的创新性在于将噪声参数估计与融合权重优化纳入统一的概率框架中。

2.5.1. 融合权重的概率化学习机制

在EM算法的M步骤中，融合权重的更新基于粒子滤波与卡尔曼滤波的相对性能指标。具体而言，权重通过以下公式计算：

$\omega =\frac{{ℒ}_{PF}}{{ℒ}_{PF}+{ℒ}_{KF}}$

其中${ℒ}_{PF}$ 和${ℒ}_{KF}$ 分别表示两个滤波器的边缘似然度，计算公式为：

${ℒ}_{PF}={\displaystyle {\prod }_{t=1}^{T}p\left(y_t|{\widehat{x}}_t^{PF},R\right)}$ , ${ℒ}_{KF}={\displaystyle {\prod }_{t=1}^{T}p\left(y_t|{\widehat{x}}_t^{KF},R\right)}$

为增强数值稳定性，实际计算采用对数域处理：

$\omega =\frac{\exp \left(\log {ℒ}_{PF}-\eta \right)}{\exp \left(\log {ℒ}_{PF}-\eta \right)+\exp \left(\log {ℒ}_{KF}-\eta \right)}$

其中$\eta =\max \left(\log {ℒ}_{PF},\log {ℒ}_{KF}\right)$ 为归一化因子。

2.5.2. 参数反馈与实时调节机制

学习到的系统参数通过闭环机制实时反馈至滤波器运行过程中：过程噪声协方差Q的更新直接影响状态转移模型。在粒子滤波中，它决定了粒子的传播范围：

$x_t^{\left(i\right)}~\mathcal{N}\left(F{x}_{t-1}^{\left(i\right)},Q\right)$

在卡尔曼滤波中，它调整预测步的不确定性：

$P_t^{-}=F{P}_{t-1}{F}^{\text{T}}+Q$

观测噪声协方差R的更新同步影响两个滤波器。在粒子滤波中，它调节权重的似然计算：

${\omega }_t^{\left(i\right)}\propto \mathcal{N}\left(y_t;H{x}_t^{\left(i\right)},R\right)$

在卡尔曼滤波中，它修正卡尔曼增益：

$K_t=P_t^{-}{H}^{\text{T}}{\left(H{P}_t^{-}{H}^{\text{T}}+R\right)}^{-1}$

参数更新采用立即反馈策略：在每个EM周期结束后，新估计的Q和R立即应用于下一时间步的滤波过程，确保系统能够快速适应动态变化的环境。

2.5.3. 创新性分析

本方法的创新性体现在三个方面：首先，提出了基于概率似然的权重学习机制，使融合决策具有严格的统计基础；其次，建立了参数估计与滤波器运行的闭环交互，实现端到端的自适应优化；最后，通过周期性的EM更新平衡了计算复杂度与系统性能。相较于传统的固定权重融合方法，本框架能够根据信号特性动态调整融合策略，在保证估计精度的同时增强系统的鲁棒性。

3. 模拟数据及结果分析

3.1. 期望最大化算法优化粒子滤波和卡尔曼滤波压制随机噪声

在状态估计中，有效地压制随机观测噪声与过程噪声，从而还原出平滑、高保真的真实状态信号，是评价滤波算法性能的关键指标。在本研究提出的EM-KF-PF融合框架中，期望最大化算法通过扰动粒子进行搜索，其中参数δ作为高斯扰动的标准差，其取值直接决定了算法在探索能力与噪声压制能力之间的权衡。

本研究所用的真实状态信号为：

$x\left(t\right)=0.5\sin \left(2\pi \cdot 0.5\cdot t\right)+0.3\tanh \left(3\cdot \left(t-5\right)\right)$

该信号兼具低频光滑特性和陡变的非线性过渡，为检验算法在平滑段压制噪声与在突变段保持响应性的双重能力提供了理想的测试基准。

为评估δ对噪声压制效果的影响，本研究对比了δ = 0.001, 0.005, 0.01, 0.05, 0.1, 0.2六组参数下的滤波输出，如表1。仿真结果表明，δ的取值对估计结果的平滑度与噪声水平有决定性影响：

① 当δ取值较大(如0.1，0.2)时，强劲的高斯扰动虽有利于全局探索，但本身已成为一个显著的噪声源。其估计结果在整个时间域内都伴随着高频抖动，严重劣化了输出信号的信噪比。即便在信号平稳的正弦阶段，也无法有效平滑观测噪声，估计方差大，噪声压制能力最弱。

② 当δ取值过小(如<0.01)时，算法虽能输出极其平滑的轨迹，但其代价是牺牲了必要的探索性。这种过度的“平滑化”会导致算法无法有效响应突变点附近的真实状态突变，表现为估计轨迹变化迟缓，将真实的动态变化误当作噪声进行了压制，从而引入较大的滞后偏差。

Table 1. Impact of different types of noise on noise suppression effectiveness

表1. 不同噪声对噪声压制效果的影响

δ	0.001	0.005	0.01	0.05	0.1	0.2
去噪前SNR (dB)	52.81	38.51	32.31	18.44	12.24	6.53
去噪后SNR (dB)	33.54	33.06	34.04	31.37	26.94	21.83
SNR变化(dB)	−19.27	−5.45	1.73	12.93	14.71	15.30
MSE	0.000078	0.000088	0.00007	0.000129	0.000358	0.001163

因此，δ = 0.05被证实是赋予EM-KF-PF算法最强噪声压制能力的关键参数。它确保算法能够甄别并滤除随机噪声，同时保留并快速响应真实的状态跃迁，最终实现全局最优的滤波平滑效果。后续所有对比实验均基于此最优参数配置展开。

3.2. 蚁群算法优化粒子滤波和卡尔曼滤波性能分析

本小节将比较不同去噪方法对单个雷克子波的处理结果，采用的子波主频为30 Hz，中心位置为0.5 s，截取的部分为0~1 s，图1直观地展示了各方法对含噪雷克子波的处理结果。可以看出：

KF结果曲线较为平滑，但在波峰和波谷等信号变化剧烈处出现了明显的平滑效应，导致波形幅值被削弱，细节特征丢失，存在过平滑现象。

PF算法在波形主体部分跟踪良好，但在信号平稳区(如波峰前后)其估计曲线存在细微的抖动，这是由于粒子退化与重采样带来的估计方差。

KF-PF固定权重融合方法综合了KF的平滑性与PF的跟踪能力，其效果优于前两种单一方法，但在参数固定不变的情况下，其性能提升有限，无法自适应信号不同阶段的特性。

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1. Processing results of different denoising methods on a single Ricker wavelet. (a) Original signal; (b) Signal after noise addition; (c) Result after KF denoising; (d) Result after PF denoising; (e) Result after KF-PF denoising; (f) Result after EM-KF-PF denoising

图1. 不同去噪方法对单个雷克子波的处理结果。(a) 原始信号；(b) 加入噪声后信号；(c) KF去噪后结果；(d) PF去噪后结果；(e) KF-PF去噪后结果；(f) EM-KF-PF去噪后结果

EM-KF-PF动态融合方法的去噪信号与原始真实信号吻合度最高。它不仅有效地抑制了背景噪声，而且最为完美地保留了波峰的幅值特征与波形的细节信息，显著优于其他三种对比方法。这证明了EM算法通过动态优化关键参数，成功地引导融合系统在“平滑噪声”与“跟踪动态”之间找到了最佳平衡点。

4. 实际数据及结果分析

实验采用地震监测数据，采样率为50 Hz。

图2展示了EM-KF-PF混合方案、PF与KF对人工添加的高斯噪声(σ = 1.0)的去噪效果对比。KF处理结果图2(a)显示出典型的线性滤波特性：信号细节过度平滑，并且明显去除了部分原始数据的细节，反映KF对非平稳噪声适应性不足。PF处理结果图2(b)在细节保留上优于KF，微幅震荡可见度提升，但粒子采样噪声导致输出信号存在随机毛刺。图2(c)为KF-PF的处理结果，相比KF和PF单一处理效果好一些，但仍然存在随机毛刺。图2(d)为EM-KF-PF的处理结果，可以看出该方法能去噪后的图像干净、平滑，并且保留了图像的细节。

5. 评价标准

5.1. 均方根误差(RMSE)

RMSE衡量估计值$x_t$ 与真实状态$x_t^{true}$ 之间的平均平方根偏差，值越小表示估计越精确。

计算公式：$\text{RMSE}=\sqrt{\frac{1}{T}{\displaystyle {\sum }_{t=1}^T{\left(x_t-x_t^{true}\right)}^2}}$ 。

5.2. 信噪比(SNR)

SNR衡量信号(真实状态)与噪声(估计误差)的能量比，单位为分贝(dB)，值越大表示噪声影响越小。

计算公式：$\text{SNR}=10\cdot {\log }_{10}\left(\frac{P_{signal}}{P_{noise}}\right)$ 。

其中：$P_{signal}=\frac{1}{T}{\displaystyle {\sum }_{t=1}^T{\left(x_t^{true}\right)}^2}$ 是信号功率(真实状态的均方值)。

$P_{noise}=\frac{1}{T}{\displaystyle {\sum }_{t=1}^T{\left(x_t^{true}-x_t\right)}^2}$ 是噪声功率(MSE本身)。

Figure 2. Denoising effects of different methods. (a) Signal processed by KF; (b) Signal processed by PF; (c) Signal processed by the KF-PF method; (d) Signal processed by the EM-KF-PF method

图2. 不同方法去噪效果。(a) KF处理后信号；(b) PF处理后信号；(c) KF-PF方法处理后信号；(d) EM-KF-PF方法处理后信号

5.3. 相关系数

$\rho =\frac{{{\displaystyle \sum }}_{i=1}^N\left(x_i-\overline{x}\right)\left({\widehat{x}}_i-\overline{\widehat{x}}\right)}{\sqrt{{{\displaystyle \sum }}_{I=1}^N{\left(x_i-\overline{x}\right)}^2}\cdot \sqrt{{{\displaystyle \sum }}_{I=1}^N{\left({\widehat{x}}_i-\overline{\widehat{x}}\right)}^2}}$

5.4. 峰值信噪比(PSNR)

$\text{PSNR}=20\cdot {\log }_{10}\left(\frac{{\text{MAX}}_{\text{Z}}}{\sqrt{\text{MSE}}}\right)$

6. 结果与讨论

表2展示了四种去噪算法在复杂噪声环境下的性能对比结果。原始信号受到低频基线漂移与非平稳噪声的严重干扰，时域波形存在明显畸变。从整体性能来看，EM-KF-PF算法在均方根误差(RMSE)与信号保真度(相关系数)两个核心指标上均表现最佳，验证了其综合优势。

具体而言，KF凭借其最优估计特性，取得了RMSE为0.105994、相关系数为0.8976的成绩，实现了高效的背景平滑，但其强烈的线性高斯假设导致对信号中突发峰值的过度抑制，造成动态特征损失。PF的性能与KF相近(RMSE：0.106128，相关系数：0.8973)，其粒子集合虽具备理论上的非线性逼近能力，但固有的粒子退化问题增大了后验分布的估计方差，在信号平稳时段引入了不应有的细微抖动。

值得深入讨论的是，EM-KF-PF方案取得了最低的RMSE(0.106006)与并列最高的相关系数(0.8976)，但其SNR改善值(1.94 dB)略低于KF-PF方案(1.97 dB)。这一表面上的不一致性，恰恰揭示了不同评价指标的物理内涵差异：RMSE与相关系数更侧重于衡量重构信号与真实信号在波形形态上的全局一致性，而SNR改善则更敏感于全频带内信号与噪声的能量比例。分析认为，EM-KF-PF框架通过EM算法在线精确估计并动态修正噪声统计特性，其优势主要体现在对信号关键动态特征(如瞬态峰值)的更精确重构上，这种局部的、时域上的精度提升，对全局RMSE和相关系数贡献显著，但由于其策略性地保留了更多具有诊断价值的瞬变细节(这些细节在严格意义上也包含一定的能量)，可能导致其在全局能量视角的SNR指标上提升不显著。这恰恰说明，EM-KF-PF方案的成功之处在于其从追求单一化的噪声抑制，转向了在抑制噪声与保留关键特征之间寻求有效平衡。

Table 2. Comparison of denoising metrics for different methods

表2. 不同方法去噪指标比较

方法	RMSE	相关系数	PSNR (dB)	SNR改善(dB)
KF	0.105994	0.8976	16.58	1.96
PF	0.106128	0.8973	16.57	1.95
KF-PF	0.106060	0.8975	16.58	1.97
EM-KF-PF	0.106006	0.8976	16.58	1.94

7. 结论

本论文以提升复杂信号去噪性能为核心目标，通过系统的理论分析和方法创新，设计并实现了一种基于EM算法与KF-PF深度融合的自适应去噪框架。现将主要研究结论归纳如下：

1) 通过对比实验验证了单一滤波器在特定场景下的局限性：标准卡尔曼滤波(KF)在线性高斯条件下具有良好的去噪效率，但在面对复杂非高斯噪声与非线性动态变化时，其固定参数模型适应性不足，导致在信号峰值与瞬变区域性能下降。粒子滤波(PF)在处理非高斯问题时表现出较好的灵活性，但其粒子退化现象及较高的计算复杂度，限制了在实际系统中的应用。

2) 模型融合策略展现出明显优势：KF-PF融合方法有效结合了KF的稳定性与PF的灵活性，整体性能优于单一滤波方法。该策略通过分工协作的方式，分别利用KF处理信号背景平滑部分与PF捕捉瞬态特征，为提升去噪效果提供了有效路径。

3) 动态权重机制有助于保留信号细节：基于峰值检测的动态权重分配策略，在信号峰值区域适当增强PF的参与度，在抑制噪声的同时，较好地保留了信号中的瞬态特征与峰值信息，从而缓解了传统去噪方法中常见的“过度平滑”现象。