3、数据预处理:去噪滤波算法与数据质量控制

激光雷达原始回波信号中不可避免地混杂着大气湍流、硬件噪声及背景光干扰。数据预处理的核心目标是在保留有效风场信息的前提下,最大限度抑制噪声,并通过质量控制剔除不可靠数据。本章将重点讲解两种经典滤波算法与两项关键质量控制手段。

3.1 去噪滤波算法

3.1.1 中值滤波

原理:中值滤波是一种非线性平滑技术,通过将每个数据点替换为其邻域窗口内所有点的中位数来消除孤立噪声点。对于激光雷达径向速度谱或风廓线中的“椒盐噪声”(如飞点、瞬时强干扰),中值滤波效果显著,且能较好地保留信号边缘(如风切变界面)。

关键参数:窗口长度 L(通常取奇数,如3、5、7)。窗口越大,平滑越强,但可能模糊细节。

Python实现示例(一维径向速度序列)

import numpy as np
from scipy.signal import medfilt

def median_filter_velocity(v_r, kernel_size=5):
    """
    对径向速度序列进行中值滤波
    :param v_r: 一维数组,径向速度 (m/s)
    :param kernel_size: 窗口大小,必须为奇数
    :return: 滤波后的速度序列
    """
    # 确保kernel_size为奇数
    if kernel_size % 2 == 0:
        kernel_size += 1
    v_filtered = medfilt(v_r, kernel_size=kernel_size)
    return v_filtered

3.1.2 小波去噪

原理:小波变换将信号分解到不同频率尺度上。噪声通常集中在高频细节分量中,而风场的大尺度结构(如平均风、低频波动)分布在低频近似分量中。通过阈值处理高频小波系数,再重构信号,即可实现去噪。

常用步骤

  1. 选择小波基(如 db4sym8)和分解层数(通常3~5层)。
  2. 对信号进行离散小波变换(DWT)。
  3. 对每一层的高频细节系数应用软阈值或硬阈值函数。
  4. 进行逆小波变换(IDWT)重构信号。

Python实现示例

import pywt
import numpy as np

def wavelet_denoise(signal, wavelet='db4', level=4, mode='soft'):
    """
    小波去噪
    :param signal: 一维信号数组
    :param wavelet: 小波基名称
    :param level: 分解层数
    :param mode: 阈值模式 'soft' 或 'hard'
    :return: 去噪后的信号
    """
    # 分解
    coeffs = pywt.wavedec(signal, wavelet, level=level)
    # 估计噪声标准差(使用第一层细节系数)
    sigma = np.median(np.abs(coeffs[-1])) / 0.6745
    # 计算通用阈值
    threshold = sigma * np.sqrt(2 * np.log(len(signal)))
    # 对每一层细节系数进行阈值处理(保留近似系数不变)
    coeffs_thresh = [coeffs[0]]  # 近似系数不处理
    for i in range(1, len(coeffs)):
        coeffs_thresh.append(pywt.threshold(coeffs[i], threshold, mode=mode))
    # 重构
    signal_denoised = pywt.waverec(coeffs_thresh, wavelet)
    return signal_denoised

对比建议:中值滤波计算快,适合实时处理;小波去噪对非平稳噪声(如湍流引起的随机波动)效果更优,但需谨慎选择小波基与阈值。

3.2 数据质量控制

3.2.1 信噪比阈值

定义:信噪比(SNR)是信号功率与噪声功率的比值,通常以dB表示。激光雷达每个距离门的SNR直接反映了该处回波的可信度。

阈值设定原则

  • 低SNR(如 < -5 dB):数据受噪声主导,应直接剔除。
  • 中等SNR(-5 ~ 0 dB):可保留但需结合其他指标(如谱宽)综合判断。
  • 高SNR(> 0 dB):数据可靠,优先采用。

实际应用:在风廓线反演前,对每个距离门的径向速度进行SNR掩膜:

def snr_mask(v_r, snr, snr_threshold=-3.0):
    """
    根据SNR阈值掩膜无效数据
    :param v_r: 径向速度数组
    :param snr: 对应的信噪比数组 (dB)
    :param snr_threshold: 信噪比阈值 (dB)
    :return: 掩膜后的速度数组(无效点设为NaN)
    """
    v_masked = np.where(snr >= snr_threshold, v_r, np.nan)
    return v_masked

3.2.2 一致性检验

目的:检测并剔除因硬件异常、强湍流或硬目标干扰导致的“野点”。

常用方法

  1. 时间连续性检验:同一距离门在相邻扫描周期内的速度变化不应超过物理极限(如最大风速变化率)。

    def temporal_consistency(v_series, max_delta=5.0):
        """
        时间一致性检验:相邻时刻速度差超过max_delta则标记为无效
        :param v_series: 同一距离门的时间序列
        :param max_delta: 最大允许变化 (m/s)
        :return: 布尔掩膜,True表示有效
        """
        diff = np.abs(np.diff(v_series, prepend=v_series[0]))
        mask = diff <= max_delta
        return mask
    
  2. 空间一致性检验:相邻距离门之间的速度梯度不应过大(如超过20 m/s/km)。

    def spatial_consistency(v_profile, range_gate_length=30.0, max_gradient=0.02):
        """
        空间一致性检验:相邻距离门速度梯度
        :param v_profile: 一维径向速度廓线
        :param range_gate_length: 距离门长度 (m)
        :param max_gradient: 最大允许梯度 (m/s/m)
        :return: 布尔掩膜
        """
        gradient = np.abs(np.gradient(v_profile, range_gate_length))
        mask = gradient <= max_gradient
        return mask
    
  3. 谱宽约束:若激光雷达同时输出谱宽数据,谱宽过大(如 > 4 m/s)通常意味着湍流极强或信号不可靠,应予以剔除。

3.3 预处理流程整合

在实际数据处理中,建议按以下顺序执行:

  1. 原始数据读取:获取径向速度、SNR、谱宽等。
  2. SNR阈值掩膜:剔除低SNR数据。
  3. 去噪滤波:对保留的数据进行中值滤波或小波去噪。
  4. 一致性检验:依次进行时间连续性、空间连续性及谱宽检验。
  5. 插值填补(可选):对少量缺失点进行线性插值,但避免过度平滑。

下表总结了各方法的适用场景:

方法主要作用适用噪声类型计算开销
中值滤波消除孤立野点椒盐噪声、飞点
小波去噪抑制随机噪声高斯噪声、湍流波动
SNR阈值剔除弱信号低信噪比区域极低
一致性检验剔除异常突变硬件故障、强干扰

通过上述预处理,可显著提升后续风场反演(如VAD、DBS)的准确性与鲁棒性。下一章将基于预处理后的数据,讲解水平风场与垂直风廓线的反演算法。