CN106909736B - 一种扰动运动统计特性的解析计算方法 - Google Patents
一种扰动运动统计特性的解析计算方法 Download PDFInfo
- Publication number
- CN106909736B CN106909736B CN201710103753.9A CN201710103753A CN106909736B CN 106909736 B CN106909736 B CN 106909736B CN 201710103753 A CN201710103753 A CN 201710103753A CN 106909736 B CN106909736 B CN 106909736B
- Authority
- CN
- China
- Prior art keywords
- autocorrelation
- sequence
- disturbance
- residual
- fitted
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Complex Calculations (AREA)
Abstract
本发明涉及一种扰动运动统计特性的解析计算方法,该方法经过扰动信息提取、自相关序列计算、自相关函数拟合、功率谱密度计算等步骤,定量给出载体运动的自相关和功率谱密度函数的具体参数,可以直接计算外部扰动频率点的功率谱密度,克服了以前功率谱密度尖峰与跳变严重,对扰动特性描述不精确的缺点,实现了对载体运动扰动特性的量化描述,精度高、数据平滑;本发明方法能够通过数值计算方法,准确、简单、定量计算出载体运动的扰动特性。
Description
技术领域
本发明涉及一种扰动运动统计特性的解析计算方法,属于力学环境试验技术领域。
背景技术
扰动运动的谱分析可以通过对扰动信号直接进行傅里叶变换实现,根据幅值谱大小及其对应频率得到扰动特性。但是这种方法得到的幅值谱曲线存在较大的尖峰与跳变(如图2所示),使得扰动特性结果不够准确。此后人们通过扰动信号的功率谱密度分析来获得扰动特性,即通过求取信号的功率谱密度,根据功率谱密度幅值及其对应频率得到扰动特性。这一方法在工程中获得了广泛的应用。
功率谱密度求取方法的主要包括周期图法、相关函数法和ARMA模型。
(1)由周期图法计算功率谱密度
周期图法首先将扰动信号进行分段,之后对各段扰动信号进行傅里叶变换,得到其幅值谱密度X(jω),最后通过公式:
得到信号的功率谱密度Φ(ω)。其中,N为信号分段总数,T为每段信号的总时长。但是,该方法的缺点是受噪声影响谱密度序列为离散点,由其构成的曲线不平滑。
(2)由自相关序列计算功率谱密度
相关函数法首先求取信号的自相关序列,之后对其进行傅里叶变换得到信号的功率谱密度。由于实际中信号的采样总是离散且时间有限的,因此直接使用公式:
得到的无偏自相关序列与信号自相关函数存在较大偏差,且与周期图法计算的功率谱基本相同,都存在波动量偏大的问题。
(3)基于ARMA模型计算功率谱密度
为减小功率谱的波动量和离散点问题,在现代随机信号处理中采用ARMA模型来求解功率谱密度,但其缺点是ARMA模型的阶次非常高,不利于数据分析。
为此,需要研究一种既具有平滑的优点,又具有简化的优点的扰动运动统计特性的计算方法。
发明内容
本发明的目的在于克服现有技术的上述缺陷,提供一种扰动运动统计特性的解析计算方法,该方法能够通过数值计算方法,准确、简单、定量计算出载体运动的扰动特性。
本发明的上述目的主要是通过如下技术方案予以实现的:
一种扰动运动统计特性的解析计算方法,包括
提取载体运动信息的等间隔周期的时间序列,根据所述时间序列和期望时间序列求出扰动时间序列;
根据所述扰动时间序列求出自相关序列;
将所述自相关序列与自相关解析模型进行拟合,得到自相关解析模型的参数;
根据所述自相关解析模型的参数,求出任意角频率处的功率谱密度。
在上述解析计算方法中,所述提取载体运动信息的等间隔周期的时间序列为x(k),k=0,1,…,N-1,N为正整数;
在上述解析计算方法中,根据所述扰动时间序列求出自相关序列r(k),具体计算公式如下:
其中:l为非负整数。
在上述解析计算方法中,将所述自相关序列与自相关解析模型进行拟合,得到自相关解析模型的参数μ、λ、α,其中μ为衰减系数,λ为周期、α为初始相位;
所述自相关解析模型为:
其中:t为时间;r(0)为自相关序列中第一个值。
在上述解析计算方法中,将所述自相关序列与自相关解析模型进行拟合,得到自相关解析模型的参数μ、λ、α,包括:
(3)、计算当前迭代次数为l时的自相关序列的待拟合残差bk,并建立向量矩阵B,判断向量矩阵B是否满足要求,或者当前迭代次数l是否为迭代次数上限,若满足,则μ、λ、α的取值分别为当前迭代次数为l时的μl、λl、αl;若不满足,l加1,进入步骤(4);
(4)、计算上一次待拟合参数的残差,包括衰减系数残差Δμ、周期残差Δλ、初始相位残差Δα;
(5)、根据衰减系数残差Δμ、周期残差Δλ、初始相位残差Δα对上一次待拟合参数进行补偿;
(6)、重复(3)~(5),利用补偿后的待拟合参数计算自相关序列的待拟合残差bk,直至满足步骤(3)中的要求。
在上述解析计算方法中,所述步骤(2)中从自相关序列中截取部分序列作为拟合的观测数据,其中k′=0,1,…,N2-1,N2=[(3t1)/Ts],[]表示取整;t1为自相关序列随时间第一次负向时刻,Ts为自相关序列的采样周期。
其中:k′=0,1,…,N2-1,l为当前迭代次数;
记向量B为:
若向量B满足BTB<10-2,或者满足l=D,D为迭代次数上限,则:
μ=μl
λ=λl
α=αl
否则,l自加1,进入步骤(4);
在上述解析计算方法中,所述步骤(3)中拟合初值的选择如下:
在上述解析计算方法中,所述步骤(4)中计算上一次待拟合参数的残差,包括衰减系数残差Δμ、周期残差Δλ、初始相位残差Δα具体方法如下;
(4.1)计算自相关序列残差的偏导参数:
(4.2)记矩阵A为:
(4.3)采用如下公式计算Δ:
Δ=(ATA)-1ATB
记向量Δ为:
Δ=[Δμ Δλ Δα]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所示为本发明解析计算方法的流程图,本发明提出的惯导设备载体运动扰动特性确定方法,具体步骤如下:
(2)、对扰动序列计算其自相关序列r(k),计算公式为
其中:l为非负整数。
(3)、自相关解析模型如下:
其中:t为时间;r(0)为自相关序列中第一个值。
将自相关序列与自相关解析模型进行拟合,得到自相关解析模型的参数μ、λ、α,其中μ为衰减系数,λ为周期、α为初始相位;自相关序列的拟合方法如下:
(3.3)、选择拟合初值,具体方法为:
其中:t1为自相关序列随时间第一次负向时刻。
其中:k′=0,1,…,N2-1,l为当前迭代次数;
记向量B为:
若向量B满足BTB<10-2,或者满足l=D,D为迭代次数上限,则:
μ=μl
λ=λl
α=αl
停止自相关函数拟合。否则,l自加1,进入步骤(3.5)。
(3.5)、计算上一次待拟合参数的残差,包括衰减系数残差Δμ、周期残差Δλ、初始相位残差Δα,具体方法如下:
对截断后的自相关序列中的每一个点,计算以下三个参数(自相关序列残差的偏导参数),
记矩阵A为:
采用如下公式计算Δ:
Δ=(ATA)-1ATB
记向量Δ为:
Δ=[Δμ Δλ Δα]T。
(3.6)、根据衰减系数残差Δμ、周期残差Δλ、初始相位残差Δα对上一次待拟合参数进行补偿,具体方法为:
μl=μl-1+Δμ
λl=λl-1+Δλ
αl=αl-1+Δα。
(3.7)、重复(3.4)~(3.6)进行迭代计算,即利用补偿后的待拟合参数计算自相关序列的待拟合残差bk,直至步骤(3)中的任意条件满足。
若步骤(3.4)中选择初值μ1=0,α1=0,首次进行计算时即满足BTB<10-2,或者满足l=D,则μ、λ、α的取值为:μ=μ1、λ=λ1、α=α1停止自相关函数拟合。若不满足条件,则l加1,进入下一步,计算μ1、λ1、α1待拟合参数的残差,即衰减系数残差Δμ、周期残差Δλ、初始相位残差Δα,并对μ1、λ1、α1进行补偿,即:
μ2=μ1+Δμ
λ2=λ1+Δλ
α2=α1+Δα
将补偿后的μ2、λ2、α2代入步骤(3.4)中公式,计算待拟合残差bk,并进行判断,若满足要求,则停止自相关函数拟合,μ、λ、α的取值为:μ=μ2、λ=λ2、α=α2。若不满足要求,则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为正整数;
根据所述扰动时间序列求出自相关序列r(k),具体计算公式如下:
其中:l为非负整数;
将所述自相关序列与自相关解析模型进行拟合,得到自相关解析模型的参数μ、λ、α,其中μ为衰减系数,λ为周期、α为初始相位;
所述自相关解析模型为:
其中:t为时间;r(0)为自相关序列中第一个值;
根据所述自相关解析模型的参数,求出任意角频率处的功率谱密度的具体方法如下:
其中:ω为角频率,ω=2πf,f为频率。
2.根据权利要求1所述的解析计算方法,其特征在于:将所述自相关序列与自相关解析模型进行拟合,得到自相关解析模型的参数μ、λ、α,包括:
(3)、计算当前迭代次数为l时的自相关序列的待拟合残差bk,并建立向量矩阵B,判断向量矩阵B是否满足要求,或者当前迭代次数l是否为迭代次数上限,若满足,则μ、λ、α的取值分别为当前迭代次数为l时的μl、λl、αl;若不满足,l加1,进入步骤(4);
(4)、计算上一次待拟合参数的残差,包括衰减系数残差Δμ、周期残差Δλ、初始相位残差Δα;
(5)、根据衰减系数残差Δμ、周期残差Δλ、初始相位残差Δα对上一次待拟合参数进行补偿;
(6)、重复(3)~(5),利用补偿后的待拟合参数计算自相关序列的待拟合残差bk,直至满足步骤(3)中的要求。
8.根据权利要求7所述的解析计算方法,其特征在于:所述步骤(5)中根据衰减系数残差Δμ、周期残差Δλ、初始相位残差Δα对上一次待拟合参数进行补偿的方法如下:
μl=μl-1+Δμ
λl=λl-1+Δλ
αl=αl-1+Δα;
其中:μl-1、λl-1、αl-1分别为上一次的衰减系数、上一次的周期和上一次的初始相位。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710103753.9A CN106909736B (zh) | 2017-02-24 | 2017-02-24 | 一种扰动运动统计特性的解析计算方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710103753.9A CN106909736B (zh) | 2017-02-24 | 2017-02-24 | 一种扰动运动统计特性的解析计算方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN106909736A CN106909736A (zh) | 2017-06-30 |
CN106909736B true CN106909736B (zh) | 2020-07-14 |
Family
ID=59208962
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201710103753.9A Active CN106909736B (zh) | 2017-02-24 | 2017-02-24 | 一种扰动运动统计特性的解析计算方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN106909736B (zh) |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US5818969A (en) * | 1995-05-12 | 1998-10-06 | Intel Corporation | Intelligent start for motion estimation search |
CN103245830A (zh) * | 2013-04-03 | 2013-08-14 | 云南电力试验研究院(集团)有限公司电力研究院 | 一种结合ar谱估计与非线性优化的间谐波检测方法 |
CN103578121A (zh) * | 2013-11-22 | 2014-02-12 | 南京信大气象装备有限公司 | 干扰运动环境下基于共享高斯模型的运动检测方法 |
-
2017
- 2017-02-24 CN CN201710103753.9A patent/CN106909736B/zh active Active
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US5818969A (en) * | 1995-05-12 | 1998-10-06 | Intel Corporation | Intelligent start for motion estimation search |
CN103245830A (zh) * | 2013-04-03 | 2013-08-14 | 云南电力试验研究院(集团)有限公司电力研究院 | 一种结合ar谱估计与非线性优化的间谐波检测方法 |
CN103578121A (zh) * | 2013-11-22 | 2014-02-12 | 南京信大气象装备有限公司 | 干扰运动环境下基于共享高斯模型的运动检测方法 |
Non-Patent Citations (3)
Title |
---|
城市电网规划辅助决策系统负荷预测子系统的开发和改进;杜柏均;《 中国优秀硕士学位论文全文数据库 工程科技Ⅱ辑》;20120215(第02期);论文第15-16页 * |
机载气象雷达信号处理若干关键技术研究;卢晓光;《中国博士学位论文全文数据库 信息科技辑》;20151115(第11期);论文第67-85页 * |
高精度航天器微振动力学环境分析;张振华 等;《航天器环境工程》;20091231;第26卷(第6期);第528-531页 * |
Also Published As
Publication number | Publication date |
---|---|
CN106909736A (zh) | 2017-06-30 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Kawahara et al. | Tandem-STRAIGHT: A temporally stable power spectral representation for periodic signals and applications to interference-free spectrum, F0, and aperiodicity estimation | |
CN107688554B (zh) | 基于自适应傅里叶分解的频域辨识方法 | |
CN107167306B (zh) | 基于阶次提取的旋转机械转子运行状态模态分析方法 | |
Greś et al. | Kalman filter-based subspace identification for operational modal analysis under unmeasured periodic excitation | |
CN103018555B (zh) | 一种高精度的电力参数软件同步采样方法 | |
CN110441018B (zh) | 火工冲击响应数据时频分析方法 | |
CN110530509B (zh) | 基于最大熵谱分析的高速动车组轴箱振动主频预测方法 | |
CN105307095B (zh) | 一种基于fft的高分辨率音频频率测量方法 | |
CN103218482B (zh) | 一种动力学系统中不确定参数的估计方法 | |
CN105137373A (zh) | 一种指数信号的去噪方法 | |
CN109374966A (zh) | 一种电网频率估计方法 | |
CN105021210A (zh) | Mems陀螺仪随机漂移误差的处理方法 | |
CN114329347B (zh) | 一种电能表计量误差预测方法、装置及存储介质 | |
CN106909736B (zh) | 一种扰动运动统计特性的解析计算方法 | |
Guschina et al. | The Exact Frequency Domain Solution for the Periodic Synchronous Averaging Performed in Discrete-Time | |
CN104535082A (zh) | 一种基于飞行试验和理论计算判断惯导元件性能的方法 | |
CN112595889B (zh) | 非理想多指数衰减正弦信号的欠Nyquist采样与参数测量方法 | |
CN112541157B (zh) | 一种信号频率精确估计方法 | |
CN108181650B (zh) | 一种数据处理方法和装置 | |
CN109711036B (zh) | 飞行控制系统试验结果的评估方法 | |
RU2306592C1 (ru) | Способ активной идентификации линейных объектов управления | |
Yang et al. | Empirical analysis of the MIMU random-drift error based on wavelet soft-threshold and ARIMA model | |
CN111781643B (zh) | 地震层速度确定方法及装置 | |
Wang et al. | Correlation coefficient stationary series method for gyroscope random drift | |
CN105372493A (zh) | 基于三条dft复数谱线的信号幅值和相位测量方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |