师学明，王家映

师学明，王家映_大人地球物理科学

师学明，王家映

（中国地质大学地球物理与空间信息学院，武汉，430074）

摘　要：详细研究了一般地球物理反演问题的迭代优化求解过程与物理学中原子跃迁过程的对应关系，建立了反演问题中模型空间、初始模型、局部极值模型、最优化模型等与原子的态空间、定态、激发态、基态等的对应关系。在此基础上，模拟了物理学中原子从激发态向基态跃迁的物理过程，建立了一种与原子跃迁过程相对应的非线性随机跃迁数学模型和模型解跃迁搜索准则，导出了适用于一般地球物理资料的模拟原子跃迁的非线性反演算法。用理论测试函数和一维层状大地电磁资料对这种新的反演方法进行了数值试验，结果表明该方法具有解不依赖于初始模型、收敛速度快等优点。

关键词：地球物理，非线性反演，全局最优化，局部最优化，模拟原子跃迁反演法

1　引言

20世纪60年代末开始，由于Backus和Gilbert等地球物理学家^［1～6］的不懈努力，地球物理的反演理论特别是广义线性反演理论已趋成熟，并在地球物理各种勘探方法中得到了广泛应用^［7～11］。然而，由于大多数的地球物理反演问题都是非线性的，因此严重地限制了线性反演理论的应用范围。为此，一些非线性法逐步被引入到地球物理反演中来。

非线性最优化反演法基本上可以分为两类：一类是线性化迭代方法^［12］，如最速下降法、牛顿法、共轭梯度法和非线性最小二乘法，以及近年来发展起来的用于解决三维反演问题的不完全分解共轭梯度法^{［13，14］}。由于梯度的局部特性，使得这类方法在目标函数多极值的情况下，反演容易陷入局部极值，其解严重依赖于初始模型，所以它和线性反演方法一样仍然是一种局部化反演方法。另一类非线性反演法是模拟自然界某种物理过程或某种生物现象的“启发式”反演法，例如模拟“赌博”的蒙特卡洛法^［15］、模拟物体从高温逐渐冷却结晶过程的模拟退火法^{［16～19］}、模拟生物从低级到高级进化过程的遗传算法^{［20～23］}和模拟人脑神经元工作原理的人工神经网络方法技术^{［24，25］}等。这类非线性反演法与局部最优化反演方法依赖于目标函数的梯度信息不同，它是在模型空间中随机地全局搜索，因而不会陷入到局部极值。另外这类方法还具有易于施加先验约束条件、适于并行计算等优点，因而获得了地球物理学家的青睐^{［26，27］}。

通过对理论物理学中原子跃迁过程^{［28，29］}与地球物理反演中的反演过程的深入对比分析与思考，王家映^［30］首次提出可以模拟原子跃迁过程来建立反演算法的思想。在此基础上，我们对理论物理学中原子跃迁的物理过程与地球物理反演过程的对应关系进行了分析，研究了地球物理反演问题的定态，研究了基于原子跃迁理论的非线性随机跃迁最优化数学模型和模型解随机跃迁搜索准则，建立了一种新的适用于一般地球物理资料反演问题的非线性反演新方法——模拟原子跃迁反演法，并用理论模型和实际资料的反演进行了数值试验。

2　模拟原子跃迁反演算法的基本原理

2.1　基本思想

当原子处于激发态时，原子有较高的能级，电子在离核较远的轨道上运动。原子从较高能级的激发态向较低能级的激发态或基态跃迁的过程，是辐射能量的过程。这时，电子从离原子核较远的轨道（激发态）回到离原子核较近的轨道（激发态）或最近的轨道（基态）上运动，这种跃迁要释放能量。显然，当电子跃迁到能量较低的激发态时，原子仍不稳定，只有跃迁到能量最低的基态时，才是最稳定的^［28，29］。地球物理资料的反演，是一目标函数（或方差函数）最优化为最小（或极小）的过程。当方差函数达到全局极小时，反演就终止。

原子跃迁的过程就是原子从能量较高的激发态向能量最低的基态的跃迁过程，它是一个原子自身寻求能量最低分布的最优化过程，此过程与地球物理反演的最优化过程基本相似。因此，可以模拟原子跃迁，建立一种新的非线性最优化反演方法，用于地球物理资料的处理与解释。

2.2　地球物理反演过程与原子跃迁过程的对应分析

根据模拟原子跃迁反演算法的基本思想，对地球物理反演过程和原子跃迁过程的一些基本概念、基本过程进行类比分析，可以建立原子跃迁过程的一些基本量和地球物理反演过程的基本量之间的相互关系，如表1所示。

