CN103973263B - 一种逼近滤波方法 - Google Patents

Abstract

本发明涉及一种新的逼近滤波方法，其步骤包括：1、建立广义延拓逼近多项式，2、构造广义延拓逼近最优化求解模型，3、求解广义延拓逼近最优化模型，4、广义延拓逼近滤波方法最优估计值的求解，5、残余误差σ大小的求解，6、求解在t_n+2时刻的最优状态值。本发明所述方法的优越效果在于，所述逼近滤波方法在递推逼近时，采用了广义延拓逼近方法，它是非线性逼近模型,兼有插值和拟合两种功能。与线性逼近方法和最小二乘逼近方法相比较，具有逼近精度高、方法灵活方便、适应性强等特点。

Description

一种逼近滤波方法

技术领域

本发明涉及数据处理和信号处理领域，是一种基于广义延拓逼近方法基础上新的高精度滤波方法，称为SHL广义逼近滤波方法。

技术背景

从包含着误差的数据中提取出需要的数据信息；从包含着干扰及噪声的信号中提取需要的信号，这种数据和信号的处理方法称为滤波方法。其中最有影响的滤波方法是卡尔曼滤波方法。卡尔曼滤波方法创立于1960年，由卡尔曼把状态量的概念引入到最小均方误差估计中，建立的一种线性、无偏、以均方误差最小为准则的最优估计方法。所述方法采用递推形式，根据前一时刻的状态量的最佳估计值，按照低阶状态量的递推公式，获得低阶状态外推后状态量新的预估值,这时再通过与观测值的加权组合,就可以实时地计算出所需状态量的实时最优估计值。卡尔曼滤波方法计算量少、计算速度快、逼近精度好等特点,所以广泛应用于众多领域。

卡尔曼滤波方法与其它滤波方法一样，其滤波的最基本准则是最小二乘逼近准则。申请人在前发明了广义延拓最小二乘逼近方法，并把广义延拓最小二乘逼近方法成功应用于滤波逼近，创立了一种新的逼近滤波方法。

发明内容

本发明提供一种逼近滤波方法，它把广义延拓最小二乘逼近模型与方法应用于滤波过程之中，构筑出的新的滤波方法。

本发明所述的逼近滤波方法的技术解决方案是：

1、建立多约束非线性方程优化求解模型，选用的目标函数是待求解代入测量方程后残差的平方值函数，通过求残差平方值的极小化实现最优逼近；选用的约束条件为位置约束、高度约束、速度及加速度约束和行进方向角约束等多类约束条件，所述约束条件能够改善解域，减少优化求解时间，提高测量解的精度。所述广义延拓逼近方法采用二次或高次非线性函数作为逼近多项式的非线性模型。非线性模型含有低阶状态量的变化情况，逼近精度高于线性模型，逼近效果更好。

2、本发明所述的滤波方法采用广义状态量和多重联立组合形式，所述广义状态量包括状态量及其多阶导数值；多重联立组合形式包括通过状态量多阶导数的积分递推，及广义测量值的并列组合形式,获得状态方程或状态方程组。由于广义状态量和多重联立组合形式引入了运动惯性，引入了数据间的联系,增加了冗余度，提高了前后解之间的相关性，能压低随机误差，获得关联度高的轨迹解。

3、本发明对状态方程或状态方程组进行融合，或进行几类方程的融合求解。求解时，采用广义延拓逼近方法，建立组合当前时刻测量值与预估值,求最优估计值的广义延拓逼近模型。在模型中，把待求时刻的最优估计值点作为插值锁定点；把一段先验逼近点上的最优估计值，或测量值，或预估值作为拟合数据值。这样，可以把逼近多项式系数和权系数同时作为优化变量，通过求目标函数逼近残差平方和的极小化，优化得到逼近多项式的系数，同时解决了权系数的优选问题,实现逼近滤波,最后,获得融合的最优估计值解。

本发明所述的逼近滤波方法主要实施步骤是：

1、建立广义延拓逼近多项式

在广义延拓逼近多项式中，假设已知t_n时刻的最优状态值及以前的状态测量值x_{i i＝m,.....,n-1}，则可以建立下述广义延拓逼近多项式：

上式(1)中，x(t_i)为广义延拓逼近多项式；t_i为时序变量；a₁,a₂,a₃为广义延拓逼近多项式的待求系数；m为分段逼近时的起始点；n为分段逼近时的终止点。

所述广义延拓逼近多项式作为t_n+1时刻最优状态估计值的插值约束，

即：

上式(2)中为t_n+1时刻状态的最佳逼近估计值。

