利用改进模拟退火算法反演甘肃兰州地区地幔转换带导电模型

翁爱华，郭峻豪

(吉林大学 地球探测科学与技术学院，吉林 长春 130026)

0 引 言

兰州地区地处欧亚大陆中部，其西部可能存在塔里木地幔柱[1-2]，而南部可能存在峨眉地幔柱[3]。最新的地震成像发现欧亚大陆中北部下地幔中存在低速异常[4]，动力学模拟推测其可能是太平洋深部下地幔中大低速体的一个分支[5]。地幔柱的存在将导致局部高温异常，在转换带中引起电导率的增高。而兰州地区位于上述可能的地幔柱的中间位置，其深部电导率可能对这些地幔柱产生响应。因此，研究兰州地区地幔转换带的电导率模型，能够从电性结构角度间接给出地幔柱可能存在的相关信息。

地球深部，尤其是在410～1 600 km深度的岩石，其电导率主要依靠地磁测深C-响应的变化进行有效表达[6]。地磁测深C-响应的估计使用地磁台站记录的水平磁场和垂直磁场数据[7-8]，它们受到地核中的导电对流流体以及地幔、地壳、电离层、磁层和洋流体系的影响而发生磁场强度随时间的改变[9-10]。因此，通过研究地表磁场的变化，可以获得地球内部电导率，进而研究地球内部与岩石电性特征密切相关的参数和地球内部动力学特征[11-12]。

Kelbert等采用非线性共轭梯度法[13]，反演获得地幔的三维电性结构[14]。Kuvshinov等采用有限内存拟牛顿法进行地磁测深反演[15]，发现欧洲、非洲地区地幔转换带具有高阻特征，而中国东北地区地幔转换带具有低阻特征[16]。Munch等采用随机优化和模型探索技术进行地磁测深一维反演，结合温度、压力和含水量实验数据，获得欧洲、非洲等地区地幔转换带的温度信息[17]。张艳辉等基于全局光滑约束技术将得到的地磁数据进行了光滑模型一维反演，获得中国东部地区地幔转换带的含水量分布[18]。这些不同反演方法获得的结果既有相似性，也有互异性。

结果上的差异有可能是反演方法或反演时初始模型不同造成的，而非启发式寻优算法应用在地球物理反演工作中能很好地摆脱初始模型对反演结果的影响[19-21]，另外其寻得的最优解是在整个约束范围搜索，可有效避免解陷入局部最小的问题。因此，本文首次尝试将模拟退火(Simulated Annealing,SA)非启发式算法应用于地磁测深反演工作中，为研究地球深部结构提供新的参考依据[22-24]。

模拟退火全局优化算法目前已经被应用于地球物理反演中[25-27]。Ingber等针对模拟退火算法的效率问题，提出了非常快速模拟算法(Very Fast Simulated Annealing，VFSA)[28-29]。师学明等应用模拟退火算法进行一维大地电磁测深反演，获得了与地震剖面相一致的解释结果[30]。Sharma用模拟退火算法对印度东部甘加盆地进行一维直流电阻率测深反演，获得准确的地下含水层信息[31]。据此可以看到，模拟退火算法在电磁探测中都取得了较好的应用效果。

模拟退火算法反演一般用最后一次迭代的模型作为反演结果，而本文的不同之处在于对所有有效的随机反演模型进行统计分析，将统计模型作为地磁测深数据一维随机反演结果。为此，本文首先简要介绍地磁测深的基本原理；接着给出结合统计分析的随机反演技术；在此基础上，利用理论模型讨论影响反演结果的因素和效果；最后，对甘肃兰州台站实测地磁C-响应数据进行反演，并讨论该地区地幔转换带的性质。

1 地磁测深原理

地球磁层中的电流可以激发产生磁场，在地磁测深中常被称为外源场，其能在地球内部导电介质中引起感应场。因此，在地表观测到的感应场可携带地幔导电信息[32]。一般通过定义参数C-响应建立地表观测数据和地球内部电性结构的关系[33]。其表达式为

(1)

式中：tanθ为源空间的补偿项；a0为地球半径，取值6 370 km；在地磁测深研究中，采用球坐标系，原点定义在地心处，则任意点的磁场H均包含指向地心的分量(Hr)、水平北向分量(Hθ)和水平东向分量(Hφ)，其中，φ为经度，θ为纬度，r为指向地心的距离。

