CN108801131B - 北斗高频变形监测数据的滤波方法及系统 - Google Patents

Abstract

本发明公开了北斗高频变形监测数据的滤波方法及系统，包括步骤：S100利用样本数据初始化量测扩增卡尔曼滤波法的状态方程和量测方程；S200利用双重粗差判别统计量对待滤波观测值进行粗差类型判断，根据粗差类型进行粗差修正；S300对粗差修正后的待滤波观测值进行噪声滤波，获得二维状态向量滤波值；S400根据二维状态向量滤波值计算观测值滤波值，将观测值滤波值反馈到S200，作为上一历元观测值的滤波值；S500重复S200‑S400，直至全部历元的待滤波观测值滤波完毕。本发明采用双重粗差判别统计量实现粗差的有效识别，依据粗差类型进行个性化修正，极大的保留了监测数据中的有效信息，从而保证滤波结果的可靠性。

Description

北斗高频变形监测数据的滤波方法及系统

技术领域

本发明属于大地测量技术领域，特别是一种北斗高频变形监测数据的滤波方法及系统。

背景技术

全球卫星导航系统GNSS(Global Navigation Satellite System)具有速度快、全天候、自动化程度高等特点，已经广泛应用于各类工程的变形监测，如大坝变形监测、桥梁变形监测、滑坡变形监测以及矿区变形监测等。我国自主研发的北斗定位系统作为全球定位系统的重要组成部分，正以快速的发展速度，逐步实现全球范围内的导航定位，将会在国内外导航定位应用领域同其他的卫星定位系统产生更激烈的竞争。但是，作为一种新型的卫星导航定位系统，其在高频变形监测数据的误差特征和滤波方法仍然是空白。北斗卫星定位系统监测结果的准确性，与众多内外因素相关，如卫星端的卫星轨道信息与卫星钟、接收机端的接收机天线相位中心偏差位置与接收机钟、传播过程中的大气中的电离层、对流层以及工程环境等，它们都是导致北斗卫星定位系统的高频变形监测数据产生多种误差的原因。针对这些因素的影响，目前主要从北斗定位系统基线高精度解算的角度提出相应的处理算法；但工程环境的复杂性和不确定性，使得基线解算处理算法难以顾及，导致实时解算的北斗高频变形监测数据序列被粗差和噪声信息所污染，大大降低了利用北斗卫星定位系统进行工程变形监测的可靠性和准确性。因此，针对实时解算得到的北斗高频变形监测数据，研究相应的数据滤波方法是非常必要的。

卫星导航定位系统的变形监测数据的滤波方法研究，长期以来一直是卫星导航定位系统在变形监测领域的热点。自R.E.Kalman于1960年提出卡尔曼滤波以来，其作为一种动态数据滤波方法被广泛应用于卫星导航定位系统的变形监测数据处理中，从而在一定程度上提高了卫星导航定位系统的可靠性。但标准的卡尔曼滤波是以状态噪声和量测噪声为互不相关的白噪声为前提，其对于量测数据中的有色噪声和粗差处理能力有限。对于量测数据中的有色噪声问题，学者们通常以一阶AR模型作为有色噪声的模型，对标准卡尔曼滤波进行改进，形成量测扩增卡尔曼滤波方法，从而实现兼顾有色观测噪声的卡尔曼滤波；对于量测数据的粗差问题，学者们利用稳健估计等价权原理，在对卡尔曼滤波环节中增益矩阵评价的基础上，通过选择合适的权函数来代替量测噪声的协方差阵，形成抗差卡尔曼滤波方法，进而达到对量测数据中粗差值的有效解决。然而，北斗卫星高频变形监测数据中会同时出现有色噪声和粗差的情况，而上述方法只能单一削弱其中一个方面的误差。由此，提出一个能够同时消除北斗高频变形监测数据中的多样性噪声和多样性粗差的方法，从而切实提高北斗卫星定位系统变形监测结果的可靠性，显得尤为迫切。