2、构造广义延拓逼近最优化求解模型

构造广义延拓逼近最优化求解模型，求解广义延拓逼近多项式的系数a₁,a₂,a₃，求解需满足：

式(3)中，I(a₁,a₂,a₃)为优化的目标函数；x_i为t_i时刻的状态测量值，其中，为t_n+1时刻状态的最佳逼近估计值，在模型中把作为插值点；x_{i i＝m,m+1,.Λn}为t_m时刻至t_n时刻的状态测量值，模型中把它们作为拟合处理点；a₁,a₂,a₃为广义延拓逼近多项式的待求系数。

式(3)即为带插值点约束的非线性最优化模型，求解上述非线性最优化模型即可得到广义延拓最优状态逼近多项式。

3、求解广义延拓逼近最优化模型

(1)非线性测量方程的直接求解方法

若测量方程为非线性测量方程，可直接对非线性测量方程采用单纯形法、复合形法等直接法进行求解，即：

F(x_i)＝p_{i i＝1,2,...,n}……(4)

式(4)中，F(x_i)为非线性函数；x_i为函数变量；p_i为测量量；i为序号；n为序号总数。

(2)递推所求时元状态预估量

在滤波中，利用状态量外推方程可以改善解的精度及轨迹解的关联性。状态量可以是与测量相关的导数值或微分值；也可以是其它差分信息量。在滤波中，需要预估下一时元的状态预估量。本方法采用在当前最优状态估计值的基础上，加上可观测的状态量的一阶导数值及高阶导数值与时间间隔的乘积，作为下一时刻的状态预估量，采用的公式如下：

${\tilde{x}}_{n+1}={\hat{x}}_n+\frac{d(x_{n+1}-x_n)}{d(t_{n+1}-t_n)}\mathrm{\Δt}+\frac{1}{2}\•\frac{d\″(x_{n+1}-x_n)}{d\″(t_{n+1}-t_n)}\mathrm{\Δ}{t}^2\mathrm{\Λ}$

通过上式(5)可以求得所求时元的状态预估量，初始状态估计量可用前几个时元的状态测量最小二乘拟合得到。

(3)优化求解权系数k_x(n+1)

在选择权系数时(在卡尔曼滤波方法里，权系数的选取需要迭代求解，本发明所述方法则不需要进行迭代求解)，首先建立广义延拓逼近模型，把待求的t_n+1时刻的最优位置坐标估计值作为插值锁定点，把其它先验点上的状态测量值做为最小二乘拟合逼近点，建立逼近残差平方和的极小化模型，通过直接法求解非线性目标函数逼近残差的极小化值，获得权系数k_x(n+1)。

本发明所述方法为权系数的获得开辟了一条新的选择途径，具体实施步骤如下：

把t_n+1时刻的状态预估量和t_n+1时刻的状态测量量x_n+1组合起来，求t_n+1时刻的状态最优估计量

其中，x_n+1为实际状态测量值。

在上述式(6)中，权系数k_x(n+1)的求解通过求解下述模型，即式(7)来实现，

在模型式(7)中，是可行区间变量k_x(n+1)的约束区间数， 是可行区间变量a₁,a₂,a₃的约束区间数，

通过对上述模型(7)进行直接法等优化算法求解，得到广义延拓逼近多项式的系数a₁,a₂,a₃，状态量组合时的权系数k_x(n+1)以及目标函数的极小化值F(I)。

4、广义延拓逼近滤波方法最优估计值的求解

把所求时元的状态预估值和权值k_x(n+1)带入到上式(6)中，即可得到广义延拓滤波方法的所求时刻状态最优估计量

5、残余误差σ大小的求解

把a₁,a₂,a₃,k_x(n+1)这些系数代入模型(7)中，可以得到广义延拓多项式逼近时目标函数的极小化值F(I)，若需要表征滤波的逼近程度，可以把优化时最后获得的目标函数的极小化值FI＝minI(a₁,a₂,a₃,k_xn))带入如下公式(8)，得到逼近误差的均方根值，即可判断滤波效果。

6、求解在t_n+2时刻的最优状态值和

用n+1代替n，用t_n+2代替t_n+1，重复上述步骤1-5，迭代渐进，便可以求得在t_n+2时刻的最优状态值

本发明所述逼近滤波方法的优越效果在于：

1、所述逼近滤波方法在递推逼近时，采用了广义延拓逼近方法，它是非线性逼近模型,兼有插值和拟合两种功能。与线性逼近方法和最小二乘逼近方法相比较，具有逼近精度高、方法灵活方便、适应性强等特点。

