CN103163554A - 利用零偏vsp资料估计速度和q值的自适应波形的反演方法 - Google Patents

Abstract

本发明公开了一种利用零偏VSP资料估计速度和Q值的自适应波形的反演方法，该方法可根据实际VSP资料的信噪比及直达波受上行波干涉的程度，自适应最大限度地选取未受干扰的初至波片段，充分地利用有效信息来构造目标函数；针对Q值比其它参数更难估计的特点，通过恰当地构造数据加权矩阵提高了目标函数对Q值变化的敏感度；采用近来发展的自适应乘性正则化方法改善了非线性反演的病态问题，并通过非线性变换限定了待求参数的取值范围；推导出了雅可比矩阵各元素的解析表达式，从而减小了反问题的计算量。该方法与谱比值法和子波包络峰值瞬时频率法等相比较，该方法法受上行波影响相对较小、抗噪性能更强，反演的参数精度高。

Description

利用零偏VSP资料估计速度和Q值的自适应波形的反演方法

技术领域

本发明属于地球物理勘探领域，涉及一种地球模型物性参数反演方法，特别涉及一种基于高斯牛顿法利用零偏VSP地震数据反演介质速度和品质因子Q的方法。

背景技术

地震波在地层中传播时会逐渐衰减，衰减因素包括几何扩散、透射/反射损失、介质的非完全弹性等。其中由介质的非完全弹性引起的衰减是介质的固有特性，这种特性常用品质因子Q来描述。非完全弹性引起的衰减会影响地震波的振幅和相位：如果已知Q值，进行反Q滤波，则可以使地震信号的高频成分得到增强，从而提高地震资料的纵向分辨率；另外，Q值也能提供岩性、饱和度、孔隙率等信息，可用于储层识别和烃类检测。所以，反演介质的Q值具有重要的实际意义。

国内外学者提出了多种Q值直接估计的方法。时间域方法有子波模拟法^[1]、解析信号法^[2]和脉冲上升时间法等，这些方法都需要振幅保真度高的数据，而实际记录常受到几何扩散，反/透射损失等因素的影响，致使这些时间域直接法估计的Q值精度降低。频率域常用的方法有对数谱比值法^[3](LSR,log-spectrum ratio)、谱匹配法、中心频率偏移法^[4]等，这些方法首先要用一个时窗截取一段地震记录，再对截出的地震记录进行分析。这类方法中采用不同时窗对估计结果影响较大，如何最佳地选择窗的类型及长度是一个难题。利用脉冲瞬时属性的变化来估计Q值也是重要的途径之一。Mathneey和Nowack^[5]提出瞬时频率匹配法，高静怀和杨森林等人^[6]提出利用小波域地震子波包络峰值处瞬时频率估计Q值的方法WEPIF(wavelet envelope peakinstantaneous frequency)。这类方法的优点是时窗易于选取或不需要加时窗，估计的Q值分辨率较高，但由于只利用了脉冲包络峰值处单个或相邻几个数据值，所以估计的Q值受随机噪声的影响很大。Tonn^[7]对近十种直接估计Q值的方法进行了比较，结果表明每种方法都有一定的适用条件。

与反射地震资料相比,VSP资料可对穿过介质的波进行直接测量,这为利用波形反演估计介质Q值及其它参数提供了条件。在这方面国内外学者开展了大量的研究工作,如：Stewart^[8]提出同时利用相邻四道记录的全波形在频率域反演速度、衰减因子、上行波和下行波等；Amundsen和Mittet^[9]只利用初至下行波和一次反射波在频率域反演相速度和品质因子Q，该方法考虑了反/透射系数随频率的变化，但要求界面位置和介质密度是已知的，这在有些实际问题中难以满足。

目前常用的各种Q值估计(或反演)方法都有其使用条件，用于具体实际问题时仍有一些问题需解决，如：在加窗截取初至下行波或一次反射波时如何选取最佳窗，在频率域如何选择最佳频带等。另外，要获得精确的Q值，还需知道准确的速度等参数。

发明内容

本发明的目的在于克服上述现有技术的缺点，提供一种利用零偏VSP资料估计速度和Q值的自适应波形的反演方法，该方法利用零偏VSP资料初至波时域波形同时反演介质速度、品质因子Q值和透射系数的自适应方法，基于高斯牛顿发，通过恰当地构造目标函数和巧妙地利用约束条件，从而使得算法抗噪能力强、反演的参数精度高。

本发明的目的是通过以下技术方案来解决的：

