CN111766628A - 一种预条件的时间域弹性介质多参数全波形反演方法 - Google Patents

Abstract

本发明公开了一种预条件的时间域弹性介质多参数全波形反演方法，本发明包括如下步骤：S1、选择多个尺度频组对观测数据和地震子波进行预处理；S2、通过交错网格有限差分正演对目标震源子波和多个参数初始模型计算获得模拟观测记录，并通过将各尺度下的观测记录与对应尺度下的模拟观测记录做差获取波场残差；S3、通过伴随状态法求得正演算子的共轭算子，并通过波场残差得到残差波场；S4、通过残差波场和模拟观测记录计算参数初始模型的参数梯度，并对所有参数梯度进行预处理；S5、通过抛物线拟合法对初始迭代步长进行优化；S6、通过LBFGS法模型更新公式以及所有参数模型，得到多参数反演结果。本发明有效提高多参数全波形反演的效率、效果及稳定性。

Description

一种预条件的时间域弹性介质多参数全波形反演方法

技术领域

本发明涉及油气物探工程领域，具体地说是一种预条件的时间域弹性介质多参数全波形反演方法。

背景技术

现阶段油气勘探的重点正从构造勘探向岩性勘探转变，相比于传统偏移成像方法，多参数弹性波全波形反演在提供更为准确的纵、横波速度的同时，还可以提供更丰富的弹性参数等信息，有利于对地震资料解释和油气开发。

由于全波形反演理论基于数据拟合，严重依赖初始模型，低精度模型会导致方法不收敛。同时，由于地震资料采集技术等方面的限制，实际数据中的低频信息缺失会导致基于数值拟合的梯度导引类全波形反演方法难以收敛至全局极小值，另一方面，由于全波形反演方法计算量巨大，导致方法需要很高的计算成本。因此，如何提高多参数全波形反演方法的稳定性及反演效率是本领域技术人员需要解决的技术问题。

发明内容

本发明的目的是针对以上不足，提供一种预条件的时间域介质多参数全波形反演方法。

本发明所采用技术方案是：

一种预条件的时间域弹性介质多参数全波形反演方法，应用于对各向同性均匀弹性介质假设下的地下构造进行高精度的多参数建模，所述方法包括如下步骤：

S1、选择多个尺度频组对观测数据和地震子波进行预处理，获得不同尺度频组下的目标震源子波和观测记录；

S2、通过交错网格有限差分正演对目标震源子波和多个参数初始模型计算获得模拟观测记录，并通过将各尺度下的观测记录与对应尺度下的模拟观测记录做差获取波场残差；

S3、通过伴随状态法求得正演算子的共轭算子，并通过波场残差得到残差波场；

S4、通过残差波场和模拟观测记录计算参数初始模型的参数梯度，并对所有参数梯度进行预处理；

S5、通过抛物线拟合法对初始迭代步长进行优化；

S6、通过LBFGS法模型更新公式以及所有参数模型，得到多参数反演结果。

作为对本发明方法的进一步优化，本发明步骤S1中，对输入数据的预处理过程为对地震子波和观测记录进行滤波得到低频目标震源子波和滤波后的观测记录，且该过程中对地震子波和观测记录进行滤波所采用的滤波器的计算式为：

式中，f_target为低频目标震源子波，f_original为原始地震子波，ε为控制数值溢出的因子。

作为对本发明方法的进一步优化，本发明步骤S2中，多个所述参数初始模型包括纵波速度初始模型、横波速度初始模型和密度初始模型；

作为对本发明方法的进一步优化，本发明步骤S4中，通过步骤S3中将波场残差反向传播得到的残差波场与步骤S2获得模拟观测记录计算参数梯度所用的梯度，在此过程中，所采用的梯度公式为：

其中：λ和μ为拉梅常数，ρ为密度。

作为对本发明方法的进一步优化，本发明步骤S4中，将基于拉梅常数及密度进行参数化时的参数梯度转换得到纵波速度初始模型、横波速度初始模型和密度初始模型的参数，其转换公式为：

作为对本发明方法的进一步优化，本发明步骤S4中，通过对计算的对角Hessian算子对各参数梯度进行预处理，得到预处理后的参数梯度，所述对角Hessian算子的计算方法为：