2、所述逼近滤波方法采用了多约束最优化模型，它以逼近偏差平方和函数，或以残差的平方值函数为最优化目标函数，以逼近多项式系数和权系数，或以状态量等为优化变量，通过最优化算法便可以获得优化结果，具有概念清晰、方法直观、求解范围大、计算机运算速度快等特点。

3、所述逼近滤波方法能够直接利用低阶状态量等作最佳预估值的递推；也可以用方向角等其它状态量作为约束条件。所以本方法具有适应面宽、扩展性强等特点，有进一步研究与探索的价值。

附图说明

图1是本发明所述用matlab模拟车辆沿道路正弦行进时的跟踪轨迹图；

图2是本发明所述用matlab模拟车辆沿道路正弦行进时的平滑轨迹图；

图3是本发明所述车辆实测数据x(东)向跟踪轨迹图；

图4是本发明所述车辆实测数据y(北)向跟踪轨迹图；

图5是本发明所述车辆实测数据xy(东北)向跟踪轨迹图；

图6是本发明所述车辆实测数据xy(东北)向平滑轨迹图。

具体实施方式

本发明以卫星定位导航为例说明具体处理时采用的滤波方法的主要求解步骤如下：

1、初始最优位置状态量

现假设不考虑加速度测量值，仅考虑获得了位置及速度的序列测量值的情况，若令已获得的位置及速度的序列测量值为:

{x_i,y_i,v_xi,v_yi i＝1,2,.....,m,.....,n} ……(9)

其中,x_i和y_i分别是沿x方向和y方向的位置坐标；v_xi和v_yi分别为沿x方向和y方向的速度分量；i为自变量，当i为时间变量时，i＝1,2,.....,m,.....,n即表示t₁,t₂,.....,t_m,.....,t_n时刻。若从已获得的位置及速度的测量值序列中,选取了一段最邻近时刻的位置及速度测量值数据序列,即从t_m时刻至t_n时刻的一段位置及速度测量值数据序列:

{x _i,y _i,v _xi,v _yi i＝k,.....,n} ……(10)

通过对t_n时刻以前的位置状态测量数据进行最小二乘拟合的方法进行已知数据的处理(长度选取10)，用处理后逼近多项式上t_n时刻的值，近似作为t_n时刻的最优初始位置状态量

2、初始最优速度量的确定

由于在导航定位中速度的测量值精确度较高，这里初始时刻t_n的最优速度值及之后各时刻的速度的最优速度值都直接选用对应时刻的速度测量值v_xi,v_yi。

(1)所求时元状态预估量的确定

1)递推所用速度量即一阶状态量的确定

将t_n+1时刻速度的测量值v_x(n+1),v_y(n+1)和t_n时刻速度的测量值v_xn,v_yn取平均作为递推t_n+1时刻位置状态预估值的一阶状态值，如下式(11)、(12)

2)t_n+1时刻位置状态量的最佳预估值和的确定

得到了平均速度值之后，再利用位置和速度运动状态量之间的关系，在t_n时刻最优位置状态量的基础上，加上平均速度值与时间间隔Δt的乘积，通过公式(13)、(14)推导得到时刻t_n+1位置状态量的最佳预估值和

(2)t_n+1时刻最优权系数k_x(n+1),k_y(n+1)的求解

1)t_n+1时刻的最优状态估计量的表示

将t_n+1时刻的位置状态测量值x_n+1,y_n+1、步骤2)求解得到的t_n+1时刻状态预估量和待求的最优权系数k_x(n+1),k_y(n+1)代入式(15)、(16)中可得最优状态估计量的表示量

2)t_n+1时刻最优权系数k_x(n+1),k_y(n+1)的求解

将t_n+1时刻及之前的位置状态测量值v_xi,v_yi和上式(15)、(16)表示的t_n+1时刻的最优状态估计量代入下面的广义延拓逼近模型(17)、(18)中得到模型式(19)、(20)

及

得

及

在模型式(19)和(20)中，是可行区间变量k_x(n+1)的约束区间数， 是可行区间变量a₁,a₂,a₃的约束区间数， 是可行区间变量k_y(n+1)的约束区间数， 是可行区间变量b₁,b₂,b₃的约束区间数，

通过直接法等最优化算法，求解模型式(19)和模型式(20)，可得广义延拓逼近多项式的系数(a₁,a₂,a₃)、(b₁,b₂,b₃)、权系数k_x(n+1)和k_y(n+1)，以及目标函数最优化极小值FI₁,FI₂的大小。

