DolphinDB 物联网滤波算法库

在工业物联网 (IIoT) 和车联网场景中，海量传感器数据的清洗与实时分析是企业所关注的重点。例如，在成千上万个温度传感器记录的数据中，往往会由于异常原因出现极端值，需要精准定位并修复。同时，对于某些特定的分析，诸如光谱分析、波形分析等，需要在平滑噪声的同时，尽可能保留信号的峰值和宽度等高阶矩信息。

传统的解决方法通常采用诸如 "数据库存储 + Python/Spark 计算"等的模式，其基本处理流程如下所示：

DB -> Network (Serialize) -> Python/Spark Memory (Deserialize) -> Compute -> Network -> DB

由于传统处理方式存在着链路长，延迟高，内存消耗大等难以忽视的缺陷，在业务处理中往往面临着如下痛点：

• 数据搬运昂贵：将每天 TB 级的传感器数据导出到 Python 环境进行清洗，网络带宽和 I/O 成为绝对瓶颈。

• 系统架构臃肿：维护两套异构系统（存储集群 + 计算集群），运维成本和硬件投入双倍增加，需要面对较为繁琐的环境和依赖配置，多语言异构系统维护的复杂度较高。

• 计算性能受限：以 Python 为主的计算方式难以逃脱其固有的解释型语言弊端，导致其在计算性能上有着显著的瓶颈。

DolphinDB 的架构及处理流程设计（DB Storage -> Memory -> Compute），提供了一套同时支撑高吞吐写入和复杂信号处理的系统，全流程在同一进程空间内完成，并利用列式存储直接进行向量化计算。

基于此，DolphinDB 推出了 IoT Filtering 模块作为海量传感器数据的清洗与实时分析这一业务场景的解决方案，相较于传统方式，该方案具有如下基本优势：

• 库内计算：滤波算法直接运行在数据存储节点。数据“零移动”，计算“就地发生”，解决了海量时序数据分析中的网络传输瓶颈。

• 运行、维护等成本降低：无需单独用于数据搬运（ETL） 的服务器集群，无需维护复杂的 Python 依赖环境，降低了多语言异构系统的维护复杂度，真正实现“即插即用”。

• 出色的运算能力：强大的向量化计算能力，以及“流批一体”的计算框架，在具备了极致性能的同时适用于多种业务场景。

限幅滤波 amplitudeLimitingFilter

• 算法含义

  限幅滤波用于消除传感器信号中的偶然大幅度突变。通过比较相邻采样值的差值与预设阈值，若差值超过阈值则舍弃当前值，保持上一次的输出值，从而滤除由干扰引起的异常跳变。

• 公式

  

  其中 Δ 为偏差阈值。

• 代码实现

  def amplitudeLimitingFilter(x, threshold){
      n = size(x)
      if(n < 2) return x
      res = array(DOUBLE, n)
      res[0] = x[0]
      lastVal = double(x[0])
      for(i in 1:n){
          currVal = x[i]
          if(abs(currVal - lastVal) > threshold){
              res[i] = lastVal
          }else{
              res[i] = currVal
              lastVal = currVal
          }
      }
      return res
  }

中位值滤波 medianFilter

• 算法含义

  中位值滤波取滑动窗口内所有采样值的中位数作为输出，能有效抑制脉冲噪声（如椒盐噪声），同时较好地保留信号的边缘信息。

• 公式

  
• 代码实现

  def medianFilter(x, window){
      return mmed(x, window)
  }

移动平均滤波 movingAverageFilter

• 算法含义

  移动平均滤波计算滑动窗口内所有采样值的算术平均值，用于平滑高频随机噪声，但会引入相位滞后。

• 公式

  

• N 为窗口大小

• 代码实现

  def movingAverageFilter(x, window){
      return mavg(x, window)
  }

Hampel 滤波 hampelFilter

• 算法含义

  Hampel 滤波是一种基于中位数和中位绝对偏差（MAD）的鲁棒离群点检测与替换方法。它识别并替换超出动态阈值的异常值，对数据中的极端值不敏感。