作为对本发明方法的进一步优化，本发明步骤S4中，利用求取的对角Hessian算子对各参数梯度进行预处理，得到预处理后的参数梯度，其预处理后的梯度的计算方法为：

式中，

为目标函数关于λ,μ,ρ的梯度，P_λ、P_μ和P_ρ为的

的预处理算子。

P_i的构造方法如下：

式中，

为对角近似的伪Hessian矩阵，κ为用来保证算法稳定性而加入的阻尼因子。

作为对本发明方法的进一步优化，本发明步骤S5，通过抛物线拟合法对初始迭代步长进行优化的计算公式为：

式中，α₀,α₁,α₂为给定初始步长，E(α_i)为由三个步长得到的对应的目标函数值，并且满足如下关系：

作为对本发明方法的进一步优化，本发明步骤S6中，通过LBFGS法模型更新公式以及所有参数模型的计算过程为：

式中，α_k为每次迭代时的求得的优化步长，H^-1为误差函数的Hessian矩阵的逆，

为误差函数关于模型参数的梯度；

其中，所用更新公式中的H^-1的求取方法，其计算方法为：

s_k＝δm_k+1-δm_k

式中，m表示计算过程中使用了最近m次迭代的模型更新信息。

作为对本发明方法的进一步优化，本发明通过LBFGS法模型更新公式以及所有参数模型的过程中，若不满足迭代终止条件，则将得到的模型反演结果作为初始模型继续进行更新直到满足迭代终止条件；若满足迭代终止条件，则接着判断所用多尺度频率是否等于最大频率，若不满足，则提高所用多尺度频率，直到其与最大频率相等并满足迭代终止条件时反演结束。

本发明具有以下优点：

本发明通过应用多尺度反演策略和梯度预处理算子作为预算条件，提高了多参数全波形反演的稳定性；可使模型深度能量得到很好的补偿，提高深处的反演质量，并且能够加速目标函数的收敛效率。

附图说明

为了更清楚地说明本发明实施例中的技术方案，下面将对实施例中描述中所需要使用的附图作简要介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域的普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

下面结合附图对本发明进一步说明：

图1为本发明的流程示意图。

图2为洼陷模型的纵波模型；

图3为等梯度模型反演的初始模型；

图4为纵波速度在8hz频率下迭代25次的反演结果；

图5为纵波速度在8hz频率下迭代50次的反演结果；

图6为纵波速度在15hz频率下迭代25次的反演结果；

图7为纵波速度在15hz频率下迭代50次的反演结果；

图8为横波速度在8hz频率下迭代25次的反演结果；

图9为横波速度在8hz频率下迭代50次的反演结果；

图10为横波速度在15hz频率下迭代25次的反演结果；

图11为横波速度在15hz频率下迭代50次的反演结果；

图12为密度在8hz频率下迭代25次的反演结果；

图13为密度在8hz频率下迭代50次的反演结果；

图14为密度在15hz频率下迭代25次的反演结果；

图15为密度在15hz频率下迭代50次的反演结果；

图16为纵波速度的梯度；

图17为纵波速度的预处理后的梯度；

图18为横波速度的梯度；

图19为横波速度的预处理后的梯度；

图20为密度的梯度；

图21为密度的预处理后的梯度；

图22为深度为1000米的纵波速度反演曲线图；

图23为深度为1000米的纵波速度反演曲线图；

图24为深度为1000米的密度反演曲线图；

图25为归一化的误差能量随迭代次数的变化曲线图。

具体实施方式

下面结合附图和具体实施例对本发明作进一步说明，以使本领域的技术人员可以更好地理解本发明并能予以实施，但所举实施例不作为对本发明的限定，在不冲突的情况下，本发明实施例以及实施例中的技术特征可以相互结合。

需要理解的是，在本发明实施例的描述中，“第一”、“第二”等词汇，仅用于区分描述的目的，而不能理解为指示或暗示相对重要性，也不能理解为指示或暗示顺序。在本发明实施例中的“多个”，是指两个或两个以上。

本发明实施例中的属于“和/或”，仅仅是一种描述关联对象的关联关系，表示可以存在三种关系，例如，A和/或B，可以表示：单独存在A，单独存在B，同时存在A和B这三种情况。另外，本文中字符“/”一般表示前后关联对象是一种“或”关系。

