CN110542406B - 基于emd-mpf改进的陀螺仪信号去噪方法 - Google Patents
基于emd-mpf改进的陀螺仪信号去噪方法 Download PDFInfo
- Publication number
- CN110542406B CN110542406B CN201910826514.5A CN201910826514A CN110542406B CN 110542406 B CN110542406 B CN 110542406B CN 201910826514 A CN201910826514 A CN 201910826514A CN 110542406 B CN110542406 B CN 110542406B
- Authority
- CN
- China
- Prior art keywords
- function
- intrinsic mode
- signal
- noise
- mode function
- 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
- 238000000034 method Methods 0.000 title claims abstract description 135
- 230000009467 reduction Effects 0.000 claims abstract description 13
- 238000012545 processing Methods 0.000 claims abstract description 10
- 230000006870 function Effects 0.000 claims description 131
- 239000002245 particle Substances 0.000 claims description 92
- 238000005311 autocorrelation function Methods 0.000 claims description 28
- 230000008569 process Effects 0.000 claims description 23
- 238000012952 Resampling Methods 0.000 claims description 20
- 238000000354 decomposition reaction Methods 0.000 claims description 15
- 238000009826 distribution Methods 0.000 claims description 14
- 238000005070 sampling Methods 0.000 claims description 9
- 238000005314 correlation function Methods 0.000 claims description 6
- 230000003595 spectral effect Effects 0.000 claims description 6
- 238000004364 calculation method Methods 0.000 claims description 5
- 238000010606 normalization Methods 0.000 claims description 4
- 230000000717 retained effect Effects 0.000 claims description 4
- 230000002194 synthesizing effect Effects 0.000 claims description 3
- 230000008901 benefit Effects 0.000 abstract description 6
- 238000001914 filtration Methods 0.000 description 25
- 238000004458 analytical method Methods 0.000 description 11
- 230000002829 reductive effect Effects 0.000 description 10
- 238000005259 measurement Methods 0.000 description 9
- 230000000694 effects Effects 0.000 description 6
- 230000003068 static effect Effects 0.000 description 6
- 238000004422 calculation algorithm Methods 0.000 description 5
- 230000000875 corresponding effect Effects 0.000 description 5
- 230000006872 improvement Effects 0.000 description 5
- 230000003044 adaptive effect Effects 0.000 description 4
- 238000012360 testing method Methods 0.000 description 4
- 238000000342 Monte Carlo simulation Methods 0.000 description 3
- 238000012216 screening Methods 0.000 description 3
- 238000000540 analysis of variance Methods 0.000 description 2
- 230000002238 attenuated effect Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 238000004519 manufacturing process Methods 0.000 description 2
- 238000002156 mixing Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 238000010187 selection method Methods 0.000 description 2
- 230000007704 transition Effects 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 241001123248 Arma Species 0.000 description 1
- 230000001133 acceleration Effects 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000015556 catabolic process Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000002596 correlated effect Effects 0.000 description 1
- 238000006731 degradation reaction Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000004069 differentiation Effects 0.000 description 1
- 230000004927 fusion Effects 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 230000002401 inhibitory effect Effects 0.000 description 1
- 238000009434 installation Methods 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 230000004807 localization Effects 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 230000036961 partial effect Effects 0.000 description 1
- 238000002360 preparation method Methods 0.000 description 1
- 238000004445 quantitative analysis Methods 0.000 description 1
- 238000011867 re-evaluation Methods 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 238000012731 temporal analysis Methods 0.000 description 1
- 238000011426 transformation method Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C25/00—Manufacturing, calibrating, cleaning, or repairing instruments or devices referred to in the other groups of this subclass
- G01C25/005—Manufacturing, calibrating, cleaning, or repairing instruments or devices referred to in the other groups of this subclass initial alignment, calibration or starting-up of inertial devices
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C9/00—Measuring inclination, e.g. by clinometers, by levels
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/02—Preprocessing
- G06F2218/04—Denoising
- G06F2218/06—Denoising by applying a scale-space analysis, e.g. using wavelet analysis
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Manufacturing & Machinery (AREA)
- Gyroscopes (AREA)
Abstract
本发明公开了一种基于EMD‑MPF改进的陀螺仪信号去噪方法,包括:将陀螺仪的有噪信号分解为本征模态函数和残差信号;通过确定的两个标识参数对本征模态函数进行阶次选择,将本征模态函数分为噪声本征模态函数、混合本征模态函数和信息本征模态函数;舍弃所述噪声本征模态函数,保留所述信息本征模态函数,并对所述混合本征模态函数进行降噪处理;对降噪处理后的混合本征模态函数和保留的所述信息本征模态函数进行信号重构,从而得到去噪的陀螺仪信号。实现精确去除噪声信号,从而提高导航精确度的优点。
Description
技术领域
本发明涉及信号去噪领域,具体地,涉及一种基于EMD-MPF改进的陀螺仪信号去噪方法。
背景技术
由于微机电***(Micro-Electro-Mechanical System,MEMS)的发展,可以制造出更小尺寸和更廉价的加速度计和陀螺仪,以此组成成本更低惯性测量单元(InertialMeasurement Unit,IMU),由于其具有体积小、重量轻、功耗低、可大批量生产等优势,使得惯性导航或结合惯性测量单元的定位、导航技术迅速发展,并在交通、测量、航空航天等众多领域得到广泛应用。MEMS陀螺仪的误差大小在很大程度上会决定***的实际性能,因此对陀螺仪各类误差的抑制具有重要的实际意义。陀螺仪可能产生的误差包括了安装误差、杆臂误差、刻度因子误差等;而噪声则主要包括了***噪声和量测噪声两大方面。其中,MEMS陀螺仪的随机漂移则是噪声部分的重要来源,随机漂移不论对于陀螺仪自身还是对于整个导航***的精度影响都是很大,由于设备的长时间工作,轻微的漂移在累积或是***消息融合之后都会导致较大的误差。所以,如何降低噪声已成为MEMS陀螺信号处理中最重要的问题之一。
MEMS陀螺仪的信号的去噪方法主要分为建模去噪方法和非建模去噪方法两大类:第一,广泛采用成熟的随机漂移模型或是采用ARMA进行建模,然后利用卡尔曼滤波或各类改进的卡尔曼滤波技术对信号中的噪声进行过滤;第二,直接使用小波分解或是经验模态分解(Empirical Mode Decomposition,EMD)方法对信号进行分析、去除其中的噪声,以此提取有效信号。
传统的去噪方法一般使用低通、高通或带通滤波器来消除噪声。例如使用数字低通滤波器消除捷联惯性导航***(Inertial Navigation System, INS)中的高频噪声;使用FIR低通滤波器过滤和抑制MEMS传感器的量测噪声;应用高通滤波器针对陀螺仪抖动频带中具有非常高的衰减并且在通带中具有线性相位响应的现象,对相应的信号抖动噪声进行过滤。但是传统的去噪方法是基于经典滤波理论的,其只适合用于信号和噪声没有重叠的情况,所以对于实际信号滤波过后的频率成分中仍然是含有白噪声的。对于现代滤波理论,高斯滤波、维纳滤波、卡尔曼滤波等是其的代表方法,这些方法均需要利用信号和噪声所统计的先验知识,例如使用了自适应鲁棒卡尔曼滤波方法,并使用创新序列的加权协方差来调整测量噪声协方差矩阵,以此来减小陀螺仪中的随机噪声;应用一种随机加权方法下的抗野值递推最小二乘法(RLS)自适应滤波算法来降低MEMS陀螺仪输出当中的噪声分量;使用自适应扩展卡尔曼滤波方法补偿侧滑角估计,以此降低由惯性传感器偏移引起的潜在大漂移估计误差。这些方法克服了经典滤波理论的缺点,但是在使用过程中需要给出一定的最有标准,并需要知道噪声的先验统计信息,而这些在实际中通常很难或不可能获得,一般只能考一个最优估计进行代替。同时,这些方法需要精准建立模型,模型的偏差同样会影响到最后的去噪效果。
在此基础上继续发展出了非建模的去噪方法。小波变换(wavelet transform,WT)是一种新的变换分析方法,它继承和发展了短时傅立叶变换局部化的思想,同时又克服了窗口大小不随频率变化等缺点,所以具有良好的多分辨率分析特性,可以对高频处时间细分,低频处频率细分,以此达到对不同频段和时间的信号进行有效分析。利用改进的小波变换方法对陀螺仪误差特性进行分析,并建立新的软阈值函数以此抑制陀螺仪的测量噪声;应用了第二代小波变换针对MEMS陀螺进行去噪。但是,小波变换分析方法需要针对不同情况选择不同的分解层数、小波基,去噪还要选择合适的阈值,这对使用十分不便;同时,WT方法更适用于线性、稳态信号。从而存在对噪声信息去除不精细,从而造成导航精确度低的问题。
发明内容
本发明的目的在于,针对上述问题,提出一种基于EMD-MPF改进的陀螺仪信号去噪方法,以实现精确去除噪声信号,从而提高导航精确度的优点。
为实现上述目的,本发明实施例采用的技术方案是:
一种基于EMD-MPF改进的陀螺仪信号去噪方法,包括:
将陀螺仪的有噪信号分解为本征模态函数和残差信号;
通过确定的两个标识参数对本征模态函数进行阶次选择,将本征模态函数分为噪声本征模态函数、混合本征模态函数和信息本征模态函数;
舍弃所述噪声本征模态函数,保留所述信息本征模态函数,并对所述混合本征模态函数进行降噪处理;
对降噪处理后的混合本征模态函数和保留的所述信息本征模态函数进行信号重构,从而得到去噪的陀螺仪信号。
作为本发明实施例的一种具体实现方式,所述将陀螺仪的有噪信号分解为本征模态函数和残差信号为:
将陀螺仪的有噪信号分解为多个本征模态函数和一个残差信号。
作为本发明实施例的一种具体实现方式,所述对降噪处理后的混合本征模态函数和保留的所述信息本征模态函数进行信号重构中的重构公式为:
作为本发明实施例的一种具体实现方式,所述MPF为改进的PF方法,包括:
把预测粒子样本的概率分布记作p(xi),1≤i≤N,把权重wi>w0的粒子样本xT的概率分布记作p1(xi),1≤i≤T,即为p,把剩余粒子样本xN-T的概率分布记作q(xi),N-T<i≤N,即为q,wi为权重,w0为权值阈值,N和T均为常数;
计算不同粒子样本之间的K-L散度如下:
对于小权重的粒子样本,产生N-T个随机数uj,其范围是uj∈(w0,1],对该粒子样本的权重进行更新,DKL为K-L散度:
作为本发明实施例的一种具体实现方式,所述确定的两个标识参数中,所述两个标识参数分别为M1和M2,选定确定标识参数M1,包括:使用自相关函数方差确定标识参数M1。
作为本发明实施例的一种具体实现方式,所述使用自相关函数方差确定标识参数M1,包括:
根据不同阶次的本征模态函数信号,计算各自对应的自相关函数;
根据连续本征模态函数的自相关函数计算两个自相关的函数之间的欧几里得距离均方差,得到计算结果MSED;
根据MSED确定标识参数M1。
作为本发明实施例的一种具体实现方式,所述标识参数M1与数值k和数值L相关联,
所述1≤k≤[2L/3],[2L/3]表示不超过2L/3的最大整数,从而保证M1的取值过大至L,k和L均为自然数。
作为本发明实施例的一种具体实现方式,确定标识参数M2,包括:
计算原始信号的概率密度函数与每个本征模态函数的概率密度函数之间的巴塔恰里雅距离,其中0≤BC≤1和0≤BD≤∞,BD为巴塔恰里雅距离,BC为巴塔恰里雅系数;
基于功率谱密度函数计算所述原始信号的概率密度函数和每个所述本征模态函数的概率密度函数之间的巴塔恰里雅距离数据在不同特征下的可分离性;
基于所述可分离性确定标识参数M2。
作为本发明实施例的一种具体实现方式,所述通过确定的两个标识参数对本征模态函数进行阶次选择,将本征模态函数分为噪声本征模态函数、混合本征模态函数和信息本征模态函数,包括:
1<i<M1时,为噪声本征模态函数;
M1<i<M2时,为混合本征模态函数;
M2<i<L时,为信息本征模态函数,i表示本征态函数的阶次
L表示本征态函数的总阶数。
本发明的实施例具有以下有益效果:
下面通过附图和实施例,通过将有噪信号分解为本征模态函数和残差信号,然后通过两个标识参数本征模态函数分为噪声本征模态函数、混合本征模态函数和信息本征模态函数;然后将舍弃掉噪声本征模态函数,降噪处理后的混合本征模态函数与信息本征模态函数进行信号重构,精确去除了噪声信号,达到提高导航精确度的优点。
附图说明
图1为本发明实施例所述的基于EMD-MPF改进的陀螺仪信号去噪方法的流程图;
图2为本发明实施例所述的EMD-MPF去噪方法框图;
图3为本发明实施例所述的MEMS陀螺仪X轴静态数据;
图4为本发明实施例所述的X轴陀螺数据的EMD分解图;
图5为本发明实施例所述的各阶分解信号的自相关函数;
图6为本发明实施例所述的MSED分布情况;
图7为本发明实施例所述的分解信号的概率密度函数;
图8为本发明实施例所述的巴氏距离分布情况;
图9为本发明实施例所述的混合IMF利用MPF的滤波结果;
图10为本发明实施例所述的X轴静态试验结果对比图;
图11为本发明实施例所述的Allen方差分析结果。
具体实施方式
以下结合附图对本发明的优选实施例进行说明,应当理解,此处所描述的优选实施例仅用于说明和解释本发明,并不用于限定本发明。
本实施例的用于MEMS陀螺信号的EMD-MPF去噪方法。通过经验模态分解将原始信号分解为若干个本征模态函数,然后通过分析和阶次选择把IMF分为噪声IMF、混合IMF和信息IMF,最后通过消除噪声和改进粒子滤波的方法得到最终的去噪信号。
EMD方法可以针对非线性、非平稳信号的自适应信号分解算法,其不仅突破了Fourier变换的局限性,而且不存在如小波变换一样需要预选小波基函数的问题,具有良好的时频分辨率和自适应性。基于经验模态分解,使用完全数据驱动的方法对信号部分重建,以此达到信号滤波的效果;应用基于经验模态分解和Allan方差的时频域FOG信号去噪方法,有效地提高了导航精度。
如图1所示,一种基于EMD-MPF改进的陀螺仪信号去噪方法,包括:
S101:将陀螺仪的有噪信号分解为本征模态函数和残差信号;
S102:通过确定的两个标识参数对本征模态函数(IMF)进行阶次选择,将本征模态函数分为噪声本征模态函数、混合本征模态函数和信息本征模态函数;
S103:舍弃所述噪声本征模态函数,保留所述信息本征模态函数,并对所述混合本征模态函数进行降噪处理;
S104:对降噪处理后的混合本征模态函数和保留的所述信息本征模态函数进行信号重构,从而得到去噪的陀螺仪信号。
作为本发明实施例可选的一种实现方式,所述将陀螺仪的有噪信号分解为本征模态函数和残差信号为:
将陀螺仪的有噪信号分解为多个本征模态函数和一个残差信号。
作为本发明实施例的一种具体实现方式,所述对降噪处理后的混合本征模态函数和保留的所述信息本征模态函数进行信号重构中的重构公式为:
作为本发明实施例可选的一种实现方式,所述MPF为改进的PF方法,包括:
把预测粒子样本的概率分布记作p(xi),1≤i≤N,把权重wi>w0的粒子样本xT的概率分布记作p1(xi),1≤i≤T,即为p,把剩余粒子样本xN-T的概率分布记作q(xi),N-T<i≤N,即为q,wi为权重,w0为权值阈值,N和T均为常数;
计算不同粒子样本之间的K-L散度如下:
引入K-L散度可以较好的保证粒子样本集的有效性和多样性;
对于小权重的粒子样本,产生N-T个随机数uj,其范围是uj∈(w0,1],uj是利用程序中的随机函数random,在范围(w0,1]中产生,DKL为K-L散度,对该粒子样本的权重进行更新:
作为本发明实施例可选的一种实现方式,所述确定的两个标识参数中,所述两个标识参数分别为M1和M2,选定确定标识参数M1,包括:使用自相关函数方差确定标识参数M1。
作为本发明实施例的可选的一种实现方式,所述使用自相关函数方差确定标识参数M1,包括:
根据不同阶次的本征模态函数信号,计算各自对应的自相关函数;
根据连续本征模态函数的自相关函数计算两个自相关的函数之间的欧几里得距离均方差,得到计算结果MSED;
根据MSED确定标识参数M1。
作为本发明实施例可选的一种实现方式,所述标识参数M1与数值k和数值L相关联,
所述1≤k≤[2L/3],[2L/3]表示不超过2L/3的最大整数,从而保证M1的取值过大至L,k和L均为自然数。
作为本发明实施例可选的一种实现方式,确定标识参数M2,包括:
计算原始信号的概率密度函数与每个本征模态函数的概率密度函数之间的巴塔恰里雅距离,其中0≤BC≤1和0≤BD≤∞,BD为巴塔恰里雅距离, BC为巴塔恰里雅系数;
基于功率谱密度函数计算所述原始信号的概率密度函数和每个所述本征模态函数的概率密度函数之间的巴塔恰里雅距离数据在不同特征下的可分离性;
基于所述可分离性确定标识参数M2。
作为本发明实施例可选的一种实现方式,所述通过确定的两个标识参数对本征模态函数进行阶次选择,将本征模态函数分为噪声本征模态函数、混合本征模态函数和信息本征模态函数,包括:
1<i<M1时,为噪声本征模态函数;
M1<i<M2时,为混合本征模态函数;
M2<i<L时,为信息本征模态函数,i表示本征态函数的阶次
L表示本征态函数的总阶数。
在一个具体的应用场景中,
1、改进的EMD-MPF方法:
EMD是一种适用于处理非平稳非线性序列的自适应的时空分析方法。该方法通过筛选的迭代过程,将信号分解为一系列不偏离时间域的本征模态函数(IMF),因此具有良好的时频分辨率和自适应性,能够完美地重构原始信号,同时具有突出信号中可能被忽视的精细地质构造的潜能。目前EMD方法已经成功的应用于降低信号噪声的领域当中。其中,本征模态函数须满足以下两个要求:(1)极值数和零交叉数最多相等或相差一个;(2)在任何一点,由局部最大值定义的包络的平均值和由局部最小值定义的包络为零。对于有噪信号x(t)分解得到的数量为L的IMF信号h(i)(t) 和一个残差信号rL(t),则可以对其进行如下形式的信号重构:
粒子滤波(Particle Filter,PF)的思想基于蒙特卡洛方法(Monte Carlomethods),它是利用粒子集来表示概率,可以用在任何形式的状态空间模型上。该方法是一种通过寻找一组在状态空间中传播的随机样本来近似的表示概率密度函数,用样本均值代替积分运算,进而获得***状态的最小方差估计的过程。因为该方法在蒙特卡洛方法的基础上又加了一层重要性采样思想,所以可以用一组样本近似表示***的后验概率分布,进而可以使用这一近似结果估计非线性***的状态。对于PF方法的一般步骤如下:
步骤1,初始化状态量:由大量粒子按照先验概率分布在空间中模拟初始状态。
步骤2,预测阶段:根据状态转移方程,每一个粒子得到一个预测粒子。
步骤3,校正阶段:对预测粒子进行评价,状态越接近真实状态的粒子,则赋予其较大的权重。
步骤4,重采样阶段:根据粒子的权重对粒子进行筛选,在筛选过程中既要保留权重较大的粒子,又要保留一部分权重较小的粒子。
步骤5,滤波阶段:将重采样后的粒子代入状态转移方程的新的预测粒子,即返回步骤2进行循环。
传统的EMD方法将原始信号分解后得到一系列频率由高到低分布的IMF,一般认为噪声基本存在于高频段,所以会直接删除部分高频段的IMF 然后进行重组得到去噪信号。PF方法虽然很好地摆脱了解决非线性滤波问题时随机量必须满足高斯分布的制约,但是,其存在的问题是需要用大量的样本数量才能很好地近似***的后验概率密度,即需要大量具有多样性的有效粒子,而实际环境越复杂就需要越多的粒子,也就导致算法的复杂程度越高。
2、改进的IMF阶次选择方法:
EMD方法能够将含噪声的信号分解成若干个不同尺度的IMF分量和一个残留剩余。通常,已知噪声信号主要集中在高频部分,同时IMF的组成一般可以分为三个组成部分,即噪声IMF、混合IMF和信息IMF。现在,为了信号去噪而不是单纯的完全重建,因此主要问题是如何区分相关的 IMF和不相关的IMF,即如何找到不同IMF组成部分的分界点是该方法去噪的关键问题。根据上述内容,重新将噪声信号描述如下:
2.1、关于分界点M1的确定:
目前的问题是如何确定两个分界点M1和M2。对于第一个参数,使用自相关函数方差来确定。在实际应用当中,噪声信号中的实际信号分量必然不会被知道。换句话说,永远无法在理论上或实践中完全得知并去除噪声信号中的噪声。因此,既不能通过最小化目标信号和噪声信号的均方误差,也不能得到每个分量各自的信噪比来确定参数。IMF的重建信号如下:
自相关函数可以确定信号在不同时域下的相关程度,根据不同阶次的 IMF信号,各自对应的自相关函数计算如下:
然后,根据连续IMF的自相关函数计算两者之间的欧几里得距离均方差,其中Pea为相关系数,并将最后(5)式的计算结果简称为MSED,计算如下:
根据随机噪声的统计特性,只有在零点处随机噪声的自相关函数最大,在其他时刻其自相关函数则会快速衰减到零值附近。使用欧几里得距离是一种常用的分类依据,而引入自相关函数的方差,则是不仅考虑到了自相关函数下的无量纲化情况,同时依靠自相关函数考虑了不同分量之间时域上的相关性,一定程度上减少了欧几里得距离在单一特征下受不同分量的干扰。
最后,得到第一个指标参数M1,确定如下:
其中,k的值设定在1和[2L/3]之间,[2L/3]表示不超过2L/3的最大整数,如此设定是为了保证M1的取值过大至L。
通过上述分析,根据自相关特性,如果是噪声主体的IMF分量,则随机噪声在零点处有最大值,在其他时域会近似的快速下降为零;如果是混合信息的IMF分量,则随机噪声在其他时域处还会出现较大波动。同时,依据MSED计算、描述连续时域的IMF分量之间的差异,MSED值越小,则对应的h(i)(t)分量比重越小;而无论是混合IMF还是信息IMF,由于它们包含MEMS陀螺的大量有用信息,所以它们中对应的h(t)具有很高的能量。综上分析,MSED的计算值可以作为第一个表示参数。
2.2、关于分界点M2的确定:
分界点M2主要用于标识、区分信息IMF和混合IMF。为了确定第二个标识参数,引入了原始信号的概率密度函数(power spectral density function,PDF)与每个IMF的PDF之间的巴塔恰里雅距离(Bhattacharyya distance,BD)。Bhattacharyya距离和Bhattacharyya系数以20世纪30年代曾在印度统计研究所工作的一个统计学家A.Bhattacharya命名。有学者提出使用欧几里得距离作为分离标识,但是欧几里得距离会明显受样本维数的影响;马哈拉诺比斯距离的方法考虑了各种特征之间的关系,消除变量之间的相关干扰,是一种间接的区分方法,同时,马氏距离有可能会夸大变化微小的变量的作用,因此降低分界点的准确性。而巴氏距离常用于两离散样本之间的概率分布,以此判断测量类之间的可分离性,不仅排除了变量之间的相关性的干扰,而且是一种更直接的区分方法。针对离散概率分布p和q之间的巴塔恰里雅距离定义如下:
BD(p,q)=-ln[BC(p,q)] (9),
其中,BC被称为Bhattacharyya系数;并且满足0≤BC≤1和0≤BD≤∞。
功率谱密度函数定义了信号或者时间序列的功率是如何随频率分布的,是对数据特征的很好描述。因此,原始信号和每个IMF的PDF之间的巴塔恰里雅距离数据在不同特征下的可分离性,计算如下:
分析上式,在IMF中混合的噪声越多,则和原始信号相比,可分离性就越低;如果是“纯粹”的信息IMF,和原始信号相比,则可分离性就会越高。所以,M2的选取方法如下:
有一个观点,如果IMF分量和原始信号两者PDF之间的巴塔恰里雅距离越小,则可分离性越小,即IMF分量中包含的噪声越多;相反,两者 PDF之间的巴塔恰里雅距离越大,则相似性越小、越可以分离,即IMF 分量中包含的噪声越少。所以将巴塔恰里雅距离最大值的下一个点视为分界点,可以跟有效地区分混合IMF分量和信息IMF分量。通常,分界点M2所确定的IMF阶数要高于分界点M1所确定的阶数,但是有可能出现一种情况,即有噪信号所分解的IMF大部分阶次都是含有噪声的,从而导致M1在M2之后,所以对(10)式改进如下,以此保证标识参数M2所确定的IMF 阶数要高于M1所确定的阶数:
综上所述,M1和M2两个标识参数将EMD方法分解得到IMF分量分为了三个部分,其中,1<i<M1时,为噪声IMF分量,进行舍弃;M1<i<M2时,认为是混合IMF分量,使用下面改进的粒子滤波作进一步降噪处理; M2<i<L时,认为时信息IMF分量,进行保留。
3、改进的PF方法:
传统的重采样过程即在算法的每次迭代过程中,根据粒子权值对粒子进行重采样,舍弃了权值较小的粒子,代之以权值较大的粒子。这种重采样过程较好的解决了粒子滤波中原有的退化现象,但是也带来了样本枯竭的问题。所以结合K-L散度进行改进,提出了一种创新的重采样过程。设定重采样过程的权值阈值为w0,对应高于权重阈值的粒子数目为T,粒子总数为N,则传统重采样过程为保留权值高于w0的粒子,然后根据权重重采样生成N-T个粒子。这样就破坏了粒子样本集的有效性和多样性。现在,改进的重采样过程具体步骤如下:
1)把预测粒子样本的概率分布记作p(xi),1≤i≤N;把权重wi>w0的粒子样本xT的概率分布记作p1(xi),1≤i≤T;把剩余粒子样本xN-T的概率分布记作 q(xi),N-T<i≤N。
2)计算不同粒子样本之间的K-L散度如下:
3)针对小权重的粒子样本,产生N-T个随机数uj,其范围是uj∈(w0,1],对该粒子样本的权重进行更新:
5)综合两部分粒子,进行权值归一化,完成PF方法的继续迭代操作。
改进的重采样过程不仅对部分粒子重新评估了权重,而且利用K-L散度增加重新评估过程的有序性,让重采样过程中随机分布差别尽可能减小,同时也尽可能地保持了粒子的有效性和多样性,以此提高了重采样过程的效率。
4、改进方法的降噪***流程:
基于上述的改进过程,本实施例提出的EMD-MPF去噪方法具体流程在图2中进行描述,执行过程分为以下几个步骤:
步骤1:将MEMS陀螺仪的有噪信号x(t)根据EMD方法分解得到若干个IMF和一个残差信号。
步骤2:通过M1和M2两个标识参数对IMF进行阶次选择,将IMF分为噪声IMF、混合IMF和信息IMF三个部分。
步骤3:舍弃噪声IMF;保留信息IMF;对混合IMF使用改进的PF 方法进行降噪处理。
步骤4:对降噪处理后的混合IMF和信息IMF进行信号重构,得到最后完整的去噪结果,信号重构公式如下:
在基于EMD-MPF改进的陀螺信号去噪方法中,因为无法在理论上或实践中完全得知并去除噪声信号中的噪声。所以,既不能通过最小化目标信号和噪声信号的均方误差,也不能得到每个分量各自的信噪比来确定参数。因此,根据随机噪声的统计特性,只有在零点处随机噪声的自相关函数最大,在其他时刻其自相关函数则会快速衰减到零值附近。使用欧几里得距离是一种常用的分类依据,而引入自相关函数的方差,则是不仅考虑到了自相关函数下的无量纲化情况,同时依靠自相关函数考虑了不同分量之间时域上的相关性,一定程度上减少了欧几里得距离在单一特征下受不同分量的干扰;根据自相关特性,如果是噪声主体的IMF分量,则随机噪声在零点处有最大值,在其他时域会近似的快速下降为零;如果是混合信息的IMF分量,则随机噪声在其他时域处还会出现较大波动。在上述相关内容的基础上进行修改后,可以不用知道任何实际信号分量,快速区分噪声主要存在的IMF分量和信息IMF以此确定两者之间的指标参数。
为了确定第二个标识参数,引入了原始信号的概率密度函数(power spectraldensity function,PDF)与每个IMF的PDF之间的巴塔恰里雅距离 (Bhattacharyyadistance,BD)。巴氏距离常用于两离散样本之间的概率分布,以此判断测量类之间的可分离性,不仅排除了变量之间的相关性的干扰,而且是一种更直接的区分方法。如果IMF分量和原始信号两者PDF 之间的巴塔恰里雅距离越小,则可分离性越小,即IMF分量中包含的噪声越多;相反,两者PDF之间的巴塔恰里雅距离越大,则相似性越小、越可以分离,即IMF分量中包含的噪声越少。所以将巴塔恰里雅距离最大值的下一个点视为分界点,可以跟有效地区分混合IMF分量和信息IMF分量。
结合K-L散度进行改进,提出了一种创新的重采样过程。设定重采样过程的权值阈值为w0,对应高于权重阈值的粒子数目为T,粒子总数为N,则传统重采样过程为保留权值高于w0的粒子,然后根据权重重采样生成 N-T个粒子。这样就破坏了粒子样本集的有效性和多样性。
通过实验对本实施验证如下:
实验环境是惯性导航转动平台,该平台的惯性测量单元是一款全姿态测量传感装置IMU-300A,该装置由三轴的MEMS陀螺仪,三轴的MEMS 加速度计和三轴的磁阻型磁强计三部分构成。该装置固定于一个双轴电动转台之上,***可分别测出三个轴X、Y、Z上的绝对角速率和加速度值。
实验中数据的有效位通过实验平台标定标准进行确定,实验***采用的是双轴转台,***包括:一个航姿参考***和一个角度传感器,能够精确实现控制转台角度转动。***能够自动标定,通过航姿参考***输出姿态角提供相对准确的参考姿态角的信息。
为了较为清楚、全面的验证本实施例所提出方法的有效性,后续设置了两部分实验内容,第一部分采用MEMS陀螺仪的静态数据;第二部分采用动态数据。为了体现改进方法的信号去噪效果,分别使用了以下4种方法进行结果的分析、对比,方法一:使用传统粒子滤波;方法二:使用传统EMD方法;方法三:使用改进的PF方法;方法四:使用本实施例提出的基于MED-MPF改进的陀螺仪去噪方法。
3.1、转台静态数据测试
试验准备阶段,上电后,首先将装有陀螺的实验***进行回零操作,其次预热等待5min左右,数据的采样频率为100Hz,待***稳定后开始采集陀螺输出的原始数据;如图3所示,采集得到陀螺的静止数据(以X 轴的信号数据为例),数据长度为49428,持续时间约为8min左右。
首先,使用EMD方法对该噪声信号进行分解,得到如图4所示的15 阶IMF信号和1个信号残差。
根据分解得到的15阶IMF信号和一阶残差信号,分别计算各阶IMF 信号的自相关函数,结果如图5所示。简单分析一下自相关函数,自相关函数是描述随机信号在任意不同时刻的取值之间的相关程度,是对信号自身的互相关;由于随机噪声在零点处的自相关系数最大,同时陀螺信号的前后时刻有一定的相关性,导致随机噪声的关联不会立刻中断,前后时刻的自相关性也一定程度上反映了随机噪声在时域上的相关性,所以,使用连续IMF信号的自相关函数之间的距离均方差,可以在一定程度上判断随机噪声和有用信号的相关性,以此达到区分主要含随机噪声的IMF阶次的目的。
根据图5所得结果以及M1参数确定方法,结果如图6所示。
概率密度函数给出了信号取不同幅值大小的概率,很好的描述了信号幅值大小的分布情况,根据每阶IMF信号的PDF和原始信号PDF之间的巴氏距离,可确定出第二分界点的位置,以此区分混合IMF和信息IMF。每阶IMF和原始信号PDF之间的概率密度函数关系如图7所示,M2标识点的确定结果如图8所示。
根据图7和图9的分界点选取结果所示,可以得知在第5阶IMF信号之前,基本均是噪声IMF,可以直接舍弃;在第10阶IMF信号之后基本为大量无噪声信号;其中,在第5阶至第10阶之间,即为选出的混合IMF 部分,对该部分信号使用本实施例改进的粒子滤波方法,得到如图10所示的对比结果,最后将滤波后的结果和信息主导的IMF信号进行重构,得到最终的处理结果。
将滤波结果和信息主导IMF信号进行重构,并按之前提到的4种方法进行对比分析,结果如图11所示,可以看到本实施例提出的方法效果较为明显,相比原始数据和其他3种方法,降噪幅度较大,为了更直观的对比实验结果,后续进行了Allen方差分析,并分别计算了各自方法的信噪比和均方误差。
Allen方差是一种常用的分析陀螺噪声的方法,该方法可以较容易的对各种误差源的统计特性进行时域上细致的表征和辨识,因此通过Allen 方差的定量分析,可以进一步验证本实施方法对于陀螺去噪的有效性和适用性。通过图11结果对比,可以发现使用本实施例的方法,Allen方差分析所得到的各项误差均有不同程度的减小。
RMSE可以很好的反映样本数据的波动情况,数值越小说明波动越小、所含噪声越少。本实施例所提方法的RMSE值较原始数据降低了约91%,较方法1降低了约89%,较方法2降低了约81%,较方法3降低了约63%,效果十分明显。信噪比则更直接的反映了信号中的噪声含量,一般信噪比越大,说明混在信号中的噪声越小。本实施例所提方法的信噪比较原始信号提高了20.96dB,较方法1提高了19.01dB,较方法2提高了14.62dB,较方法3提高了8.67dB,表明该方法针对MEMS陀螺信号的去噪效果较为明显。分析结果如表1所示。
表1实验结果误差分析。
本实施例针对MEMS陀螺仪输出中随机噪声的问题,提出了一种新的结合EMD和改进PF的去噪方法。对EMD分解得到的多个IMF分量进行分析,通过两个分解点的计算得到噪声IMF、混合IMF和信息IMF,再对混合IMF使用改进的PF方法进行去噪,最后将去噪的混合IMF和信息 IMF进行重构,得到最终的去噪信号。在实验部分,分别使用了静态数据和动态数据,针对原始数据和4种去噪方法进行了结果对比,发现本实施例提出的方法有一定的参考价值,可以有效的去除MEMS陀螺仪输出信号中的噪声,较大的改善了信号质量。
本实施例具有以下优势:
1、相对于传统的EMD方法而言,提出了更为细致的分界点选择方式,将阶层分为三个部分,对于随机噪声更有针对性;
2、优势之二在于结合K-L散度针对PF算法进行了改进,提出了一种创新的重采样过程,使得粒子样本集保持了有效性和多样性;
3、优势之三经验模态分解将原始信号分解为若干个本征模态函数,然后通过分析和阶次选择把IMF分为噪声IMF、混合IMF和信息IMF,最后通过消除噪声和改进粒子滤波的方法得到最终的去噪信号。
最后应说明的是:以上所述仅为本发明的优选实施例而已,并不用于限制本发明,尽管参照前述实施例对本发明进行了详细的说明,对于本领域的技术人员来说,其依然可以对前述各实施例所记载的技术方案进行修改,或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内,所作的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。
Claims (3)
1.一种基于EMD-MPF改进的陀螺仪信号去噪方法,其特征在于,包括:
将陀螺仪的有噪信号分解为本征模态函数和残差信号;
通过确定的两个标识参数对本征模态函数进行阶次选择,将本征模态函数分为噪声本征模态函数、混合本征模态函数和信息本征模态函数;
舍弃所述噪声本征模态函数,保留所述信息本征模态函数,并对所述混合本征模态函数进行降噪处理;
对降噪处理后的混合本征模态函数和保留的所述信息本征模态函数进行信号重构,从而得到去噪的陀螺仪信号;
所述对降噪处理后的混合本征模态函数和保留的所述信息本征模态函数进行信号重构中的重构公式为:
所述MPF为改进的PF方法,包括:
把预测粒子样本的概率分布记作p(xi),1≤i≤N,把权重wi>w0的粒子样本xT的概率分布记作p1(xi),1≤i≤T,即为p,把剩余粒子样本xN-T的概率分布记作q(xi),N-T<i≤N,即为q,wi为权重,w0为权值阈值,N和T均为常数;
计算不同粒子样本之间的K-L散度如下:
对于小权重的粒子样本,产生N-T个随机数uj,其范围是uj∈(w0,1],对该粒子样本的权重进行更新,DKL为K-L散度:
所述确定的两个标识参数中,所述两个标识参数分别为M1和M2,选定确定标识参数M1,包括:使用自相关函数方差确定标识参数M1;
所述使用自相关函数方差确定标识参数M1,包括:
根据不同阶次的本征模态函数信号,计算各自对应的自相关函数;
根据连续本征模态函数的自相关函数计算两个自相关的函数之间的欧几里得距离均方差,得到计算结果MSED;
根据MSED确定标识参数M1;
所述标识参数M1与数值k和数值L相关联,
所述1≤k≤[2L/3],[2L/3]表示不超过2L/3的最大整数,从而保证M1的取值过大至L,k和L均为自然数;
确定标识参数M2,包括:
计算原始信号的概率密度函数与每个本征模态函数的概率密度函数之间的巴塔恰里雅距离,其中0≤BC≤1和0≤BD≤∞,BD为巴塔恰里雅距离,BC为巴塔恰里雅系数;
基于功率谱密度函数计算所述原始信号的概率密度函数和每个所述本征模态函数的概率密度函数之间的巴塔恰里雅距离数据在不同特征下的可分离性;
基于所述可分离性确定标识参数M2。
2.根据权利要求1所述的基于EMD-MPF改进的陀螺仪信号去噪方法,其特征在于,所述将陀螺仪的有噪信号分解为本征模态函数和残差信号为:
将陀螺仪的有噪信号分解为多个本征模态函数和一个残差信号。
3.根据权利要求1所述的基于EMD-MPF改进的陀螺仪信号去噪方法,其特征在于,所述通过确定的两个标识参数对本征模态函数进行阶次选择,将本征模态函数分为噪声本征模态函数、混合本征模态函数和信息本征模态函数,包括:
1<i<M1时,为噪声本征模态函数;
M1<i<M2时,为混合本征模态函数;
M2<i<L时,为信息本征模态函数,i表示本征态函数的阶次
L表示本征态函数的总阶数。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910826514.5A CN110542406B (zh) | 2019-09-03 | 2019-09-03 | 基于emd-mpf改进的陀螺仪信号去噪方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910826514.5A CN110542406B (zh) | 2019-09-03 | 2019-09-03 | 基于emd-mpf改进的陀螺仪信号去噪方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN110542406A CN110542406A (zh) | 2019-12-06 |
CN110542406B true CN110542406B (zh) | 2022-05-27 |
Family
ID=68711083
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201910826514.5A Active CN110542406B (zh) | 2019-09-03 | 2019-09-03 | 基于emd-mpf改进的陀螺仪信号去噪方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN110542406B (zh) |
Families Citing this family (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111157019A (zh) * | 2020-01-06 | 2020-05-15 | 金陵科技学院 | 一种基于EMD-Allan微机械陀螺信号分析的方法 |
CN111523076B (zh) * | 2020-03-24 | 2021-04-02 | 中国人民解放军军事科学院评估论证研究中心 | 基于Fal函数计算角加速度的方法、装置及*** |
CN112633225B (zh) * | 2020-12-31 | 2023-07-18 | 矿冶科技集团有限公司 | 矿用微震信号滤波方法 |
CN113537012A (zh) * | 2021-07-06 | 2021-10-22 | 国网江苏省电力有限公司常州供电分公司 | 接地网干扰信号的去噪方法、装置及计算机设备 |
CN114088077B (zh) * | 2021-12-10 | 2023-03-21 | 哈尔滨工业大学 | 一种改进的半球谐振陀螺信号去噪方法 |
CN116299562B (zh) * | 2023-05-26 | 2023-08-04 | 中国海洋大学 | 一种高度计测距电离层误差校正滤波处理方法 |
CN117433519A (zh) * | 2023-12-21 | 2024-01-23 | 武汉优米捷光电子制造有限责任公司 | 一种mems惯性测量组件的高精度温度补偿方法 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107766793A (zh) * | 2017-09-20 | 2018-03-06 | 天津大学 | 基于混合型方法的mems陀螺仪信号去噪处理方法 |
CN108338784A (zh) * | 2017-01-25 | 2018-07-31 | 中国科学院半导体研究所 | 基于eemd的小波熵阈值的心电信号去噪方法 |
CN110070031A (zh) * | 2019-04-18 | 2019-07-30 | 哈尔滨工程大学 | 一种基于emd和随机森林的海底底质声呐回波特征提取融合方法 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CA2676125A1 (en) * | 2007-02-09 | 2008-08-21 | Ge Healthcare Bio-Sciences Corp. | System and method for tracking the movement of biological materials |
-
2019
- 2019-09-03 CN CN201910826514.5A patent/CN110542406B/zh active Active
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108338784A (zh) * | 2017-01-25 | 2018-07-31 | 中国科学院半导体研究所 | 基于eemd的小波熵阈值的心电信号去噪方法 |
CN107766793A (zh) * | 2017-09-20 | 2018-03-06 | 天津大学 | 基于混合型方法的mems陀螺仪信号去噪处理方法 |
CN110070031A (zh) * | 2019-04-18 | 2019-07-30 | 哈尔滨工程大学 | 一种基于emd和随机森林的海底底质声呐回波特征提取融合方法 |
Non-Patent Citations (3)
Title |
---|
EMD阈值滤波在光纤陀螺漂移信号去噪中的应用;崔冰波,陈熙源,宋锐;《光学学报》;20150228;第35卷(第2期);第0207001-1-0207001-6页 * |
基于改进自适应粒子滤波的无线传感网船舶追踪;梅骁峻,吴华锋,陈彦臻,蒋恩青;《上海海事大学学报》;20180630;第39卷(第2期);第12-16页 * |
粒子滤波算法综述;胡士强, 敬忠良;《控制与决策》;20050430;第20卷(第4期);第361-366页 * |
Also Published As
Publication number | Publication date |
---|---|
CN110542406A (zh) | 2019-12-06 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN110542406B (zh) | 基于emd-mpf改进的陀螺仪信号去噪方法 | |
Narasimhappa et al. | A modified Sage-Husa adaptive Kalman filter for denoising Fiber Optic Gyroscope signal | |
Hušková et al. | Bootstrapping sequential change-point tests for linear regression | |
Kopparapu et al. | Identifying optimal Gaussian filter for Gaussian noise removal | |
Karaim et al. | Low-cost IMU data denoising using Savitzky-Golay filters | |
CN111260131B (zh) | 一种短期交通流的预测方法及装置 | |
CN112163536B (zh) | 基于粒子群算法改进的小波阈值函数去噪方法 | |
CN109798920A (zh) | 基于改进emd的mems陀螺随机误差建模滤波方法 | |
CN108469609B (zh) | 一种用于雷达目标跟踪的检测信息滤波方法 | |
Ramalingam et al. | Microelectromechnical systems inertial measurement unit error modelling and error analysis for low-cost strapdown inertial navigation system | |
CN112766224B (zh) | 畸变信号中提取真实信号的方法、装置、设备和存储介质 | |
Wang et al. | Radar HRRP target recognition in frequency domain based on autoregressive model | |
CN116541668B (zh) | 一种游泳划水次数确定方法、装置、设备及存储介质 | |
CN104111109B (zh) | 一种基于不同阶次统计量及支持向量机的机械振动状态识别方法 | |
Mukherjee et al. | New method for enhanced efficiency in detection of gravitational waves from supernovae using coherent network of detectors | |
CN111289800B (zh) | 一种基于广义回归神经网络的小电阻振动监测方法 | |
CN109443393B (zh) | 一种基于盲分离算法的捷联惯导信号提取方法及*** | |
Hayama et al. | Coherent network analysis for triggered gravitational wave burst searches | |
CN109344678B (zh) | 基于小波阈值的mems陀螺去噪方法 | |
Zhao et al. | An optimized heading de-noising method of the polarization compass for improving orientation accuracy | |
Ma et al. | Outlier correction method of telemetry data based on wavelet transformation and Wright criterion | |
CN110858309B (zh) | 一种多基准时钟加权合成方法 | |
Mahata et al. | Direct identification of continuous-time errors-in-variables models | |
Abdel-Hamid et al. | Improving the performance of MEMS-based inertial sensors by removing short-term errors utilizing wavelet multi-resolution analysis | |
Brogna et al. | Modelling vehicles NVH performances: a probabilistic approach. |
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 |