(2)

(3)

在球坐标系中，对式(3)进行变量分离，并取时谐因子eiωt，得到的磁场分量[35-36]为

(4)

(5)

而每层中Rn(r)与层导电性相关，满足

[n(n+1)-(iωμσ)2r2]Rn(r)=0

(6)

考虑层界面上磁场连续性，相邻的l-1和l层的r=rl边界上Rn(r)满足边界条件

Rn,l-1(rl)=Rn,l(rl)

(7)

(8)

2 反演算法

2.1 算法流程

(9)

式中：λ为正则化因子；Nω为频率个数；Cobs(ω)为观测数据计算得到C-响应；Cmod(ω)为模型理论C-响应；fm为模型约束项，用以描述模型粗糙度；σ为电导率模型向量，是需要反演的参数。

需要说明的是，组成地球的各薄球层厚度固定，不参与反演。目标函数最小化的过程与固体物质冷却退火达到能量最小稳定状态的过程相似。因此，模拟退火全局优化算法应该可以实现一维地磁测深反演。为此，将式(9)定义的目标函数作为退火过程中的能量函数，通过合理降低温度，当目标函数相应的能量达到最小时各个分子对应的状态就是反演问题的解。

本文提出的改进模拟退火(Improved Simulated Annealing，ISA)算法是在传统模拟退火算法反演的基础上，结合统计分析改进了反演的效果。其主要步骤包括：

步骤1，初始化T0、M、σmax、σmin、σ0、L、Rmax，j=0。

步骤2，置l=1，执行链内循环L次。

②当R

④更新l=l+1。

步骤3，降低温度(T)；j=j+1。

步骤4，当j>M时，退火结束；否则，返回步骤2。

步骤5，统计最优解。对拟合差小于Rmax的所有有效模型进行统计分析用以确定最终解。

在上述迭代过程中，T0为初始温度；j为当前迭代次数；M为迭代次数上限；l为链内当前循环次数，受马尔可夫链的链长(L)控制，链长越大则同一温度下全局搜索次数越多，更有机会得到全局最优解[37]，考虑时间成本，结合实验验证，本文采用L=10可以达到搜索效果；σmax和σmin分别为电导率参数取值上、下限，它们可以根据已知的岩石物理模型[38]或者前人经验[13]进行估计；σ0为模型初始电导率，因为模拟退火算法反演不需要固定的初始模型，所以本文采用有效范围内的随机数作为初始值；Rmax表示拟合差阈值，其中数据拟合差(R)的定义为

(10)

拟合差阈值的合理选择对反演至关重要。合理的阈值选择根据C-响应的数据质量、迭代次数等综合考虑。根本原则是要求该阈值能在退火过程中筛选出足够大且符合对数正态分布的样本量，否则需要在反演过程中调整拟合差阈值。

需要指出的是，在上述步骤中，传统模拟退火算法执行步骤1～4，而目标函数的最优解由最后一次退火结果确定。本文增加步骤5，使得反演的结果不单由最后一次退火结果确定，而是由所有符合约束条件解的统计期望值确定，从而更能反映随机反演的特点。

2.2 模型参数随机生成

随机电导率参数的产生基于温度的似柯西分布法，具体公式为

(11)

(12)

抽样和反演的结果受模型参数的范围约束。如果模型参数区间选择不合理，则反演结果集中在区间边界处(图1)[31]，在区间不合理时应适当扩大参数区间的上、下限。

2.3 Metropolis准则

当前被随机修改后的模型需要依据Metropolis准则去判断是否被接受。该准则需要计算新状态的接纳概率(P)，其表达式为

(13)

Δ=(σj+1)-(σj)

(14)

阴影部分为预置参数变化区间；虚线代表理论真实值(黑色或者红色)；实线代表反演结果图1 模型约束范围对反演结果的影响Fig.1 Influence of Constraint Range of Model on Inversion Results

2.4 温度降低机制

合理的退火温度降低机制可以高效控制迭代过程的收敛，从而获得近似最优解。反演过程中控制温度的公式为

T(j)=T0exp(-cj1/m)

(15)

式中：T为当前温度；c为常数。

过高的初始温度会使解在迭代初期收敛很慢，因此，选择过高的初始温度是没有必要的。Sharma在直流电测深反演实验中发现，c取值为1，T0取值为0.01，m取值为2.5，可以保证反演的效率，并且在反演过程中它们不依赖于模型参数[31]。本文通过实验验证了这组参数对地磁测深也较为适用。