本发明提供一种预条件的时间域介质多参数全波形反演方法的实施例，应用于对各向同性均匀弹性介质假设下的地下构造进行高精度的多参数建模，其特征在于：如图1所示，所述方法包括如下步骤：

S1、选择多个尺度频组对观测数据和地震子波进行预处理，获得不同尺度频组下的目标震源子波和观测记录；所述观测数据即为通过检波器实际采集到的数据，这里的观测数据包括但不限于起伏地表高程、观测系统、野外观测记录(即炮点数据)和可以预设的计算参数等，其中偏移参数包括偏移速度场横向采样点及纵向采样点、空间采样间隔、时间采样间隔、时间采样点数和主频。

通过选择多个尺度频组对地震子波和观测记录进行维纳滤波，得到低频目标震源子波和滤波后的观测记录，其主要目的是通过低通滤波把观测数据和地震子波分解成不同频带的观测数据和子波成分。然后由低频带到高频带逐频带进行反演，并把低频带反演结果作为高频带反演时输入的初始模型，从低频信息到高频信息分别重建模型的低、高波数成分，逐步收敛得到全局极值点，增加了反演过程的稳定性。在此滤波过程中，采用Wiener滤波器，其计算共算公式如下：

式中，f_target为低频目标震源子波，f_original为原始地震子波，ε为控制数值溢出的因子

S2、通过交错网格有限差分正演对目标震源子波和多个参数初始模型计算获得模拟观测记录，并通过将各尺度下的观测记录与对应尺度下的模拟观测记录做差获取波场残差；

S3、通过伴随状态法求得正演算子的共轭算子，并通过波场残差反向传播得到残差波场；

S4、通过残差波场和模拟观测记录计算参数初始模型的参数梯度，并对所有参数梯度进行预处理；

S401、参数梯度的求取过程如下：

在二维情况下，一阶速度-应力波动方程：

其中，v_x和v_z为质点在不同方向的速度分量，τ_ij为应力，λ和μ为拉梅常数，ρ为密度。

可将上式写为

其中，u＝(v_x,v_z,τ_xx,τ_zz,τ_xz)^T为波场向量；f＝(f_x,f_z,f_xx,f_xz,f_zz)^T为震源项；矩阵A的具体形式为

定义最小二乘意义下的误差函数为：

其中，L是正演算子，m是模型参数，d_obs是观测得到的地震数据。基于伴随状态法进行梯度的推导：

伴随算子满足关系式：

将正演波动方程代入式(6)并整理得到伴随波动方程：

从式(6)可以看出，基于一阶速度-应力波动方程推导得到的伴随方程，与公式(1)的原始波动方程有很大的差别，为了简化后续计算，定义线性变换算子T：

T＝T'Λ(T')^-1 (7)

其中