表1　地球物理反演过程与原子跃迁过程的对应关系表

2.3　地球物理反演问题的定态

从表1可以看出，模拟原子跃迁反演算法的思想，是将地球物理反演问题的目标函数值等价于原子的能量，将目标函数朝减小的方向移动，等价于原子从较高能级的激发态向较低能级的激发态和基态跃迁。众所周知，在原子跃迁的量子理论中，原子只能处在一系列不连续的、独立的能级上。因此，我们首先需要采用某种算法，将地球物理反演问题的目标函数转换为相应于原子的离散的能级。

事实上，地球物理的局部最优化方法例如线性化迭代方法，是在给定一个初始模型的情况下，根据目标函数的梯度逐步向局部极小值的迭代优化过程。随着迭代的进行，初始模型将逐步向局部极值对应的模型逼近。当解陷入局部极值时，目标函数值将稳定不变。此时，可将该目标函数值的大小看作是“原子的能级”，将该局部极值所对应的模型看作是“原子的定态”^［31］。

图1（a）是典型的地球物理反演问题的目标函数示意图。由图可见，地球物理反演问题的目标函数在模型空间存在许多局部极值。因此，用基于梯度的线性化迭代反演方法，如果给定的初始值在局部极小附近，解就会逐步迭代收敛到当前的局部极小值；如果给定的初始模型在全局极小附近，就可以很快地收敛到全局极小，如图1（b）所示。

对于一个地球物理反演问题，假定目标函数的局部极值点的数目是有限个，假定在目标函数每个局部极值点附近有一个初始模型，对这些初始模型都采用线性迭代反演算法并使其收敛，将得到与目标函数相对应的目标能级曲线图。换句话说，也就可以得到该地球物理反演问题的“能级”，如图1（c）中粗黑线所表示。

图1　地球物理反演问题的目标函数（a）、线性化反演过程（b）及目标函数“能级化”（c）示意图

E（x）为目标函数；x为模型空间；E₁，E₂，…，E₅为能级化目标函数记地球物理反演问题的目标函数为E（x），其中x为模型参数，那么目标函数的能级化数学模型为^［32］：

式中：^～E（x^＊）称为目标函数的能级值或称为“目标能级”，x为任意模型值，x^＊为局部极小对应的模型值，也即是所谓的“定态”。

显然，目标能级具有下列性质。

（1）目标能级小于或等于反演问题的目标函数值。

（2）目标能级没有改变目标函数局部极小值和全局极小值的大小。

假定地球物理反演问题的局部极值点是有限个，那么对地球物理反演问题的目标函数进行能级化的结果是：目标能级是一系列独立的、离散的数值，这可以类比为原子跃迁中原子的能级。目标能级所对应的模型空间中的模型，进行局部最优化（线性化迭代反演法或非线性方法）反演，得到局部最优化模型或全局最优化模型。这些局部最优化模型可类比为原子的激发态，全局最优化模型可类比为原子的基态；地球物理反演问题的目标函数局部最优化模型和全局最优化模型所张的空间，可类比为原子的态空间。地球物理反演问题的一些基本概念与原子跃迁的一些基本概念的对应关系见表2。

表2　地球物理反演的一些基本概念与原子跃迁的一些基本概念的对应关系

2.4　模拟原子跃迁反演算法的设计

从原理上讲，只要初始模型群体的数目足够多，目标函数的局部极小值和全局极小值所对应的模型空间至少有一个初始模型，则采用局部线性化反演算法可以一次求得该反演问题的“激发态”和“基态”，采用排序或查找的方法就可以找到反演问题的“基态”。

但是，由于受到计算机的内存容量与速度的限制，对于一般的地球物理反演问题，初始模型群体的数目不可能太大，一般在几十个到几百个。这样，初始模型往往不能遍布整个目标函数的局部极小值和全局极小值区域，一次迭代后反演的地球物理参数可能还位于局部极小值，即还处在能级较高的激发态，还需要继续在模型空间里寻找目标函数的基态。因此，必须采用模拟原子跃迁的方法，让原子继续在态空间进行演化，寻找对应反演问题的“基态”。

在原子跃迁过程中，既有原子从较高能级的激发态向较低能级的激发态和基态的跃迁过程，也有从较低能级的基态或者激发态向更高能级的激发态的跃迁过程。根据原子物理学，在平衡态，各能级的原子数的分布服从玻尔兹曼（Boltzmannn）分布^［29］，即某个能级E_i原子跃迁到另外一个目标能级Ej的跃迁几率可以采用下列公式进行计算：

式中：k为Boltzmann常数，T为体系的绝对温度。