根据对北斗高频变形监测数据的滤波方法研究现状分析，可以发现：虽然国内外已有众多学者分别针对北斗高频变形监测数据中的噪声和粗差提出了很多滤波方法，但是由于没有同时考虑到北斗卫星高频变形监测数据中的噪声和粗差的双重影响，导致这些滤波算法所得到的滤波结果准确性有限。

发明内容

本发明的目的是提供一种考虑到噪声和粗差的双重影响的、北斗高频变形监测数据的滤波方法及系统。

本发明是为弥补传统北斗高频变形监测数据滤波方法中对误差信息处理的片面性，导致无法准确地消除北斗高频变形监测数据中误差。本发明提出了一种改进量测扩增卡尔曼滤波法，并将其用于北斗高频变形监测数据的滤波，可同时滤除北斗高频变形监测数据中的噪声和粗差，从而有效提高北斗卫星定位系统在变形监测应用中的精度。

本发明北斗高频变形监测数据的滤波方法，包括步骤：

S100利用样本数据初始化量测扩增卡尔曼滤波法的状态方程和量测方程，所述样本数据为预先采集的北斗高频变形监测数据序列；

S200对待滤波观测值进行粗差修正，本步骤具体为：

利用双重粗差判别统计量δ₁和δ₂对待滤波观测值进行粗差类型判断，根据粗差类型，利用公式

进行粗差修正；

其中，双重粗差判别统计量

Z(k+1)表示待滤波观测值；Z′(k+1)表示待滤波观测值的粗差修正值；

表示待滤波观测值的上一历元观测值的滤波值；r(k)和r(k-1)分别为待滤波观测值所在历元和上一历元对应的一阶差分值；

S300对粗差修正后的待滤波观测值进行噪声滤波，获得二维状态向量滤波值；

S400根据二维状态向量滤波值计算观测值滤波值，将观测值滤波值反馈到步骤S200，作为上一历元观测值的滤波值；

S500重复S200-S400，直至全部历元的待滤波观测值滤波完毕。

进一步的，步骤S100具体为：

利用样本数据估算状态方程和量测方程中二维状态向量的初值

及其对应的误差协方差矩阵初值P(0)、状态噪声的协方差矩阵初值Q(0)、量测噪声中高斯白噪声的协方差矩阵初值R(0)、以及量测噪声相关系数矩阵初值ρ(0)；

估算样本数据一阶差分数据序列的平均值m和标准差σ′。

所述估算样本数据一阶差分数据序列的平均值m和标准差σ′，具体为：

首先，利用四分位间距准则对样本数据的一阶差分数据序列进行粗差剔除，即剔除一阶差分数据序列中所对应统计量绝对值|U|＞3的一阶差分数据；

接着，计算粗差剔除后的一阶差分数据序列的平均值m和标准差σ′。

本发明北斗高频变形监测数据的滤波系统，包括：

第一模块，用来利用样本数据初始化量测扩增卡尔曼滤波法的状态方程和量测方程，所述样本数据为预先采集的北斗高频变形监测数据序列；

第二模块，用来对待滤波观测值进行粗差修正；

具体为：用来利用双重粗差判别统计量δ₁和δ₂对待滤波观测值进行粗差类型判断，根据粗差类型，利用公式

进行粗差修正；

其中，双重粗差判别统计量Z(k+1)表示待滤波观测值；Z′(k+1)表示待滤波观测值的粗差修正值；

表示待滤波观测值的上一历元观测值的滤波值；r(k)和r(k-1)分别为待滤波观测值所在历元和上一历元对应的一阶差分值；

第三模块，用来对粗差修正后的待滤波观测值进行噪声滤波，获得二维状态向量滤波值；

第四模块，用来根据二维状态向量滤波值计算观测值滤波值，将观测值滤波值反馈到第二模块，作为上一历元观测值的滤波值。

和现有技术相比，本发明具有如下优点和有益效果：

(1)相比传统的北斗高频变形监测数据的滤波方法，本发明不仅可消除噪声信息的干扰，还可消除粗差信息的干扰，使得滤波后的变形监测数据精度得到显著提升。

