不等式约束的非线性共轭梯度法研究范文

《石油物探杂志》2016年第3期

摘要：

井中激电三维反演存在多解性且计算时间长，制约了其由理论走向实际应用。针对这一问题，采用非线性共轭梯度（NLCG）反演方法，在反演目标函数中直接施加约束条件压制其解的非唯一性，也就是将介质电阻率取值范围作为先验信息和约束条件以惩罚函数法的方式引入到反演目标函数中。与常规三维电阻率反演目标函数相比，增加了不等式约束项的目标函数理论上可以压制反演的多解性。理论模型的反演结果表明，基于不等式约束的井中激电三维反演方法有效地改善了反演结果的精度，以惩罚函数法施加不等式约束条件的方式现实可行。

关键词：

井中激电；反演；非线性共轭梯度法；不等式约束；惩罚函数法

由SEIGEL［1］理论可知，时间域激发极化的正反演理论建立在电阻率法正反演基础上，直流电阻率法正反演的发展与激发极化法正反演的发展有着紧密联系，因此直流电阻率法的反演方法可以方便地应用到激发极化数据反演中。早期的直流电阻率法反演由于计算速度及内存的限制，大多数反演为一维、二维的最小二乘反演［2－4］；随着计算机的快速发展及数值技术的进步，越来越多的研究针对三维反演［5－10］。反演中网格化的方式、偏导数矩阵的处理方法、反演方法技术的选择是反演问题关键，直接决定着反演的效率和精度。SASAKI［11］利用互换定理计算偏导数矩阵对复杂地电模型进行了反演，但是计算时间依然较长；共轭梯度法能够避免矩阵的直接求解，只需求解矩阵的向量积，大大节省了对计算时间和内存的需求，基于这个优势，ZHANG等［12］使用线性共轭梯度法完成了直流电阻率法三维模型的反演，成为了直流电阻率法三维反演的主要方法；在国内，三维电阻率反演研究也取得了较多的研究成果［13－20］。但反演的计算速度、非唯一性以及稳定性仍然是三维反演的主要难题。SASAKI［11］提出的最小构造约束使相邻网格之间电阻率差异极小，已经成为电阻率反演中常用的约束形式；黄俊革等［18］将体积因子作为光滑约束，使得三维电阻率反演的精度和深部网格的分辨率得到显著提高；宛新林等［21］对粗糙度矩阵进行改进，使之适用于各种情况下光滑度矩阵的求取，将二阶等间距光滑因子进行加权处理，改变光滑因子值大小以控制模型的拟合程度，并用最小二乘正交分解法求解反演方程，改善了三维反演的效果；刘斌等［22］引入自适应加权光滑约束控制深部分辨率，改善了深部网格的反演效果，随后，讨论了不等式约束条件和先验结构约束条件的加入方式［23－24］，有效地去除了反演成像中的假异常和多余构造。非线性问题的线性化［25］往往会造成反问题的不适定性，因此迭代拟合法正逐渐发展为非线性反演方法，其中非线性共轭梯度法（NLCG）不仅减小了方法的不适定性，而且由于其将偏导数矩阵的计算融入梯度计算，因而提高了计算效率。RODI等［26］和NEWMAN等［27］利用非线性共轭梯度法分别解决了大地电磁二维和三维反演问题；胡祖志等［28］使用非线性共轭梯度法实现了大地电磁拟三维反演，计算过程中采用一维灵敏度矩阵来代替三维灵敏度矩阵，提高了反演速度；翁爱华等［29］和林昌洪等［30］分别使用非线性共轭梯度法解决了可控源频率测深的三维反演问题；刘云鹤等［31－32］使用非线性共轭梯度法研究了三维频率域航空电磁法反演；董浩等［33－34］利用非线性共轭梯度法对带地形的三维大地电磁数据反演进行了研究。我们结合非线性共轭梯度反演方法以及ZHANG等［12］提出的偏导数矩阵计算方法，将介质电阻率取值范围作为不等式约束以惩罚函数法的方式加入到反演目标函数中进行井中激电三维反演。利用理论模型验证了基于不等式约束的三维电阻率反演方法的可行性与有效性。

1.基于不等式约束的三维电阻率反演方法