2.5 解分布统计分析

由于电导率参数不能为负值且为随机数，所以迭代过程中所有的有效模型参数可以用对数正态分布来描述，且有lnσ服从期望值为μ、方差为ε2的正态分布。就某一层而言，反演最终的电导率等于该层有效电导率参数集合的期望值E(σ)，其表达式为

E(σ)=eμ+ε2/2

(16)

其中，ε作为标准差可以用来衡量结果的不确定性。当ε越小，表示估计值越可靠。距期望值一个标准差以内的范围所包含的数据占总体数据的68%，本文将此范围作为可靠区间。

估计最终解的可靠性也要通过检验样本数据是否符合对数正态分布加以验证，并采用对数正态概率分布P-P图来判断[39]。P-P图是根据变量的累积概率对应于自然对数正态理论分布累积概率所绘制的。其横坐标表示电导率样本累积概率密度；纵坐标为理论累积概率密度p(σ)，由概率密度函数计算得到。其表达式为

(17)

当被检测的数据符合对数正态分布时，在P-P图中样本累积概率和理论累积概率近似相等，从而在图中表现为45°的似直线分布。从图2可见，反演的电导率σ4样本分布基本符合对数正态分布规律。

图2 σ4样本分布直方图和对数正态概率分布P-P图Fig.2 Histogram and P-P Graph of Lognormal Probability Distribution of Sample σ4

3 结果分析与讨论

3.1 模型和理论C-响应

为验证模拟退火算法应用于地磁测深一维反演中的可行性及其准确性，对合成模型进行反演测试并对结果进行讨论。理论模型(表1)以Kelbert等的一维正演层状模型为基础[13]，其C-响应(图3)作为观测数据进行反演。反演过程中模型变化范围根据Püthe等应用蒙特卡罗抽样方法获得的400个独立一维全球平均模型的最大变化区间确定[21]。

表1 理论模型电导率参数Tab.1 Electric Parameters of Theoretical Model

图3 理论模型的C-响应Fig.3 C-response of Theoretical Model

图(b)中层电导率是反演结束时的模型参数图4 理论模型的地磁响应模拟退火算法反演结果Fig.4 Results of SA Algorithm Inversion of Geomagnetic Response of Theoretical Model

图5 反演结果稳定性对比Fig.5 Stability Comparisons of Inversion Result

3.2 模拟响应数据反演

首先，对图3中的纯理论数据进行反演测试。模拟退火时，迭代次数上限设置为20 000，链长为1，其他参数同Sharma在直流电测深反演实验中的参数设置[31]。图4给出反演过程中拟合差和模型电导率变化及迭代结束时各层电导率。由图4可以看出，拟合差在反演过程中整体为降低趋势。由于理论C-响应不含噪声且足够精确，所以根据拟合差计算公式，其值可以降到接近0。当迭代次数超过3 000后，各个层的电导率变化集中在理论真实值的附近；当迭代次数大于10 000，反演得到的各层电导率与真实值相差极小，但其中依然存在经过Metropolis准则的判断，较差模型被接受，而拟合差升高的情况出现。

由于噪声在实际地磁数据中普遍存在[18]，一般认为，受过大噪声影响的C-响应不能用于反演研究，所以本文对图3中理论C-响应添加可接受的8%高斯噪声进行反演方法测试。反演采用与理论数据反演相同的初始条件，分别进行了3次独立模拟退火算法反演，反演结束时的电导率模型如图5所示。从图5(a)可见，在反演结束时，3次反演数据拟合效果基本一样；但从图5(b)可见，受噪声影响，反演结束时，传统模拟退火算法获得的最终解(3条蓝色虚线)存在较大差异。例如，σ3分别为0.160、0.168和0.367 S·m-1，都与0.200 S·m-1的理论电导率差异较大。这表明反演结果除了受方法本身的模型灵敏度差异影响外，即使采用全局反演方法，受噪声等因素的影响，依然存在明显多解性，不容易确定最合理的模型。

3.3 改进模拟退火算法反演结果

