CN103576167A - 基于hht和支持向量机的周跳探测与修复方法 - Google Patents

Abstract

本发明提供了一种基于HHT和支持向量机的周跳探测与修复方法，首先采用EMD算法对码相二次差时间序列进行多时间尺度分解，进而对分解得到的IMF分量作Hilbert变换；其次根据IMF分量的Hilbert幅值是否同时发生突变来探测周跳；然后对Hilbert幅值发生突变的IMF分量分别建立支持向量机预测模型，并预测周跳发生时刻的值；最后比较预测值与实际值的大小来确定周跳大小。本发明可用于单频或多频载波相位观测值的周跳探测与修复，提高了周跳探测的灵敏度以及周跳修复的准确性。

Description

基于HHT和支持向量机的周跳探测与修复方法

技术领域

本发明涉及卫星导航精密定位技术领域，具体涉及一种周跳探测与修复方法。

背景技术

卫星导航定位的原理是，当卫星位置已知时，通过测量四颗以上的卫星到同一接收机的距离可以解算出接收机的三维坐标。当前用于定位解算的测量值主要有伪距测量值和载波相位测量值这两种，其中伪距测量值的误差大，定位精度不高；同伪距测量值相比，载波相位测量值的误差小，用于定位时，绝对定位的精度可达厘米级，相对定位精度可达毫米级。因此在高精度定位接收机中，普遍采用载波相位测量值来进行定位。

载波相位测量值由初始整周未知数、载波相位整周数和不足整数的小数三部分组成，其中，初始整周未知数无法事先知道，所以又称整周模糊度；小数部分由接收机鉴相器获得；整数部分由接收机电子计数器获得。由于山脉或建筑物的遮挡、恶劣的电离层条件和接收机内置软件错误等原因，导致信号失锁，从而引起整周计数错误，产生所谓的周跳现象。周跳是载波相位测量值的特有问题，它给整周模糊度引入了整周数的偏差，进而影响定位的精度以及可靠性，仅一周的周跳就会导致二十厘米的定位误差。因此，在高精度载波相位定位中，必须能及时发现周跳并正确的给予处理，这就需要一种有效的周跳探测与修复方法。

目前，已有多种方法用于探测与修复周跳，这些方法大致可分为两类：

第一类方法是既可以用于单频接收机，又可用于多频接收机，如高次差法、多项式拟合法、多普勒频移法、Kalman滤波法和码相二次差等方法。高次差法只能探测大周跳，且难以修复周跳；多项式拟合法探测周跳的能力取决于数据采样间隔，难以探测小周跳；多普勒频移法要求接收机能观测多普勒频移量，其测周跳的能力与数据采样间隔有关；Kalman滤波方法仅适用于载体作匀速运动且信号失锁较短的情况；码相二次差方法易受伪距测量噪声的影响，探测效果不理想。

第二类方法是仅适用于多频接收机，如电离层残差法、宽巷组合法、伪距/载波相位组合法等。电离层残差法能探测小周跳，但是当双频载波相位同时发生周跳且周跳大小之比接近频率比时，该方法失效；宽巷组合法探测周跳的能力强，但要求接收机能接收双频精密P码数据，这限制了其应用；伪距/载波相位组合法探测周跳的能力依赖于伪距测量精度，通常难以探测小周跳。

可见上述方法都存在一定的局限性，尤其对于单频接收机而言，普遍难以探测和修复小周跳。

发明内容

为了克服现有技术的不足，本发明提供一种载波相位测量中基于HHT和支持向量机的周跳探测与修复方法，对周跳检验量的时间序列经EMD分解后，对IMF分量作Hilbert变换，根据高频IMF分量的Hilbert幅值探测周跳，再结合支持向量机回归技术修复周跳，提高了周跳探测的灵敏度以及周跳修复的准确性。

本发明解决其技术问题所采用的技术方案包括以下步骤：

步骤1、根据伪距和载波相位观测数据构成码相二次差时间序列D(t)，历元t＝1～L,其中L为大于1的自然数；

步骤2、对时间序列D(t)通过EMD分解算法进行多时间尺度分解，得到n个IMF分量和一个趋势分量，其中n为大于1的自然数；

步骤3、对各IMF分量进行Hilbert变换，得到各分量的Hilbert幅值；

步骤4、判断各分量的Hilbert幅值在t历元是否发生突变；

如果各分量imf_i的Hilbert幅值在t历元同时发生突变，则认为在t历元发生周跳，记录在该历元幅值发生突变的分量，执行步骤5；

如果各分量imf_i的Hilbert幅值在t历元没有同时发生突变，则认为没有发生周跳，对历元t加1，继续执行步骤4；