类似地，模拟原子跃迁的反演方法，在求出地球物理反演问题的定态之后，就可以根据定态对应的目标能级值大小进行模拟原子跃迁。具体地，首先在局部最优解中选择一个模型作为m_i，然后随机地在局部最优解中选择另外一个模型作为m_j。如果E（m_i）＞E（mj），表示m_j模型对应的目标能级小，此时模型m_i跃迁到模型m_j。如果E（m_i）≤E（m_j），表示m_i模型对应的目标能级小，此时，要根据两个模型的能级差来计算跃迁的概率，即概率p＝exp｛-［E（m_j）-E（m_i）］/kT｝。然后随机给定一个概率p_r（0＜p_r＜1），如果p＞p_r，则模型m_i跃迁到模型m_j；否则，不发生跃迁。上面计算中，玻耳兹曼常数一般设为1，温度T为控制参数。

地球物理反问题的定态模型，按照上面的方式进行模拟原子跃迁后，不同能级上模型个数分布开始发生改变，模型群体开始向能级较低的激发态和基态演化。模拟原子跃迁之后，模型群可能还处在能级较高的激发态上。为了寻找“基态”，一种可行的方法是重新给定初始模型，重复上述的过程，直至找到反演问题的基态。但这种重新启动的迭代反演法的缺点是本次迭代与上一次迭代之间没有反演信息的继承。为了利用上一次迭代反演得到的信息，第二次迭代的初始模型是在第一次迭代反演的定态模型基础上进行某种“粒子轰击”得到的新模型。所谓粒子轰击，就是对原子系统进行高能照射，使原子从较低的能级跃迁到较高的能级的物理过程。对处于定态的地球物理反演问题进行“粒子轰击”之后^［28］，模型参数就从当前的定态值跃迁到了另外一个值。进行“粒子轰击”的机制有很多，其目的是使原子从某个定态跃迁到另外一个定态。因此，模型参数将发生一个随机的扰动，即

式中：x^＊为“定态”解，δx为一个随机数。

扰动的新模型对应的是哪一个定态，还必须通过第二次迭代过程中的局部最优化方法求取才能知道。通过这样的一个“粒子轰击”的物理过程，就可以得到第二次迭代的初始模型。第二次迭代的初始模型，既利用了上一次迭代的反演参数值与目标能级，又创造出了新的反演参数值，新的反演参数值可能对应了目标函数中新的目标能级。有了第二次迭代的初始模型，就可以仿照第一次迭代的过程进行模拟原子跃迁反演。如此反复迭代，直至迭代达到设定的条件。

综上所述，模拟原子跃迁反演的算法步骤为：

第一步：随机给定一初始模型群体；

第二步：对每个模型进行线性化迭代反演；求取目标能级（局部极值）及其对应的定态（反演参数）。如果目标能级小于给定的阈值，转到步骤五；

第三步：根据目标能级，进行模拟原子跃迁；

第四步：对原子跃迁之后的定态，进行“粒子轰击”，产生新的模型群体，转步骤二；

第五步：打印目标能级和反演参数；

第六步：结束。

2.5　算法分析

在目标能级空间里，采用模拟原子跃迁进行反演的突出优点是：模型参数可以很容易地穿越目标函数空间中的“能量壁垒”。如图2所示，在目标函数空间里，假设目标能级值E₃＞E₁。当初始模型在目标能级为E₃的模型空间内，采用线性算法所求得的解都在定态x^＊₃，而很难穿越E₃和E₁之间的能量壁垒而到达定态x^＊₁。但是采用模拟原子跃迁算法，由于目标能级值E₃＞E₁，因此采用模拟原子跃迁算法，高能级的激发态能自发地跃迁到低能级的激发态，因此定态解x^＊₃很容易跃迁到定态解x^＊₁。即在目标能级空间里，反演的模型参数在模型空间里可以很容易从能量值较高的E₃跃迁到能量值较低的E₁态，即可以很容易穿越目标函数极小值与极小值之间的“能量壁垒”，迅速达到另一个能级更低的局部极小值和全局极小值。

同时，能级跃迁还有原子从低能级的激发态向高能级的激发态的跃迁过程，即原子吸收能量而发生跃迁。例如图2中地球物理反演参数从目标能级E₃→E₁的跃迁过程就相当于原子的自发跃迁过程，而反演参数从目标能级E₂→E₃的跃迁过程就相当于原子吸收能量发生跃迁的过程。这个过程体现了模拟原子跃迁算法具有非线性寻优的特点。

