CN112083494B - 一种基于改进萤火虫步长因子的多级粒子滤波算法 - Google Patents

一种基于改进萤火虫步长因子的多级粒子滤波算法 Download PDF

Info

Publication number
CN112083494B
CN112083494B CN202011012994.0A CN202011012994A CN112083494B CN 112083494 B CN112083494 B CN 112083494B CN 202011012994 A CN202011012994 A CN 202011012994A CN 112083494 B CN112083494 B CN 112083494B
Authority
CN
China
Prior art keywords
particle
particles
weight
time
firefly
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
Application number
CN202011012994.0A
Other languages
English (en)
Other versions
CN112083494A (zh
Inventor
于涵
王典
张敬权
滕明
刘财
赵军
李鹏
刘元骐
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Jilin University
Original Assignee
Jilin University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Jilin University filed Critical Jilin University
Priority to CN202011012994.0A priority Critical patent/CN112083494B/zh
Publication of CN112083494A publication Critical patent/CN112083494A/zh
Application granted granted Critical
Publication of CN112083494B publication Critical patent/CN112083494B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/364—Seismic filtering
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/004—Artificial life, i.e. computing arrangements simulating life
    • G06N3/006—Artificial life, i.e. computing arrangements simulating life based on simulated virtual individual or collective life forms, e.g. social simulations or particle swarm optimisation [PSO]
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/30—Noise handling
    • G01V2210/32—Noise reduction
    • G01V2210/324—Filtering

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Remote Sensing (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Biophysics (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Evolutionary Computation (AREA)
  • General Health & Medical Sciences (AREA)
  • Molecular Biology (AREA)
  • Computing Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Biomedical Technology (AREA)
  • Computational Linguistics (AREA)
  • Software Systems (AREA)
  • Artificial Intelligence (AREA)
  • Health & Medical Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本申请公开了一种基于改进萤火虫步长因子的多级粒子滤波算法，该方法对萤火虫粒子滤波算法中的亮度公式以及位置更新公式进行修改，采用自适应步长萤火虫算法进行迭代更新，并对结果进行多级滤波。本申请通过自适应步长萤火虫算法，驱使权值小的粒子向高似然区域移动，减弱了粒子退化问题，可以在使用较少粒子的情况下达到较高的滤波精度，节约了计算资源。

Description

一种基于改进萤火虫步长因子的多级粒子滤波算法
技术领域
本发明涉及地震勘测领域，具体为一种基于改进萤火虫步长因子的多级粒子滤波算法。
背景技术
地震勘探是以地下岩层的结构物性差异为基础，在地表利用人工方法产生地震波，通过观测分析地震波向地下传播的响应，从而推断地下地层结构，物质，组成及演化等的探测方法，是勘测地下石油，天然气等资源的重要手段，在油田和工程地质勘探，区域地质及地壳研究等方面起着至关重要的作用。地震勘探过程包括野外地震数据采集，室内资料处理以及地震记录解释分析等工作。其中第一阶段野外地震数据采集，工作人员首先选定指定区域人工激发震源，利用阵列排布的检波器接收地下反射到地表的地震波，由于地表条件，震源激发，地震波传播及检波器接收等诸多方面的原因，随机噪声通常叠加在地震资料中，当其能量过大时，有用信号就被淹没，严重影响最终成像效果，不利于反射信号的追踪和识别，对后续的资料分析和读取造成很大困扰，因此，如何更好地从含噪记录中提取有用信号，消除随机干扰，最大限度地还原有效信号的能量和信息，达到地震勘探的三高要求，即“高信噪比，高分辨率，高保真度”，一直是地震勘探工作者们深入研究的课题。因此，压制地震勘探数据中的随机噪声是地震勘探过程中至关重要的环节之一。
应用于地震勘探随机噪声压制的方法有很多，包括：叠加处理，多项式拟合，经验模态分解，卡尔曼滤波，小波变换等，并且与之相应的改进算法也相继被提出，用于完善地震资料降噪过程中的不足。但是这些降噪方法，在较高信噪比的情况下可以得到比较理想的效果。如果地震资料的信噪比较低，数据受到强随机噪声的干扰，一些常规的方法就很难达到预期的去噪效果，甚至导致有效信号的损失及衰减。因此，在低信噪比下，如何对强随机噪声进行有效地压制，还原信号的信息与波形，更好得提取出有效的地震信号仍是亟待解决的课题。
粒子滤波技术在非线性，非高斯系统中具有优越性，并且应用领域十分广泛，如视觉跟踪；目标定位，导航，跟踪领域；通信与信号处理领域；图像处理；卫星遥感等。本发明主要应用于通信与信号处理领域，对地震勘探随机噪声压制领域具有重要的现实意义，可以为地震信号处理领域提供新的思路，为满足地震探勘的精度要求提供新的方法。
粒子滤波是蒙特卡洛思想结合递归贝叶斯估计的非线性滤波方法，可以对非线性非高斯系统的参数或者状态进行估计，估计精度可以逼近最优。粒子是指采样样本，定性地表示状态的空间信息，滤波是指通过粒子统计方式获取当前系统的准确状态，粒子滤波通常理解为利用当前有效的观测值和预测的系统状态值估计出当前系统的有效状态值，这样可以从带噪声的数据中提取出系统有效信号的状态，可以用于有效数据的估计，在处理非高斯，非线性时变系统的参数估计和状态滤波问题方面具有独特的优势和广阔的前景。
近年来，随着计算机水平的不断提高，粒子滤波的研究取得了很大的发展。许多改进的粒子滤波方法被用于解决不同的实际问题。例如，Pitt和Shephard提出的辅助粒子滤波器（Auxiliary Particle Filter, APF）以解决粒子多样性丧失的正则化粒子滤波器(Regularized Particle Filter, RPF)。自适应粒子滤波（Adaptive PF, APF）主要用于需要在线实时处理的应用，使粒子随环境的变化而改变数量，避免使冗余的粒子使算法的复杂度和运算量增加。Kotecha等人提出高斯粒子滤波(Gaussian PF, GPF)解决了重采样引起样本枯竭的问题，并改善算法的收敛性和并行运算性。
现有技术中，针对萤火虫算法以及粒子滤波的相关研究有：1.田梦楚,薄煜明,陈志敏,吴盘龙,赵高鹏.“萤火虫算法智能优化粒子滤波”[J].《自动化学报》,2016,42(01):89-97. DOI 10.16383/j.aas.2016.c150221，该文献中，针对粒子滤波重采样导致的粒子贫化以及需要大量粒子才能进行状态估计的问题，该方法结合粒子滤波的运行机制，引入萤火虫群体的优胜劣汰机制及萤火虫个体的吸引和移动的行为，使粒子群智能地向高似然区域移动，提高粒子群的整体质量。2.方正,佟国峰,徐心和.“粒子群优化粒子滤波方法”[J].《控制与决策》,2007(03):273-277，其针对粒子滤波方法存在粒子贫乏及初始状态未知时需要大量粒子才能进行鲁棒状态预估问题，将粒子群优化思想引入粒子滤波中。该方法将最新观测值融合到采样过程中，并对采样过程利用粒子群优化算法进行优化，将粒子集合朝后验概率密度分布较大的区域运动，从而克服粒子贫乏问题，降低精确预估所需的粒子数。3.乔美玉. “基于EMD的粒子滤波在地震勘探随机噪声压制中的应用”[D].吉林大学,2014，由于粒子滤波算法所要求的状态空间模型受到观测噪声和状态噪声共同作用，当信噪比较低时，强随机噪声的干扰会影响粒子序列的生成，给粒子分布带来偏差。该方法结合时，频域均有较好分辨率的经验模态分解方法，修正噪声对粒子序列带来的偏差，使粒子分布更接近信号值真实后验分布，在低信噪比条件下突出有效信号。
现有技术中主要存在如下缺陷：
1. 在标准的萤火虫算法改进粒子滤波过程中，利用全局最优值指导粒子群体的运动，使所有的粒子都朝着全局最优粒子方向移动，虽然克服了陷入局部极值的现象，但是随着迭代次数的增加，粒子都朝着全局最优值附近移动，会造成最优值附近粒子密度过大，降低了粒子的多样性，会影响下一时刻粒子的寻优能力。
2.由于地震资料较为复杂，而粒子滤波的滤波结果很大程度依赖于所建立的模型，因此得到适合建立模型的地震信号非常重要。由于采用原始带噪声的数据所建立的模型是较为粗糙的，受到噪声的影响较大。
发明内容
由于萤火虫算法自身运行机制的特殊性，直接将萤火虫算法与粒子滤波进行直接融合会导致许多问题，例如粒子交互会导致算法的复杂性大幅提高，出现局部最优现象，因此需要对萤火虫算法进行修正。在标准的萤火虫算法中，每只萤火虫需要相对于其他所有萤火虫进行一次位置更新，即对于种群数量为N的萤火虫算法，一次迭代中的时间复杂度为
Figure DEST_PATH_IMAGE002
，若算法需要T次迭代，则总的时间复杂度为
Figure DEST_PATH_IMAGE004
，将其直接应用到粒子滤波中会影响粒子滤波的实时性。
本发明对标准的萤火虫算法进行改进，提出一种新的步长因子，让其随着迭代的增加，自适应调整步长的大小。在迭代初期，使用较大的步长来提高粒子的探索能力，以保证粒子获得较为宽广的分布，避免陷入局部最优解；在迭代后期，粒子已经处于最优值周围因此采用小步长。这一改进在保证算法寻优能力的同时，提高了算法的效率。同时，本发明在对数据进行去噪的过程中，提出多级滤波的概念，将滤波的结果进行多次建模并滤波。通过多级滤波提高建模的精度，得到与纯净数据更为接近的模型再进行滤波处理，从而提高滤波的精度。
本发明包括以下内容：
1.修改萤火虫亮度公式
在萤火虫算法中，萤火虫种群将会向着亮度高的地方移动，因此将使萤火虫的荧光亮度等于粒子的归一化权值，即
Figure DEST_PATH_IMAGE006
，其中
Figure DEST_PATH_IMAGE008
表示 t时刻的粒子
Figure DEST_PATH_IMAGE010
，
Figure DEST_PATH_IMAGE012
表示t时刻粒子
Figure 19634DEST_PATH_IMAGE010
的归一化权值。将发光亮度定义为归一化权值，将使粒子向权值大的区域移动，因此，能够解决权值退化问题。选取合适的迭代次数或者设置阈值，可以使粒子在大权值附近有较好的覆盖，同时增大粒子的多样性。
由于粒子滤波问题中往往没有局部最优的需求，因此，在每次迭代寻优过程中，只需选出当前的全局最优值，让萤火虫向当前状态的全局最优值移动即可。另外，通过定义一个阈值threshold及最大迭代次数，当粒子权值
Figure DEST_PATH_IMAGE014
大于threshold或达到最大迭代次数时，粒子停止迭代更新，因此只对分布在概率尾部的粒子进行调整。
2.修改位置更新公式
在标准萤火虫算法中，步长因子α是一个固定值，这会导致在迭代次数达到一定程度时，粒子会在最优值附近反复震荡。在传统的自适应步长萤火虫算法中，步长因子与迭代次数成反比，往往是一个以迭代次数为自变量的单调递减函数，因此，本申请设计了一种新的步长因子α的自适应调整策略，即
Figure DEST_PATH_IMAGE016
，
Figure DEST_PATH_IMAGE018
表示t时刻粒子
Figure DEST_PATH_IMAGE020
的权值，即荧光亮度，
Figure DEST_PATH_IMAGE022
表示t时刻粒子集的平均权值，
Figure DEST_PATH_IMAGE024
表示t时刻粒子集的最大权值。这使得α与粒子的归一化权重成反比，权值大的说明距离最优解近，因此采用小步长，权值小的距离最优解远，采用大步长。修正的萤火虫算法的位置更新公式为：
Figure DEST_PATH_IMAGE026
，
其中，g代表当前迭代次数，d代表第d维度，
Figure DEST_PATH_IMAGE028
代表t时刻第g次迭代中的全局最优值，
Figure DEST_PATH_IMAGE030
为服从均匀分布的一个随机数。改进后的算法FAPF，虽然在一定程度上增加了算法的复杂程度，但是在算法性能上得到了提高。
一种基于改进萤火虫步长因子的多级粒子滤波算法，其特征在于包括以下步骤：
S1：使用时变自回归模型（AR）对地震数据进行分段建模，得到系统的状态方程，误差及数据长度T；
S2：从初始概率密度
Figure DEST_PATH_IMAGE032
采样N个粒子，每个粒子的权值为1/N，得到粒子集合
Figure DEST_PATH_IMAGE034
,。
S3：根据建模所得的状态方程：
Figure DEST_PATH_IMAGE036
进行采样更新，得到新的粒子集合：
Figure DEST_PATH_IMAGE038
。其中，
Figure 822593DEST_PATH_IMAGE008
表示t时刻的粒子
Figure 8855DEST_PATH_IMAGE020
， y表示观测值，q为参考分布即建模所得的状态方程。
S4:根据似然函数
Figure DEST_PATH_IMAGE040
计算所有粒子的权值
Figure DEST_PATH_IMAGE042
并归一化处理，
Figure DEST_PATH_IMAGE044
，使萤火虫的荧光亮度等于粒子的归一化权值，即
Figure 773373DEST_PATH_IMAGE006
，其中
Figure 728691DEST_PATH_IMAGE008
表示 t时刻的粒子
Figure 816732DEST_PATH_IMAGE020
，
Figure 755869DEST_PATH_IMAGE012
表示t时刻粒子
Figure 233731DEST_PATH_IMAGE020
的归一化权值；
S5: 设置阈值threshold及最大迭代次数。采用自适应步长萤火虫算法进行迭代更新，根据修改后的位置更新公式
Figure DEST_PATH_IMAGE046
将小权值的粒子向高似然区域移动并计算更新后粒子的归一化权值
Figure 574713DEST_PATH_IMAGE012
。其中，
Figure DEST_PATH_IMAGE048
代表t时刻第g次迭代中的全局最优值，β(r)为吸引度函数，定义为：
Figure DEST_PATH_IMAGE050
,
Figure DEST_PATH_IMAGE052
为最大吸引度，γ为光强吸收系数，
Figure DEST_PATH_IMAGE054
为粒子
Figure DEST_PATH_IMAGE056
与全局最优值
Figure 879662DEST_PATH_IMAGE048
之间的空间距离，
Figure DEST_PATH_IMAGE058
为服从均匀分布的一个随机数；其中，
Figure DEST_PATH_IMAGE060
，
Figure 444111DEST_PATH_IMAGE012
表示t时刻粒子
Figure DEST_PATH_IMAGE062
的归一化权值，即荧光亮度，，
Figure DEST_PATH_IMAGE064
表示t时刻粒子集的平均权值，
Figure DEST_PATH_IMAGE066
表示t时刻粒子集的最大权值。
S6:随着粒子的迭代，不断计算并更新全局最优值，当粒子的最小权值
Figure DEST_PATH_IMAGE068
时停止迭代，否则继续迭代至最大次数并实施重采样技术，选择并复制权重较大的粒子。
S7：状态估计：
Figure DEST_PATH_IMAGE070
；
S8：令t=t+1，进行下一时刻的滤波。当t=T时，完成滤波;
S9：将S8的结果重复步骤S1-S8，重复次数≥1。
优选地，步骤S7中，根据噪声的大小选择合适的重复次数。
优选地，步骤S1中，将建模产生的误差设置为采样方程中的过程噪声。
有益效果
1. 在传统的自适应步长萤火虫算法中，步长因子与迭代次数成反比，一般是以迭代次数为自变量的单调递减函数，而本发明提出的新的自适应步长因子
Figure 770005DEST_PATH_IMAGE016
，与粒子的权值成反比，随着迭代的进行，根据粒子的权值大小自适应调整步长因子的大小，因此能够通过较小的迭代次数来达到较好的收敛效果。在基本的粒子滤波算法中，由于粒子权值的方差会随着时间迭代而不断增加，普遍存在退化问题。本方案通过自适应步长萤火虫算法，驱使权值小的粒子向高似然区域移动，减弱了粒子退化问题，可以在使用较少粒子的情况下达到较高的滤波精度，节约了计算资源。
2. 在对地震数据进行噪声压制时，采用多级滤波的策略，将滤波后的数据进行多次建模并滤波。随着滤波级数的增加，可以明显提高滤波的精度，提高数据的信噪比。本方案对地震数据进行降噪处理时，对地震数据不需要进行较多的预处理，通过多级滤波策略提高地震数据的信噪比。
3. 由于地震数据资料较为复杂，对其观测噪声的定量较为困难，本方法在使用时变自回归（AR）算法对数据进行建模时，由于建模存在一定的误差，将此误差设定为采样方程中的过程噪声，降低了人为计算观测方差的难度，并且操作简便，减少了人为设定参数对结果的影响。
附图说明
图1为本发明技术方案的流程图
图2为Ricker子波示意图
图3为加噪信号示意图
图4为PF，FAPF滤波对比图
图5为PF，FAPF多级滤波对比图
图6为PF，FAPF信噪比对比图
图7为频谱对比
图8为地质模型
图9为理论模型
图10为带入噪声的含噪记录
图11为FAPF滤波结果
图12为单道记录时域对比图
图13为单道记录频域对比图
具体实施方式
本发明技术方案的流程图的如图1所示，一种基于改进萤火虫步长因子的多级粒子滤波算法，其特征在于包括以下步骤：
S1：使用时变自回归模型（AR）对地震数据进行分段建模，得到系统的状态方程，误差及数据长度T，将建模产生的误差设置为采样方程中的过程噪声；
S2：从初始概率密度
Figure 699915DEST_PATH_IMAGE032
采样N个粒子，每个粒子的权值为1/N，得到粒子集合
Figure 660918DEST_PATH_IMAGE034
,。
S3：根据建模所得的状态方程：
Figure 309068DEST_PATH_IMAGE036
进行采样更新，得到新的粒子集合：
Figure 863153DEST_PATH_IMAGE038
。其中，
Figure 913148DEST_PATH_IMAGE008
表示t时刻的粒子
Figure 779473DEST_PATH_IMAGE020
， y表示观测值，q为参考分布即建模所得的状态方程。
S4:根据似然函数
Figure 383761DEST_PATH_IMAGE040
计算所有粒子的权值
Figure 869100DEST_PATH_IMAGE042
并归一化处理，
Figure 646038DEST_PATH_IMAGE044
，使萤火虫的荧光亮度等于粒子的归一化权值，即
Figure 824210DEST_PATH_IMAGE006
，其中
Figure 306007DEST_PATH_IMAGE008
表示 t时刻的粒子
Figure 595037DEST_PATH_IMAGE020
，，
Figure 963833DEST_PATH_IMAGE012
表示t时刻粒子
Figure 309976DEST_PATH_IMAGE020
的归一化权值；
S5: 设置阈值threshold及最大迭代次数。采用自适应步长萤火虫算法进行迭代更新，根据修改后的位置更新公式
Figure 888856DEST_PATH_IMAGE046
，将小权值的粒子向高似然区域移动并计算更新后粒子的归一化权值
Figure 106210DEST_PATH_IMAGE012
。其中，
Figure 719726DEST_PATH_IMAGE048
代表t时刻第g次迭代中的全局最优值，β(r)为吸引度函数，定义为：
Figure 115066DEST_PATH_IMAGE050
,
Figure 837034DEST_PATH_IMAGE052
为最大吸引度，γ为光强吸收系数，
Figure 211077DEST_PATH_IMAGE054
为粒子
Figure 679098DEST_PATH_IMAGE056
与全局最优值
Figure 245340DEST_PATH_IMAGE048
之间的空间距离，
Figure 64391DEST_PATH_IMAGE058
为服从均匀分布的一个随机数；其中，
Figure 623549DEST_PATH_IMAGE060
，
Figure 208726DEST_PATH_IMAGE012
表示t时刻粒子
Figure 680290DEST_PATH_IMAGE062
的归一化权值，即荧光亮度，
Figure 721058DEST_PATH_IMAGE064
表示t时刻粒子集的平均权值，
Figure 428114DEST_PATH_IMAGE066
表示t时刻粒子集的最大权值。
S6:随着粒子的迭代，不断计算并更新全局最优值，当粒子的最小权值
Figure DEST_PATH_IMAGE072
＞threshold时停止迭代，否则继续迭代至最大次数并实施重采样技术，选择并复制权重较大的粒子。
S7：状态估计：
Figure 664536DEST_PATH_IMAGE070
；
S8：令t=t+1，进行下一时刻的滤波。当t=T时，完成滤波;
S9：将S8的结果重复步骤S1-S8，重复次数≥1。
采用本申请的方法进行正演实验，本次正演模拟采用主频为30Hz的Ricker子波进行实验，Ricker子波的波函数为：
Figure DEST_PATH_IMAGE074
，其中，
Figure DEST_PATH_IMAGE076
为中心频率。给Ricker子波加入3db的高斯白噪噪声进行对比试验，结果如图2-4所示，从图2~4可以看出，对带噪声的Ricker子波信号采用PF和FAPF算法进行噪声压制处理都可以取得一定的效果，在一定程度上提高信噪比。由于去噪结果图并不能非常直观得体现出某一种算法的效果更好，因此利用绝对误差平均值，绝对误差方差及均方根误差计算其统计特性来体现去噪能力，结果如下表所示：
表1 Ricker子波去噪结果对比
算法
Figure DEST_PATH_IMAGE078
（绝对误差平均值）
Figure DEST_PATH_IMAGE080
（绝对误差方差）		 REMS（均方根误差）
PF 0.0531 0.0016 0.0668
FAPF 0.0507 0.0015 0.0638
表2 Ricker子波去噪信噪比对比
噪声大小(db) 噪声类型 PF滤波(db) FAPF滤波（db）
3 高斯白噪声 7.3567 8.1384
从表1-2中可以看出，FAPF较PF，绝对误差均值
Figure DEST_PATH_IMAGE082
，方差
Figure DEST_PATH_IMAGE084
，均方根误差REMS都要稍小,且FAPF滤波后的信噪比更高。由于一次滤波得到的结果提高了信噪比，而粒子滤波的结果很大程度上依赖于所建立的模型精度，与建模过程中的过程噪声及观测噪声的大小密切相关，利用原始含噪数据只进行一次滤波，得到的结果是较为粗糙的，因此提出多级滤波的概念，将一次滤波后的结果进行多次建模并滤波，实验结果如图4所示。通过观察图5可以看出，多级FAPF滤波后的结果可以更好地还原雷克子波的波形，整体波形比较干净，很好地跟踪了雷克子波的时变特征，基本恢复了纯净雷克子波的形态，高频噪声被有效抑制，低频信号得到了较好得还原，并且信噪比得到了明显提高。根据图6可以看出，FAPF滤波结果的信噪比明显增加，但增加的幅度不断减少，因此，为了得到最好的滤波效果并且兼顾节约计算资源，需要根据噪声的大小选择合适的滤波级数。
地震记录的合成是地震子波与反射系数相互褶积的结果，将合成的地震记录中加入-1db的高斯白噪声，去噪结果如图8-13所示。图8为地质模型，从图9~图11可以发现，FAPF滤波后，结果整体上较为平滑，可以有效地压制噪声，突出有效信号。为了更好地体现出FAPF滤波的效果，选取其中第十道记录，进行时域和频域分析，结果如图12-13所示，从图12~13中可以看出，在时域上，滤波结果已经较好地拟合了纯净信号，在地层界面处的子波波形还原较好；在频域上，高频信号在一定程度在被压制，突出了低频有效信号，与原始纯净记录相比，虽然仍还有部分高频信号，但是在时域波形中不影响子波的识别及同相轴的追踪。
以上为本发明较佳的实施方式，本发明所属领域的技术人员还能够对上述实施方式进行变更和修改，因此，本发明并不局限于上述的具体实施方式，凡是本领域技术人员在本发明的基础上所作的任何显而易见的改进、替换或变型均属于本发明的保护范围。

Claims (3)

1.一种基于改进萤火虫步长因子的多级粒子滤波算法，其特征在于包括以下步骤：
S1：使用时变自回归模型(AR)对地震数据进行分段建模，得到系统的状态方程，误差及数据长度T；
S2：从初始概率密度p(x0)采样N个粒子，每个粒子的权值为1/N，得到粒子集合
Figure FDA0003025589330000011
S3：根据建模所得的状态方程：
Figure FDA0003025589330000012
进行采样更新，得到新的粒子集合：
Figure FDA0003025589330000013
其中，
Figure FDA0003025589330000014
表示t时刻的粒子xi，y表示观测值，q为参考分布即建模所得的状态方程；
S4:根据似然函数p(yt|xt)计算所有粒子的权值
Figure FDA0003025589330000015
并归一化处理，
Figure FDA0003025589330000016
使萤火虫的荧光亮度等于粒子的归一化权值，即
Figure FDA0003025589330000017
其中
Figure FDA0003025589330000018
表示t时刻的粒子xi，
Figure FDA0003025589330000019
表示t时刻粒子xi的归一化权值；
S5:设置阈值threshold及最大迭代次数；采用自适应步长萤火虫算法进行迭代更新，根据修改后的位置更新公式
Figure FDA00030255893300000110
Figure FDA00030255893300000111
将小权值的粒子向高似然区域移动并计算更新后粒子的归一化权值
Figure FDA00030255893300000112
其中，
Figure FDA00030255893300000113
代表t时刻第g次迭代中的全局最优值，β(r)为吸引度函数，定义为：
Figure FDA00030255893300000114
β0为最大吸引度，γ为光强吸收系数，ri为粒子
Figure FDA00030255893300000115
与全局最优值
Figure FDA00030255893300000116
之间的空间距离，
Figure FDA00030255893300000117
为服从均匀分布的一个随机数；其中，
Figure FDA00030255893300000118
Figure FDA00030255893300000119
表示t时刻粒子xi的归一化权值，即荧光亮度，
Figure FDA00030255893300000120
表示t时刻粒子集的平均权值，
Figure FDA00030255893300000121
表示t时刻粒子集的最大权值；
S6:随着粒子的迭代，不断计算并更新全局最优值，当粒子的最小权值
Figure FDA00030255893300000122
时停止迭代，否则继续迭代至最大次数并实施重采样技术，选择并复制权重较大的粒子；
S7：状态估计：
Figure FDA0003025589330000021
S8：令t＝t+1，进行下一时刻的滤波；当t＝T时，完成滤波；
S9：将S8的结果重复步骤S1-S8，重复次数≥1。
2.根据权利要求1中所述的一种基于改进萤火虫步长因子的多级粒子滤波算法，其特征在于：步骤S7中，根据噪声的大小选择合适的重复次数。
3.根据权利要求1中所述的一种基于改进萤火虫步长因子的多级粒子滤波算法，其特征在于：步骤S1中，将建模产生的误差设置为采样方程中的过程噪声。
CN202011012994.0A 2020-09-23 2020-09-23 一种基于改进萤火虫步长因子的多级粒子滤波算法 Active CN112083494B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011012994.0A CN112083494B (zh) 2020-09-23 2020-09-23 一种基于改进萤火虫步长因子的多级粒子滤波算法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011012994.0A CN112083494B (zh) 2020-09-23 2020-09-23 一种基于改进萤火虫步长因子的多级粒子滤波算法

Publications (2)

Publication Number Publication Date
CN112083494A CN112083494A (zh) 2020-12-15
CN112083494B true CN112083494B (zh) 2021-07-06

Family

ID=73739694

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011012994.0A Active CN112083494B (zh) 2020-09-23 2020-09-23 一种基于改进萤火虫步长因子的多级粒子滤波算法

Country Status (1)

Country Link
CN (1) CN112083494B (zh)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102779336A (zh) * 2012-07-09 2012-11-14 重庆大学 基于混沌的光电图像海杂波抑制方法
CN102932870A (zh) * 2012-10-30 2013-02-13 河南科技大学 一种无线传感器网络节点部署方法
WO2018182765A1 (en) * 2017-03-28 2018-10-04 Saudi Arabian Oil Company A seismic processing workflow for broadband single-sensor single-source land seismic data
CN110782815A (zh) * 2019-11-13 2020-02-11 吉林大学 一种全息立体探测系统及其方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102323617B (zh) * 2011-06-13 2014-03-12 中国石油化工股份有限公司 一种复杂地表的二维地震资料连片处理方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102779336A (zh) * 2012-07-09 2012-11-14 重庆大学 基于混沌的光电图像海杂波抑制方法
CN102932870A (zh) * 2012-10-30 2013-02-13 河南科技大学 一种无线传感器网络节点部署方法
WO2018182765A1 (en) * 2017-03-28 2018-10-04 Saudi Arabian Oil Company A seismic processing workflow for broadband single-sensor single-source land seismic data
CN110782815A (zh) * 2019-11-13 2020-02-11 吉林大学 一种全息立体探测系统及其方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于萤火虫群优化的空间非合作目标相对导航粒子滤波算法;张大力，等;《中国惯性技术学报》;20170430;第25卷(第2期);269-274 *

Also Published As

Publication number Publication date
CN112083494A (zh) 2020-12-15

Similar Documents

Publication Publication Date Title
CN104614718B (zh) 基于粒子群算法的激光雷达波形数据分解的方法
CN110850482B (zh) 一种基于变分模态分解原理的瞬变电磁信噪分离方法
CN112946749B (zh) 基于数据增广训练深度神经网络压制地震多次波的方法
CN111190227B (zh) 基于残差卷积生成对抗模型的低信噪比地震数据消噪方法
CN102819043B (zh) 阵列信号随机噪声自适应模型去噪方法
CN109031422A (zh) 一种基于CEEMDAN与Savitzky-Golay滤波的地震信号噪声抑制方法
CN108550116A (zh) 低信噪比下的硅单晶生长图像的自适应随机共振去噪方法
CN112598593B (zh) 基于非均衡深度期望块对数似然网络的地震噪声压制方法
CN113887398A (zh) 一种基于变分模态分解和奇异谱分析的gpr信号去噪方法
Lv et al. Noise removal for semi-airborne data using wavelet threshold and singular value decomposition
CN112882115A (zh) 基于gwo优化小波阈值的大地电磁信号去噪方法及系统
CN112083494B (zh) 一种基于改进萤火虫步长因子的多级粒子滤波算法
CN111624558A (zh) 一种基于De-Chirp技术的SAR干扰抑制方法及装置
CN111856590B (zh) 海洋大地电磁探测的海浪磁干扰压制方法
Hao et al. Denoising method based on spectral subtraction in time-frequency domain
CN111767887B (zh) 基于小波分解与ime频率估计的瞬变电磁数据处理方法
CN110596756B (zh) 基于自适应混合复扩散模型的沙漠地震勘探噪声压制方法
CN112363217A (zh) 一种地震数据随机噪声压制方法及系统
Chang et al. Random noise suppression for seismic data using a non-local Bayes algorithm
CN112200069A (zh) 时频域谱减法和经验模态分解联合的隧道滤波方法及系统
CN110687605A (zh) 基于改进的k-svd算法在地震信号处理的算法分析应用
CN111898567B (zh) 基于经验小波变换的弹丸激波信号噪声抑制方法
Yang et al. Sidelobe suppression of SAR images by spectrum shaping
CN112711074B (zh) 一种地震初至波的去噪方法及装置
CN116091882B (zh) 一种基于自适应双通道pcnn的偏振图像融合方法

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