其中，判断各分量的Hilbert幅值在t历元发生突变的条件为：

$a_{\mathrm{imfi}}\left(t\right)-\stackrel{\‾}{a_{\mathrm{imfi}}}>{3\mathrm{\σ}}_{\mathrm{imfi}}$

其中，a_imfi(t)为分量imf_i在t历元的Hilbert幅值；其均值

其方差 ${\mathrm{\σ}}_{\mathrm{imfi}}=\sqrt{\frac{\underset{i=1}{\overset{t}{\mathrm{\Σ}}}{[a_{\mathrm{imfi}}\left(i\right)-\stackrel{\‾}{a_{\mathrm{imfi}}}]}^2}{t}};$

步骤5、对步骤4所记录的幅值发生突变的各分量的时间序列imf_i(j)进行相空间重构，j＝1～t，将一维的时间序列转化成矩阵形式，构成训练样本集；

步骤6、采用步骤5所述的各分量训练样本集，建立与之对应的支持向量机回归预测模型，对每个支持向量机回归预测模型进行预测，得到幅值发生突变的各分量在t历元的预测值；

步骤7、通过对比预测值与实际值的不同来确定周跳的大小ΔN，其中ΔN为整数；

周跳的大小ΔN根据下式计算：

ΔN＝Int[∑imf_i(t)-∑imf_i′(t)]

其中，Int表示取最接近的整数；∑imf_i′(t)表示幅值发生突变的各分量在t历元的预测值之和；∑imf_i(t)表示幅值发生突变的各分量在t历元的实际值之和；

步骤8、对周跳进行修复，重复步骤1～步骤7，直至t＞L。

本发明的有益效果是：本发明提出的基于HHT和支持向量机的周跳探测与修复方法，适用于单频和多频接收机，可以探测到一周的周跳，在确定周跳的大小时所受伪距测量噪声的影响小，因此计算出的周跳较为准确。本文明提高了周跳探测与修复的精度和可靠性。

附图说明

图1是基于HHT和支持向量机的周跳探测与修复方法的流程图；

图2是对码相二次差时间序列进行EMD分解的具体处理流程图。

具体实施方式

下面结合附图和实施例对本发明进一步说明，本发明包括但不仅限于下述实施例。

基于HHT和支持向量机的周跳探测与修复方法，包括以下步骤：

步骤1、根据伪距和载波相位观测数据构成码相二次差时间序列D(t)，t＝1～L,其中L为大于1的自然数；

步骤2、对步骤1所述的时间序列D(t)，通过EMD分解算法进行多时间尺度分解，得到n个IMF分量和一个趋势分量，其中n为大于1的自然数；

步骤3、对步骤2所述的各分量进行Hilbert变换，得到各分量的Hilbert幅值；

步骤4、判断步骤3所述的各分量的Hilbert幅值在t历元是否发生突变；

如果分量imf₁和imf₂的Hilbert幅值在t历元同时发生突变，则认为在t历元发生周跳，记录在该历元幅值发生突变的分量，执行步骤5；

如果分量imf₁和imf₂的Hilbert幅值在t历元没有同时发生突变，则认为没有发生周跳，令t＝t+1，继续执行步骤4；

其中，判断各分量的Hilbert幅值在t历元发生突变的条件为：

$a_{\mathrm{imfi}}\left(t\right)-\stackrel{\‾}{a_{\mathrm{imfi}}}>{3\mathrm{\σ}}_{\mathrm{imfi}}$

其中，a_imfi(t)为分量

在t历元的Hilbert幅值； ${\mathrm{\σ}}_{\mathrm{imfi}}=\sqrt{\frac{\underset{i=1}{\overset{t}{\mathrm{\Σ}}}{[a_{\mathrm{imfi}}\left(i\right)-\stackrel{\‾}{a_{\mathrm{imfi}}}]}^2}{t}}.$

步骤5、对步骤4所记录的幅值发生突变的各分量的时间序列imf_i(j)，j＝1～t，进行相空间重构，将一维的时间序列转化成矩阵形式，构成训练样本集；

步骤6、采用步骤5所述的各分量训练样本集，建立与之对应的支持向量机回归预测模型，对每个模型进行预测，得到幅值发生突变的各分量在t历元的预测值；

步骤7、通过对比预测值与实际值的不同来确定周跳的大小ΔN，其中ΔN为整数；

周跳的大小ΔN根据下式计算：

ΔN＝Int[∑imf_i(t)-∑imf_i′(t)]