井中激电反演包括电阻率反演和极化率反演两个阶段，极化率反演建立在电阻率反演基础之上，因此激发极化反演首先需要对电阻率进行反演计算，完成电阻率反演后再进行极化率反演。

1．1目标函数依据正则化理论，最小二乘意义下的电阻率反演的目标函数为：Φ（m）＝Φd（m）＋λΦm（m）＝‖Wd［dobs－F（m）］‖22＋λ‖Wmm‖22（1）式中：F（m）为正演响应函数；m（mi，i＝1，2，…，m）为模型参数；dobs（dj，j＝1，2，…，n）为观测数据；Wd为N×N阶加权矩阵；Wd＝diag（1／σ1，1／σ2，…，1／σn），σi为第i个数据的标准差；λ为正则化参数；Wm为M×M阶光滑度矩阵。设反演的观测参数为单极－单极电位φobs，三维反演参数为各网格单元的电导率值σ，由于这两个参数变化范围很大，为了反演的稳定性，通常采用对数来标定反演数据及模型参数，即d＝lnφobs，m＝lnσ。

1．2不等式约束三维电阻率反演中的解通常存在不适定性问题，例如在源点附近会出现奇异点，或者出现负值（采用取对数的方式消除）以及很大的电阻率异常值，为了降低这种不适定性，往往会在程序中人为设定某些参数的上限值，甚至下限值：ρmini≤mi≤ρmaxi（2）式中：mi为第i个网格单元的电阻率；ρmini和ρmaxi分别为第i个网格单元的电阻率下限和上限。这种经验性的在程序中加入上、下限约束的过程和反演计算本身没有任何关系，这种方式也不适于复杂的实际反演情况，所以我们以惩罚函数的形式将不等式约束引入到目标函数中很好地解决了三维电阻率反演的不适定问题。将不等式约束（（2）式）加入到常规的电阻率反演目标函数，构造出带有不等式约束的反演目标函数（不等式约束施加的方式详见附录A）：Z（m）＝Φd（m）＋λΦm（m）＋μP（m）＝［d－F（m）］TC－1d［d－F（m）］＋λmTCTmCmm＋μ［min（m－mmin，0）Tmin（m－mmin，0）＋min（mmax－m，0）Tmin（mmax－m，0）］（3）式中：mmin，mmax为每个单元网格电阻率的上、下限。公式（3）即为本文带不等式约束条件的三维电阻率反演的目标函数，与常规三维电阻率反演目标函数相比，该方法增加了不等式约束项，但表达式简单、容易理解、易于实现。理论上讲，不等式约束对抑制多解性有重要作用，能有效地改善反演效果。

1．3非线性共轭梯度反演非线性共轭梯度反演需要解决的关键问题有3个：①目标函数梯度的计算；②最优步长的计算；③预条件因子的选取与计算。1．3．1梯度的求解从NLCG算法的计算过程来看，需要计算目标函数的梯度，目标函数（（3）式）的梯度表达式为：g＝－2JTC－1d［d－F（m）］＋2λCTmCmm＋2μ［min（m－mmin，0）－min（mmax－m，0）］（4）从（4）式中可以看出，梯度g的计算关键在于雅克比矩阵J的计算，而且只需要计算雅克比矩阵的转置与任意一维向量的乘积JTy，因此无需考虑雅克比矩阵的存储问题，并且每次JTy的计算均可在每次反演计算中的正演计算后一起得到，从而大大加快了反演的计算速度。ZHANG等［12］对JTy计算进行了详细的分析。

1．3．2步长的确定在极小化目标函数Φ（m（i）＋αu（i））时，需要确定搜索步长α，因为只包含一个变量α，所以可以用全局最优化的一维线性搜索方法，如精确线性搜索法、黄金分割法、插值法和非精确线性搜索法等求取。而通常精确线性搜索花费的时间较长，因此一般采用非精确的一维线性搜索。我们采用充分下降条件和二次插值法来求取反演步长，具体的搜索方式见参考文献［31］。非精确线性搜索中的Wofe－Powell准则为：Φ（mk＋αkpk）≤Φ（mk）＋c1αkΔΦ（mk）TpkΔΦ（mk＋1＋αkpk）Tpk≥c2ΔΦ（mk）Tp烅烄烆k（5）式中：Φ为正演算子；mk为第k次迭代的模型参数向量；αk为迭代步长；pk为搜索方向；c1，c2为常数，满足0＜c1＜c2＜1；k为迭代次数。在实际应用中，一般取c1＝0．0001，c2＝0．1000。