(2)针对北斗高频变形监测数据中粗差的多样性特点，例如离散型和区域型粗差，本发明采用双重粗差判别式实现粗差的有效识别，并在修正方法上依据粗差类型特点进行个性化修正，极大的保留了原始北斗高频变形监测数据中的有效信息，从而保证了滤波结果的可靠性。

附图说明

图1是本发明滤波方法的具体流程图；

图2是实施例中原始信号S1的示意图；

图3是实施例中滤波信号S2的示意图；

图4是实施例中采用传统滤波方法处理后的滤波信号S3示意图。

图5是实施例中真实参照信号S4的示意图。

具体实施方式

为了更清楚地说明本发明和/或现有技术中的技术方案，下面将对照附图说明本发明的具体实施方式。显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图，并获得其他的实施方式。

应当理解，此处所描述的实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的实施例中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

实施例

本实施例选择了一组包含600个数据点的北斗变形监测数据序列作为样本数据，样本数据为中国湖北省武汉市白沙洲长江大桥于2016年建立的北斗卫星定位连续监测网中监测点BD43，于2017年8月22日12时0分至12时10分采集的北斗高频变形监测数据序列。该监测网共包含60个监测点，每个监测点每隔1s采集一次三维坐标。由于变形监测环境的复杂性和不确定性，使得实时解算的桥梁北斗高频变形监测数据被噪声和粗差信息所污染，导致桥梁管理机构无法实时准确掌握桥梁的安全状态。图2给出的是监测点BD43于2017年8月22日12时0分至12时10分采集的Z方向的位移变形监测数据，即原始信号S1。显然，原始信号S1包含着严重的误差，以至于变形监测点的变形趋势无法准确确定。

参见图1，本实施例中滤波方法的具体步骤如下：

步骤1，利用样本数据估算量测扩增卡尔曼滤波的初始参数，所述样本数据即预先采集的北斗高频变形监测数据序列。

本步骤的具体实施过程如下：

利用量测扩增卡尔曼滤波法的位置模型和速度模型，构建目标对应的状态方程和量测方程，所述目标指监测对象，本实施例中，目标即白沙洲长江大桥。所构建的状态方程和量测方程见式(1)：

式(1)中：

X(k)和X(k+1)分别为目标在时刻k和k+1的二维状态向量；

φ(k)为状态转移矩阵，取值为

G(k)为状态噪声系数，取值为

w(k)为时刻k的状态噪声；

Z(k+1)为目标在时刻k+1的观测值；

H(k+1)为量测转移矩阵，取值为(1 0)；

v(k)和v(k+1)分别为目标在时刻k和k+1的量测噪声；

ρ(k)为时刻k+1与时刻k的量测噪声相关系数矩阵；

ξ(k)为时刻k的量测噪声中的高斯白噪声。

所述的初始参数包括二维状态向量的初值

及其对应的误差协方差矩阵初值P(0)、状态噪声的协方差矩阵初值Q(0)、量测噪声中高斯白噪声的协方差矩阵初值R(0)、量测噪声相关系数矩阵初值ρ(0)、以及样本数据一阶差分数据序列的平均值m和标准差σ′。一阶差分数据序列的平均值m和标准差σ′用于粗差处理环节。

本实施例中，根据样本数据前5个北斗变形监测数据序列估计目标二维状态向量的初值

其中，位置状态量初值取前5个北斗变形监测数据序列中目标的位置均值，为-0.0033m；速度状态量初值为0.0014m/s；故

对应的误差协方差矩阵P(0)取样本数据第一组观测量所对应的方差阵，即

量测噪声v(k+1)对应的相关系数矩阵ρ(k)和高斯白噪声协方差矩阵R(k)，可通过尤尔-沃克(Yule-Walker)方程对其初值进行确值，本实施例中，ρ(0)＝0.0264，R(0)＝0.00001348m²。

