Sage-Husa自适应卡尔曼滤波：海洋磁测海浪磁场噪声抑制实战

简介：这是一份基于Sage-Husa自适应卡尔曼滤波器的海浪磁场噪声抑制Matlab项目源码，面向从事海洋电磁探测、水下目标磁异常检测或信号处理方向的研究者，也适合有一定Matlab基础的新手学习。代码完整覆盖从海浪磁场噪声产生、PSD功率谱分析到自适应卡尔曼滤波Q/R矩阵修正与抑噪效果对比的全流程，可帮助读者快速理解Sage-Husa自适应算法在海浪磁场干扰场景下的参数整定与实现技巧。压缩包共6个文件，全部为.m脚本，体积仅8KB，结构紧凑，便于逐段研读和二次开发。目前已有829人学习下载，该校正版本经过开发者亲测可运行，下载后若遇问题可联系作者指导，适合作为算法验证与论文复现的参考实现。

1. 海浪磁场噪声抑制：为什么是 Sage-Husa 自适应卡尔曼滤波

海洋磁测和水中磁性目标探测有个反直觉现象：目标信号还没出现，磁传感器输出已经在缓慢起伏，幅度可从 0.1 nT 到十几 nT，周期从几秒到二十几秒。第一反应通常是高通滤波，但目标磁异常信号和海浪磁场噪声同处于 0.05～1 Hz 频段，频域滤波只会把两者一起削平。真正可行的是自适应滤波：把海浪磁场噪声的统计特性纳入状态空间模型，用观测数据在线修正噪声方差，再让滤波器自动调节增益。Sage-Husa 自适应卡尔曼滤波是这一类估计器里工程落地最常见的方案，它不需要额外传感器测量浪高或波向，只靠磁力仪自身的测量新息就能跟踪噪声变化。这篇文章按「噪声怎么产生—算法原理—可复现代码—参数整定—工程增强」的顺序展开，适合做海洋磁测数据处理、水下平台磁探测和信号降噪的工程师参考。

2. 海浪磁场噪声的产生机理与频谱分布

2.1 海水切割地磁力线：感应电流如何变成磁场噪声

海水是良导体，电导率在 3～5 S/m 之间。地磁场在海水中并不是静态背景场，波浪运动会推动导电海水以轨道速度切割地磁力线，形成感应电场和感应电流，这个电流体系会在海面附近产生二次磁场。磁传感器感受到的"海浪磁场噪声"，本质就是这种感应磁场的短周期起伏。

量级估算可以用一个粗粒度关系：感应磁场幅度正比于 μ0、海水电导率 σ、波浪轨道速度 u、地磁场强度 B 以及电流环的特征尺度 H 的乘积。简单代入开阔海域的典型参数，感应磁场的幅度通常落在 0.1 到几 nT 的范围。浪高 2～3 m、波周期 8～10 s 的涌浪条件下，轨道速度增大，感应磁场可以抬升到 10 nT 以上。这片海域如果再有明显的地磁异常背景场，二次感应的幅度还会进一步放大。

2.2 海浪谱的主频：噪声能量集中在哪个频段

海浪不是单一频率的简谐波，而是由风浪和涌浪叠加成的随机过程。描述海浪能量分布的常用工具是海浪谱，Pierson-Moskowitz 谱足够说明问题：谱峰值频率随风速下降而降低，风速 10 m/s 时，充分成长风浪的峰值周期大约在 9～11 s，对应的主频在 0.1 Hz 附近。涌浪周期更长，可以到 15～20 s，对应的频率低到 0.05 Hz 左右。

海洋磁测受这个频段的直接影响：低通滤波能压掉高频电磁干扰，但对海浪感应磁场无能为力，因为它的能量分布恰好落在磁异常目标信号所在的位置。

2.3 为什么常规滤波无法分离：目标信号与噪声的频谱混叠

磁力仪测量水中磁性目标时，目标在传感器前方通过，信号包络持续时间由最近接近距离和平台速度共同决定。拖体以 2 m/s 速度航行、与目标最近距离约 30 m 时，信号包络持续约 15～20 s，主瓣能量折算到频率域就是 0.05～0.2 Hz，和中浪到大浪情况下的海浪磁场噪声主频带几乎完全重叠。

这种情况下，固定截止频率的滤波器无法在保留目标信号的同时抑制噪声，经典的固定参数卡尔曼滤波也只是把 R 矩阵设成一个估计值，海况一旦变化就失效。真正需要的是能在线感知噪声方差变化的估计器，这正是 Sage-Husa 自适应卡尔曼滤波器的切入点。

3. Sage-Husa 自适应卡尔曼滤波器的原理与可估性分析

3.1 标准卡尔曼滤波在海上磁测中的失效模式

离散线性卡尔曼滤波的状态空间模型写作：

