车辆轨迹数据的小波去噪

        小波去噪具有受噪声影响小，相较于卡尔曼滤波、滑动平均等方法而言，对非单一频率、非周期变化的信号具有更好的效果，因而，很适用于对车辆轨迹数据进行去噪。小波去噪基于小波变换实现，小波变换解决了傅里叶变换无法描述时间和短时傅里叶变换时间窗大小严重影响信号分析的缺点，在进行小波变换时，需要明确四个参数：小波基函数、分解层数、阈值、阈值函数。四个参数的选取均会影响到小波去噪的效果，特别是分解层数，分解层数选择不当可能导致去噪不彻底或信号失真。

• 小波基函数

        小波基函数有常用六项指标：支撑长度、对称性、消失矩、正则性、相似性

        常用的小波基函数族有symN、dbN、cofiN等，其中N为小波分解的级数。        

• 正交性

        一个正交小波基含三重定义：同一尺度下小波函数与尺度函数内积为0；不同尺度之间的小波函数内积为0；同一尺度下时间轴上任意两个小波函数内积为0。上述常用的三类小波基函数均有着正交性。

• 支撑长度

        支撑长度太长，会耗费过多的计算时间，产生更多高幅值的小波系数，产生边界问题，支撑长度太短，消失矩太低，不利于信号的集中。大部分应用选择支撑长度为5~9之间的小波。symN小波族和dbN小波族的支撑长度为2N-1，coifN小波族的支撑长度为6N-1。

• 对称性

       具有对称性的小波，因其对应的滤波器具有线型相位，在进行数据处理时可有效避免相位畸变。 dbN小波族中除db1外均不具有对称性，在对信号进行分析和重构时会产生一定的相位失真，symN小波族、coifN小波族比dbN小波族具有更好的对称性，在对信号进行分解和重构时能一定程度上减少信号的相位失真。

• 消失矩

       消失矩用于使尽量多的小波系数为0或产生尽量少的非0小波系数，有利于数据压缩和噪声的消除，消失矩越大，越多的小波系数为0，但是在一般情况下，消失矩越高，支撑长度也越长，需综合考量。symN小波族、dbN小波族的消失矩阶数为N，coifN小波族的消失矩阶数为2N。

• 正则性

        正则性好的小波，能在信号或图像的重构中获取较好的平滑效果，减小量化和舍入误差的视觉影响，但是在一般情况下，正则性越好，支撑长度也越长，计算时间也更长，因此需综合考量。通常情况下，正则性随消失矩的增加而加强。

• 相似性

        在进行数据处理时，使用与信号波形相似的小波基函数，在进行信号的压缩和去噪时具有更好的效果。

• 分解层数

        小波去噪基本思想是对原始信号进行多层小波分解，生成细节系数和近似系数，而后使用阈值函数和阈值规则对系数进行处理，再进行重构，进而得到去噪后信号。小波分解和重构大多使用Mallat算法。Mallat算法于1986年提出，从空间的概念上形象的说明了小波的多分辨率特性，随着尺度由大到小变化，在各尺度上可以由粗到细的观察图像的不同特征的一种算法，算法示意如下图所示：

        小波分解和重构的的层数并不是越大越好，也不是越小越好，根据多分辨率分析理论，高层分解的小波系数对应信号低频部分，低频部分主要由信号构成，因此分解层数越高，去掉的低频成分越多，去噪效果越明显，但信号失真也相应增大，为保守起见，分阶层数不宜太高，通常最大不超过5层，对于那些本身波动性强的序列，其信号本身的高频成分较多，更不能取太高的层次，一般不超过3层。现有针对车辆轨迹数据去噪的研究中，小波去噪的分解层数大多基于经验或直接给定得到。