状态噪声w(k)的协方差矩阵是在前面取值的基础上，通过经验进行试算确定，本实施例中，状态噪声协方差矩阵初值Q(0)＝0.00000337m²。

一阶差分数据的平均值m和标准差σ′的确值方法如下：

(1)利用四分位间距(IQR，inter-quartile range)准则，对样本数据的一阶差分数据序列进行粗差剔除，以获取更加接近真实情况的平均值m和标准差σ′。

更具体的，采用公式(2)进行粗差剔除：

式(2)中：

dZ(i)表示一阶差分数据序列中第i个一阶差分数据；

median(dZ(i))表示一阶差分数据序列的中位数；

i表示一阶差分数据序列中一阶差分数据的编号，i＝1,2，......P，P为一阶差分数据序列长度；

s表示标准化处理后的IQR，即一阶差分数据变异性的量度，本实施例中，s取0.7413IQR。

剔除一阶差分数据序列中统计量|U|＞3的一阶差分数据，即完成粗差剔除。

(2)计算粗差剔除后的一阶差分数据序列的平均值m和标准差σ′。

平均值m和标准差σ′的计算公式见式(3)-(4)：

式(3)-(4)中：

dZ＇(i)表示粗差剔除后的一阶差分数据序列中第i个一阶差分数据；

Q为粗差剔除后的一阶差分数据序列长度，i＝1,2，......Q。

本实施例中，计算得到m＝0.000080898m，σ′＝0.0049m。

步骤2，对待滤波观测值进行粗差识别与修正。

实时滤波中，待滤波观测值为当前历元的观测值；事后滤波中，待滤波观测值则为选定的观测值。

本步骤进一步包括：

(1)利用一阶差分数据序列离散度稳定性构建双重粗差判别统计量δ₁和δ₂，实现离散型粗差和区域型粗差的同时识别。

所构建的双重粗差判别统计量δ₁和δ₂如下：

式(5)中：

Z(k)和Z(k+1)分别为时刻k和k+1的观测值，其中，Z(k+1)为当前历元的观测值；

为时刻k观测值的滤波值，即前一历元观测值的滤波值，的初值

m为一阶差分数据序列的平均值，由步骤1获得；

σ′为一阶差分数据序列的标准差，由步骤1获得。

(2)根据待滤波观测值的粗差判别结果进行粗差修正。

当δ₁＞3且δ₂＞3时，则判别为离散型粗差，对离散型粗差的观测值，修正方法是利用上一时刻的一阶差分值对当前一阶差分值进行替换。当δ₁＞3且δ₂≤3，则判别为区域型粗差，对区域性粗差的观测值，观测值间的一阶差分的内在联系仍然存在，故修正方式是保留当前一阶差分值。

粗差判别及粗差修正见式(6)：

式(6)中：

Z(k+1)为时刻k+1的观测值；

Z′(k+1)为时刻k+1观测值的粗差修正值；

为时刻k观测值的滤波值；

r(k)和r(k-1)分别为时刻k+1和k对应的一阶差分值。

步骤3，对粗差修正后的待滤波观测值进行噪声滤波，获得二维状态向量滤波值。

本步骤可采用量测扩增卡尔曼滤波法中已有的噪声滤波模型，为便于理解，下面将对噪声滤波模型及其滤波过程进行详细说明。

本步骤进一步包括：

(1)计算标准化后的量测转移矩阵H^*(k)、准化后的高斯白噪声协方差矩阵R^*(k)以及状态噪声和量测噪声间的协方差矩阵S(k)，如下：

式(7)中：

H(k)和H(k+1)分别为时刻k和k+1下标准化前的量测转移矩阵；

φ(k)为状态转移矩阵；

ρ(k)为时刻k+1与时刻k的量测噪声相关系数矩阵；

R(k)为时刻k下标准化前的高斯白噪声协方差矩阵；

G(k)为状态噪声系数；

Q(k)为状态噪声的协方差矩阵；

S(k)为状态噪声和量测噪声间的协方差矩阵。