3、广义延拓逼近滤波方法最优估计值的求解

将求解式(15)、(16)得到的t_n+1时刻最优权系数k_x(n+1),k_y(n+1)代入到式(19)、(20)中，即可得到t_n+1时刻广义延拓逼近滤波方法最优估计值4、残余误差σ大小的求解

通过求解模型(19)和模型(20)，可以得到广义延拓多项式逼近式目标函数的极小化值FI₁,FI₂。若需要表征滤波的逼近程度，可以把优化时最后获得的目标函数的极小化值FI₁,FI₂分别代入如下公式(21)，得到逼近误差的均方根值，即可判断滤波效果

用n+1代替n，用t_n+2代替t_n+1，重复上述步骤2)-5)，迭代渐进，便可以求得在t_n+2时刻的最优位置坐标值和这样逐次递进便可以得到一系列的经过滤波处理后的位置坐标值，实现位置坐标数据的实时滤波。

下面结合说明书附图详细介绍本发明所述逼近滤波方法的实例。

实例一

本发明的仿真试验是在matlab上仿真模拟了车辆沿道路正弦行进时的轨迹情况。原始测量值用表示，经本滤波方法滤波后的数据用表示，理论值用表示。

假设位置测量值的测量误差为±1m/s之间，k的速度为(s_k+1-s_k)/(t_k+1-t_k)(无误差理论值)。当做广义延拓逼近时，可以获得新的逼近滤波仿真的效果(如图1所示)。当广义延拓逼近时插值点前移后，则可以起到一定的平滑作用，获得的新的平滑滤波仿真的效果，如图2所示。

实例二

若考虑在二维平面上数据的滤波，现以GPS终端在平面地图上的导航应用为例进行说明，如图3--5所示。

在导航应用时，要在地图上表示出终端的位置坐标，描绘出终端的运动轨迹。因为x方向和y方向的位置坐标分别是时间的函数，所以采用分别计算终端在x方向及y方向的位置坐标状态量的变化来实现，最后组合起来获得x方向及y方向的位置坐标状态同时变化的轨迹。

申请人采用中科微公司的导航模块在国家天文台园区实测的一组位置测量值数据。图3为X方向(东向)的滤波数据处理(横坐标为时间t)图。图4为Y方向(北向)的滤波数据处理(横坐标为时间t)图。图5为最终XY方向(东北方向)的滤波数据处理图(横坐标为东向坐标X,纵坐标为北向坐标Y)。当广义延拓逼近时插值点前移后，则可以起到一定的平滑作用，获得的新的平滑滤波仿真的效果图，如图6所示，由于实测范围不大，Z方向(天方向)的坐标设置为常数，图中未予表示。由图3、图4、图5中滤波前后数据对比可以发现，经过新的滤波后的数据要比原始测量值数据整体上光滑，在有数据突变出的峰值处的也要比原始的测量数据幅度低，平滑和跟随效果较好，且不是单纯的平滑拟合，反映了滤波效果。

通过仿真和应用例的试验，可以看到逼近滤波数据处理方法确实能从包含着误差的数据中提取出需要的信息，能从包含着噪声的信号中提取出需要的信号，而且不论信号是有起伏变化，或有方向变化的，都能提取出来，不会完全丢失。它也不会像数据平滑处理方法那样，对观测数据进行修匀时，把外界的干扰和影响消除了，让修匀后的观测值落在一条较光滑的线上，但有时这样做时，往往把变化的信号也平滑掉了，这是不希望的。所以说逼近滤波数据处理方法不是一种简单的数据平滑方法，确实是一种新的、有广泛应用价值的滤波方法。

Claims (1)

1.一种逼近滤波方法，其步骤如下：

步骤1、建立广义延拓逼近多项式

在广义延拓逼近多项式(1)中，假设已知t_n时刻的最优状态值及以前的状态测量值x_{i i＝m,.....,n-1}，则可以建立下述广义延拓逼近多项式(1)：

$x\left(t_i\right)=a_1+a_2{t}_i+a_3{t}_i^2,i=m,m+1,......n-1......\left(1\right);$

上式(1)中，x(t_i)为广义延拓逼近多项式；t_i为时序变量；a₁,a₂,a₃为广义延拓逼近多项式的待求系数；m为分段逼近时的起始点；n为分段逼近时的终止点，所述广义延拓逼近多项式作为t_n+1时刻最优状态估计值的插值约束，

即：

${\hat{x}}_{n+1}=a_1+a_2{t}_{n+1}+a_3{t}_{n+1}^2......\left(2\right);$