不同的原子结构，其激发态和基态的能级与结构是不同的。对应地，不同的地球物理反演问题，由于反演问题的性质和观测误差等因素的影响，目标函数的局部极小值的分布即反演问题“定态”的分布是很不一样的。因此，采用模拟原子跃迁算法，追踪具体地球物理反演问题的模拟原子跃迁过程，可定性地分析模型空间的结构。

另外，与模拟退火、遗传算法等完全非线性算法相比，模拟原子跃迁算法的优点是收敛速度相对较快。这是由于在模拟原子跃迁算法中采用了线性优化算法寻找局部极值和全局极值。当解在全局极小值附近时，模拟原子跃迁算法能很快地找到真解，而模拟退火、遗传算法由于搜索的随机性而收敛较慢。

图2　模拟原子的目标能级跃迁示意图

（光滑实线为目标函数曲线；台阶实线为能级化的目标函数曲线；虚线为能级跃迁示意线）

3　数值试验

为了检验模拟原子跃迁反演算法的有效性，用MATLAB语言进行了程序设计，并进行了大量的数值模拟试验，均取得了较好的结果。下面给出几个典型反演问题的试验过程和结果。

3.1　测试函数一：单极值最优化模型

简单的平方和函数：

该函数是一条抛物线，只有一个全局极小值点，f（2.5×10^-4）＝0。

模拟原子跃迁算法的初始模型数目是10个，线性化迭代反演算法采用牛顿法。由于在模型空间范围内只有一个全局极值点，10个初始模型线性化反演的结果都收敛于全局极小解，线性迭代后只有一个“定态”，即x＝0.000 25，对应的目标函数值f（0.000 25）＝0。由于该值为全局极小值，因此该“定态”是该反演问题所对应的“基态”，即这个问题只需要进行一次模拟原子跃迁反演，就可以搜索到反演问题所对应的“基态”。

3.2　测试函数二：多极值最优化模型

多极值函数：

该函数只有一个全局极小值点，f（0）＝0，有10个局部极值点，函数曲线如图3（a）所示。

采用模拟原子跃迁反演算法，初始模型取10个；线性迭代方法采用共轭梯度算法，最大线性迭代次数设为30；模拟原子跃迁反演的最大迭代次数设为10。模拟原子跃迁反演法的计算结果如表3所示，误差下降曲线如图3（b），其中的迭代次数和计算时间是经过10次运行的平均值。为了对比其他非线性反演的计算效率，对该函数采用了蒙特卡洛和模拟退火反演法，见表3。从运行时间上看，蒙特卡洛平均需要2.454s，模拟退火需要0.636s，而模拟原子跃迁算法只需要0.040s，所需时间明显低于上述两种非线性算法；从运行结果来看，蒙特卡洛和模拟退火算法的最优解已经达10-6量级，非常接近于全局最优解，而模拟原子跃迁算法由于采用了共轭梯度算法，可以准确地求取全局极小值，结果要优于其他两种非线性算法。

图3　测试函数二的曲线图（a）及模拟原子跃迁的误差下降曲线（b）

表3　测试函数二几种算法结果的比较

3.3　三层（H型）大地电磁测深数据的模拟原子跃迁反演

表4　三层（H型）模型地电参数和模拟原子跃迁算法反演结果

表4给出了三层（H型）地电断面的参数及模拟原子跃迁反演后的结果。首先对没有加随机噪声的数据进行了反演。为了模拟实际的大地电磁观测数据，还分别对添加了1%、5%和10%的随机误差的数据进行了反演。模型群体取10个，模拟原子跃迁迭代次数取5。线性反演采用阻尼最小二乘法，线性反演的最大迭代次数为100次，误差阈值为1.0。不同噪声情况下，采用模拟原子跃迁算法均得到了较好的反演结果。图4是加10%随机噪声情况下的理论曲线和模拟原子跃迁反演法得到的最优模型的拟合曲线。由图可以看出，两者的拟合程度非常好。

图4　模拟原子跃迁反演法对三层（H型）地电模型反演结果对比图

（a）视电阻率曲线对比图；（b）阻抗相位对比图

4　结论与讨论

笔者提出了一种新的地球物理反演方法——模拟原子跃迁反演法，并给出了模拟原子跃迁的算法步骤，同时用几个理论测试函数进行了数值模拟。结果表明，模拟原子跃迁反演法是一种计算上比较现实的地球物理反演方法，用它可以较快地逼近全局最优解，从而既优于传统的线性或者近线性的局部最优化方法，避免了选取初始模型、陷入局部极值等缺点；又克服了蒙特卡洛、模拟退火反演等方法模型搜索时间长的弱点。模拟原子跃迁的全局收敛性在理论上还有待于进一步的推导和证明。