• 公式

  

  其中系数 1.4826 的作用是使 MAD 与标准差尺度一致。

• 代码实现

• def hampelFilter(x, window, k){
      med = mmed(x, window)
      absDev = abs(x - med)
      mad = mmed(absDev, window)
      sigma = 1.4826 * mad
      threshold = k * sigma
      isOutlier = absDev > threshold
      res = iif(nullFill(isOutlier, false), med, x)
      return res
  }

双指数平滑 doubleExponentialSmoothing

• 算法含义

  双指数平滑适用于具有趋势的时间序列。它通过两个平滑参数 α 和 β 分别估计当前水平（Level）和趋势（Trend），能跟踪数据的线性变化趋势。

• 公式

  

• 其中，初始值： L_1 = x_1,\quad T_1 = x_2 - x_1 。

• 代码实现

  def doubleExponentialSmoothing(x, alpha, beta){
      n = size(x)
      if(n < 2) return x
      a1 = 2.0 - alpha - alpha * beta
      a2 = alpha - 1.0
      F = alpha * x - (alpha * (1.0 - beta)) * prev(x)
      F[0] = x[0]
      F[1] = x[1] - a1 * x[0]
      return iterate([0.0, 0.0], [a2, a1], F)
  }

FIR 滤波器 firFilter

• 算法含义

  有限脉冲响应滤波器，通过卷积实现线性滤波，其输出仅由当前及过去输入的加权和决定，系统总是稳定且易于设计线性相位。

• 公式

  

  其中 b_k为滤波器系数（脉冲响应）。

• 代码实现

  def firFilter(x, b){
      M = size(b)
      n = size(x)
      
      if(n == 0 or M == 0) return array(DOUBLE, 0)
      if(M == 1) return x*b[0]
      
      //调用mavg函数，要保证weights参数，也就是b的长度在2-1024范围内
      if(M >= 2 and M <= 1024) return mavg(take(0.0, (M-1)).append!(x*sum(b)), reverse(b))[(M-1):]
      
      //当b的长度大于1024时，调用矩阵运算
      if(M > 1024) return dot(nullFill(matrix(loop(move{x, }, 0..(M-1))), 0.0), b)[0]
  }

一阶滞后滤波 firstOrderLagFilter

• 算法含义

  一阶滞后滤波是模拟 RC 低通电路的数字实现，也是最简单的 IIR 滤波器。

• 公式

  

  其中 α 为滤波系数。

• 代码实现

  def firstOrderLagFilter(x, alpha){
      return ewmMean(x, alpha=alpha, adjust=false)
  }

IIR 滤波器 iirFilter

• 算法含义

  无限脉冲响应滤波器利用了过去的输出值作为反馈，因此可以用比 FIR 低得多的阶数实现陡峭的过渡带。

• 公式：

  

  通常归一化使 α_0=1。

• 代码实现

  def iirFilter(x, b, a) {
      n = size(x)
      M = size(b)
      N = size(a)
      
      if(n == 0 or N < 1 or M < 1) return array(DOUBLE, 0)
      
      // 使用局部变量存储系数 (Use local variables for coefficients)
      b_local = b
      a_local = a
      
      // 向量化系数归一化（避免循环内重复计算）
      if(a[0] != 1.0){
          b_local = b / a[0]
          a_local = a / a[0]
      }
      
      if(M == 1){
          feedforward = x*b_local[0]
      }
      //调用mavg函数，要保证weights参数，也就是b的长度在2-1024范围内
      if(M >= 2 and M <= 1024){
          feedforward = mavg(take(0.0, (M-1)).append!(x*sum(b_local)), reverse(b_local))[(M-1):]
      }
      
      //当b的长度大于1024时，调用矩阵运算
      if(M > 1024){
          feedforward = dot(nullFill(matrix(loop(move{x, }, 0..(M-1))), 0.0), b_local)[0]
      }
      
      if((N-1) == 0) {
          return feedforward
      }
       
      return iterate(take(0.0, (N-1)), -reverse(a_local[1:]), feedforward)
  }