上式(2)中为t_n+1时刻状态的最佳逼近估计值；

步骤2、构造广义延拓逼近最优化求解模型

构造广义延拓逼近最优化求解模型，求解广义延拓逼近多项式的系数a₁,a₂,a₃，求解需满足：

$\left\{\begin{array}{cc}\min & I(a_1,a_2,a_3)=\underset{i=m}{\overset{n}{\Σ}}{\[a_1+a_2{t}_i+a_3{t}_i^2-x_i\]}^2\\ S.T.& a_1+a_2{t}_{n+1}+a_3{t}_{n+1}^2={\hat{x}}_{n+1}\end{array}\right.......\left(3\right);$

式(3)中，I(a₁,a₂,a₃)为优化的目标函数；x_i为t_i时刻的状态测量值，其中，为t_n+1时刻状态的最佳逼近估计值，在模型中把作为插值点；x_{i i＝m,m+1,.Λn}为t_m时刻至t_n时刻的状态测量值，模型中把它们作为拟合处理点；a₁,a₂,a₃为广义延拓逼近多项式的待求系数，式(3)即为带插值点约束的非线性最优化模型，求解上述非线性最优化模型即可得到广义延拓最优状态逼近多项式；

步骤3、求解广义延拓逼近最优化模型

步骤3.1非线性测量方程的直接求解方法

若测量方程为非线性测量方程，可直接对非线性测量方程采用单纯形法、复合形法直接法进行求解，即：

F(x_i)＝p_{i i＝1,2,...,n}……(4)；

式(4)中，F(x_i)为非线性函数；x_i为函数变量；p_i为测量量；i为序号；n为序号总数，

步骤3.2递推所求时元状态预估量

在滤波中，利用状态量外推方程可以改善解的精度及轨迹解的关联性，状态量可以是与测量相关的导数值或微分值，也可以是其它差分信息量，在滤波中，需要预估下一时元的状态预估量，本方法采用在当前最优状态估计值的基础上，加上可观测的状态量的一阶导数值及高阶导数值与时间间隔的乘积，作为下一时刻的状态预估量，采用的公式如下：

${\tilde{x}}_{n+1}={\hat{x}}_n+\frac{d(x_{n+1}-x_n)}{d(t_{n+1}-t_n)}\mathrm{\Δ}t+\frac{1}{2}\·\frac{d^{\′\′}(x_{n+1}-x_n)}{d^{\′\′}(t_{n+1}-t_n)}{\mathrm{\Δt}}^2+\mathrm{\Λ}......\left(5\right);$

通过上式(5)可以求得所求时元的状态预估量，初始状态估计量可用前几个时元的状态测量最小二乘拟合得到，

步骤3.3优化求解权系数k_x(n+1)

把t_n+1时刻的状态预估量和t_n+1时刻的状态测量量x_n+1组合起来，求t_n+1时刻的状态最优估计量

${\hat{x}}_{n+1}={\tilde{x}}_{n+1}+k_{x(n+1)}(x_{n+1}-{\tilde{x}}_{n+1})......\left(6\right);$

其中，x_n+1为实际状态测量值，在式(6)中，权系数k_x(n+1)的求解通过求解下述模型，即求解以下模型即式(7)来实现，

在模型式(7)中，是可行区间变量k_x(n+1)的约束区间数， 是可行区间变量a₁,a₂,a₃的约束区间数， 通过对上述模型(7)进行直接法优化算法求解，得到广义延拓逼近多项式的系数a₁,a₂,a₃，状态量组合时的权系数k_x(n+1)以及目标函数的极小化值F(I)；

步骤4、广义延拓逼近滤波方法最优估计值的求解

把所求时元的状态预估值和权值k_x(n+1)带入到上式(6)中，即可得到广义延拓滤波方法的所求时刻状态最优估计量

步骤5、残余误差σ大小的求解

把a₁,a₂,a₃,k_x(n+1)系数代入模型(7)中，得到广义延拓多项式逼近时目标函数的极小化值F(I)，若需要表征滤波的逼近程度，可以把优化时最后获得的目标函数的极小化值FI＝minI(a₁,a₂,a₃,k_xn))带入如下公式(8)，得到逼近误差的均方根值：

$\mathrm{\σ}=\sqrt{\frac{FI}{n-m+1}}......\left(8\right);$

步骤6、求解在t_n+2时刻的最优状态值和

用n+1代替n，用t_n+2代替t_n+1，重复上述步骤1-5，迭代渐进，便可以求得在t_n+2时刻的最优状态值