其中，Int表示取最接近的整数；∑imf_i′(t)表示幅值发生突变的各分量在t历元的预测值之和；∑imf_i(t)表示幅值发生突变的各分量在t历元的实际值之和。

步骤8、对周跳进行修复，重复步骤1～步骤7，直至t＞L。

实施例：

图1是基于HHT和支持向量机的周跳探测与修复方法的流程图，包括如下步骤：

步骤1、根据伪距和载波相位观测值构成码相二次差时间序列D(t)，t＝1～L,其中L为大于1的自然数，表示序列长度，时间序列D(t)如图2所示；

步骤2、对步骤1所述的时间序列D(t)，通过EMD分解算法进行多时间尺度分解，得到若干IMF分量和一个趋势分量；

图2提供了对时间序列D(t)进行EMD分解的处理流程，具体包括以下步骤：

步骤21、判断D(t)的极值点个数是否大于2。如果是则说明序列有平稳化的必要，执行步骤22，否则说明序列是平稳的，无需分解或分解已结束，执行步骤28；

步骤22、令H(t)＝D(t)，开始在H(t)中进行提取IMF分量的迭代操作；

步骤23、找出序列H(t)中所有的局部极大值点和局部极小值点；

步骤24、对所有的局部极大值点和局部极小值点，通过三次样条插值法形成上包络线H_max(t)和下包络线H_min(t)，使H(t)满足H_min(t)＜H(t)＜H_max(t)；

步骤25、计算上下包络线的算术平均值并令H(t)＝H(t)-m(t)；

步骤25a、判断H(t)是否已满足成为IMF分量的判定条件，包括以下三点：

①极值点个数和过零点个数一致或至多相差一个；

②上下包络线的算术平均值接近0；

③两次连续处理结果之间的标准差SD的取值在0.2～0.3之间，SD定义为：

$\mathrm{SD}=\underset{t=1}{\overset{L}{\mathrm{\Σ}}}\left[\frac{{|H_{k-1}\left(t\right)-H_k\left(t\right)|}^2}{{H^2}_{k-1}\left(t\right)}\right]$

其中，H_k-1(t)和H_k(t)是在步骤23至步骤25中连续两次的处理结果，k代表这一过程的处理次数。

若同时满足上述三个条件，执行步骤26，否则返回，继续执行步骤23；

步骤26、成功提取一个IMF分量，即imf_i(t)＝H_k(t)；

步骤27、从D(t)减去提取的IMF分量，即D(t)＝D(t)-imf_i(t)，并执行步骤21判断D(t)的极值点个数；

步骤28、此时所有的IMF分量都被提取出来，剩余的D(t)则称为趋势分量r_n(t)；

步骤29、经EMD分解的时间序列D(t)可表示为：

$D\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{imf}}_i\left(t\right)+r_n\left(t\right)$

步骤3、对步骤2所述的各分量进行Hilbert变换，得到各分量的Hilbert幅值；

其中，各分量的Hilbert幅值的计算公式为：

$a_{\mathrm{imfi}}\left(t\right)=\sqrt{{{\mathrm{imf}}_i}^2\left(t\right)+{\left[\frac{1}{\mathrm{\π}}{\∫}_{-\∞}^{\∞}\frac{{\mathrm{imf}}_i\left(\mathrm{\τ}\right)}{t-\mathrm{\τ}}\mathrm{d\τ}\right]}^2}$

步骤4、判断步骤3所述的各分量的Hilbert幅值在t历元是否发生突变；

如果分量imf₁和imf₂的Hilbert幅值在t历元同时发生突变，则认为在t历元发生周跳，记录在该历元幅值发生突变的分量，执行步骤5；

如果分量imf₁和imf₂的Hilbert幅值在t历元没有同时发生突变，则认为没有发生周跳，执行步骤4a；

步骤4a、令t＝t+1，执行步骤4b；

步骤4b、判断t是否大于L；

如果t≤L，继续执行步骤4；

如果t＞L，退出程序；

其中，判断各分量的Hilbert幅值在t历元发生突变的条件为：

$a_{\mathrm{imfi}}\left(t\right)-\stackrel{\‾}{a_{\mathrm{imfi}}}>{3\mathrm{\σ}}_{\mathrm{imfi}}$

其中， $\stackrel{\‾}{a_{\mathrm{imfi}}}=\frac{\underset{i=1}{\overset{t}{\mathrm{\Σ}}}{a}_{\mathrm{imfi}}\left(i\right)}{t};{\mathrm{\σ}}_{\mathrm{imfi}}=\sqrt{\frac{\underset{i=1}{\overset{t}{\mathrm{\Σ}}}{[a_{\mathrm{imfi}}\left(i\right)-\stackrel{\‾}{a_{\mathrm{imfi}}}]}^2}{t}}.$

