CN106909736B - 一种扰动运动统计特性的解析计算方法 - Google Patents

Abstract

本发明涉及一种扰动运动统计特性的解析计算方法，该方法经过扰动信息提取、自相关序列计算、自相关函数拟合、功率谱密度计算等步骤，定量给出载体运动的自相关和功率谱密度函数的具体参数，可以直接计算外部扰动频率点的功率谱密度，克服了以前功率谱密度尖峰与跳变严重，对扰动特性描述不精确的缺点，实现了对载体运动扰动特性的量化描述，精度高、数据平滑；本发明方法能够通过数值计算方法，准确、简单、定量计算出载体运动的扰动特性。

Description

一种扰动运动统计特性的解析计算方法

技术领域

本发明涉及一种扰动运动统计特性的解析计算方法，属于力学环境试验技术领域。

背景技术

扰动运动的谱分析可以通过对扰动信号直接进行傅里叶变换实现，根据幅值谱大小及其对应频率得到扰动特性。但是这种方法得到的幅值谱曲线存在较大的尖峰与跳变(如图2所示)，使得扰动特性结果不够准确。此后人们通过扰动信号的功率谱密度分析来获得扰动特性，即通过求取信号的功率谱密度，根据功率谱密度幅值及其对应频率得到扰动特性。这一方法在工程中获得了广泛的应用。

功率谱密度求取方法的主要包括周期图法、相关函数法和ARMA模型。

(1)由周期图法计算功率谱密度

周期图法首先将扰动信号进行分段，之后对各段扰动信号进行傅里叶变换，得到其幅值谱密度X(jω)，最后通过公式：

得到信号的功率谱密度Φ(ω)。其中，N为信号分段总数，T为每段信号的总时长。但是，该方法的缺点是受噪声影响谱密度序列为离散点，由其构成的曲线不平滑。

(2)由自相关序列计算功率谱密度

相关函数法首先求取信号的自相关序列，之后对其进行傅里叶变换得到信号的功率谱密度。由于实际中信号的采样总是离散且时间有限的，因此直接使用公式：

得到的无偏自相关序列与信号自相关函数存在较大偏差，且与周期图法计算的功率谱基本相同，都存在波动量偏大的问题。

(3)基于ARMA模型计算功率谱密度

为减小功率谱的波动量和离散点问题，在现代随机信号处理中采用ARMA模型来求解功率谱密度，但其缺点是ARMA模型的阶次非常高，不利于数据分析。

为此，需要研究一种既具有平滑的优点，又具有简化的优点的扰动运动统计特性的计算方法。

发明内容

本发明的目的在于克服现有技术的上述缺陷，提供一种扰动运动统计特性的解析计算方法，该方法能够通过数值计算方法，准确、简单、定量计算出载体运动的扰动特性。

本发明的上述目的主要是通过如下技术方案予以实现的：

一种扰动运动统计特性的解析计算方法，包括

提取载体运动信息的等间隔周期的时间序列，根据所述时间序列和期望时间序列求出扰动时间序列；

根据所述扰动时间序列求出自相关序列；

将所述自相关序列与自相关解析模型进行拟合，得到自相关解析模型的参数；

根据所述自相关解析模型的参数，求出任意角频率处的功率谱密度。

在上述解析计算方法中，所述提取载体运动信息的等间隔周期的时间序列为x(k)，k＝0,1,…,N-1，N为正整数；

所述期望时间序列为

k＝0,1,…,N-1，N为正整数；

所述扰动时间序列为

计算公式如下：

在上述解析计算方法中，根据所述扰动时间序列求出自相关序列r(k)，具体计算公式如下：

其中：l为非负整数。

在上述解析计算方法中，将所述自相关序列与自相关解析模型进行拟合，得到自相关解析模型的参数μ、λ、α，其中μ为衰减系数，λ为周期、α为初始相位；

所述自相关解析模型为：

其中：t为时间；r(0)为自相关序列中第一个值。

在上述解析计算方法中，将所述自相关序列与自相关解析模型进行拟合，得到自相关解析模型的参数μ、λ、α，包括：

(1)、对自相关序列r(k)进行归一化处理，得到归一化处理后的自相关序列

(2)、从自相关序列

中截取部分自相关序列

作为拟合的观测数据；

(3)、计算当前迭代次数为l时的自相关序列

的待拟合残差b_k，并建立向量矩阵B，判断向量矩阵B是否满足要求，或者当前迭代次数l是否为迭代次数上限，若满足，则μ、λ、α的取值分别为当前迭代次数为l时的μ_l、λ_l、α_l；若不满足，l加1，进入步骤(4)；

(4)、计算上一次待拟合参数的残差，包括衰减系数残差Δμ、周期残差Δλ、初始相位残差Δα；