根据图5可知，传统模拟退火算法反演只取最后一次迭代值作为结果，显然忽略电导率参数在迭代过程中变化的分布特征，因此，如果结合统计分析方法可能会更准确地计算模型电导率的估计值。本文采用上述解分布特征分析的方法后，再对3次反演结果进行统计分析，得到的结果模型如图5(b)中3条红色实线所示。图5(b)中，红色虚线表示估计值上、下限，可以作为可靠区间。

对比图5(a)、(b)可知，虽然传统模拟退火算法反演结果相差较大，但改进模拟退火算法反演结果不受随机性影响，统计估计的解基本一样，反映了解收敛的方向(图中红色实线基本重叠)。因此，改进模拟退火算法反演具有更强的稳定性和准确性，受多解性影响小。3次反演过程中保留的数据集不同，但改进模拟退火算法反演得到3个最终模型几乎一致，且与理论模型高度近似，如σ3分别为0.214、0.218和0.220 S·m-1，接近0.200 S·m-1的理论电导率。虽然估计值与理论值有一定的差异，但每层电导率真实值都在可靠区间内。需要指出的是，经过实验，更多次的改进模拟退火算法反演结果依旧保持高度近似。因此，在实际反演中，2、3次改进模拟退火算法反演就可以实现电导率的可靠估计。

3.4 抗噪能力测试

为了进一步考察改进模拟退火算法反演技术抗噪能力，对图3中理论C-响应分别加上5%、10%、15%的高斯噪声后进行单次模拟退火算法反演，之后再进行统计分析。除拟合差阈值外(因为该参数的选择需要考虑噪声水平)，反演的其他参数都与图5一致。随着噪声影响的增加，拟合差阈值受噪声影响分别设置为：R5%,max=0.070、R10%,max=0.145、R15%,max=0.190。图6展示了不同噪声影响下C-响应反演拟合情况。最终反演得到的拟合差依次为0.068、0.141、0.187，指示噪声越大，拟合差越大。但从图6可见，不同噪声影响都得到与理论数据较好的拟合曲线。

图6 不同噪声影响下的C-响应反演拟合Fig.6 Fitting Diagrams of C-response Inversion at Different Noises

图(b)中不同颜色对应的虚线表示可靠区间图7 不同噪声影响下的理论模型C-响应反演结果Fig.7 Inversion Results of C-response of Theoretical Model at Different Noises

图7(a)为4个样本在3种噪声影响下反演结果的P-P图。从图7(a)可见，反演得到的有效模型参数基本符合对数正态分布，且在当前的噪声水平条件下，基本不受噪声的影响。图7(b)给出了最终电导率模型和可靠区间。从图7(b)可以看到，虽然受到不同噪声影响，但真实值都在可靠区间内。另外，反演电导率模型在410 km以浅和900 km以深基本重合。一方面可能是这些区间处于地磁测深的弱敏感区域，另一方面也可能反映本文方法压制噪声影响的能力。在地幔转换带深度附近，与理论模型比较，反演结果较好地反映理论模型参数。值得指出的是，根据噪声水平对反演结果的影响，噪声的存在导致反演的电导率稍微偏小，且存在噪声越大，反演的导电性越差，可靠区间范围越大。因此，当数据存在小于15%的噪声时，可以认为本文的反演方法依然有效，能保证该方法对受噪声影响的实测C-响应数据反演的准确性。

3.5 兰州地区地磁测深数据反演

本次研究将改进模拟退火算法应用在兰州台站的实测地磁测深数据的一维反演中。兰州台站地处内陆，受海洋效应小[40]，因此，直接利用基于BIRRP方法[41]处理得到的C-响应[图8(a)]进行反演。反演时根据深部地幔导电性随着主要矿物相变而跃变的规律，将地核以上划分为9层[13]。各层电导率区间的取值依据Semenov等反演结果[16]估计。模拟退火算法反演控制参数同图5。参考图6，估计图8(a)的噪声水平约为10%，设置Rmax为0.145。从图8(a)可以看出，最终反演得到与原始数据相同的变化趋势且光滑的C-响应曲线，计算得到的最终模型对应的拟合差为0.135。