步骤5、对步骤4所记录的幅值发生突变的各分量的时间序列imf_i(j)，j＝1～t，进行相空间重构，将一维的时间序列转化成矩阵形式，构成训练样本集{(x_imfi,j,y_imfi,j)|j＝1,2,…,t-2-(m-1)τ}：

$\left\{\begin{array}{c}{x}_{\mathrm{imfi},j}=[{\mathrm{imf}}_i\left(j\right),{\mathrm{imf}}_i(j+\mathrm{\τ}),\·\·\·,{\mathrm{imf}}_i(j+(m-1)\mathrm{\τ})]\\ y_{\mathrm{imfi},j}={\mathrm{imf}}_i(j+(m-1)\mathrm{\τ}+1)\end{array}\right.j=\mathrm{1,2},\·\·\·,t-2-(m-1)\mathrm{\τ}$

其中，m、τ分别为相空间重构的嵌入维数和时间延迟，采用C-C法选取。

步骤6、采用步骤5所述的各分量训练样本集，建立与之对应的支持向量机回归预测模型，对每个模型进行预测，得到幅值发生突变的各分量在t历元的预测值；

其中，所述的支持向量机回归模型采用最小二乘支持向量机回归模型，核函数选用径向基函数，采用基于网格搜索的交叉验证法分别确定径向基核宽度与正则化参数。

其中，支持向量机回归模型的输入向量为x_imfi＝[imf_i(t-1-(m-1)τ),imf_i(t-1-(m-1)τ+τ),…,imf_i(t-1)]，输出为imf_i′(t)。

步骤7、通过对比预测值与实际值的不同来确定周跳的大小ΔN，其中ΔN为整数，并对周跳进行修复；

其中，周跳的大小ΔN根据下式计算：

ΔN＝Int[∑imf_i(t)-∑imf_i′(t)]

其中，Int表示取最接近的整数；∑imf_i′(t)表示幅值发生突变的各分量在t历元的预测值之和；∑imf_i(t)表示幅值发生突变的各分量在t历元的实际值之和。

步骤8、对周跳进行修复，重复步骤1～步骤7，直至t＞L。

Claims (1)

1.一种基于HHT和支持向量机的周跳探测与修复方法，其特征在于包括下述步骤：

步骤1、根据伪距和载波相位观测数据构成码相二次差时间序列D(t)，历元t＝1～L,其中L为大于1的自然数；

步骤2、对时间序列D(t)通过EMD分解算法进行多时间尺度分解，得到n个IMF分量和一个趋势分量，其中n为大于1的自然数；

步骤3、对各IMF分量进行Hilbert变换，得到各分量的Hilbert幅值；

步骤4、判断各分量的Hilbert幅值在t历元是否发生突变；

如果各分量imf_i的Hilbert幅值在t历元同时发生突变，则认为在t历元发生周跳，记录在该历元幅值发生突变的分量，执行步骤5；

如果各分量imf_i的Hilbert幅值在t历元没有同时发生突变，则认为没有发生周跳，对历元t加1，继续执行步骤4；

其中，判断各分量的Hilbert幅值在t历元发生突变的条件为

其中，a_imfi(t)为分量imf_i在t历元的Hilbert幅值；其均值

其方差 ${\mathrm{\σ}}_{\mathrm{imfi}}=\sqrt{\frac{\underset{i=1}{\overset{t}{\mathrm{\Σ}}}{[a_{\mathrm{imfi}}\left(i\right)-\stackrel{\‾}{a_{\mathrm{imfi}}}]}^2}{t}};$

步骤5、对步骤4所记录的幅值发生突变的各分量的时间序列imf_i(j)进行相空间重构，j＝1～t，将一维的时间序列转化成矩阵形式，构成训练样本集；

步骤6、采用步骤5所述的各分量训练样本集，建立与之对应的支持向量机回归预测模型，对每个支持向量机回归预测模型进行预测，得到幅值发生突变的各分量在t历元的预测值；

步骤7、通过对比预测值与实际值的不同来确定周跳的大小ΔN，其中ΔN为整数；

周跳的大小ΔN＝Int[∑imf_i(t)-∑imf_i′(t)]，

其中，Int表示取最接近的整数；∑imf_i′(t)表示幅值发生突变的各分量在t历元的预测值之和；∑imf_i(t)表示幅值发生突变的各分量在t历元的实际值之和；

步骤8、对周跳进行修复，重复步骤1～步骤7，直至t＞L。

Images