(5)、根据衰减系数残差Δμ、周期残差Δλ、初始相位残差Δα对上一次待拟合参数进行补偿；

(6)、重复(3)～(5)，利用补偿后的待拟合参数计算自相关序列的待拟合残差b_k，直至满足步骤(3)中的要求。

在上述解析计算方法中，所述步骤(1)中对自相关序列r(k)进行归一化处理，得到归一化处理后的自相关序列

具体公式如下：

在上述解析计算方法中，所述步骤(2)中从自相关序列

中截取部分序列

作为拟合的观测数据，其中k′＝0,1,…,N₂-1，N₂＝[(3t₁)/T_s]，[]表示取整；t₁为自相关序列随时间第一次负向时刻，T_s为自相关序列的采样周期。

在上述解析计算方法中，所述步骤(3)中计算当前迭代次数为l时的自相关序列

的待拟合残差b_k，具体公式如下：

其中：k′＝0,1,…,N₂-1，l为当前迭代次数；

记向量B为：

若向量B满足B^TB<10^-2，或者满足l＝D，D为迭代次数上限，则：

μ＝μ_l

λ＝λ_l

α＝α_l

否则，l自加1，进入步骤(4)；

在上述解析计算方法中，所述步骤(3)中拟合初值的选择如下：

在上述解析计算方法中，所述步骤(4)中计算上一次待拟合参数的残差，包括衰减系数残差Δμ、周期残差Δλ、初始相位残差Δα具体方法如下；

(4.1)计算自相关序列残差的偏导参数：

(4.2)记矩阵A为：

(4.3)采用如下公式计算Δ：

Δ＝(A^TA)^-1A^TB

记向量Δ为：

Δ＝[Δμ Δλ Δα]^T。

在上述解析计算方法中，所述步骤(5)中根据衰减系数残差Δμ、周期残差Δλ、初始相位残差Δα对上一次待拟合参数进行补偿的方法如下：

μ_l＝μ_l-1+Δμ

λ_l＝λ_l-1+Δλ

α_l＝α_l-1+Δα；

其中：μ_l-1、λ_l-1、α_l-1分别为上一次的衰减系数、上一次的周期和上一次的初始相位。

在上述解析计算方法中，根据所述自相关解析模型的参数，求出任意角频率处的功率谱密度的具体方法如下：

其中：ω为角频率，ω＝2πf，f为频率。

本发明与现有技术相比的有益效果是：

(1)本发明给出了一种扰动运动统计特性的解析计算方法，该方法经过扰动信息提取、自相关序列计算、自相关函数拟合、功率谱密度计算等步骤，定量给出载体运动的自相关和功率谱密度函数的具体参数，可以直接计算外部扰动频率点的功率谱密度，克服了以前功率谱密度尖峰与跳变严重，对扰动特性描述不精确的缺点，实现了对载体运动扰动特性的量化描述，精度高、数据平滑；

(2)本发明采用的自相关序列拟合方法对自相关序列的拟合结果更加精确，对载体运动过程的量化描述更加精确，得到的功率谱密度曲线更加平滑；

(3)本发明给出了一种扰动运动统计特性的解析计算方法，具有运算量小，便于编程与工程实现的优点。

附图说明

图1为某次试验东向速度的扰动曲线；

图2为某次试验东向速度扰动的线谱；

图3为采用周期图法计算的扰动功率谱密度；

图4为某次试验东向速度扰动的自相关序列；

图5为本发明解析计算方法的流程图；

图6为本发明某次试验自相关序列的拟合过程图；

图7为本发明自相关拟合结果对比图；

图8为本发明方法与周期图法得到的功率谱密度曲线对比图。

具体实施方式

下面结合附图和具体实施例对本发明作进一步详细的描述：

本发明扰动运动统计特性的解析计算方法，经过扰动信息提取、自相关序列计算、自相关函数拟合、功率谱密度计算等步骤，定量给出载体运动的自相关和功率谱密度函数的具体参数。在此基础上，就可直接计算外部扰动频率点的功率谱密度。

如图5所示为本发明解析计算方法的流程图，本发明提出的惯导设备载体运动扰动特性确定方法，具体步骤如下：

(1)、设采样频率为f_s，采样周期为T_s＝1/f_s，提取载体运动信息的等间隔周期的时间序列值x(k)，k＝0,1,…,N-1；根据x(k)及其期望值序列

可求得扰动序列

具体表达式为

(2)、对扰动序列计算其自相关序列r(k)，计算公式为

其中：l为非负整数。

(3)、自相关解析模型如下：

其中：t为时间；r(0)为自相关序列中第一个值。