x_k = F x_{k-1} + w_k z_k = H x_k + v_k

其中 w_k 和 v_k 分别是过程噪声和量测噪声，协方差矩阵记为 Q 和 R。标准卡尔曼滤波把 Q、R 当作已知常数，在磁测场景中这就是最大的失配来源：海浪磁场噪声的方差随海况变化，R 在 0.0025 nT² 这样一个量级的水平上完全无法表达"上午还是轻浪、下午变成大浪"的真实过程。

R 设小了，滤波器过度信任当前量测，输出会呈现出和海浪噪声几乎一致的起伏；R 设大了，滤波器对量测响应迟钝，小幅度磁异常目标信号被平滑到检不出。实际处理时很多团队会拿一段数据试出"合适"的 R，但这一套在海况稳定时可以工作，遇到阵风过境海况突变，输出的信噪比立刻退化。

3.2 新息驱动噪声统计估计：Sage-Husa 的递推公式

Sage-Husa 的核心思想是用量测新息在线估计噪声统计量。定义新息为 y_k = z_k - H x_{k|k-1}，它反映的是实际量测偏离预测的程度。引入遗忘因子 b（0 < b < 1），定义权重 d_k = (1-b) / (1-b^(k+1))，量测噪声方差的递推式为：

R_k = (1 - d_k) R_{k-1} + d_k (y_k y_k^T - H P_{k|k-1} H^T)

这个式子的物理意义很清楚：y_k y_k^T 是当前量测方差的样本估计，H P_{k|k-1} H^T 是滤波器预测的不确定度，两者相减才是真正由量测噪声贡献的部分。如果新息突然变大，R_k 会跟着变大，滤波器自动降低对量测的信任，增益矩阵 K_k 收小，从而抑制海浪噪声的冲击。

需要强调的是，这个估计器还远不是被认为可以随意使用的工程工具，Sage-Husa 的完整形式还包含过程噪声 Q 的估计式，但工程上实际采用时通常要放弃同时估计 Q 和 R。

3.3 Q 与 R 的可辨识性：固定 Q 估 R 是工程默认

Q 和 R 都从新息序列中提取信息，在标量或低阶系统中同时估计两者，会出现所谓的"竞争"现象：新息变大时，R 估计器和 Q 估计器各自把变化归因到自己的头上，估计结果来回跳动，甚至导致协方差矩阵失去非负定性，滤波器发散。这是个理论上有解、工程上不实用的典型例子。

海洋磁测场景有一个有利条件：状态方程描述的是传感器平台附近背景磁场的缓慢变化，模型相对可信，Q 可以固定成一个小值；而变化剧烈的海浪感应磁场，作用在量测方程上，正好对应 R 的时变。固定 Q 估 R，既规避了同时估计的发散风险，又命中了问题的主要矛盾。

3.4 遗忘因子和新息权重的配合

遗忘因子 b 控制记忆长度。令 L = 1/(1-b)，b = 0.975 时约等效保留最近 40 个采样点的信息，采样率 10 Hz 下就是 4 秒左右。这个时间尺度和海浪主周期在同一量级，R 估计既不会对单次噪声波动过度反应，也不会在海况突变时反应太慢。b 的具体选择在第 5 章给出，原理上它是"跟踪速度"和"估计方差"之间的折中，这一点在下文仿真中可以直接看到。

4. 用 Python 仿真海浪磁场噪声并实现 Sage-Husa 抑制算法

4.1 仿真信号设计：AR(1)海浪噪声与过顶磁异常脉冲

海浪磁场噪声的频率集中在低频段，工程上常用一阶自回归模型 AR(1) 近似它的谱特性。下面的代码生成 10 分钟仿真数据，先产生平稳噪声，再在 400 s 处把噪声幅度放大到 4 倍，模拟海况恶化。

import numpy as np fs = 10.0 # 采样率 10 Hz N = 6000 # 总样本数 a = 0.92 # AR(1) 系数，决定噪声谱峰位置 sigma_n = 0.05 # 平稳段噪声标准差，单位 nT rng = np.random.default_rng(42) noise = np.zeros(N) innov = rng.standard_normal(N) noise[0] = innov[0] * sigma_n for k in range(1, N): # AR(1) 递推：当前噪声由上一时刻乘以系数 a，叠加上新息 noise[k] = a * noise[k-1] + sigma_n * np.sqrt(1 - a*a) * innov[k] noise[4000:] *= 4.0 # 400 s 后幅度放大到 4 倍，方差变为 16 倍

AR(1) 系数 a 控制自相关的衰减速度。10 Hz 采样率下，a 取 0.92 时噪声谱峰约在 0.12 Hz，适合模拟中浪海况；如果目标海域以长周期涌浪为主，可以把 a 提到 0.96～0.98，能量会更靠近 0.05 Hz。

