简介这是一份基于Sage-Husa自适应卡尔曼滤波器的海浪磁场噪声抑制Matlab项目源码面向从事海洋电磁探测、水下目标磁异常检测或信号处理方向的研究者也适合有一定Matlab基础的新手学习。代码完整覆盖从海浪磁场噪声产生、PSD功率谱分析到自适应卡尔曼滤波Q/R矩阵修正与抑噪效果对比的全流程可帮助读者快速理解Sage-Husa自适应算法在海浪磁场干扰场景下的参数整定与实现技巧。压缩包共6个文件全部为.m脚本体积仅8KB结构紧凑便于逐段研读和二次开发。目前已有829人学习下载该校正版本经过开发者亲测可运行下载后若遇问题可联系作者指导适合作为算法验证与论文复现的参考实现。1. 海浪磁场噪声抑制为什么是 Sage-Husa 自适应卡尔曼滤波海洋磁测和水中磁性目标探测有个反直觉现象目标信号还没出现磁传感器输出已经在缓慢起伏幅度可从 0.1 nT 到十几 nT周期从几秒到二十几秒。第一反应通常是高通滤波但目标磁异常信号和海浪磁场噪声同处于 0.051 Hz 频段频域滤波只会把两者一起削平。真正可行的是自适应滤波把海浪磁场噪声的统计特性纳入状态空间模型用观测数据在线修正噪声方差再让滤波器自动调节增益。Sage-Husa 自适应卡尔曼滤波是这一类估计器里工程落地最常见的方案它不需要额外传感器测量浪高或波向只靠磁力仪自身的测量新息就能跟踪噪声变化。这篇文章按「噪声怎么产生—算法原理—可复现代码—参数整定—工程增强」的顺序展开适合做海洋磁测数据处理、水下平台磁探测和信号降噪的工程师参考。2. 海浪磁场噪声的产生机理与频谱分布2.1 海水切割地磁力线感应电流如何变成磁场噪声海水是良导体电导率在 35 S/m 之间。地磁场在海水中并不是静态背景场波浪运动会推动导电海水以轨道速度切割地磁力线形成感应电场和感应电流这个电流体系会在海面附近产生二次磁场。磁传感器感受到的海浪磁场噪声本质就是这种感应磁场的短周期起伏。量级估算可以用一个粗粒度关系感应磁场幅度正比于 μ0、海水电导率 σ、波浪轨道速度 u、地磁场强度 B 以及电流环的特征尺度 H 的乘积。简单代入开阔海域的典型参数感应磁场的幅度通常落在 0.1 到几 nT 的范围。浪高 23 m、波周期 810 s 的涌浪条件下轨道速度增大感应磁场可以抬升到 10 nT 以上。这片海域如果再有明显的地磁异常背景场二次感应的幅度还会进一步放大。2.2 海浪谱的主频噪声能量集中在哪个频段海浪不是单一频率的简谐波而是由风浪和涌浪叠加成的随机过程。描述海浪能量分布的常用工具是海浪谱Pierson-Moskowitz 谱足够说明问题谱峰值频率随风速下降而降低风速 10 m/s 时充分成长风浪的峰值周期大约在 911 s对应的主频在 0.1 Hz 附近。涌浪周期更长可以到 1520 s对应的频率低到 0.05 Hz 左右。海况有效波高主周期海浪磁场噪声主频带轻浪风浪初生0.20.8 m36 s0.20.5 Hz中浪常见作业海况12 m710 s0.10.2 Hz大浪涌浪为主2.54 m1116 s0.060.12 Hz海洋磁测受这个频段的直接影响低通滤波能压掉高频电磁干扰但对海浪感应磁场无能为力因为它的能量分布恰好落在磁异常目标信号所在的位置。2.3 为什么常规滤波无法分离目标信号与噪声的频谱混叠磁力仪测量水中磁性目标时目标在传感器前方通过信号包络持续时间由最近接近距离和平台速度共同决定。拖体以 2 m/s 速度航行、与目标最近距离约 30 m 时信号包络持续约 1520 s主瓣能量折算到频率域就是 0.050.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}它反映的是实际量测偏离预测的程度。引入遗忘因子 b0 b 1定义权重 d_k (1-b) / (1-b^(k1))量测噪声方差的递推式为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.960.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.0H 1.0Q 固定为一个极小值让滤波器主要跟随背景磁场的缓慢变化海浪噪声的方差变化交给 R 估计器。实现用标量运算比矩阵版本短得多逻辑也更直观。class SageHusaKF: 单通道 Sage-Husa 自适应卡尔曼滤波固定 Q 估计 R def __init__(self, F, H, Q, R0, P0, x0, b0.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) 个采样点不同取值对应完全不同的跟踪行为。b等效记忆采样点10 Hz 下的时间尺度适用场景0.9010 点1 s海况突变极快但 R 估计抖动大0.9520 点2 s阵风导致的快速海况变化0.97540 点4 s常规风浪变化跟踪与平滑的折中0.99100 点10 s涌浪缓慢变化R 估计平稳经验法则是把记忆长度设在海浪主周期的 0.51 倍之间。主周期 8 s、采样率 10 Hz 时主周期对应 80 个采样点b 取 0.9750.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.8499% 分位点约 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 点超限这个条件把磁异常目标脉冲和海况突变区分开目标过顶信号只持续 13 个采样点不足以连续触发 5 次海况突变导致噪声方差上升后新息会在相当长的时间里持续超限。应用时还需要维护一个长度为 W 的观测滑窗W 取记忆长度的 23 倍z_window更新方式很简单用collections.deque(maxlenW)即可。这层保护的意义在于它和 R 估计器是并行关系R 递推负责跟踪缓慢变化卡方检核负责识别快速跳变并把 R 直接拉到样本方差水平。两者配合海况突变后滤波器不需要等几十秒而是可以在一两个记忆周期内恢复平滑输出目标信号被误吞的概率也会明显降低。本文还有配套的精品资源点击获取
Sage-Husa自适应卡尔曼滤波：海洋磁测海浪磁场噪声抑制实战
简介这是一份基于Sage-Husa自适应卡尔曼滤波器的海浪磁场噪声抑制Matlab项目源码面向从事海洋电磁探测、水下目标磁异常检测或信号处理方向的研究者也适合有一定Matlab基础的新手学习。代码完整覆盖从海浪磁场噪声产生、PSD功率谱分析到自适应卡尔曼滤波Q/R矩阵修正与抑噪效果对比的全流程可帮助读者快速理解Sage-Husa自适应算法在海浪磁场干扰场景下的参数整定与实现技巧。压缩包共6个文件全部为.m脚本体积仅8KB结构紧凑便于逐段研读和二次开发。目前已有829人学习下载该校正版本经过开发者亲测可运行下载后若遇问题可联系作者指导适合作为算法验证与论文复现的参考实现。1. 海浪磁场噪声抑制为什么是 Sage-Husa 自适应卡尔曼滤波海洋磁测和水中磁性目标探测有个反直觉现象目标信号还没出现磁传感器输出已经在缓慢起伏幅度可从 0.1 nT 到十几 nT周期从几秒到二十几秒。第一反应通常是高通滤波但目标磁异常信号和海浪磁场噪声同处于 0.051 Hz 频段频域滤波只会把两者一起削平。真正可行的是自适应滤波把海浪磁场噪声的统计特性纳入状态空间模型用观测数据在线修正噪声方差再让滤波器自动调节增益。Sage-Husa 自适应卡尔曼滤波是这一类估计器里工程落地最常见的方案它不需要额外传感器测量浪高或波向只靠磁力仪自身的测量新息就能跟踪噪声变化。这篇文章按「噪声怎么产生—算法原理—可复现代码—参数整定—工程增强」的顺序展开适合做海洋磁测数据处理、水下平台磁探测和信号降噪的工程师参考。2. 海浪磁场噪声的产生机理与频谱分布2.1 海水切割地磁力线感应电流如何变成磁场噪声海水是良导体电导率在 35 S/m 之间。地磁场在海水中并不是静态背景场波浪运动会推动导电海水以轨道速度切割地磁力线形成感应电场和感应电流这个电流体系会在海面附近产生二次磁场。磁传感器感受到的海浪磁场噪声本质就是这种感应磁场的短周期起伏。量级估算可以用一个粗粒度关系感应磁场幅度正比于 μ0、海水电导率 σ、波浪轨道速度 u、地磁场强度 B 以及电流环的特征尺度 H 的乘积。简单代入开阔海域的典型参数感应磁场的幅度通常落在 0.1 到几 nT 的范围。浪高 23 m、波周期 810 s 的涌浪条件下轨道速度增大感应磁场可以抬升到 10 nT 以上。这片海域如果再有明显的地磁异常背景场二次感应的幅度还会进一步放大。2.2 海浪谱的主频噪声能量集中在哪个频段海浪不是单一频率的简谐波而是由风浪和涌浪叠加成的随机过程。描述海浪能量分布的常用工具是海浪谱Pierson-Moskowitz 谱足够说明问题谱峰值频率随风速下降而降低风速 10 m/s 时充分成长风浪的峰值周期大约在 911 s对应的主频在 0.1 Hz 附近。涌浪周期更长可以到 1520 s对应的频率低到 0.05 Hz 左右。海况有效波高主周期海浪磁场噪声主频带轻浪风浪初生0.20.8 m36 s0.20.5 Hz中浪常见作业海况12 m710 s0.10.2 Hz大浪涌浪为主2.54 m1116 s0.060.12 Hz海洋磁测受这个频段的直接影响低通滤波能压掉高频电磁干扰但对海浪感应磁场无能为力因为它的能量分布恰好落在磁异常目标信号所在的位置。2.3 为什么常规滤波无法分离目标信号与噪声的频谱混叠磁力仪测量水中磁性目标时目标在传感器前方通过信号包络持续时间由最近接近距离和平台速度共同决定。拖体以 2 m/s 速度航行、与目标最近距离约 30 m 时信号包络持续约 1520 s主瓣能量折算到频率域就是 0.050.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}它反映的是实际量测偏离预测的程度。引入遗忘因子 b0 b 1定义权重 d_k (1-b) / (1-b^(k1))量测噪声方差的递推式为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.960.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.0H 1.0Q 固定为一个极小值让滤波器主要跟随背景磁场的缓慢变化海浪噪声的方差变化交给 R 估计器。实现用标量运算比矩阵版本短得多逻辑也更直观。class SageHusaKF: 单通道 Sage-Husa 自适应卡尔曼滤波固定 Q 估计 R def __init__(self, F, H, Q, R0, P0, x0, b0.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) 个采样点不同取值对应完全不同的跟踪行为。b等效记忆采样点10 Hz 下的时间尺度适用场景0.9010 点1 s海况突变极快但 R 估计抖动大0.9520 点2 s阵风导致的快速海况变化0.97540 点4 s常规风浪变化跟踪与平滑的折中0.99100 点10 s涌浪缓慢变化R 估计平稳经验法则是把记忆长度设在海浪主周期的 0.51 倍之间。主周期 8 s、采样率 10 Hz 时主周期对应 80 个采样点b 取 0.9750.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.8499% 分位点约 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 点超限这个条件把磁异常目标脉冲和海况突变区分开目标过顶信号只持续 13 个采样点不足以连续触发 5 次海况突变导致噪声方差上升后新息会在相当长的时间里持续超限。应用时还需要维护一个长度为 W 的观测滑窗W 取记忆长度的 23 倍z_window更新方式很简单用collections.deque(maxlenW)即可。这层保护的意义在于它和 R 估计器是并行关系R 递推负责跟踪缓慢变化卡方检核负责识别快速跳变并把 R 直接拉到样本方差水平。两者配合海况突变后滤波器不需要等几十秒而是可以在一两个记忆周期内恢复平滑输出目标信号被误吞的概率也会明显降低。本文还有配套的精品资源点击获取