(2)计算标准化后的卡尔曼增益、状态向量滤波值及其对应的误差协方差矩阵，如下：

式(8)中：

为时刻k+1下标准化后的卡尔曼增益；

Z′(k)和Z′(k+1)分别为时刻k和k+1的观测值的粗差修正值；

和

分别为时刻k和时刻k+1的二维状态向量滤波值；

P(k)和P(k+1)分别为二维状态向量时刻k和时刻k+1的误差协方差矩阵。

利用模型(7)-(8)对粗差修正后的待滤波观测值进行噪声滤波，获得当前历元的二维状态向量滤波值

步骤4，根据步骤3输出的二维状态向量滤波值，计算最终滤波值，即观测值滤波值。

本步骤具体包括：

(1)利用公式(9)计算最终滤波值

式(9)中：

H(k+1)为时刻k+1下标准化前的量测转移矩阵；

为时刻k+1的二维状态向量滤波值。

(2)将最终滤波值

反馈到步骤2，作为上一时刻观测值的滤波值，继续下一历元观测值的粗差修正和噪声滤除，直至全部历元观测值处理完毕。

本实施例滤波结果的对比分析：

采用本发明滤波方法对原始信号S1进行滤波，主要包括四个步骤：(1)设置初始参数；(2)读取最新观测值，进行粗差的识别与修正；(3)对粗差处理后的观测值进行噪声滤波；(4)输出当前历元观测值的最终滤波值，将最终滤波值反馈至步骤(2)，并进行下一轮滤波，直到全部历元观测值处理完毕结束滤波。

循环600次后结束，得到原始信号S1相应的滤波信号S2，见图3。通过对比图2和图3，可以明显看出，滤波信号S2在粗差和噪声滤除上都具有很好的效果。对原始信号S1中24、133、148、239、285、357、358、431、443、473、507、569等位置的离散型粗差，本发明滤波方法可以利用前一差分值修正实现平稳过渡；对于97-114等位置的区域态粗差，本发明滤波方法保留了其差分数据的相关性，也实现了平稳过渡。对于原始信号S1中的噪声，本发明滤波方法能够有效对其进行滤除。相比原始信号S1，滤波信号S2变形信息更加直观，有效地提高了北斗高频变形监测数据的精度和可靠性。

为了更加全面的证明本发明滤波方法的优越性，选择一传统滤波方法的滤波结果进行对比。由于本发明滤波方法是基于量测扩增卡尔曼滤波方法，并对其进行改进所形成，为了使对比结果更加具有说服力，本实施例选择传统的量测扩增卡尔曼滤波结果进行对比。图4是运用传统的量测扩增卡尔曼滤波方法滤波后得到的滤波信号S3。通过对比图3和图4，可以毫无疑义的看出，滤波信号S3的滤波效果不如滤波信号S2。传统滤波方法的劣势主要表现在粗差值的滤波效果，其对于离散型的粗差值具有一定的抵抗作用；但对于区域型粗差，其抵抗能力显然不足。但在滤除北斗高频变形监测数据中误差的过程中，同时考虑噪声和粗差的存在是非常必要的。

上述仅提供了对原始信号S1的直接滤波效果的定性对比分析，下面将通过滤波前后的均方根误差RMSE(root mean squared error)和信噪比SNR(signal-to-noise ratio)来进行定量对比分析。由于缺乏相对真实信号的对照，本实施例中把奇异谱分析(singularspectrum analysis,SSA)引入到数据事后处理中，利用奇异谱分析结果作为真实参照信号S4，见图5所示。以此分别计算不同滤波方法滤波前后的均方根误差和信噪比，对比结果见表1所示。

表1 信号滤波前后均方根误差和信噪比的对比

通过表1可以看出，原始信号S1的均方根误差为0.0066m，信噪比为6.16dB。经过传统滤波方法和本发明滤波方法处理后，均在一定程度上降低了均方根误差并提高了信噪比，其中：传统滤波方法可将滤波信号均方根误差降低至0.0039m，将信噪比提升至8.9573dB，滤波后信号质量虽有提升，但提升程度有限。本发明滤波方法可将滤波信号均方根误差降至0.0022m，精度提升约3倍，信噪比提升至16.1196dB，信噪比提升约2倍。通过均方根误差和信噪比的定量指标的分析与比较，进一步表明本发明滤波方法的有效性和优越性。