模拟原子跃迁反演优化算法，是从一个初始模型群出发，寻找反演问题的激发态和基态的过程，在算法的搜索过程中具有非线性特点，可以适用于非线性地球物理反演问题。因此，开展对模拟原子跃迁反演法的进一步研究和对其他地球物理资料进行模拟原子跃迁反演，具有一定的现实意义。

参考文献

［1］Backus G E，Gilbert E.Numerical applications of a formalism for geophysical inverse problems［J］.Geophys.J.R.Astr.Soc.，1967，13（1）：247-267.

［2］Jackson D D.Interpretation of inaccurate，insufficient and inconsistent data［J］.Geophys J.R.Astr Soc.，1972，28：97-109.

［3］Parker R L.Understanding inverse theory［J］.Ann.Rev.Earth.Plan.Phys.，1977，5（1）：35-64.

［4］Tikhonov A V，Arsenin V Y.Solution of ill-posed problems［M］.New York：Wiley，1977.

［5］Oldenburg D W.Funnel functions in linear and non-linear appraisal［J］.J.Geophys.Res.，1983，88：7 387-7 398.

［6］Tarantola A.Inverse problem theory：methods for data fitting and model parameter estimation［M］.New York：Elsevier Science Publishing Company Inc.，1987.

［7］何宝侃，周熙襄，钟本善.地球物理反问题中的最优化方法［M］.北京：地质出版社，1980.

［8］栾文贵.地球物理中的反问题［M］.北京：科学出版社，1989.

［9］杨文采.地球物理反演的理论与方法［M］.北京：地质出版社，1997.

［10］王家映.地球物理反演理论［M］.北京：高等教育出版社，2002.

［11］李世雄，刘家琦.小波变换和反演数学基础［M］.北京：地质出版社，1994.

［12］杨文采.非线性地球物理反演讲座之四——非线性地震反演方法的补充及比较［J］.石油物探，1995，43（4）：110-116.

［13］Mackie R L，Madden T R.Congjugate gradient relaxation solutions for three-dimensional magnetotel-luric modeling［J］.Geophysics，1993，58：1 052-1 057.

［14］吴小平，汪彤彤.利用共轭梯度算法的电阻率三维有限元正演［J］.地球物理学报，2003，46（3）：428-432.

［15］姚姚.蒙特卡洛非线性反演方法及应用［M］.北京：地质出版社，1997.

［16］Kirkpatrick S，Gelatt C D，Vecchi M P.Optimization by simulated annealing［J］.Science，1983，220（4598）：671-680.

［17］Sen M K，Stoffa P L.Nonlinear one-dimensional seismic waveform inversion using simulated annea-ling［J］.Geophysics，1991，56：1 624-1 638.

［18］姚姚.地球物理非线性反演模拟退火法的改进［J］.地球物理学报，1995，38（5）：643-650.

［19］师学明，王家映.一维层状介质大地电磁模拟退火反演法［J］.地球科学——中国地质大学学报，1998，23（5）：542-545.

［20］Holland J H.Adaptation in natural and artificial systems［M］.Michigan：University Michigan Press，1975.

［21］Stoffa P L，Sen M K.Non-linear multiparameter optimization using genetic algorithm：Inversion of plane-wave seismograms［J］.Geophysics，1991，56：1 794-1 810.

［22］师学明，王家映，张胜业，等.多尺度逐次逼近遗传算法反演大地电磁资料［J］.地球物理学报，2000，43（1）：122-130.

［23］姚姚.用非线性遗传反演算法进行自动静校正［J］.地球科学——中国地质大学学报，1996，21（1）：90-92.

［24］Hopfield J J，Tank D W.Neural computation of decision in optimization problems［J］.Biological cybernetics，1985，52（3）：141-152.

［25］Zhang Y，Paulson K V.Magnetotelluric inversion using regularized hopfield neural networks［J］.Geophysical prospecting，1997，45：725-743.

［26］Chunduru R K，Sen M K，Stoffa P L.Hybrid optimization methods for geophysical inversion［J］.Geophysics，1997，62：1 196-1 207.

［27］杨文采.评地球物理反演的发展趋向［J］.地学前缘，2002，9（4）：389-396.

［28］关洪.量子力学基础［M］.北京：高等教育出版社，1999.

［29］赵凯华，罗蔚茵.量子物理［M］.北京：高等教育出版社，2001.

［30］王家映.地球物理资料的原子跃迁反演法（手稿）.2000.

［31］Wales D J，Scheraga H A.Global optimization of clusters，crystals and biomolecules［J］.Science，1999，285：1 368-1 372.

［32］陈乐寿，刘任，王天生.大地电磁测深资料处理与解释［M］.北京：石油工业出版社，1988.