该种利用零偏VSP资料估计速度和Q值的自适应波形的反演方法，包括以下步骤：

1）在井中布置检波器采集原始地震数据，然后对采集到的地震数据进行校正处理、噪音压制处理、真振幅恢复处理；处理后得到零偏VSP记录，记为Seis(n,t)，n=1,2,…,N.，其中N表示总的检波器个数，t表示时间轴；

2）计算VSP地震波形Seis(n,t)的振幅包络，将每道记录的初至波振幅包络峰值进行归一化处理，记归一化后的振幅包络为IA(n,t)；并分别拾取每道数据直达波的瞬时振幅峰值到达时刻t_n，用于在步骤4）中计算速度初值；

3）以检波器的位置为界面，将介质划分为N-1个薄层，第n个薄层的厚度为第n个检波器和第n+1个检波器之间的距离，记为Δz_n；则第n个薄层上方的地震记录为Seis(n,t)，下方的地震记录为Seis(n+1,t)，波由界面上方传播到界面下方的频率域正演算子为：

$\mathrm{Seis}(n+1,\mathrm{\ω})=T_n\·\mathrm{Seis}(n,\mathrm{\ω})\·\exp [\mathrm{i\ω}(\frac{{\mathrm{\Δz}}_n}{V_n}+\frac{{\mathrm{\Δz}}_n}{\mathrm{\π}{V}_n{Q}_n}\ln \frac{{\mathrm{\ω}}_r}{\mathrm{\ω}})-\frac{\mathrm{\Δ}{z}_n\mathrm{\ω}}{2V_n{Q}_n}$

其中：ω为角频率，ω_r为参考频率，T_n为透射系数，V_n为速度，Q_n为品质因子；

4）对于任意第n个薄层，选定待反演的模型参数为

并给定误差精度ε；比较该层上下两道记录初至波的归一化振幅包络IA(n,t)和IA(n+1,t)，以层下方初至波起始点作为时窗的起始点t_s，以两振幅包络误差大于给定值ε时所对应的时刻作为时窗的终点t_e，截取Seis(n+1,t)中位于[t_s t_e]之间的波形作为观测数据d_o，d_o为含N个采样值的向量；

5）给定一个初始模型参数m⁰，以该模型参数将第n层上方Seis(n,t)的初至波传至层下方，截取[t_s t_e]之间的波形作为计算数据d_c，d_c为含N个采样值的向量；

6）以计算数据d_c与观测数据d_o的误差的模平方为目标函数，并采用乘性自适应正则化方法以改善非线性反演的病态问题，采用数据加权矩阵以提高目标函数对Q值变化的敏感度；

采用乘性正则化方法后的目标函数为：

$E\left(m\right)=\frac{1}{2}{\left|\right|W_{d}e\left(m\right)\left|\right|}^2\·{\mathrm{\η}}_k({\left|\right|W_m(m-m^p)\left|\right|}^2+{\mathrm{\δ}}_k^2)$

其中：

$e\left(m\right)=\frac{d_c-d_o}{\left|\right|d_o\left|\right|},$

${\mathrm{\η}}_k=\frac{1}{{\left|\right|W_m(m^k-m^p)\left|\right|}^2+{\mathrm{\δ}}_k^2},$

式中k表示迭代次数，W_d是数据加权矩阵，W_m是模型参数加权矩阵，m^p是根据先验信息人为给定的期望模型参数向量，该正则化方法对δ_k值不敏感，所以δ_k取为一个常数，并由数值试验加以确定；

由Q值误差引起的数据误差与观测数据振幅成正比，而由速度误差引起的数据误差与观测数据的振幅成反比。所以，为了提高目标函数对Q的敏感度，本发明采用如下的数据加权矩阵：

$\left(W_d\right)=\mathrm{diag}\{\frac{\left|d_i\right|}{\left|\right|d_o\left|\right|},d_i\∈d_o,i=1...N\}.;$

7）利用高斯牛顿法进行迭代反演求取目标函数的极小值，该极小值对应的模型参数即为第n个薄层的最优解。

进一步的，以上步骤7）中，迭代过程中采用非线性变换策略将反演的模型参数约束在一定的范围内；高斯牛顿法的迭代公式为：

$m^{k+1}=m^k-H_a^{-1}G,$

式中，k表示迭代次数，G和H_a分别表示目标函数的梯度矩阵和近似Hessian矩阵，表达式分别为：

$G=J^T{W}_d^T{W}_{d}e\left(m^k\right)+{\mathrm{\η}}_{k}F\left(m^k\right)W_m^T{W}_m(m^k-m^p),$

$H_a=J^T{W}_d^T{W}_{d}J+{\mathrm{\η}}_{k}F\left(m^k\right)W_m^T{W}_m$

$+{\mathrm{\η}}_k[W_m^T{W}_m(m^k-m^p){]}^T{J}^T{W}_d^T{W}_{d}e\left(m^k\right)$

$+{\mathrm{\η}}_k[J^T{W}_d^T{W}_{d}e\left(m^k\right){]}^T{W}_m^T{W}_m(m^k-m^p),$

$F\left(m^k\right)=\frac{1}{2}{\left|\right|W_{d}e\left(m^k\right)\left|\right|}^2$

其中雅可比矩阵J是正演模型对m的偏导数，每次迭代必须重新计算。

与现有技术相比，本发明具有以下有益效果：

本发明利用零偏VSP资料估计速度和Q值的自适应波形的反演方法可以根据实际VSP资料的信噪比及直达波受上行波干涉的程度，自适应最大限度地选取未受干扰的初至波片段，充分地利用有效信息来构造目标函数；针对Q值比其它参数更难估计的特点，通过恰当地构造数据加权矩阵提高了目标函数对Q值变化的敏感度；采用近来发展的自适应乘性正则化方法改善了非线性反演的病态问题，并通过非线性变换限定了待求参数的取值范围；推导出了雅可比矩阵各元素的解析表达式，从而减小了反问题的计算量。该方法与谱比值法和子波包络峰值瞬时频率法等相比较，其受上行波影响相对较小、抗噪性能更强，反演的参数精度高。

附图说明

图1是本发明中零偏VSP地震资料采集示意图和按检波器将介质划分为薄层的示意图。

具体实施方式

下面结合附图对本发明做进一步详细描述：

本发明基于单程波传播理论，采用Q值与频率无关的假设，提出利用零偏VSP资料初至下行波的时域波形同时反演介质速度和Q值的方法。

本发明的物质基础是在井中采集的垂直地震剖面（VSP）数据，即在地表设置震源，在垂直井内安置检波器接收的地震波。本发明基于高斯牛顿法利用零偏VSP数据反演介质参数的具体步骤分别为：

1）通过在井口放炮、井中布置检波器采集原始地震数据（如图1所示），然后对采集到的地震资料做常规预处理，如校正处理、噪音压制处理、真振幅恢复处理等。处理后得到零偏VSP记录，记为Seis(n,t)，n=1,2,…,N.，其中N表示总的检波器个数，t表示时间轴。

2）计算VSP地震波形Seis(n,t)的振幅包络，将每道记录的初至波振幅包络峰值进行归一化处理，记归一化后的振幅包络为IA(n,t)，并分别拾取每道数据直达波的瞬时振幅峰值到达时刻t_n。

3）以检波器的位置为界面，将介质划分为N-1个薄层（如图1所示），第n个薄层的厚度为第n个检波器和第n+1个检波器之间的距离，记为Δz_n。则第n个薄层上方的地震记录为Seis(n,t)，下方的地震记录为Seis(n+1,t)，波由界面上方传播到界面下方的频率域正演算子为：

Seis(n+1,ω)=T_n·Seis(n,ω)·P(Δz_n,ω)

$P({\mathrm{\Δz}}_n,\mathrm{\ω})=\exp [\mathrm{i\ω}(\frac{{\mathrm{\Δz}}_n}{V_n}+\frac{{\mathrm{\Δz}}_n}{{\mathrm{\πV}}_n{Q}_n}\ln \frac{{\mathrm{\ω}}_r}{\mathrm{\ω}})-\frac{{\mathrm{\Δz}}_n\mathrm{\ω}}{{2V}_n{Q}_n}]$

其中：ω为角频率，ω_r为参考频率，T_n为透射系数，V_n为速度，Q_n为品质因子。

4）对于任意第n个薄层，选定待反演的模型参数为

并给定误差精度ε，比较该层上下两道记录初至波的归一化振幅包络IA(n,t)和IA(n+1,t)，选取振幅包络误差小于ε的时段内对应的初至波片段作为观测数据d_o。

比较时，需将层上方初至波的振幅包络峰值处平移至层下方初至波的振幅包络峰值处。然后以层下方初至波起始点作为时窗的起始点t_s，以两振幅包络误差大于给定值ε时所对应的时刻作为时窗的终点t_e，截取Seis(n+1,t)中位于[t_s t_e]之间的波形作为观测数据d_o，d_o为含N个采样值的向量。

5）给定一个初始模型参数m⁰，以该模型参数将第n层上方Seis(n，t)的初至波传至层下方，截取[t_s t_e]之间的波形作为计算数据d_c，d_c为含N个采样值的向量。

考虑到现在的检波器响应一致性较好，透射系数一般在1左右，所以T_n的初值可在1左右选取。Q的初值根据已有的研究和测量可在10－200之间选取，也可用其他方法（如谱比值法）估计初值。V_n的初值可由初至波到达时差来估计，即V_n=Δz_n/(t_n+1t_n)。

6）以计算数据d_c与观测数据d_o的误差的模平方为目标函数，并采用乘性自适应正则化方法以改善非线性反演的病态问题，采用特殊的数据加权矩阵以提高目标函数对Q值变化的敏感度。

乘性正则化方法是将正则化因式作为一个乘性因子引入到目标函数中，正则化后的目标函数可写为：

$E\left(m\right)=\frac{1}{2}{\left|\right|W_{d}e\left(m\right)\left|\right|}^2\·{\mathrm{\η}}_k({\left|\right|W_m(m-m^p)\left|\right|}^2+{\mathrm{\δ}}_k^2)$

其中：

$e\left(m\right)=\frac{d_c-d_o}{\left|\right|d_o\left|\right|},$

${\mathrm{\η}}_k=\frac{1}{{\left|\right|W_m(m^k-m^p)\left|\right|}^2+{\mathrm{\δ}}_k^2},$

式中k表示迭代次数，W_d是数据加权矩阵，W_m是模型参数加权矩阵，m^p是根据先验信息人为给定的期望模型参数向量，该正则化方法对δ_k值不敏感，所以δ_k可取为一个常数，并由数值试验加以确定。

由Q值误差引起的数据误差与观测数据振幅成正比，而由速度误差引起的数据误差与观测数据的振幅成反比。所以，为了提高目标函数对Q的敏感度，本发明采用如下的数据加权矩阵：

$\left(W_d\right)=\mathrm{diag}\{\frac{\left|d_i\right|}{\left|\right|d_o\left|\right|},d_i\∈d_o,i=1...N\}.$

7）利用高斯牛顿法进行迭代反演求取目标函数的极小值，该极小值对应的模型参数即为第n个薄层的最优解。迭代过程中采用非线性变换策略将反演的模型参数约束在一定的范围内。

高斯牛顿法的迭代公式为：

$m^{k+1}=m^k-H_a^{-1}G,$

这里k表示迭代次数，G和H_a分别表示目标函数的梯度矩阵和近似Hessian矩阵，表达式分别为：

$G=J^T{W}_d^T{W}_{d}e\left(m^k\right)+{\mathrm{\η}}_{k}F\left(m^k\right)W_m^T{W}_m(m^k-m^p),$

$H_a=J^T{W}_d^T{W}_{d}J+{\mathrm{\η}}_{k}F\left(m^k\right)W_m^T{W}_m$

$+{\mathrm{\η}}_k[W_m^T{W}_m(m^k-m^p){]}^T{J}^T{W}_d^T{W}_{d}e\left(m^k\right)$

$+{\mathrm{\η}}_k[J^T{W}_d^T{W}_{d}e\left(m^k\right){]}^T{W}_m^T{W}_m(m^k-m^p),$

$F\left(m^k\right)=\frac{1}{2}{\left|\right|W_{d}e\left(m^k\right)\left|\right|}^2$

其中雅可比矩阵J是正演模型对m的偏导数，每次迭代必须重新计算，为减小计算量，本发明推导出其解析表达式。雅可比矩阵元为：

$J=\left[\begin{array}{ccc}\frac{{\∂d}_c\left(t_1\right)}{{\∂m}_1}& \frac{{\∂d}_c\left(t_1\right)}{{\∂m}_2}& \frac{{\∂d}_c\left(t_1\right)}{{\∂m}_3}\\ \·& \·& \·\\ \·& \·& \·\\ \·& \·& \·\\ \frac{{\∂d}_c\left(t_N\right)}{{\∂m}_1}& \frac{{\∂d}_c\left(t_N\right)}{{\∂m}_2}& \frac{{\∂d}_c\left(t_N\right)}{{\∂m}_3}\end{array}\right],$

其中：

m₁=T_n, $m_2=Q_n^{-1},$ $m_3=V_n^{-1},$ $d_c\left(t\right)=T_n\mathrm{Seis}(n,t)\⊗\mathrm{IFFT}\left[P({\mathrm{\Δz}}_n,\mathrm{\ω})\right]$

$\frac{{\∂d}_c\left(t\right)}{{\∂m}_1}=\mathrm{Seis}(n,t)\⊗\mathrm{IFFT}\left[P({\mathrm{\Δz}}_n,\mathrm{\ω})\right],$

$\frac{{\∂d}_c\left(t\right)}{{\∂m}_2}=m_1\mathrm{Seis}(n,t)\⊗\mathrm{IFFT}\left[P({\mathrm{\Δz}}_n,\mathrm{\ω}){\mathrm{\ω\Δz}}_n(\frac{{\mathrm{im}}_3}{\mathrm{\π}}\ln \frac{{\mathrm{\ω}}_r}{\mathrm{\ω}}-\frac{{\mathrm{\ω}}_3}{2})\right],$

$\frac{{\∂d}_c\left(t\right)}{{\∂m}_3}=m_1\mathrm{Seis}(n,t)\⊗\mathrm{IFFT}\left[P({\mathrm{\Δz}}_n,\mathrm{\ω}){\mathrm{\ω\Δz}}_n(\frac{{\mathrm{im}}_2}{\mathrm{\π}}\ln \frac{{\mathrm{\ω}}_r}{\mathrm{\ω}}-\frac{{\mathrm{\ω}}_2}{2})\right],$

这里表示褶积，IFFT表示反Fourier变换，本文中反Fourier变换定义为：

$f\left(t\right)=\frac{1}{2\mathrm{\π}}{\∫}_{-\∞}^{\∞}\hat{f}\left(\mathrm{\ω}\right)\exp (-\mathrm{i\ωt})\mathrm{d\ω}.$

多数方法估计的Q值都会出现负值，而负Q是没有物理意义的。如果已知模型参数的上下限，可以通过非线性变换来限制反演参数的取值范围，使其不超过上下限。这样的变换有很多种，本发明采用如下变换：

$m=\frac{m_{\max }+m_{\min }}{2}+\frac{m_{\max }-m_{\min }}{2}\sin \mathrm{\&theta;}-\&infin;<\mathrm{\&theta;}<+\&infin;,$

其中m_max、m_min分别是m的上限和下限。用变量θ代替模型参数m后，雅克比矩阵中相应的偏导数变为：

$\frac{{\&PartialD;d}_c\left(t\right)}{\&PartialD;\mathrm{\&theta;}}=\frac{\mathrm{dm}}{\mathrm{d\&theta;}}\frac{{\&PartialD;d}_c\left(t\right)}{\&PartialD;m}=\sqrt{(m_{\max }-m)(m-m_{\min })}\frac{{\&PartialD;d}_c\left(t\right)}{\&PartialD;m}$

仍由高斯牛顿法可得到每次迭代时θ的更新量Δθ，将Δθ代入非线性变换式中并整理后可得m的更新迭代公式：

$m^{k+1}=\frac{m_{\max }+m_{\min }}{2}+(m^k-\frac{m_{\max }+m_{\min }}{2})\cos \left(\mathrm{\&Delta;\&theta;}\right)+\sqrt{(m_{\max }-m^k)(m^k-m_{\min })}\sin \left(\mathrm{\&Delta;\&theta;}\right).$

按照步骤4）至7）估计出每层的模型参数，最后即可得到整个Q值随深度变化的曲线，以及整个深度的速度曲线。

Claims (2)

1.一种利用零偏VSP资料估计速度和Q值的自适应波形的反演方法，其特征在于，包括以下步骤：

1）在井中布置检波器采集原始地震数据，然后对采集到的地震数据进行校正处理、噪音压制处理、真振幅恢复处理；处理后得到零偏VSP记录，记为Seis(n,t)，n=1,2,…,N.，其中N表示总的检波器个数，t表示时间轴；

2）计算VSP地震波形Seis(n,t)的振幅包络，将每道记录的初至波振幅包络峰值进行归一化处理，记归一化后的振幅包络为IA(n,t)；并分别拾取每道数据直达波的瞬时振幅峰值到达时刻t_n，用于在步骤4）中计算速度初值；

3）以检波器的位置为界面，将介质划分为N-1个薄层，第n个薄层的厚度为第n个检波器和第n+1个检波器之间的距离，记为Δz_n；则第n个薄层上方的地震记录为Seis(n,t)，下方的地震记录为Seis(n+1,t)，波由界面上方传播到界面下方的频率域正演算子为：

$\mathrm{Seis}(n+1,\mathrm{\&omega;})=T_n\&CenterDot;\mathrm{Seis}(n,\mathrm{\&omega;})\&CenterDot;\exp [\mathrm{i\&omega;}(\frac{{\mathrm{\&Delta;z}}_n}{V_n}+\frac{{\mathrm{\&Delta;z}}_n}{\mathrm{\&pi;}{V}_n{Q}_n}\ln \frac{{\mathrm{\&omega;}}_r}{\mathrm{\&omega;}})-\frac{\mathrm{\&Delta;}{z}_n\mathrm{\&omega;}}{2V_n{Q}_n}$

其中：ω为角频率，ω_r为参考频率，T_n为透射系数，V_n为速度，Q_n为品质因子；

4）对于任意第n个薄层，选定待反演的模型参数为

并给定误差精度ε；比较该层上下两道记录初至波的归一化振幅包络IA(n,t)和IA(n+1,t)，以层下方初至波起始点作为时窗的起始点t_s，以两振幅包络误差大于给定值ε时所对应的时刻作为时窗的终点t_e，截取Seis(n+1,t)中位于[t_s t_e]之间的波形作为观测数据d_o，d_o为含N个采样值的向量；

5）给定一个初始模型参数m⁰，以该模型参数将第n层上方Seis(n,t)的初至波传至层下方，截取[t_s t_e]之间的波形作为计算数据d_c，d_c为含N个采样值的向量；

6）以计算数据d_c与观测数据d_o的误差的模平方为目标函数，并采用乘性自适应正则化方法以改善非线性反演的病态问题，采用数据加权矩阵以提高目标函数对Q值变化的敏感度；

采用乘性正则化方法后的目标函数为：

$E\left(m\right)=\frac{1}{2}{\left|\right|W_{d}e\left(m\right)\left|\right|}^2\&CenterDot;{\mathrm{\&eta;}}_k({\left|\right|W_m(m-m^p)\left|\right|}^2+{\mathrm{\&delta;}}_k^2)$

其中：

$e\left(m\right)=\frac{d_c-d_o}{\left|\right|d_o\left|\right|},$

${\mathrm{\&eta;}}_k=\frac{1}{{\left|\right|W_m(m^k-m^p)\left|\right|}^2+{\mathrm{\&delta;}}_k^2},$

式中k表示迭代次数，W_d是数据加权矩阵，W_m是模型参数加权矩阵，m^p是根据先验信息人为给定的期望模型参数向量，该正则化方法对δ_k值不敏感，所以δ_k取为一个常数，并由数值试验加以确定；

由Q值误差引起的数据误差与观测数据振幅成正比，而由速度误差引起的数据误差与观测数据的振幅成反比；所以，为了提高目标函数对Q的敏感度，本发明采用如下的数据加权矩阵：

$\left(W_d\right)=\mathrm{diag}\{\frac{\left|d_i\right|}{\left|\right|d_o\left|\right|},d_i\&Element;d_o,i=1...N\}.;$

7）利用高斯牛顿法进行迭代反演求取目标函数的极小值，该极小值对应的模型参数即为第n个薄层的最优解。

2.根据权利要求1所述的利用零偏VSP资料估计速度和Q值的自适应波形的反演方法，其特征在于，步骤7）中，迭代过程中采用非线性变换策略将反演的模型参数约束在一定的范围内；高斯牛顿法的迭代公式为：

$m^{k+1}=m^k-H_a^{-1}G,$

式中，k表示迭代次数，G和H_a分别表示目标函数的梯度矩阵和近似Hessian矩阵，表达式分别为：

$G=J^T{W}_d^T{W}_{d}e\left(m^k\right)+{\mathrm{\&eta;}}_{k}F\left(m^k\right)W_m^T{W}_m(m^k-m^p),$

$H_a=J^T{W}_d^T{W}_{d}J+{\mathrm{\&eta;}}_{k}F\left(m^k\right)W_m^T{W}_m$

$+{\mathrm{\&eta;}}_k[W_m^T{W}_m(m^k-m^p){]}^T{J}^T{W}_d^T{W}_{d}e\left(m^k\right)$

$+{\mathrm{\&eta;}}_k[J^T{W}_d^T{W}_{d}e\left(m^k\right){]}^T{W}_m^T{W}_m(m^k-m^p),$

$F\left(m^k\right)=\frac{1}{2}{\left|\right|W_{d}e\left(m^k\right)\left|\right|}^2$

其中雅可比矩阵J是正演模型对m的偏导数，每次迭代必须重新计算。

Images