将自相关序列与自相关解析模型进行拟合，得到自相关解析模型的参数μ、λ、α，其中μ为衰减系数，λ为周期、α为初始相位；自相关序列的拟合方法如下：

(3.1)、对自相关序列r(k)进行归一化处理，得到归一化处理后的自相关序列

具体方法为：

(3.2)、从自相关序列

中截取部分序列

作为拟合的观测数据，其中k′＝0,1,…,N₂-1，N₂＝[(3t₁)/T_s]，[]表示取整；t₁为自相关序列随时间第一次负向时刻，T_s为自相关序列的采样周期。

(3.3)、选择拟合初值，具体方法为：

其中：t₁为自相关序列随时间第一次负向时刻。

(3.4)、计算自相关序列

的待拟合残差b_k，具体方法为：

其中：k′＝0,1,…,N₂-1，l为当前迭代次数；

记向量B为：

若向量B满足B^TB<10^-2，或者满足l＝D，D为迭代次数上限，则：

μ＝μ_l

λ＝λ_l

α＝α_l

停止自相关函数拟合。否则，l自加1，进入步骤(3.5)。

(3.5)、计算上一次待拟合参数的残差，包括衰减系数残差Δμ、周期残差Δλ、初始相位残差Δα，具体方法如下：

对截断后的自相关序列中的每一个点，计算以下三个参数(自相关序列残差的偏导参数)，

记矩阵A为：

采用如下公式计算Δ：

Δ＝(A^TA)^-1A^TB

记向量Δ为：

Δ＝[Δμ Δλ Δα]^T。

(3.6)、根据衰减系数残差Δμ、周期残差Δλ、初始相位残差Δα对上一次待拟合参数进行补偿，具体方法为：

μ_l＝μ_l-1+Δμ

λ_l＝λ_l-1+Δλ

α_l＝α_l-1+Δα。

(3.7)、重复(3.4)～(3.6)进行迭代计算，即利用补偿后的待拟合参数计算自相关序列的待拟合残差b_k，直至步骤(3)中的任意条件满足。

若步骤(3.4)中选择初值μ₁＝0，

α₁＝0，首次进行计算时即满足B^TB<10^-2，或者满足l＝D，则μ、λ、α的取值为：μ＝μ₁、λ＝λ₁、α＝α₁停止自相关函数拟合。若不满足条件，则l加1，进入下一步，计算μ₁、λ₁、α₁待拟合参数的残差，即衰减系数残差Δμ、周期残差Δλ、初始相位残差Δα，并对μ₁、λ₁、α₁进行补偿，即：

μ₂＝μ₁+Δμ

λ₂＝λ₁+Δλ

α₂＝α₁+Δα

将补偿后的μ₂、λ₂、α₂代入步骤(3.4)中公式，计算待拟合残差b_k，并进行判断，若满足要求，则停止自相关函数拟合，μ、λ、α的取值为：μ＝μ₂、λ＝λ₂、α＝α₂。若不满足要求，则l加1，进入下一步，如此循环迭代，直至满足步骤(3.4)中的要求。

(4)、利用步骤(3)中求取的三个参数μ、λ、α，计算出任意角频率ω＝2πf处的功率谱密度，具体计算方法为：

其中：ω为角频率，ω＝2πf，f为频率。

上述计算出的功率谱密度和自相关函数能够估计载体的扰动特性，能够指导控制系统的设计，并给予参考。

实施例1

对某次惯性导航系统载体环境摸底试验数据为例，分析载体的扰动特性。

从原始加速度和陀螺仪中提取载体的运动信息，并进行导航计算。惯导系统采样频率f为200Hz。从开始测试到停止的过程中，提取得到的东向速度扰动序列

如图1所示，进行傅里叶变换后的线谱如图2所示，周期图法得到的信号功率谱密度如图3所示，无偏自相关曲线如图4所示。由上述附图可以看出，功率谱密度函数尖峰与跳变严重，无法直接从曲线中获得系统运动扰动的周期和幅值。

利用本发明对无偏自相关序列进行拟合，拟合过程迭代进行，最后四次的拟合过程如图6所示，图6a、6b、6c、6d分别为最后四次的拟合过程，拟合得到的自相关曲线如图7所示，然后对自相关函数进行傅里叶变换后得到的功率谱密度曲线如图8所示。

由图8可以看出：本发明计算得到的系统运动功率谱密度曲线准确、清晰且平滑，系统运动的主周期和最大幅值也可以从曲线中精确获得，方便后续对惯导系统进行定量分析。利用本发明的拟合方法对自相关函数的拟合结果与传统的拟合方法对自相关函数的拟合结果对比如图7所示，从图中可以看出，本发明的方法弥补了现有自相关函数拟合方法误差大的缺点，可以更好地对载体运动扰动特性进行描述。