1．3．3预条件因子的选择预条件因子能够改善系数矩阵的特征值分布，减少其条件数以提高反演的稳定性，在NLCG法中，NEWMAN等［27］认为近似Hessian矩阵的逆作为预条件因子使得搜索方向具有与牛顿方向类似的性质，因此预条件因子M－1越近似Hes－sian矩阵的逆，越能加速收敛。预条件矩阵可以在每步反演迭代中固定或者更新变化，NEW－MAN等［27］认为Hessian矩阵不是一个常数，因此预条件矩阵也应该时刻更新变化；而计算全Hessian矩阵或者近似Hessian矩阵代价较大，因此只需要计算它的对角元素来作为预条件因子，NEWMAN等［27］选择拟牛顿法中的BFGS算法产生预条件因子，只需要计算并存储Hessian矩阵的对角元素，其对角元素可以依靠目标函数的梯度和上一次迭代的结果给出［27，33］。我们选取的预条件因子为：Mk＋1＝Mk＋gkgTkgTkuk＋ykyTkαkyTkuk（6）式中：gk＝ΔΦ（mk），yk＝ΔΦ（mk＋1）－ΔΦ（mk）。通常第一次迭代时，M1设为单位矩阵I。按公式（6）更新预条件因子能够确保Mk＋1为正定且NLCG法中的搜索方向是下降方向。

2.井中激电的三维反演

根据体激发极化理论，假设σi（i＝1，2，…，M）与ηi（i＝1，2，…M）分别表示地下介质的电导率和极化率，由等效电阻率法可知，激发极化总场φ（σ＊）由等效电导率σ＊＝（1－η）σ决定，且视极化率可表示为：ηs＝φ（σ＊）－φ（σ）φ（σ＊）（7）由（7）式可推导出视极化率ηs与地下介质真极化率ηi之间的近似线性关系：ηs＝－∑Mi＝1lnφlnσiηi＝－Jη（8）其中，J刚好是电阻率三维反演中偏导数矩阵的对数形式，因此根据正则化理论，得到最小二乘意义下极化率反演的目标函数为：Φ（m）＝Φd（m）＋λΦm（m）＝‖Wd（ηs＋Jη）‖22＋λ‖Wmη‖22（9）对公式（9）求极小，即令Φη＝0，有：g（η）＝Φη＝2JTWTdWd（ηs＋Jη）＋2λWTmWmη＝0（10）（JTWTdWdJ＋λWTmWm）η＝－JTWTdWdηs（11）可以看出，（11）式满足Ax＝b的形式，因此（11）式可以直接利用线性共轭梯度法求解。（11）式中各项参数和电阻率反演中的各项参数相同。

3.数值模型反演算例

为了验证方法的可靠性和计算效率，我们给出了几个理论模型来进行反演计算。三维非线性共轭梯度反演算法具体流程如图1所示。需要特别说明的是，本文在目标函数中施加不等式约束均是在电阻率三维反演中实现的，极化率的反演直接利用线性共轭梯度法反演完成。流程图中各参数含义如下：初始预处理矩阵H0一般取为单位矩阵；收敛系数ε1为正则化因子的阈值；收敛系数ε2为前、后两次相对拟合差Rms的绝对值，即目标函数值的下降量小于某个阈值时，说明目标函数值不再下降，达到极小值，因为理论数据也不能使数据拟合差降低到0；收敛系数ε3为目标函数的相对拟合差阈值；收敛系数ε4为目标函数的梯度值的阈值，如果目标函数的梯度很小，目标函数值也达到极小值。相对拟合差Rms定义如下：Rms＝∑Mj＝1dobsj－d0jdobsj2槡M（12）式中：M是观测数据的个数；dobsj和d0j分别为合成的观测数据与正演理论值。值得注意的是，在大多数情况下用公式（12）的值来表示三维电阻率反演精度［35］，而屈有恒［36］认为其值的大小只对反演的过程有意义，不能用来衡量反演精度，任何一个给定的初始模型的反演精度的估算都是相对的，因此对于给定的初始模型，相对拟合差只要能下降到初始模型误差的1／20～1／10即可认为达到了较高的精度。