高斯滤波 gaussianFilter

• 算法含义

  高斯滤波是一种线性平滑滤波器，其使用高斯核进行加权移动平均，核函数由标准差 σ 决定。离中心点越近的样本，权重越大；越远，权重越小。这使得它的平滑效果非常自然，在时域和频域都没有振铃效应。

• 公式

  核函数：

  

  归一化后与信号卷积：

  

• 代码实现

  def gaussianFilter(x, window, sigma){
      if(window % 2 == 0) throw "Window size must be odd."
      half_window = (window - 1) / 2
      r = (0-half_window) .. half_window
      numerator = pow(r, 2)
      denominator = 2 * pow(sigma, 2)
      kernel = exp(0.0 - (numerator / denominator))
      kernel = kernel / sum(kernel)
      n = size(x)
      y = array(DOUBLE, n, n, 0)
      for(j in -half_window: (half_window+1)){
          idx = j + half_window
          w = kernel[idx]
          shifted_x = move(x, -j)
          shifted_x = nullFill(shifted_x, 0.0)
          y += w * shifted_x
      }
      return y
  }

  优化后的代码实现：

  // def gaussianFilter(x, window, sigma) {
  //     // 窗口奇数检查
  //     if(window % 2 == 0) throw "Window size must be odd."
  //     // 标准差检查
  //     if(sigma <= 0) throw "sigma 必须大于 0，当前值: " + string(sigma)
      
  //     half_window = int((window - 1) \ 2)
  //     n = size(x)
      
  //     kernel = exp(-pow(seq(-half_window, half_window), 2) / (2 * pow(sigma, 2)))
  //     kernel = kernel \ sum(kernel)
      
  //     x_extended = take(0.0, half_window) join x join take(0.0, half_window)
      
  //     // 缺失值以 0 填充参与卷积，边界按原逻辑返回可能含 NULL 的序列
  //     return mavg(x_extended.nullFill(0), kernel)[(window - 1): (window - 1 + n)]
  // }
  
  def gaussianFilter(x, window, sigma) {
      // 窗口奇数检查
      if(window % 2 == 0) throw "Window size must be odd."
      // 标准差检查
      if(sigma <= 0) throw "sigma 必须大于 0，当前值: " + string(sigma)
      
      half_window = int((window - 1) \ 2)
      n = size(x)
      
      kernel = exp(-pow(seq(-half_window, half_window), 2) / (2 * pow(sigma, 2)))
      kernel = kernel \ sum(kernel)
      
      x_extended = take(0.0, half_window) join x join take(0.0, half_window)
      
      // 缺失值以 0 填充参与卷积，边界按原逻辑返回可能含 NULL 的序列
      return mavg(x_extended.nullFill(0), kernel)[(window - 1): (window - 1 + n)]
  }

卡尔曼滤波 kalmanFilter

• 算法含义

  线性离散卡尔曼滤波是一种递归的状态估计算法，在存在过程噪声和测量噪声的情况下，通过预测-更新步骤，以最小均方误差估计系统状态。适用于多传感器融合与目标追踪。

• 公式

卡尔曼滤波基于如下两个方程：

状态方程（预测）：x_k = Fx_k-1 + Bu_k + w_k-1, w ∼ N(0, Q)

观测方程（更新）：z_k = Hx_k + v_k, v ∼ N(0, R)

其中，卡尔曼滤波的五个核心公式，分别是：

状态先验估计：x̂^-_k = Fx̂_k-1

先验误差协方差：P^-_k = FP_k-1F^T + Q

卡尔曼增益计算：K_k = P^-_kH^T(HP^-_kH^T + R)^-1

状态后验估计（更新）：x̂_k = x̂^-_k + K_k(z_k − Hx̂^-_k)

后验误差协方差（更新）：P_k = (I − K_kH) P^-_k