以上所述，仅为本发明一个具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。

本发明未详细说明部分属于本领域技术人员公知常识。

以上所述，仅为本发明最佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。

本发明说明书中未作详细描述的内容属于本领域专业技术人员的公知技术。

Claims (8)

1.一种扰动运动统计特性的解析计算方法，其特征在于：包括

提取载体运动信息的等间隔周期的时间序列，根据所述时间序列和期望时间序列求出扰动时间序列；

根据所述扰动时间序列求出自相关序列；

将所述自相关序列与自相关解析模型进行拟合，得到自相关解析模型的参数；

根据所述自相关解析模型的参数，求出任意角频率处的功率谱密度；

所述提取载体运动信息的等间隔周期的时间序列为x(k)，k＝0,1,…,N-1，N为正整数；

所述期望时间序列为

N为正整数；

所述扰动时间序列为

计算公式如下：

根据所述扰动时间序列求出自相关序列r(k)，具体计算公式如下：

其中：l为非负整数；

将所述自相关序列与自相关解析模型进行拟合，得到自相关解析模型的参数μ、λ、α，其中μ为衰减系数，λ为周期、α为初始相位；

所述自相关解析模型为：

其中：t为时间；r(0)为自相关序列中第一个值；

根据所述自相关解析模型的参数，求出任意角频率处的功率谱密度的具体方法如下：

其中：ω为角频率，ω＝2πf，f为频率。

2.根据权利要求1所述的解析计算方法，其特征在于：将所述自相关序列与自相关解析模型进行拟合，得到自相关解析模型的参数μ、λ、α，包括：

(1)、对自相关序列r(k)进行归一化处理，得到归一化处理后的自相关序列

(2)、从自相关序列

中截取部分自相关序列

作为拟合的观测数据；

(3)、计算当前迭代次数为l时的自相关序列

的待拟合残差b_k，并建立向量矩阵B，判断向量矩阵B是否满足要求，或者当前迭代次数l是否为迭代次数上限，若满足，则μ、λ、α的取值分别为当前迭代次数为l时的μ_l、λ_l、α_l；若不满足，l加1，进入步骤(4)；

(4)、计算上一次待拟合参数的残差，包括衰减系数残差Δμ、周期残差Δλ、初始相位残差Δα；

(5)、根据衰减系数残差Δμ、周期残差Δλ、初始相位残差Δα对上一次待拟合参数进行补偿；

(6)、重复(3)～(5)，利用补偿后的待拟合参数计算自相关序列的待拟合残差b_k，直至满足步骤(3)中的要求。

3.根据权利要求2所述的解析计算方法，其特征在于：所述步骤(1)中对自相关序列r(k)进行归一化处理，得到归一化处理后的自相关序列

具体公式如下：

4.根据权利要求3所述的解析计算方法，其特征在于：所述步骤(2)中从自相关序列

中截取部分序列

作为拟合的观测数据，其中k′＝0,1,…,N₂-1，N₂＝[(3t₁)/T_s]，[]表示取整；t₁为自相关序列随时间第一次负向时刻，T_s为自相关序列的采样周期。

5.根据权利要求4所述的解析计算方法，其特征在于：所述步骤(3)中计算当前迭代次数为l时的自相关序列

的待拟合残差b_k，具体公式如下：

其中：k′＝0,1,…,N₂-1，l为当前迭代次数；

记向量B为：

若向量B满足B^TB<10^-2，或者满足l＝D，D为迭代次数上限，则：

μ＝μ_l

λ＝λ_l

α＝α_l

否则，l自加1，进入步骤(4)。

6.根据权利要求5所述的解析计算方法，其特征在于：所述步骤(3)中拟合初值的选择如下：

7.根据权利要求4～6之一所述的解析计算方法，其特征在于：所述步骤(4)中计算上一次待拟合参数的残差，包括衰减系数残差Δμ、周期残差Δλ、初始相位残差Δα具体方法如下；

(4.1)计算自相关序列残差的偏导参数：

(4.2)记矩阵A为：

(4.3)采用如下公式计算Δ：

Δ＝(A^TA)^-1A^TB

记向量Δ为：

Δ＝[Δμ Δλ Δα]^T。

8.根据权利要求7所述的解析计算方法，其特征在于：所述步骤(5)中根据衰减系数残差Δμ、周期残差Δλ、初始相位残差Δα对上一次待拟合参数进行补偿的方法如下：

μ_l＝μ_l-1+Δμ

λ_l＝λ_l-1+Δλ

α_l＝α_l-1+Δα；

其中：μ_l-1、λ_l-1、α_l-1分别为上一次的衰减系数、上一次的周期和上一次的初始相位。

Images