3．1算例一：2个低阻立方体模型首先建立三维反演模型，然后利用常规反演目标函数（公式（1））对电阻率进行非线性共轭梯度反演计算，再进行极化率的反演。模型如图2所示，在电阻率为ρ0＝100Ω•m，极化率为η0＝0的均匀半空间中存在2个大小为20m×20m×20m（边长a＝20m），电阻率ρ1＝ρ2＝10Ω•m，极化率η1＝η2＝0．1的低阻体，顶面埋深均为h＝5m，底面埋深25m，异常体边界离x轴与y轴距离相等且距离d＝5m，红色圆点代表井口所在位置，井口坐标为（50m，50m，0）。分别在井中不同深度布设多个源：A0（z＝0），A1（z＝－10m），A2（z＝－20m），A3（z＝－30m），A4（z＝－40m），A5（z＝－50m），在这些源处固定A极供电；观测电极距2m，在地面51条测线上逐点移动M极，通过三维有限元正演程序共得到15606个一次场电位数据和视极化率的“实测数据”。为了提高反演速度，将整个区域分为边界区域和目标区域，目标区域选为100m×100m×100m，扩展区域分别向3个方向延展200m。反演初始模型选取电阻率ρ＝100Ω•m，极化率为η＝0的均匀半空间。反演参数选择如下：正则化因子初始值l0＝0．05，减小系数q＝0．1，给定反演终止收敛系数ε2＝ε3＝1．0×10－5。图3为模型的电阻率与极化率反演结果。图4为y＝40m和y＝60m处的电阻率反演切片（xoz垂直断面）。图5为y＝40m和y＝60m处的极化率反演切片（xoz垂直断面）。低阻异常体反演的电阻率数值基本接近真实值，高阻体差异较大；高阻体的极化率数值接近真实值；反演结果基本圈定了高、低阻异常体的水平中心位置和边界，反演的高低阻异常体的中心位置都有不同程度的向上偏移现象。为了对比施加约束条件前后的反演效果，我们应用基于不等式约束的反演目标函数（（13）式）对电阻率进行非线性共轭梯度反演计算。首先，利用有限元法对模型进行正演模拟，井中点电源供电点一共11个，点源间隔为5m，观测电极距2m，在地面51条测线上逐点移动M极，通过三维有限元正演程序共得到28611个视电阻率的“实测数据”。以点电源深度为纵轴，得到视电阻率断面图（图6），图6中虚线方框分别为低阻异常体在xoz平面上的投影。由视电阻率断面图可以看到观测数据中的视电阻率基本在120Ω•m以内，因此将整个模型中网格单元的电阻率上限设为150Ω•m（有一定冗余），下限设为0。即mmin＝0，mmax＝150Ω•m。为了加快反演计算速度，提取A0（z＝0），A1（z＝－10m），A2（z＝－20m），A3（z＝－30m），A4（z＝－40m），A5（z＝－50m）6个供电点处15606个一次场电位数据进行反演。计算区域与模型大小的设置、反演参数与单个低阻立方体模型一样，其中惩罚因子μ取为常数1，反演结果见图7至图9。由未施加不等式约束条件下的反演结果（图3至图5）可看出，图4电阻率反演切片图中源点附近出现了较大的电阻率值（大于200Ω•m），导致整个三维反演成像效果较差，尤其对2个低阻异常体的成像精度较低。而在加入不等式约束条件后的反演结果（图7至图9）中，从图7中可以看出两个低阻异常体的成像效果大大提高，在规模、形状、电阻率值等特征方面与原模型基本一致。