接下来加入磁异常目标信号。磁性目标过顶时传感器测得的典型波形是双极性脉冲，用一个高斯导函数近似。

def mag_target(t, t0, A, sigma_t): # 磁性目标过顶信号：过零点前后极性相反，峰-峰间隔约 2*sigma_t return A * (t - t0) / sigma_t * np.exp(-((t - t0) ** 2) / (2 * sigma_t ** 2)) t = np.arange(N) / fs target = mag_target(t, 200.0, 0.4, 8.0) # 目标在 200 s 处出现 z = noise + target # 磁力仪观测值

t0 取 200 s，σ 取 8 s，对应平台速度 2 m/s、最近接近距离约 30 m 的目标信号。目标在噪声突变之前出现，方便对比两种滤波器在有目标和无目标时段的表现。

4.2 Sage-Husa 自适应卡尔曼滤波的最小实现

状态模型选一阶随机游走，F = 1.0，H = 1.0，Q 固定为一个极小值，让滤波器主要跟随背景磁场的缓慢变化，海浪噪声的方差变化交给 R 估计器。实现用标量运算，比矩阵版本短得多，逻辑也更直观。

class SageHusaKF: """单通道 Sage-Husa 自适应卡尔曼滤波，固定 Q 估计 R""" def __init__(self, F, H, Q, R0, P0, x0, b=0.975): self.F = F # 状态转移系数，随机游走取 1.0 self.H = H # 量测系数，直接观测状态取 1.0 self.Q = Q # 固定的过程噪声方差 self.R = R0 # 初始量测噪声方差 self.P = P0 # 初始协方差 self.x = x0 # 初始状态 self.b = b # 遗忘因子 self.k = 0 def step(self, z): # 时间更新：预测当前状态和协方差 x_pred = self.F * self.x P_pred = self.F * self.P * self.F + self.Q # 计算新息和卡尔曼增益 y = z - self.H * x_pred S = self.H * P_pred * self.H + self.R K = P_pred * self.H / S # 状态和协方差更新 self.x = x_pred + K * y self.P = (1 - K * self.H) * P_pred # R 估计更新：d_k 随步数增加而递减 self.k += 1 d = (1 - self.b) / (1 - self.b ** (self.k + 1)) self.R = (1 - d) * self.R + d * (y * y - self.H * P_pred * self.H) self.R = max(self.R, 1e-4) # 下限保护，防止 R 变负 return self.x, self.R

代码里有三个要点。第一，时间更新先于量测更新，y 用的是量测减去预测值的残差，而不是减去上一次滤波输出。第二，R 估计式中的y * y是当前时刻的样本方差，需要减去H * P_pred * H这部分预测不确定度，才是纯粹的量测噪声贡献。第三，理论上两个量相减可能得到负值，这是协方差估计中最经典的失效路径，必须加下限保护。

4.3 对比基准：固定 R 的经典卡尔曼滤波

固定 R 的卡尔曼滤波可以直接复用同一个类，把 R 的更新部分去掉：

class FixedRKF: """固定量测噪声方差的标准卡尔曼滤波，作为对比基准""" def __init__(self, F, H, Q, R, P0, x0): self.F = F; self.H = H; self.Q = Q; self.R = R self.P = P0; self.x = x0 def step(self, z): x_pred = self.F * self.x P_pred = self.F * self.P * self.F + self.Q y = z - self.H * x_pred S = self.H * P_pred * self.H + self.R K = P_pred * self.H / S self.x = x_pred + K * y self.P = (1 - K * self.H) * P_pred return self.x

注意这里 R 初始值取的是平稳段噪声方差 0.0025，这个取值在对比实验中是最优的——海况突变前两个滤波器几乎没有差别。正是因此，突变后的差异才完全归因于 R 的自适应能力。

4.4 结果怎么看：R 跟踪曲线和输出标准差

跑完两组滤波后，重点看两个指标。一是 R 估计曲线：Sage-Husa 的 R 应该在 400 s 后的几十秒内从 0.0025 爬向 0.04 附近，与噪声方差的实际变化对应；固定 R 的滤波器的 R 永远停在 0.0025。二是输出标准差：固定 R 滤波器在突变后输出会明显毛糙，因为它的增益仍然按低噪声环境计算，把海浪噪声当成真实信号跟随。

实际运行代码会看到，Sage-Husa 在海况突变后输出的幅值起伏明显小于固定 R 滤波器，但因为遗忘因子 b = 0.975 约保留 4 秒记忆，R 估计有一个爬升过程，突变头几秒仍然会漏进一部分噪声。这个滞后正是后面第 6 章要解决的问题。

5. Sage-Husa 的参数整定与工程避坑

5.1 遗忘因子 b 的物理含义和选择表