图8(b)显示模型的各个电导率概率分布P-P曲线呈对角线，表明反演中拟合差小于阈值的模型符合对数正态分布特征，从而可以计算得到可靠的近似最优解。图8(c)给出了本文估计的最终近似最优解。在0～250 km深度内，兰州台站下方0～250 km电导率高于全球平均值[13]，且可靠区间与张艳辉等研究结果的置信区间[18]有较好的重合。在地磁测深的敏感深度范围内(410～1 200 km)，反演得到了与Semenov等相似的兰州台站下方电导率结果[16,18]。在410 km的不连续面上，没有明显的电导率跃迁，这与Yoshino等得到的大陆转换带导电性特征[42]一致；670～900 km深度的电导率略高于中国地区平均水平，与李世文的研究结论[43]相一致。需要指出的是：地幔转换带上层，本文反演得到的电导率大致为0.074 S·m-1，接近全球平均电导率；但转换带下层电导率为0.147 S·m-1，相对前人的结果，本文反演得到的电导率偏低。这个特点和欧洲中部地区转换带具有相同特征[13]。另外，虽然反演结果在转换带附近偏低，但是其可靠区间很好地包含了Semenov等在该深度的反演电导率[16,18]；且从对理论数据的噪声测试来看，存在噪声时，反演得到的电导率会有偏小的可能，反演结果在图7(b)中显示会左偏。因此，可以认为本文的反演结果是合理的，而实际导电性可能更大，但不会超过全球平均模型。

为了讨论兰州地区地幔转换带下层低的导电性，本文采用Yoshino等提出的高温高压导电模型[42]对转换带下层做岩石物理计算。模型表达式为

(18)

式中：σV为体积电导率；Cw为含水量；α为几何因子；k为玻尔兹曼常数；Ttr为热力学温度；σOH、σOP为指数前因子；HH、HP分别为电子跳跃传导和质子传导的活化焓。

Yoshino等的模型参数[42]如表2所示，其中，k=1.38×10-23J·K-1。假设地幔转换带是绝热的，暂时不考虑温度异常，以全球平局进行分析，在410、520、670 km深度处对应的温度分别为1 780 K、1 850 K和1 925 K[42]。

表2 岩石物理模型参数

图(b)中电导率是最终模型参数，对应图(c)中的黑色实线。图(c)中黑色实线代表改进模拟退火算法反演结果；黑色虚线代表可靠区间；黄色实线从左至右分别代表含水量(质量分数)为0%、0.1%、0.5%、1.0%的电导率剖面图8 兰州台站地磁测深反演结果Fig.8 Inversion Results of Geomagnetic Depth Sounding at Lanzhou Station

图8(c)中黄色实线为计算的转换带电导率剖面。从图8(c)可以发现，即使采用全球平均温度，为了重现反演结果，在转换带上层也需要要求0.5%的含水量，这与全球观测结果一致。而在转换带下层，由于反演的导电性接近全球平均电导率的下限，所以要求该层基本不含水。导电性同时取决于温度和含水量，因此，如果适当降低温度，该地区的转换带下层也可以适当含水。由于兰州地区转换带具体的温度参数难以获取，不排除兰州地区转换带湿冷的可能性。在这种情况下，兰州地区地磁观测不能提供地幔柱的证据。但需要指出的是，由于C-响应噪声导致转换带下层导电性可能偏低，保守估计，其最多可以接近全球平均电导率，显然重现这个电导率也不需要用高温进行解释。

4 结 语

(1)基于模拟退火算法反演，结合统计分析，提出了地磁测深一维反演的新思路。理论数据反演测试结果表明，结合统计分析的模拟退火算法反演得到的结果稳定，且具有良好的抗噪能力，但较大噪声影响会使得反演的导电率偏低。

(2)对兰州台站的实测数据进行反演得到的地幔转换带导电模型，其上层与前人反演结果相一致,但转换带下层反演得到的电导率偏低，推测含水量非常低，接近干枯状态。兰州地区的地磁测深不容易给出欧亚大陆中部存在地幔柱的证据。

谨以此文庆祝母校七十周年华诞，祝母校桃李遍天下，绿野追唐裴！1989年金秋，我从安徽农村带着木制的行李箱，坐汽车倒火车来到西安地质学院，图书馆大厅里的新生报道开启了我四年的求学生活。图书馆厚大的枣红色阅读书桌，物探楼阶梯教室明亮的灯光，422教室老师们时而慷慨激昂、时而循循善诱的课堂，混合着蔬菜和肥肉香味的特有钢丝面食堂，还有89211班情同手足的兄弟姐妹，是我人生中最宝贵的回忆。母校是灯塔，指引着我们前行，更是我心中的牵挂!