其中：

x_k：系统状态

z_k：测量值

F：状态转移矩阵

H：观测矩阵

Q：过程噪声协方差

R：测量噪声协方差

w,v：过程噪声和测量噪声

• 代码实现

  def kalmanFilter(z, F, H, Q, R, P0, x0){
      if(form(z) == 3){ numSteps = z.rows() } else { numSteps = size(z) }
      stateDim = size(x0)
      x_est = matrix(DOUBLE, numSteps, stateDim)
      x_curr = x0
      P_curr = P0
      I = matrix(DOUBLE, stateDim, stateDim, 0)
      for(i in 0:stateDim) I[i,i] = 1.0
      for(t in 0:numSteps){
          x_pred = F ** x_curr
          F_trans = transpose(F)
          P_pred = (F ** P_curr ** F_trans) + Q
          if(form(z) == 3){ zt = z[t,:] } else { zt = z[t] }
          y = flatten(zt) - flatten(H ** x_pred)
          H_trans = transpose(H)
          S = (H ** P_pred ** H_trans) + R
          try {
              K = P_pred ** H_trans ** inv(S)
          } catch(ex) {
              K = matrix(DOUBLE, stateDim, size(y), 0.0)
          }
          x_curr = x_pred + (K ** y)
          P_curr = (I - (K ** H)) ** P_pred
          x_est[t,:] = x_curr
      }
      return x_est
  }

• 硬件环境

  CPU：12th Gen Intel(R) Core(TM) i5-12500 (3.00 GHz)

  操作系统： Windows 11 专业版 64 位

  内存：16.0 GB

• 软件环境

  Python：3.9.12（内置 pandas、numpy、scipy、dolphindb 库）

  DolphinDB：V3.00.4.1 单机版（非JIT版本），v3.00.5_JIT 单机版（JIT版本）

本次测试在两种不同实现方式的算法上使用了 1 千万条数据（除 Hampel 算法为 1 百万条以及 Kalman 算法为 1 千条外），其余参数保持一致，其中 amplitude 在 DolphinDB 的普通版本和 JIT 版本上分别进行了测试，并对各自的滤波算法实现分别在随机生成的数据上且保持其它参数一致的情况下计算 10 次，最终测试结果以 10 次计算后所消耗时间的平均值为计量，如下表所示：

通过测试可以看出：

1. 在 8/10 的算法中，DolphinDB 性能显著优于 Python，最高提升达 65 倍（Gaussian），平均提升约 12 倍（剔除 iir 后）；

2. DolphinDB 内置纯向量化函数的实现以及在增量计算等方面的优化，与常规 Python 的实现相比有极大的性能优势；

3. DolphinDB 支持即时编译（JIT）功能，可以优化路径依赖的循环实现，性能达到与静态编译语言相近的执行效率；

4. DolphinDB 的脚本实现的滤波函数运行性能与 Python 的 scipy 库性能相当，接近 C++ 实现的性能极致。

在使用滤波功能前，需要预先将模块下载并存放在指定位置（一般是 DolphinDB 对应的 modules 文件夹下），在脚本中按如下方式引用对应的模块文件：

// 加载模块 (请根据实际部署路径修改)
use IotModules::Filtering::BasicStats        // 基础统计 (均值、中值、Hampel等)
use IotModules::Filtering::ClassicSignal     // 信号处理 (FIR, IIR, SG等)
use IotModules::Filtering::OptimalEstimation // 状态估计 (Kalman)

有关 DolphinDB 模块的教程请参考：模块概述 。

场景一：传感器数据清洗 (去除毛刺与离群点)

问题：温度传感器偶尔会输出极大或极小的异常值（如100度突变为0度或1000度），影响后续计算。

方案：使用 hampelFilter，它基于中位绝对偏差 (MAD)，比标准差法更稳健。