• 阈值

        常用阈值函数有固定阈值、无偏风险估计阈值、启发式阈值三类。启发式阈值结合了固定阈值和无偏风险估计阈值。

        计算固定阈值的Python代码如下：

        def calculate_sgtwolog(signal):
            '''
            计算固定阈值
            :param signal:信号，即coeffs[i]，coeffs为[cA_n, cD_n, cD_n-1, ..., cD2, cD1]
            :return:sgtwolog阈值，lambda_s，float
            '''
            lambda_s = math.sqrt(2 * math.log(len(signal)))

            return lambda_s

        计算无偏风险估计阈值的Python代码如下：

        def calculate_rigrsure(signal):
            '''
            计算无偏风险估计阈值
            :param signal:信号，即coeffs[i]，coeffs为[cA_n, cD_n, cD_n-1, ..., cD2, cD1]
            :return:rigrsure阈值，lambda_k，float
            '''
            abs_signal = [abs(item) for item in signal]
            abs_signal_order = sorted(abs_signal)
            new_signal = [item_1 ** 2 for item_1 in abs_signal_order]
            lst_risk = []
            N = len(new_signal)
            for i_1 in enumerate(new_signal):
                risk = (N - 2 * (i_1[0] + 1) + sum([i_2 for i_2 in new_signal if i_2 <= i_1[1]]) + (N - (i_1[0] + 1)) * i_1[1]) / (N)
                lst_risk.append(risk)
            index = lst_risk.index(min(lst_risk))
            lambda_k = math.sqrt(new_signal[index])

            return lambda_k

        计算启发式阈值时，首先计算变量Eta和Crit的大小，计算公式见下：

       当Eta<Crit时，选取固定阈值，反之选取无偏风险估计阈值和固定阈值中的较小者。

• 阈值函数

        常用的阈值函数有软阈值函数、硬阈值函数等。其基本思想为将幅度小于阈值的系数都置0，区别在于对大于阈值的系数的处理方式不同。

软阈值函数

        软阈值函数把大于阈值的小波系数中减去阈值，将幅度大于阈值的小波系数进行缩小，小于阈值的小波系数置0，其式见下：

        这样处理会产生一个显著的问题，每个大于阈值的系数都减去了阈值，将改变所有小波系数值，毫无提问会在信号重构时引入误差。但因其特性，软阈值函数获取的重构信号具有更好的光滑性。

硬阈值函数

        硬阈值函数将幅度大于阈值的小波系数保持不变，仅仅将小于阈值的小波系数置0，其式见下：

         硬阈值函数将会导致信号出现跳变、不连续，产生附加震荡，因而在进行信号去噪时，往往使用较少，但硬阈值函数获得的重构信号具有更好的逼近性。

        软阈值函数得到的小波系数整体连续性较好，从而使信号不会产生附加震荡，但是由于会压缩信号，产生一定的偏差，直接影响到重构信号与真实信号的逼近成都，硬阈值函数在均方根误差意义上优于软阈值函数，但是信号会产生附加震荡，产生跳跃点，不具有原始信号的平滑性。

• 代码示例

        以下给出以NGSIM数据集为例的小波去噪Python代码，代码中小波基函数选取为db8，分解层数选取为最大分解层数，阈值为0.2，阈值函数为软阈值函数：

def Denoising_wavelet(origin_data):
    index = []
    data = []
    for i in range(len(origin_data) - 1):
        X = float(i)
        Y = float(origin_data[i])
        index.append(X)
        data.append(Y)

    w = pywt.Wavelet('db8')  # 选用Daubechies8小波
    maxlev = pywt.dwt_max_level(len(data), w.dec_len)
    print("maximum level is " + str(maxlev))

    coeffs = pywt.wavedec(data, 'db8', level=maxlev)  # 将信号进行小波分解

    for i in range(1, len(coeffs)):
        coeffs[i] = pywt.threshold(coeffs[i], 0.2, 'soft', 0)  # 将噪声滤波，默认采用软阈值，替换值为0

    wavelet_data = pywt.waverec(coeffs, 'db8')  # 将信号进行小波重构

    if len(origin_data) == len(wavelet_data):
        return list(wavelet_data)
    if len(origin_data) != len(wavelet_data):
        print('小波处理后长度发生变化'+str(len(origin_data)-len(wavelet_data)))
        mintime = 0
        maxtime = mintime + len(data) + 1
        new_wavelet_data = np.append(wavelet_data[mintime:maxtime - 1],
                                     wavelet_data[mintime:maxtime - 1][-1])  # 小波去噪后，元素少一个，添加一个末尾元素
        return list(new_wavelet_data)

---

 若在小波去噪及车辆轨迹数据领域有相关疑问，或交通数据处理需求欢迎VX探讨：A2528945820