(T')^-1为T'的逆矩阵；

将T作用到伴随波场向量得到

其中，

将波场向量w代入式(7)后整理得到新的伴随波动方程

由式(10)可以看到，通过线性变换，伴随波动方程式(10)与正演波动方程式(2)的形式是一样的。将其进一步整理得到：

为了表示方便，将误差函数写成关于波场向量u(m)与模型参数m的函数：

E(m)＝H(u,m) (11)

其中，波场向量u(m)与模型参数m的状态方程形式为：

M(u,m)＝0 (12)

在Born近似假设下，对模型参数m引入模型参数扰动δm时，波场向量u存在相应的扰动量δu。将加入扰动后的状态方程Taylor一阶展开后可以得到

其中，

为数据空间与模型空间中的任意变量，由式(14)与式(13)可得：

向式(12)中引入模型参数扰动δm可得

其中，<·>_D表示在数据空间上的内积。

若可以显式写出

将式(15)代入式(16)可得

借助算子

的伴随算子，将式(17)整理为

定义伴随波场λ，使其满足伴随波动方程：

则目标函数随模型参数的扰动可以表示为

将误差函数、公式(2)的正演波动方程与公式(10)的伴随波动方程代入式(20)得到

其中，算子A'的具体形式为

经过线性变换后的得到的新的伴随波场w，w由式(23)决定

将式(22)和式(23)代入式(21)，整理变换后可以得到误差函数关于参数λ,μ,ρ的梯度：

因为拉梅常数和纵横波速度及密度之间具有如下关系：

应用链式求导法则，可以得到误差函数关于v_p,v_s,ρ_v的梯度

S402、对参数梯度预处理的过程如下：

时间域全波形反演中，由于深部照明不足，导致深部构造反演效果不是很理想，为了增强深部能量对梯度进行预处理。Hessian矩阵的对角元素表征了对介质参数的地下照明，利用Hessian矩阵对梯度进行预处理可以补偿地下不均匀照明的现象。由于全Hessian矩阵的计算成本很高，可将其非线性部分舍去，近似为其线性部分。

基于伪Hessian矩阵的对角元素构造纵波速度和横波速度梯度预处理算子。在只考虑炮点的Green函数不考虑接收点的Green函数情况下,得到伪Hessian算子如下：

式中i为炮数，u_i为第i炮对应的波场向量。对式(26)取其对角元素，得到纵波、横波的伪Hessian矩阵表达形式

构造密度梯度的预处理算子时，采用的近似Hessian格式如下

其中，f(ω)为地震子波，G(x'；x_s,ω)是炮点的Green函数，G(x'；x_r,ω)为接收点的Green函数。G(x'；x_s,ω)可以在正演时得到，但是在时间域求解G(x'；x_r,ω)需要很大的计算量，为了进一步简化H_a，假设炮点与检波点分布相同，即令G(x；x_r,ω)＝G(x；x_s,ω)，并令x'＝x”，得到

在利用

对梯度进行预处理时,为保证算法稳定性加入阻尼因子κ，κ是一个很小的正数，得到的预处理算子为

按照公式(30)的P_i的构造方法，可分别得到P_λ、P_μ和P_ρ。

根据式(31)，利用P_ρ、P_λ和P_μ得到预处理后的误差函数关于λ,μ,ρ的梯度

然后代入式(25)，这样就得到了预处理后的误差函数关于v_p,v_s,ρ_v的梯度

S5、通过抛物线拟合法对初始迭代步长进行优化；

相比于固定步长算法来说，应用三点抛物线拟合法估算得到的优化步长，可以减少反演迭代次数，从而提高计算效率。

给定初始步长：α₀,α₁,α₂，由三个步长得到的目标函数值满足如下关系时(式(32))，可以拟合得到一条抛物线，抛物线极小值点对应的步长值即为求得的最优步长。

求取最优步长α_opt的具体公式如下

S6，最后利用模型更新公式更新各参数模型；所用更新公式如下：

其中，α_k为每次迭代时的求得的优化步长，H^-1为误差函数的Hessian矩阵的逆，

为误差函数关于模型参数的梯度。

其中，H^-1在每次迭代过程中利用最近m次的模型更新信息进行隐式计算，其更新公式如下

其中

S6、通过LBFGS法模型更新公式以及所有参数模型，得到多参数反演结果。

判断是否满足迭代终止条件：不满足迭代终止条件则将得到的模型反演结果作为初始模型继续进行更新直到满足迭代终止条件；满足迭代终止条件的情况下判断所用多尺度频率是否等于最大频率，不满足则提高所用多尺度频率，直到其与最大频率相等并满足迭代终止条件时反演结束，此时得到最终的多参数反演结果。

为了便于对本发明方法做进一步的说明，以洼陷模型为例说明本发明的方法，流程如图1所示，所用洼陷模型的纵波模型如图2所示，可以看到，在模型的650m深处有五个等距分布的绕射体，绕射体纵波速度从左到右依次为4000m/s、3700m/s、3200m/s、2000m/s、2500m/s。在泊松体假设下，通过纵波速度模型生成横波速度模型，应用Gardner公式换算得到密度模型，这样就得到了所用的弹性介质洼陷模型。以如图3所示的等梯度模型作为反演的初始模型。

模型参数为：模型网格大小为100×200，纵向网格间距为10m，横向网格间距为10m，水平坐标范围0-2000m。模型的观测系统为：震源子波采用15Hz的零相位Ricker子波，激发方式为地表激发，一共50炮，第一炮位于第40米处，接收方式采用地表全孔径接收。时间采样率0.7ms，共2000个采样点，记录时间长度为1.4s，模型四周均为PML吸收边界。多尺度情况下，使用8Hz和15Hz频率先后各迭代25次，共50次。

将纵波速度模型、横波速度模型和密度模型应用多尺度反演策略之后，其结果如图4-图15所示，图4-图7为纵波速度的演算结果，图8-图11为横波速度的演算结果，图12-图15为密度的演算结果，通过对比同频率下迭代25与迭代50次的结果可以看出，纵波速度模型、横波速度模型和密度模型的反演结果中的洼陷界面附近和绕射体周围的假象(图中箭头处)得到了明显的压制。而且应用了梯度预处理算子之后，通过对比相同迭代次数不同频率的反演结果下可以看出模型深部反演效果得到了明显的改善，650米深度的5个绕射体可以较清晰的反映在各参数的反演结果上，最深层的反射界面的位置和形态也更清晰和准确。

图16-图21显示了纵波速度模型、横波速度模型和密度模型的梯度及其预处理后的梯度，通过同一参数模型预处理前后对比可以直观看出，应用了梯度预处理算子之后的各参数梯度的中深层能量得到了明显的加强，洼陷构造和五个绕射体以及深层反射界面都能反映在梯度上。总体来看，各参数预处理后的梯度能量也更加均衡，这样可以让每次迭代时的梯度更加有用，从而使迭代更加高效，也自然可以改善反演结果并且提高计算效率。

选用简单洼陷模型进行方法测试，分别测试了常规EFWI、应用多尺度反演策略的EFWI、应用梯度预处理的EFWI以及同时应用梯度预处理和多尺度反演策略的EFWI算法，并对比分析了多尺度反演策略和预条件算子的有效性及其优点。图22-图24分别为纵波速度模型、横波速度模型和密度模型在深度为1000米处的反演曲线图，图中：点线为初始模型、点划线为真实模型、虚线为常规EFWI、实线是同时应用梯度预处理算子和多尺度反演策略的EFWI。在各层介质分界处，常规EFWI得到的纵波速度、横波速度的振幅值都与真实值有较大的偏差，其曲线形态与真实曲线形态的偏差也很大，在700米深度以下其曲线形态与初始模型十分相近，说明其对模型深部几乎没有进行有效的反演；而文中方法即同时应用梯度预处理算子和多尺度反演策略的EFWI其反演得到的纵波速度和横波速度的曲线形态和振幅值都最接近真实值，说明其对整个模型都进行了较好的反演。而由于密度差异主要控制反射波场的强弱，对波场旅行时不敏感，因此密度是三参数中最难反演的参数，密度的反演过程不如纵波速度和横波速度稳定，在反演曲线上可以看出，其振幅值及曲线形态与真实值的偏差也比较大，但是通过对比我们也能看到，通过文中方法得到的密度反演结果，相比于常规EFWI方法而言，其曲线形态和振幅值都更加接近真实值，说明该方法也能有效改善密度参数的反演结果。

图25的归一化误差能量随迭代次数变化曲线中，虚线为常规EFWI，点划线为应用多尺度反演策略的EFWI，点线为应用了梯度预处理算子后的EFWI，实线是同时应用梯度预处理算子和多尺度反演策略的EFWI。需要注意的是，图中点划线和实线的误差曲线在第25次迭代后误差能量有一个明显的突变，这是多尺度反演策略的应用导致的。通过对比虚线和点划线可以看出，梯度预处理算子的应用可以有效加快反演效率，特别是前几次迭代时这种改善尤为明显，对比点划线和实线也得到相同的结论。

通过上述洼陷模型试算及反演结果进行对比，可以得出：

(1)与常规EFWI相比，应用了多尺度反演策略之后，反演效果得到了较好的改善，模型构造形态得到了较好的恢复，同时假象更少，构造边界也更清晰；

(2)利用预处理算子对梯度进行预处理后，模型深层梯度得到了很好的增强。从反演结果上也可以看出，与常规EFWI相比，模型深层照明不足的问题得到了很好的改善。通过对比误差能量曲线，可以看出梯度预处理可以明显提高误差函数的收敛速度；

(3)文中方法即同时应用梯度预处理算子和多尺度反演策略的EFWI得到的多参数反演结果，无论是从参数曲线形态上还是参数值上都最接近真实值，是四种方法中最好的反演结果，证明了本发明所用方法的准确性和有效性。

(4)相比于常规EFWI相比，利用本发明的方法，可以更快得到更好的反演结果。

以上所述实施例仅是为充分说明本发明而所举的较佳的实施例，本发明的保护范围不限于此。本技术领域的技术人员在本发明基础上所作的等同替代或变换，均在本发明的保护范围之内。本发明的保护范围以权利要求书为准。

Claims (10)

1.一种预条件的时间域弹性介质多参数全波形反演方法，应用于对各向同性均匀弹性介质假设下的地下构造进行高精度的多参数建模，其特征在于：

所述方法包括如下步骤：

S1、选择多个尺度频组对观测数据和地震子波进行预处理，获得不同尺度频组下的目标震源子波和观测记录；

S2、通过交错网格有限差分正演对目标震源子波和多个参数初始模型计算获得模拟观测记录，并通过将各尺度下的观测记录与对应尺度下的模拟观测记录做差获取波场残差；

S3、通过伴随状态法求得正演算子的共轭算子，并通过波场残差得到残差波场；

S4、通过残差波场和模拟观测记录计算参数初始模型的参数梯度，并对所有参数梯度进行预处理；

S5、通过抛物线拟合法对初始迭代步长进行优化；

S6、通过LBFGS法模型更新公式以及所有参数模型，得到多参数反演结果。

2.根据权利要求1所述的方法，其特征在于：

步骤S1中，对输入数据的预处理过程为对地震子波和观测记录进行滤波得到低频目标震源子波和滤波后的观测记录，且该过程中对地震子波和观测记录进行滤波所采用的滤波器的计算式为：

式中，f_target为低频目标震源子波，f_original为原始地震子波，ε为控制数值溢出的因子。

3.根据权利要求2所述的方法，其特征在于：

步骤S2中，多个所述参数初始模型包括纵波速度初始模型、横波速度初始模型和密度初始模型。

4.根据权利要求3所述的方法，其特征在于：

步骤S4中，通过步骤S3中将波场残差反向传播得到的残差波场与步骤S2获得模拟观测记录计算参数梯度所用的梯度，在此过程中，所采用的梯度公式为：

其中：λ和μ为拉梅常数，ρ为密度。

5.根据权利要求4所述的方法，其特征在于：

步骤S4中，将基于拉梅常数及密度进行参数化时的参数梯度转换得到纵波速度初始模型、横波速度初始模型和密度初始模型的参数，其转换公式为：

6.根据权利要求1所述的方法，其特征在于：

步骤S4中，通过对计算的对角Hessian算子对各参数梯度进行预处理，得到预处理后的参数梯度，所述对角Hessian算子的计算方法为：

7.根据权利要求6所述的方法，其特征在于：

步骤S4中，利用求取的对角Hessian算子对各参数梯度进行预处理，得到预处理后的参数梯度，其预处理后的梯度的计算方法为：

式中，

为目标函数关于λ,μ,ρ的梯度，P_λ、P_μ和P_ρ为的

的预处理算子。

P_i的构造方法如下：

式中，

为对角近似的伪Hessian矩阵，κ为用来保证算法稳定性而加入的阻尼因子。

8.根据权利要求1所述的方法，其特征在于：

步骤S5，通过抛物线拟合法对初始迭代步长进行优化的计算公式为：

式中，α₀,α₁,α₂为给定初始步长，E(α_i)为由三个步长得到的对应的目标函数值，并且满足如下关系：

9.根据权利要求1所述的方法，其特征在于：

步骤S6中，通过LBFGS法模型更新公式以及所有参数模型的计算过程为：

式中，α_k为每次迭代时的求得的优化步长，H^-1为误差函数的Hessian矩阵的逆，▽_mE(m)为误差函数关于模型参数的梯度；

其中，所用更新公式中的H^-1的求取方法，其计算方法为：

s_k＝δm_k+1-δm_k

y_k＝(▽_mE(m))_k+1-(▽_mE(m))_k

式中，m表示计算过程中使用了最近m次迭代的模型更新信息。

10.根据权利要求9所述的方法，其特征在于：

通过LBFGS法模型更新公式以及所有参数模型的过程中，若不满足迭代终止条件，则将得到的模型反演结果作为初始模型继续进行更新直到满足迭代终止条件；若满足迭代终止条件，则接着判断所用多尺度频率是否等于最大频率，若不满足，则提高所用多尺度频率，直到其与最大频率相等并满足迭代终止条件时反演结束。

Images