业务价值：故障预警——从“事后分析”转变为“毫秒级预警”。在设备温度出现异常波动的瞬间，而不是第二天出报表时，就能触发停机保护。这意味着过去需要长时间的离线异常检测，现在可以实现实时完成。

// 模拟数据：包含噪声和两个明显的离群点
temp_raw = 25.0 + rand(1.0, 100)
temp_raw[10] = 100.0  // 异常高值
temp_raw[50] = -50.0  // 异常低值
​
// 调用 Hampel 滤波，窗口=5，阈值=3
temp_clean = IotModules::Filtering::BasicStats::hampelFilter(temp_raw, window=5, k=3)
​
// 验证结果
print("Outlier at 10: " + temp_raw[10] + " -> " + temp_clean[10])
// 输出: Outlier at 10: 100 -> 25.4 (近似局部中位数)

场景二：GPS 轨迹平滑 (卡尔曼滤波)

问题：GPS 定位存在几米到几十米的随机误差，导致轨迹呈锯齿状。

方案：使用 kalmanFilter 融合位置观测值与运动模型。

业务价值：实时轨迹修正——对于物流车辆，能够实时修正 GPS 漂移，提供精准的实时位置服务，优化调度效率。

// 1. 定义状态空间模型 (CV模型: Constant Velocity)
// 状态 x = [position, velocity]'
​
F = matrix([1.0, 1.0], [2.0, 1.0]) // 状态转移矩阵
H = matrix([1.0, 1.0], [1.0, 2.0]) // 观测矩阵 (只观测位置)
Q = matrix([0.1, 2.0], [1.0, 0.1]) // 过程噪声
R = matrix([2.0,1], [2.0,1])       // 观测噪声 (GPS误差较大)
P0 = matrix([10.0, 2.0], [2.0, 10.0])
x0 = [0.0, 0.0]
​
// 2. 模拟观测数据 (带噪声的直线运动)
true_pos = 0..99
noise = randNormal(0, 1.5, 100)
z_obs = true_pos + noise
​
// 3. 执行滤波
res = IotModules::Filtering::OptimalEstimation::kalmanFilter(z_obs, F, H, Q, R, P0, x0)
​
// res 是 n x 2 的矩阵，第一列为平滑后的位置，第二列为速度估计
smooth_pos = res[:, 0]

前两节展示了批量历史数据的滤波处理。在生产物联网场景中，传感器数据持续实时产生，需要流式处理框架对数据逐条（或微批）进行滤波。DolphinDB 提供了完善的流计算框架，本节演示将滤波模块集成到流处理中的典型方式。

场景：温度传感器持续上报数据，需实时使用 hampelFilter 清洗离群点。

方式：基于订阅回调的实时异常值清洗

实现流程：

初始化并共享两张流表，分别用于接收温度数据和存储清洗结果

share streamTable(1:0, `time`deviceID`temperature, [TIMESTAMP, SYMBOL, DOUBLE]) as sensorStream
share streamTable(10000:0, `time`deviceID`temperature`cleaned, [TIMESTAMP, SYMBOL, DOUBLE, DOUBLE]) as cleanResult

定义订阅回调函数

def cleanSensorData(msg){
    cleaned = IotModules::Filtering::BasicStats::hampelFilter(msg.temperature, 5, 2.0)
    t = table(msg.time as time, msg.deviceID as deviceID, msg.temperature as temperature, cleaned as cleaned)
    objByName("cleanResult").append!(t)
}

订阅流数据表

subscribeTable(tableName="sensorStream", actionName="cleanData", offset=-1, handler=cleanSensorData, msgAsTable=true, batchSize=50, throttle=1.0)

模拟数据输入

n = 100
timev = now() + 1..n * 1000
deviceIDs = take(`sensor01, n)
temp = 25.0 + rand(1.0, n)
temp[20] = 100.0   // 异常高值
temp[60] = -30.0   // 异常低值
insert into sensorStream values(timev, deviceIDs, temp)

查看清洗结果

select top 5 * from cleanResult where time >= timestamp(timev[18])

预期输出：