遗忘因子是 Sage-Husa 中最先要定的参数。它不直接参与滤波，却决定了 R 估计能看多远的历史。有效记忆长度约等于 L = 1/(1-b) 个采样点，不同取值对应完全不同的跟踪行为。

经验法则是把记忆长度设在海浪主周期的 0.5～1 倍之间。主周期 8 s、采样率 10 Hz 时，主周期对应 80 个采样点，b 取 0.975～0.99 都算合理。b 太小，R 估计会被单个大浪新息带着跳；b 太大，R 估计平滑但跟不上海况。

5.2 R0、Q、P0 初始化：先预热，再启动滤波

R 的初值直接给 0 或任意小值都会让第一步增益异常放大。常见做法是先取一段不含目标信号的数据做预热，用样本方差作为 R0。

# 用前 500 个点（50 秒）的观测方差初始化 R0 warm = z[:500] R0 = np.var(warm) Q = 1e-6 # 状态演化不确定度，远小于 R 的典型值 P0 = R0 # 初始协方差直接设为 R0 量级 x0 = z[0] # 初始状态用第一个观测值

Q 的经验取值是R0 的 1/1000 到 1/100，给太大会让滤波器过度相信新息，Sage-Husa 的 R 估计会失去意义。P0 表示初始状态的不确定度，设成 R0 量级是安全的，不要给零矩阵，否则滤波前期增益计算会出现数值异常。

5.3 三个高频坑：R 变负、同时估 Q/R 发散、模型阶数不足

R 变负是最常见的失效模式。y_k² 小于 H P_{k|k-1} H^T 时，R 的递推值会变成负数。这种情形在噪声方差突降或模型误差偏大时很容易触发，必须加下限保护，建议下限取 R0 的 1/100。

同时估计 Q 和 R 的完整 Sage-Husa 形式在标量系统里也容易发散，因为新息无法同时区分过程噪声和量测噪声的贡献。工程上只估 R、固定 Q，除非系统模型本身有明确的物理依据证明 Q 确实在变化。

第三个坑是状态模型阶数不足。用一阶随机游走模型时，海浪磁场噪声的有色特性会残留在新息里，导致 R 估计略偏高。这不是滤波器的错，而是白噪声假设和真实有色噪声之间的失配。要改善，可以从 AR(2) 状态模型或量测差分入手，但通常收益有限，且参数变多后调优成本明显上升，对于大部分磁测任务，固定 Q 估 R 的一阶模型已经够用。

6. 用新息卡方检核给 Sage-Husa 加一道抗海况突变保护

6.1 为什么 R 估计追不上海况突变

Sage-Husa 的 R 估计天然假设噪声统计是慢变的，遗忘因子让它对突变响应有滞后。第 4 章的 400 s 突变场景中，R 从 0.0025 爬向 0.04 需要几十秒，这段时间内滤波器处于"模型失配"状态，输出会混入明显的噪声尖峰，可能被误判为磁异常目标。

更麻烦的是磁异常目标本身也会让新息突增，如果只依赖 R 估计来吸收新息，目标信号会反过来抬高 R，滤波输出被压钝。要同时应对这两种情况，需要引入一个独立于 R 估计器的检测机制：新息卡方检核。

6.2 归一化新息平方与重初始化补丁

归一化新息平方定义为 γ = y² / (H P H^T + R)，在模型匹配的理想情况下服从自由度为 1 的卡方分布，95% 分位点约 3.84，99% 分位点约 6.63。但因为海浪磁场噪声是有色的，实际新息分布偏"重尾"，门限取 9.0 更稳。连续超限则判定海况变化，触发 R 重初始化：

# 在 SageHusaKF.step 内部增加的状态统计 gamma = y * y / (H * P_pred * H + self.R) if gamma > 9.0: alarm += 1 else: alarm = 0 if alarm >= 5: # 连续 5 点超限，判定噪声特性发生变化 self.R = max(np.var(z_window), R_min) # 用最近窗口的样本方差重置 R alarm = 0

连续 5 点超限这个条件把磁异常目标脉冲和海况突变区分开：目标过顶信号只持续 1～3 个采样点，不足以连续触发 5 次；海况突变导致噪声方差上升后，新息会在相当长的时间里持续超限。应用时还需要维护一个长度为 W 的观测滑窗，W 取记忆长度的 2～3 倍，z_window更新方式很简单，用collections.deque(maxlen=W)即可。

这层保护的意义在于它和 R 估计器是并行关系：R 递推负责跟踪缓慢变化，卡方检核负责识别快速跳变并把 R 直接拉到样本方差水平。两者配合，海况突变后滤波器不需要等几十秒，而是可以在一两个记忆周期内恢复平滑输出，目标信号被误吞的概率也会明显降低。

本文还有配套的精品资源，点击获取