3．2算例二：单个低阻立方体模型模型如图10所示，在极化率η0＝0，电阻率ρ0＝100Ω•m的均匀半空间中存在1个大小10m×10m×10m（a＝10m），电阻率ρ1＝10Ω•m，极化率η1＝0．1的低阻低极化体，顶面埋深h＝45m，底面埋深55m，立方体两边距离x轴与y轴相等，d1＝d2＝5m，红色圆点代表井口，井口坐标为（50m，50m，0）。4个方位电极点离井口20m；网格剖分单元总数为289280，节点总数为49797。采用二级装置进行地－井五方位观测，分别在地面井口方位A0（50m，50m，0），主方位A1（30m，50m，0），A2（70m，50m，0），辅助方位A3（50m，30m，0），A4（50m，70m，0）处固定A极供电进行五方位观测，井中观测电极距2m，在井中逐点移动M极，通过三维有限元正演程序共得到255个一次场电位数据和视极化率的“实测数据”。最常用的地井激化法工作方式是地井五方位观测，由于地井方式单次观测获取的观测数据量非常少，所以三维反演有很强的多解性，很难实现精确的成像反演，较常规的反演方法是通过正演拟合的方式针对某些特定的参数进行最优化反演；而井地方式中，由于测点在地面，可以得到更多的观测数据，获取更多的地下信息，能够更好地进行三维成像反演。地井方式利用单次观测的数据进行三维电阻率和极化率的反演结果多解性较为严重，本文尝试采取两种途径来减小传统反演方式带来的多解性问题：1）将地井五方位观测的数据进行同时反演，增加观测数据的个数；2）增加更多的先验信息约束。整个区域分为边界区域和目标区域，目标区域选为100m×100m×100m，扩展区域分别向3个方向延展200m。反演初始模型选取电阻率ρ＝100Ω•m，极化率η＝0的均匀半空间。反演参数选择：正则化因子初始值λ0＝0．05，减小系数q＝0．1，给定反演终止收敛系数ε2＝ε3＝1．0×10－5。首先进行地井五方位观测数据的同时反演，不施加约束条件，利用非线性共轭梯度法对模型进行电阻率反演计算，图11与图12为地井五方位观测数据的同时反演电阻率结果。从图12中可以看出，基本能反映有低阻异常体的存在，异常体的中心位置基本吻合，无多余构造产生，但是异常体的边界范围模糊，延伸较远；由于异常体埋深较大，在地面供电的情况下，地下电性异常体产生的电位异常幅值很小，因此反演出来的电阻率值与真实电阻率值相差较大。结合五方位观测数据，将不等式约束引入到反演目标函数中。观测数据中的视电阻率都在110Ω•m以内，因此将整个模型中网格单元的电阻率上限设为150Ω•m，下限设为0；在使用电阻率探测时，经常会配合其它地球物理探测方法（如浅层地震勘探法、地质雷达探测法），这些方法对异常体界面识别和定位具有较好的效果，由此获得的异常体的位置与形态的信息是一种极为重要的先验构造信息。因此将这种构造位置的先验信息作为局部约束加入到三维反演目标函数中。假设由其它地球物理探测方法可知异常体的构造位置与形态信息，由探区矿石标本物性参数可得：55m≤x≤65m，55m≤y≤65m，－55m≤z≤－45m处已知异常体的电阻率值在50Ω•m以内。因此将这处网格单元的电阻率上限设为50Ω•m，下限设为0。惩罚因子μ取为常数1。图13至图15为施加不等式约束后的电阻率与极化率的反演结果，从反演结果可以看出，异常体水平和纵向分辨率较五方位同时反演明显提高，在位置、规模、形状等特征方面均与原模型基本一致，验证了局部不等式约束能够得到更为精确的三维反演成像结果。

4.结论

1）我们以惩罚函数法的方式将介质电阻率取值范围的不等式约束施加到三维井中激电反演目标函数中，改善了反演的效果，因此本文提出的以惩罚函数法施加不等式约束条件的方式是可行的，为约束条件的施加方式提供了一条新的途径；2）在地井观测模式中，我们尝试采取两种途径来减小传统反演方法带来的多解性问题，理论模型反演试算表明，利用构造位置的先验信息作为局部约束加入到三维反演目标函数中得到的反演结果更为精确，反演的多解性显著降低；3）惩罚因子的取值是本文的难点，本文惩罚因子的取值采用的是实验法，过大或者过小的惩罚因子都会影响反演的效果，有时甚至会出现不收敛的情况，因此给出一个合理的惩罚因子非常重要；4）本文的非线性共轭梯度反演算法计算依然很耗时，在反演多个点源多次观测数据的情况下显得更为突出，下一步需要研究并行反演算法进行计算，实现反演算法的快速计算。