以上给出了本发明涉及的具体实施方式，但本发明不局限于所描述的实施方式。在本发明给出的思路下，采用对本领域技术人员而言容易想到的方式对上述实施例中的技术手段进行变换、替换、修改，并且起到的作用与本发明中的相应技术手段基本相同、实现的发明目的也基本相同，这样形成的技术方案是对上述实施例进行微调形成的，这种技术方案仍落入本发明的保护范围内。

Claims (4)

1.北斗高频变形监测数据的滤波方法，其特征是，包括：

S100利用样本数据初始化量测扩增卡尔曼滤波法的状态方程和量测方程，所述样本数据为预先采集的北斗高频变形监测数据序列；

S200对待滤波观测值进行粗差修正，本步骤具体为：

利用双重粗差判别统计量δ₁和δ₂对待滤波观测值进行粗差类型判断，根据粗差类型，利用公式

进行粗差修正；

其中，双重粗差判别统计量Z(k+1)表示待滤波观测值；Z′(k+1)表示待滤波观测值的粗差修正值；

表示待滤波观测值的上一历元观测值的滤波值；r(k)和r(k-1)分别为待滤波观测值所在历元和上一历元对应的一阶差分值；m和σ′为样本数据一阶差分数据序列的平均值和标准差；

S300对粗差修正后的待滤波观测值进行噪声滤波，获得二维状态向量滤波值；

S400根据二维状态向量滤波值计算观测值滤波值，将观测值滤波值反馈到步骤S200，作为上一历元观测值的滤波值；

S500重复S200-S400，直至全部历元的待滤波观测值滤波完毕。

2.如权利要求1所述的北斗高频变形监测数据的滤波方法，其特征是：

步骤S100具体为：

利用样本数据估算状态方程和量测方程中二维状态向量的初值

及其对应的误差协方差矩阵初值P(0)、状态噪声的协方差矩阵初值Q(0)、量测噪声中高斯白噪声的协方差矩阵初值R(0)、以及量测噪声相关系数矩阵初值ρ(0)；

估算样本数据一阶差分数据序列的平均值m和标准差σ′。

3.如权利要求2所述的北斗高频变形监测数据的滤波方法，其特征是：

所述估算样本数据一阶差分数据序列的平均值m和标准差σ′，具体为：

首先，利用四分位间距准则对样本数据的一阶差分数据序列进行粗差剔除，即剔除一阶差分数据序列中所对应统计量绝对值|U|＞3的一阶差分数据；

接着，计算粗差剔除后的一阶差分数据序列的平均值m和标准差σ′。

4.北斗高频变形监测数据的滤波系统，其特征是，包括：

第一模块，用来利用样本数据初始化量测扩增卡尔曼滤波法的状态方程和量测方程，所述样本数据为预先采集的北斗高频变形监测数据序列；

第二模块，用来对待滤波观测值进行粗差修正；

具体为：用来利用双重粗差判别统计量δ₁和δ₂对待滤波观测值进行粗差类型判断，根据粗差类型，利用公式

进行粗差修正；

其中，双重粗差判别统计量

Z(k+1)表示待滤波观测值；Z′(k+1)表示待滤波观测值的粗差修正值；表示待滤波观测值的上一历元观测值的滤波值；r(k)和r(k-1)分别为待滤波观测值所在历元和上一历元对应的一阶差分值；m和σ′为样本数据一阶差分数据序列的平均值和标准差；

第三模块，用来对粗差修正后的待滤波观测值进行噪声滤波，获得二维状态向量滤波值；

第四模块，用来根据二维状态向量滤波值计算观测值滤波值，将观测值滤波值反馈到第二模块，作为上一历元观测值的滤波值。

Images