参考文献：

［4］苏朱刘，胡文宝，严良俊．电阻率和极化率测深法的正演修正法反演［J］．石油物探，2005，44（2）：194－198

［10］张树彬，孙东，郭秀娟．航空电磁法三维电阻率反演成像技术［J］．石油物探，2006，45（1）：98－101

［13］吴小平，汪彤彤．电阻率三维反演稳定性和可靠性研究［J］．煤田地质与勘探，2002，30（1）：57－60

［14］吴小平，徐果明．利用共轭梯度法的电阻率三维反演研究［J］．地球物理学报，2000，43（3）：420－427

［15］吴小平，徐果明．电阻率三维反演中偏导数矩阵的求取与分析［J］．石油地球物理勘探，1999，34（4）：363－372

［16］吴小平，汪彤彤．利用共轭梯度方法的电阻率三维反演若干问题研究［J］．地震地质，2001，23（2）：321－327

［17］吴小平．利用共轭梯度方法的激发极化三维快速反演［J］．煤田地质与勘探，2004，32（5）：62－64

［18］黄俊革．三维电阻率／极化率有限元正演模拟与反演成像［D］．长沙：中南大学，2003

［19］黄俊革，阮百尧，鲍光淑．基于有限单元法的三维地电断面电阻率反演［J］．中南大学学报（自然科学版），2004，35（2）：295－299

［20］强建科．起伏地形三维电阻率正演模拟与反演成像研究［D］．武汉：中国地质大学（武汉），2006

［21］宛新林，席道瑛，高尔根，等．用改进的光滑约束最小二乘正交分解法实现电阻率三维反演［J］．地球物理学报，2005，48（2）：439－444

［22］刘斌，李术才，聂利超，等．基于自适应加权光滑约束与PCG算法的三维电阻率探测反演成像［J］．岩土工程学报，2012，34（9）：1646－1653

［23］刘斌，李术才，李树忱，等．基于不等式约束的最小二乘法三维电阻率反演及其算法优化［J］．地球物理学报，2012，55（1）：260－268

［24］刘斌，聂利超，李术才，等．三维电阻率空间结构约束反演成像方法［J］．岩石力学与工程学报，2012，31（11）：2258－2268

［25］吴小平，汪彤彤．电阻率三维反演方法研究进展［J］．地球物理学进展，2002，17（1）：156－162

［28］胡祖志，胡祥云，何展翔．大地电磁非线性共轭梯度拟三维反演［J］．地球物理学报，2006，49（4）：1226－1234

［29］翁爱华，刘云鹤，贾定宇，等．地面可控源频率测深三维非线性共轭梯度反演［J］．地球物理学报，2012，55（10）：3506－3515

［30］林昌洪，谭捍东，舒晴，等．可控源音频大地电磁三维共轭梯度反演研究［J］．地球物理学报，2012，55（11）：3829－3838

［31］刘云鹤，殷长春．三维频率域航空电磁反演研究［J］．地球物理学报，2013，56（12）：4278－4287

［32］刘云鹤．三维可控源电磁法非线性共轭梯度反演研究［D］．吉林长春：吉林大学，2011

［33］董浩，魏文博，叶高峰，等．基于有限差分正演的带地形三维大地电磁反演方法［J］．地球物理学报，2014，57（3）：939－952

［34］董浩．基于有限差分法正演的大地电磁测深带地形三维反演研究［D］．北京：中国地质大学（北京），2013

［35］胡朝俊．三维直流电法共轭梯度反演算法研究［D］．北京：中国地质大学（北京），2006

［36］屈有恒．井地电阻率法及双频激电三维数值模拟与反演研究［D］．北京：中国地质大学（北京）

作者：王智 潘和平 吴爱平 李刚 方思南 单位：长江大学电子信息学院 中国地质大学（武汉）地球物理与空间信息学院