地球物理学报  2013, Vol. 56 Issue (10): 3596-3606   PDF    
大地电磁资料精细处理和二维反演解释技术研究(三)--构建二维反演初始模型的印模法
叶涛1 , 陈小斌1 , 严良俊2     
1. 中国地震局地质研究所, 北京 100029;
2. 长江大学油气资源与勘探技术教育部重点实验室, 湖北荆州 434023
摘要: 在主流的线性最优化大地电磁二维反演中, 如何合理构建初始模型是一个亟待解决的问题.常用的是采用均匀半空间或一维反演结果构建初始模型, 不易获得稳定可靠的反演效果.实践表明, 尽管基于不同初始模型的大地电磁二维反演结果差别较大, 但均较初始模型更为接近真实模型.基于这样一种认识, 经过反复的理论和实践探索, 我们提出构建大地电磁二维反演初始模型的印模法.印模法的基本思想是依据已有反演结果和均匀半空间模型之间的加权来确定下一步二维反演的初始模型, 它一方面保留了已有反演结果中关于真实模型的宏观轮廓信息, 另一方面, 保证了深部电性结构的均匀性, 从而满足大地电磁二维正演所要求的底边界条件.基于印模法, 本文进一步提出了迭代重构的反演思想.通过多个理论模型和实测数据的反演计算, 验证了上述方法可在很大程度上压制初始模型对反演结果的影响.
关键词: 大地电磁      初始模型      二维反演      印模法      迭代重构     
Refined techniques for data processing and two-dimensional inversion in magnetotelluric (Ⅲ):using the Impressing Method to construct starting model of 2D magnetotelluric inversion
YE Tao1, CHEN Xiao-Bin1, YAN Liang-Jun2     
1. Institute of Geology, China Earthquake Administration, Beijing 100029, China;
2. Key Laboratory of Exploration Technologies for Oil and Gas Resources, MOE, Yangzi University, Jingzhou Hubei 434023, China
Abstract: In the Mainstream Magnetotelluric (MT)2D inversion currently based on the linear optimization technology, the construction of starting model has a great influence on the inversion result.How to construct reasonable starting model for the 2D MT inversion is still an urgent problem.A starting model of uniform half space or from 1D inversion result may not easily obtain reliable inversion effect.Theory and long-time MT inversion experience shows that despite the difference of 2D magnetotelluric inversion results based on different starting models, the results all are closer to the true model than the starting model and obtain macro-outline characteristics of the true model.But in some details, they may be biased.Based on the understanding of such a phenomenon, and after the repeated exploration in theory and practice, we propose the Impressing method for the starting model building in MT 2D inversion, of which the basic idea is to determine the next 2D inversion's starting model by weighting models between uniform half space and the existing inversion result.The starting model built by the Impressing method contains the real information of true model and meanwhile, ensures the uniformity of deep electric structure which meets magnetotelluric 2D forward bottom boundary condition requirements.Based on the Impressing method, this paper further put forward the idea of multi-iterative inversion.Through the 2D inversions of several theoretical models and field measurement MT data, we verify the good effect of the above-mentioned methods on suppressing the non-uniqueness of inversion results resulting from starting model construction and so improving the reliability of MT 2D inversion results.This may also provide a new idea for constructing the starting model of other geophysical inversion..
Key words: Magnetotelluric      Starting model      2D inversion      The Impressing method      Multi-reconstruction     
1 引言

大地电磁反演本身是一个非线性问题.各种非线性反演方法,如蒙特卡洛法、遗传算法、模拟退火法、人工神经网络法等,由于计算能力的制约,目前在大地电磁高维反演中还没有得到实际性的应用,而以高斯-牛顿法为基础的线性化技术在大地电磁高维反演实践中依然占据着统治地位.高斯-牛顿法又称为非线性最小二乘法,其首先需要将非线性问题线性化,采用高斯最小二乘法求解线性问题,然后通过逐次迭代来求解非线性最优化问题.在非线性问题线性化过程中,需要事先给定一个初始模型.理论上,该初始模型应位于真实解的邻域内才能获得稳定可靠的反演结果.由于真实解是未知的,因而初始模型很可能在其邻域之外,导致反演结果陷入局部极值.反演实践结果也表明,采用不同的初始模型进行反演,会得到差异很大的反演结果[1-3].由于大地电磁二维反演的模型参数成千上万,反演搜索的模型空间很大,局部极值可能很多,故初始模型对二维反演结果的影响很大.

一维反演结果是对真实模型的一阶近似,在一定程度上反映了观测数据所对应的真实模型特征.由于反演参数较少,故一维反演本身的非唯一性远远弱于二维反演,似乎更适合构建初始模型.然而,实践表明,直接由一维反演结果构建二维反演的初始模型,效果不理想,在深部会出现严重的挂面条现象,与人们通常所认为的地球深部结构状况不相符,同时也不满足二维模型正演所要求的底部边界条件.目前实践中主要还是采用均匀半空间模型作为大地电磁二维反演的初始模型.但采用不同的均匀半空间模型会得到差别较大的反演结果.因此,如何有效地构建大地电磁二维反演初始模型,仍然是一个亟待解决的棘手的问题.

基于大地电磁资料处理与解释可视化集成系统MT-Pioneer[4]批量反演及多任务结果分析管理功能,本文采用不同的初始模型对多个理论模型的正演响应进行了大量的二维反演计算,并在结果对比分析的基础上,提出利用已有反演结果和均匀半空间模型的加权组合来构建大地电磁二维反演初始模型的印模法;在此基础上,进一步提出迭代重构的反演思想.经过理论和实测数据的反演结果对比实践,表明这种方法比单纯的均匀半空间或者单纯的一维反演更为有效,为合理构建地球物理反演的初始模型提供了新的思路.

2 三维模型的合成数据

为使理论反演算例更接近于野外观测真实状况,我们设计了多个二维、三维模型,获得了相应的二维、三维正演响应,然后采用不同方法构建初始模型进行二维反演计算,对反演结果进行了对比分析.限于篇幅,本文主要介绍两个三维理论模型(图 1所示)反演结果的对比情况.

图 1 内嵌三维异常体的二维阶梯模型(引自文献[5], 有修改) (a)模型1, 二维阶梯模型内嵌一个10 Ωm的3D低阻异常体; (b)模型2, 二维阶梯状模型中内嵌2500 Ωm的3D高阻异常体, 三维异常体沿x方向的构造延伸长度为4 km, 其他参数如图所示. Fig. 1 Three-dimensional anomalies in the two-dimensional ladder model (from Ref.[5], modified) The two-dimensional ladder models embedded with (a) a 10 Ωm low-resistivity 3D anomaly body (Model 1) and (b) a 2500 Ωm high-resistivity 3D anomaly body (Model 2).The extended lengths of anomaly bodies along the x direction are 4 km, other parameters as shown in the picture.

在如图 1所示的两个模型中,二维阶梯构造背景下分别包含一个三维异常体:高阻异常体和低阻异常体.三维正演采用大地电磁三维交错网格正演算法[6],正演计算网格为56×61×51,频段范围为0.001~100Hz,共计40个频点数据.在本研究系列(Ⅱ)[7]中,对这两个模型的正反演情况进行了非常详细的讨论.证实了三维模型正演响应数据的可靠性,为本项研究提供了可靠的理论合成数据.

3 反演软件、数据、网格和参数

针对图 1所示的三维理论响应数据,我们采用了不同方法构建初始模型进行二维反演.二维反演和结果对比分析采用MT-Pionner软件[4],该软件集成了Bostick[8]、曲线对比法[9]、OCCAM[10]、自适应正则化算法[11](ARIA1D)、RouPlus[12]等多种一维反演算法和带地形二维非线性共轭梯度反演方法(NLCG)[13],拥有强大的数据可视化分析和管理功能.在对比研究二维初始模型构建的过程中,采用了Bostick、OCCAM、ARIA1D多种一维反演算法,其中ARIA1D算法效果最好.

当前主流的二维反演算法有RRI[14]、OCCAM[15]、REBOCC[16]、NLCG[13]等.RRI反演速度很快,但效果一般;OCCAM反演过程稳定,结果可信度高,但速度很慢;REBOCC是OCCAM的变种,提高了OCCAM反演速度的同时,降低了其可靠性;NLCG既有较快的反演速度,又有较好的稳定性和可靠性,是当前世界上应用最多的反演算法.对于图 1所示的模型,在本研究系列(Ⅱ)[7]中做过NLCG和OCCAM反演结果的比较,NLCG反演结果优于OCCAM法.故本文二维反演采用NLCG算法.

根据已有研究结果[7],本项研究中将采用TM极化模式数据进行二维反演计算.计算选取的用于二维反演研究的测线剖面通过三维异常体的中心位置,剖面展布为y方向.为使数据更加接近于实测数据,在三维理论模型的正演响应中加入最大误差为10%的均匀随机误差.已有研究表明误差的概率分布对反演结果影响不大[17].在NLCG二维反演的参数设置中,数据门限误差为2%.

二维反演网格对于反演结果具有一定的影响.本项研究中采用矩形网格实现像素式反演成像.在MT-Pioneer中集成了测点中心网格自动生成算法[18],可自动优化生成横向网格,同时也可很方便地构建纵向网格.二维反演网格的横向网格为82个,其中不包括模型两侧自动扩展的横向边界网格,每两个测点之间至少插入一个节点,最后自动生成测点中心网格,使得横向网格既比较光滑,又满足测点基本位于网格单元中心的要求;纵向网格为114个,从地表开始往下递增,地表纵向单元尺度为6m,分几个节段设置不同的递增比例因子:<1 km段,递增比例因子为1.2;1~15 km段,递增比例因子为1.02;15~80 km段,递增比例因子为1.1;>80 km,递增比例因子为1.2.

正则化因子对于反演结果有较大的影响.为了提高反演结果的可对比性,在使用传统的方法构建初始模型时,选取1~100000之间共17个正则化因子进行二维反演计算,然后对正则化因子进行L曲线分析,选定L曲线转折点附近的正则化因子值,使得反演结果既能拟合反演数据,又能足够光滑.对图 1所示的响应数据,各种情况下L曲线的拐点位置较为一致,都在200~1000之间.经过对反演结果的对比分析后,最后确定正则化因子为700.

4 不同初始模型的反演结果对比:传统方法 4.1 初始模型为均匀半空间模型

图 2显示了初始模型分别为50 Ωm、200 Ωm、500 Ωm、1000 Ωm均匀半空间模型时的二维反演结果.可以看出,不同初始模型的大地电磁二维反演结果差别很大.相比较而言,初始模型为50 Ωm均匀半空间的反演结果与真实模型较为接近;而200 Ωm、500 Ωm、1000 Ωm均匀半空间的反演结果仅在浅部对真实模型有所反应,深部差异很大.此外,还可看到,均匀半空间电阻率越高,结果偏离越大.这主要是因为模型深部电阻率为10 Ωm,采用较低电阻率的均匀半空间模型,与真实模型更为接近一些.

图 2 不同均匀半空间构建初始模型的二维反演结果 上为三维高导异常体模型的二维反演结果, 下为高阻异常体模型的二维反演结果, 自左至右, 初始模型依次为50 Ωm、200 Ωm、500 Ωm、1000 Ωm的均匀半空间模型. Fig. 2 The 2D inversion results by using different uniform half space to construct starting model The upper inversion results are from Model 1 and the below results from Model 2, the starting model for 2D inversion are respectively 50 Ωm、200 Ωm、500 Ωm、1000 Ωm from left to right.
4.2 一维反演结果构建初始模型

由以上算例可知,采用均匀半空间法构建大地电磁二维反演的初始模型,当均匀半空间的电阻率值不合理时,反演结果不可靠.一维反演是真实模型的一阶近似,采用一维反演结果构建二维反演的初始模型可能比采用均匀半空间模型更为合理.

图 3所示为采用一维自适应正则化反演[11]结果构建初始模型的二维反演结果.从图中可以看出,二维反演结果与一维初始模型差别不大,表明二维反演本身对模型的修正能力是有限的.在模型一维性比较好的区域,一维反演很好地重构了真实模型的性质.但是,二维阶梯和三维异常体所在的区域,反演结果与真实模型的差异明显,尤其对于高导体模型,挂面条现象非常严重.许多其他算例实践也表明,大地电磁二维NLCG反演中,初始模型中阶跃式突变的电性特征在反演结果中一般都会得以保留,尤其是较深或者其他数据约束较弱的部位.由于视电阻率曲线的静位移效应广泛存在,在对实测数据的一维反演结果中,挂面条现象也就非常普遍,而这种直上直下的电性结构特征与实际的大地构造是不相符合的,而且也不满足二维正演所要求的深部为均匀半空间模型的底边界条件.这是目前大地电磁实测资料的二维反演中,大多采用均匀半空间模型而少有采用一维反演结果构建初始模型的主要原因.

图 3 基于一维反演结果构建初始模型的大地电磁二维反演 (a)三维高导异常体模型的二维反演初始模型(ARIA[11]一维反演结果); (b)三维高导异常体模型的二维反演结果; (c)三维高阻异常体模型的二维反演初始模型(ARIA一维反演结果); (d)三维高阻异常体模型的二维反演结果. Fig. 3 The 2D inversion results based on 1D inversion result to construct starting model (a) The starting model of 2D inversion for Model 1 (1D inversion results by ARIA[11]); (b) 2D inversion result of Model 1;(c) The starting model of 2D inversion for Model 2 (1D inversion results by ARIA); (d) 2D inversion result of Model 2.
5 构建初始模型的新方法:印模法 5.1 寻找新方法的基本思路和尝试

从以上算例可以看出,以均匀半空间模型作为大地电磁二维反演的初始模型,能够满足二维正演的所有边界条件,但若其电阻率值不合理,反演结果与真实模型差别会很大;而采用一维反演结果构建二维反演的初始模型,可以在某些一维性较好的地方获得比较好的反演结果,但在横向电性差异较大的地方,会出现无法修正的挂面条现象.正如以上所分析的,一维反演结果中普遍存在的挂面条现象无法满足二维反演计算中的正演底边界条件,从而影响到二维反演结果的整体特征.是否能找到一种方法,能够综合上述两种构建初始模型方法的优点呢?

为了结合二者的优势,我们尝试了多种方法.已有算例表明,无论是采用均匀半空间还是一维反演结果构建初始模型,在浅部,二维反演结果与真实模型的差别较小.这种现象并非偶然.因为一般来说浅部数据约束比较好,同时也是数据一维性比较好的地方(对应于高频部分).因此,我们尝试的第一种方法是将一维反演结果在某个深部截断,该深度以上保留原有一维反演结果,而该深度以下为均匀半空间.这种简单拼接的方法克服了一维反演结果中的挂面条问题,使得二维正演的边界条件得以满足,但却引入了一条不存在的水平向的电性突变带.通过具体的反演实例发现,这条人为引入的水平突变带在二维反演中无法消除.图 4a图 4b给出了这一方法的一个算例.

图 4 结合一维反演结果和均匀半空间模型构建二维反演初始模型的两种尝试 (a)简单拼接法构建的初始模型; (b)基于简单拼接法初始模型的二维反演结果; (c)深部加权过渡法构建的初始模型(5 km以上浅部保持原有反演结果, 5 km以下深部, 为已有反演结果与100 Ωm均匀半空间的线性加权结果); (d)基于深部加权过渡法, 初始模型的二维反演结果. Fig. 4 Two attempts to build the starting model of two-dimensional inversion by combining with the one-dimensional inversion results and homogeneous half space model (a) The starting model constructed by pure simple combination with 1D model and half space; (b) The 2D inversion results of the starting model constructed by pure simple combination with 1D model and half space; (c) The starting model constructed by deep weighted transition method (in the starting model, the model above d=5 km is the original inversion model.Below the depth, the resistivity of the starting model is the linear weighting model by the original model and ρ=100 Ωm); (d) The 2D inversion results of the starting model constructed by deep weighted transition method.

我们尝试的第二种方法是,在某个深度以上,依然保留一维反演的结果,而在该深度以下,至模型底边界,逐渐过渡到均匀半空间模型.亦即在某个深度以下,初始模型的电阻率值为一维反演结果和底边界以下的某个均匀半空间的加权和,其权重值大小依据电阻率单元的深度位置确定.显然,这种深部加权过渡的方法可以消除水平突变,又基本满足底边界条件.但反演实例表明,虽然这种方法可以减淡一维反演中的挂面条现象,但无法完全消除.对比图 4c图 3a可以看到,深部加权过渡法对挂面条模型的修正能力很有限.

最后,我们在上述思路的启发下,找到了效果良好的构建大地电磁二维反演初始模型的新方法--印模法.

5.2 印模法的基本原理

在初始模型中彻底消除挂面条现象、完全满足二维正演底边界条件,是利用已有的一维反演结果构建二维反演初始模型的第一要务.因此,在新的方法中,将初始模型的某个深度以下设为均匀半空间,在该深度以上区域(浅部)采用加权求和的方式确定该处的电阻率值.由这种方法构建的初始模型,既在浅部保留了已有反演结果的主要轮廓,又在深部消除了挂面条现象,完全满足二维正演底边界条件.经过大量的算例我们发现,用这种方式构建的初始模型较之单纯由均匀半空间或者一维反演结果构建初始模型的方法更为合理.同时,在此方法基础上,可以引入更多的二维反演技巧,为大地电磁二维反演初始模型的构建提供了新的思路.

由这种方法获得的初始模型,其浅部保留了已有反演结果的主要几何轮廓,就像一枚印章盖在一张白纸上一样,因此,我们将其形象地称作“印模法”.为了叙述方便,我们将印模法中已有反演结果模型称为“元模型”,均匀半空间模型称为“辅模型”,用于区分均匀半空间和混合加权区边界的深度被称为“印模深度”,将印模深度以上保留原有反演结果轮廓的区域称为“印模区”.形象地看,所谓“印模法”就是“元模型”在“辅模型”上留下“印模”痕迹的过程.

可以有多种加权求和的方法,如可采用线性加权和、平方加权和、开方加权和等.在本项研究中,主要采用线性加权和的方式.如图 5,假设初始模型的印模深度为d,均匀半空间(辅模型)的电阻率值可以是已有反演结果(元模型)在印模深度d处的平均电阻率值(如不特殊说明,本研究均按此处理),也可以人为设置.印模区的电阻率值为元模型与辅模型电阻率值的线性加权和:

图 5 利用印模法构建大地电磁二维反演初始模型示意图 (a)元模型; (b)辅模型; (c)初始模型.图中h为模型中某计算点的深度, d为印模深度, ρhρd分别为元模型和辅模型的电阻率值, ρ为初始模型的电阻率. Fig. 5 The sketch map of constructing 2D inversion's starting model using the Impressing method (a), (b) and (c) stand for the primary model, the auxiliary model, and starting model, respectively.Here h is depth of a calculated point in starting model and d is impressing depth.ρh and ρd are the resistivity of the primary model and the auxiliary model, respectively.ρ is resistivity of the starting model.

式中ρh为元模型在某一横向位置下深度h处的值,ρd为辅模型的电阻率值,ρ为新构建的初始模型中某一点处的电阻率值.

5.3 反演算例:基于一维反演结果的印模法

现在我们以一维反演结果为元模型构建大地电磁二维反演的初始模型.通过图 3可以看到,一维反演结果确实包含了真实模型的许多信息.大致在8 km以上,除去挂面条部分,其他位置的反演结果看起来比较合理,反演数据的约束性比较强.因此,在用印模法构建初始模型时,取印模深度d为8 km,辅模型的电阻率值为元模型在印模深度处沿剖面的平均值,由此计算得到三维高导异常体的辅模型电阻率为10.7 Ωm,三维高阻异常体的辅模型电阻率为15.1 Ωm.

图 6展示了在一维反演结果基础上利用印模法构建的初始模型以及基于此初始模型的二维反演结果.由图可知,印模法构建的初始模型中,既保留了原一维反演结果(元模型)的主要真实轮廓,又在8 km(印模深度)处逐渐转变为均匀半空间模型(辅模型),完全消除了挂面条现象.对比图 3图 6可以看出,基于印模法的二维反演较之单纯基于一维反演结果的二维反演,模型向真实模型大大地靠近了一步,而且图 6中的反演结果对初始模型的修正是显著的,不像图 3中那样,二维反演结果与初始模型的差别很小,反演基本上局限于一维反演结果做初始模型所形成的局部极值的陷阱中.

图 6 在一维反演结果基础上利用印模法构建初始模型的大地电磁二维反演 (a)印模法构建的三维高导异常体模型的二维反演初始模型(d=8 km); (b)三维高导异常体模型的二维反演结果; (c)印模法构建的三维高阻异常体模型的二维反演初始模型(d=8 km); (d)三维高阻异常体模型的二维反演结果. Fig. 6 MT-2D inversion based on the Impressing method in which 1D inversion models were used to construct the starting models (a) Starting model for Model 1 (d=8 km); (b) Corresponding 2D inversion result of starting model in (a); (c) Starting model for Model 2 (d=8 km); (d) Corresponding 2D inversion result of starting model in (c).
5.4 反演算例:基于二维反演结果的印模法

很显然,印模法的元模型绝不局限于一维反演结果,任何已有的反演结果都是可以的.图 7为三维高导异常体模型的二维反演结果及相应初始模型的对比图.在利用印模法构建初始模型时,以图 2所示的高导异常体的二维反演结果为元模型,亦即元模型分别由初始模型为50、200、500、1000 Ωm均匀半空间模型反演得到.对于不同的元模型,印模深度最终分别确定为10、7、6 km和5 km.辅模型的电阻率依据印模深度自动从初始印中求取.

图 7 三维高导异常体的二维反演结果的对比 上为反演初始模型, 采用基于二维反演结果的印模法构建; 下为二维反演结果.由左至右, 元模型的初始模型依次为50、200、500、1000 Ωm, 印模深度依次为10, 7, 6, 5 km. Fig. 7 Comparison of 2D inversion results of Model 1 The upper models are starting models constructed by the Impressing method, and the below are corresponding inversion results.From left to right, the starting model of the primary models are 50, 100, 500, 1000 Ωm, and the impressing depth d are 10, 7, 6, 5 km, respectively.

图 8为三维高阻异常体模型的二维反演结果及相应初始模型的对比图.在利用印模法构建初始模型时,以图 2所示的高阻异常体的二维反演结果为元模型,亦即元模型本身分别由初始模型为50 Ωm、200、500、1000 Ωm均匀半空间模型反演得到.基于不同元模型的印模深度分别为10、16.5、4.5 km和4 km.辅模型电阻率依据印模深度自动由元模型中求取.

图 8 三维高阻异常体的二维反演结果的对比 上为反演初始模型, 采用基于二维反演结果的印模法构建; 下为二维反演结果.由左至右, 元模型的初始模型依次为50、200、500、1000 Ωm, 印模深度依次为10、16.5、4.5、4 km. Fig. 8 Comparison of 2D inversion results of Model 2 The upper models are starting models constructed by the Impressing method, and the below are corresponding inversion results.From left to right, the starting model of the primary models are 50, 100, 500, 1000 Ωm, and the impressing depth are 10, 16.5, 4.5, 4 km, respectively

比较图 2图 7图 8可以看到,当采用印模法构建初始模型以后,无论是高导异常体模型还是高阻异常体模型,二维反演的结果都得到了明显的改善.其中,原始均匀半空间初始模型为50 Ωm的反演结果最佳,与真实模型吻合程度高.同时也要看到,200、500、1000 Ωm的反演结果与真实模型差异较大.这说明,单次印模法虽然可以使得反演结果有较大的改善,但其修正量也是有限的.

5.5 初始模型的迭代重构:基于印模法的多次反演

如上所述,采用印模法构建二维反演的初始模型,可以较大幅度改进反演结果的准确性,但如果原有初始模型与真实模型差距太大,印模法也无法使得反演结果立刻满足我们的解释要求.然而,一个很有意义的现象是,采用印模法以后的反演结果都有着明显的改进.如果我们在现有反演结果上继续采用印模法构建初始模型,做进一步的二维反演,如此迭代下去,能否得到满足要求的二维反演结果呢?这就是我们为克服初始模型的影响所提出的“迭代重构”反演的基本思想,为大地电磁二维反演提供了新的实现思路.

图 7图 8中所有的反演结果分别按照印模法迭代重构初始模型,进行了多次反演,都取得了非常理想、一致的反演结果.限于篇幅,将图 6图 7中反演结果最差的基于1000 Ωm的迭代重构反演过程展示于图 9图 10中.由图可知,从最原始的1000 Ωm均匀半空间初始模型开始,三维高导异常体模型总共经历了5次反演,1次初始反演,4次迭代重构;而三维高阻异常体模型则总共经历了6次反演,1次初始反演,5次重构迭代.其他情形的反演次数都相应地比1000 Ωm要少,而最终的反演结果间差别非常细微(见图 9图 10中的最终的结果),对结果的解释没有影响.这表明,采用基于印模法的迭代重构反演可以在很大程度上压制初始模型对反演结果影响.

图 9 高导异常体模型完整的迭代重构反演过程 上为反演初始模型, 采用基于二维反演结果的印模法构建; 下为反演结果. Fig. 9 The processing of multi-reconstructive inversion for Model 1 The upper are starting models constructed by the Impressing method where the primary model are the results of 2D inversion; the below are corresponding inversion results.
图 10 高阻异常体模型完整的迭代重构反演过程 上为反演初始模型, 采用基于二维反演结果的印模法构建; 下为反演结果. Fig. 10 The processing of multi-reconstructive inversion for Model 2 The upper are starting models constructed by the impressing method where the primary models are the results of 2D inversion; the below are corresponding inversion results.
5.6 关于印模深度d的确定

在印模法中,印模深度的确定非常重要.通过大量的算例分析研究,我们发现这是有一定规律可循的.

对比分析图 2图 3图 6-图 10可以发现,当元模型的电阻率值较小时,二维反演结果较好.这个现象不是偶然的.因为真实模型(图 1)在深部的电阻率值是10 Ωm,比较小.一般情况下,二维反演结果模型中,存在一个数据约束程度由强变弱再到无的区域,其中,数据约束程度由弱变无,在反演结果上表现为电阻率逐渐向初始模型值转变的过渡区.如果辅模型的电阻率值与真实模型的电阻率值差别很大,则过渡区的影响范围也会很大,可能影响到探测目标区,从而使得反演结果的准确性降低.

过渡区的存在为我们确定合理的印模深度提供了一个有效的手段.如图 11所示,如果元模型的深部电阻率值比较大,大于真实模型的基底电阻率值,则可以选择过渡区的极小值点作为印模深度,这样可以获得与真实模型的基底电阻率最接近的初始模电阻率值,并使得印模区电性结构的确定性较好.由图 11可知,元模型不同,印模深度存在差异.同样,如果元模型的深部电阻率值比较小,小于真实模型的基底电阻率值,则可以选择过渡区的极大值点作为印模深度.我们在前面的算例研究中实际上已经多次应用这一方法确定印模深度,经验表明,这种方法是比较有效的.

图 11 不同初始模型的反演结果的曲线对比 曲线上的圆圈为过渡区的极小值, 可取为印模深度. Fig. 11 The comparison of inversion results based on different starting models Circles on the curves stand for the minimum value of the transition area, which can be choose as the impression depth.

另一方面,印模深度对反演结果的影响也是有限的.图 12显示了三维高导体模型迭代重构Ⅲ的基础上(见图 9),采用不同印模深度进行重构IV反演的情况.由图可见,印模深度d在7~12 km之间时,反演结果的差别很小.只有当d差别较大时(5 km或20 km),反演结果存在较为明显的差别.而即使是这两个差异较大的结果,较图 9中重构Ⅲ的反演结果也更接近于真实模型,可以在此基础上再进行重构反演计算.

图 12 不同印模深度对反演结果的影响:高导异常体模型迭代重构IV的反演 Fig. 12 The impression depth′s influence on the inversion results:results of the multi-iterative inversion IV of the Model 1
6 实测资料的二维反演

为验证印模法对于实测资料的有效性,对巴西某地区大地电磁实测资料进行了二维反演解释.反演单独采用TM极化模式数据,数据门限误差为2%.利用测点中心网格自动生成法[18]构建二维反演网格,网格尺度为57×95.在0.1~10000范围内选择了16个不同正则化因子进行反演计算,利用L曲线分析方法确定折中的正则化因子为70.

直接采用不同均匀半空间模型或一维自适应正则化反演结果模型构建初始模型的二维反演结果如图 13a所示,可以看到不同初始模型的二维反演结果之间相差较大.在各反演基础上,利用印模法分别进行了3次迭代重构反演,得到了基本一致的二维反演结果模型,见图 13b所示.通过比较图 13a13b的结果可以看出,迭代重构反演结果反映了原二维反演结果中共有的电性结构特征,可以认为是可靠的.

图 13 实测资料的迭代重构反演结果与传统初始模型构建方式的二维反演结果对比 Fig. 13 Comparison of inversion results for field data by using different starting models
7 结论与讨论

如何构建大地电磁二维反演初始模型是迄今悬而未决的大地电磁二维反演过程中的问题之一.本文所给出的算例验证了传统的利用均匀半空间或者一维反演结果构建法都无法为我们提供合理的初始模型.由于长期以来未曾找到有效的突破点,关于这个问题的研究基本处于停顿状态,本文提出的印模法为这一问题的解决提供了新的思路.

印模法综合了均匀半空间法和一维反演结果法的优点.它一方面像均匀半空间方法那样,能够完全满足二维正演所需要的边界条件;另一方面,也类似于一维反演结果法,将数据约束信息带入到初始模型之中.印模法以一种模糊混杂的方式保留已有反演结果的宏观轮廓,去掉其可能因为过拟合而带来的局部结构细节,从而可以反复使用,步步递进,使得初始模型逐渐向真实模型逼近,提高了初始模型构建的合理性.

印模法中,元模型、辅模型、印模深度的选取问题需要仔细考虑.元模型可以是一维反演结果或者任何其他已有的二维剖面模型,因而比较容易确定.而辅模型的选取可以像本文那样,直接由元模型在印模深度处剖面电阻率的平均值确定,也可手动拟定一个均匀半空间的电阻率值,因而也相对简单.比较复杂的是印模深度的确定,需要具有一定的反演经验和技巧.从本文研究中可以看到,辅模型的电阻率值最好与真实模型的深部基底的值一致,但在实测数据的反演中,这个值的具体大小是未知的.本文已经给出了一种参考意见,即如果真实模型的深部基底电阻率相对于其上覆地层是低阻,那么可以根据过渡区的极小值点来确定印模深度;类似的,如果真实模型的深部基底相对于其上覆地层是高阻,那么可以根据深部过渡区的极大值来确定印模深度.由于横向不均匀性的存在,不同剖面位置下的极值深度和类型可能不同,给印模深度的确定造成一定的困难.本文在有关印模深度的不同取值对反演结果的影响方面作了一定程度的分析,发现印模深度可以在一个较大范围内变化,而反演结果变化不大.

此外,我们在多次使用迭代重构印模法反演过程中还发现,使用这一方法似乎可以减轻反演结果对正则化因子的依赖性.在对理论模型进行的所有反演计算中,正则化因子取值均为700.这个值是最开始在均匀半空间法中经过L曲线分析和反演结果的对比确定的,但在后续大量印模法及迭代重构反演中,都没有改变,均获得了理想的反演结果.为了确证这个认识,我们将正则化因子取值为30,对三维高阻异常体模型响应进行了反演实验,获得了较好的结果,但浅部与700的结果稍有不如.有关这个问题的认识可能还需要更多的反演实践的积累.

本文的算例已经证明,采用印模法及其基础上的迭代重构反演方法,可以将一个偏离真实模型很远的初始模型逐步修正到真实模型的附近,进而获得好的反演结果.很明显,印模法并不局限于大地电磁二维反演,也不局限于大地电磁法,对于需要构建初始模型的地球物理反演方法,都可以尝试应用.

致谢

中国地震局地质研究所的蔡军涛博士和长江大学地球物理与石油资源学院的谢兴兵老师分别为本文的研究提供了理论数据和实测数据,审稿专家对本文提出了宝贵的修改意见,在此谨致谢意!

参考文献
[1] 陈孝雄, 王友胜. 利用初始模型提高MT反演分辨率. 江汉石油科技 , 2006, 16(4): 19–22. Chen X X, Wang Y S. Constructing reasonable starting model to improve the resolution of MT inversion. Jianghan Petroleum Science and Technology (in Chinese) , 2006, 16(4): 19-22.
[2] 马晓冰. 青藏高原东部电性结构特征及地球动力学意义. 北京: 中国科学院地质与地球物理研究所, 2004 . Ma X B. The electrical conductivity structure of the Eastern Tibetan Plateau and its geodynamic implications (in Chinese). Beijing: Institute of Geology and Geophysics, Chinese Academy of Sciences, 2004 .
[3] Bai D H, Maxwell M J. Deep structure of the Longling-Ruili fault underneath Ruili Basin near the eastern Himalayan Syntax: insight from magnetotelluric imaging. Tectonophysics , 2003, 364: 135-146. DOI:10.1016/S0040-1951(03)00054-4
[4] 陈小斌, 赵国泽, 詹艳. MT资料处理与解释的Windows可视化集成系统. 石油地球物理勘探 , 2004, 39(Suppl.): 11–16. Chen X B, Zhao G Z, Zhan Y. A visual integrated windows system for MT data process and interpretation. Oil Geophysical Prospecting (in Chinese) , 2004, 39(Suppl.): 11-16.
[5] 蔡军涛. 大地电磁三维畸变特征及校正新技术的应用研究. 北京: 中国地震局地质研究所, 2009 . Cai J T. The feature of three-dimensional distortion in MT and research of new correction methods (in Chinese). Beijing: Institute of Geology, China Earthquake Administration, 2009 .
[6] Mackie R L, Madden T R, Wannamaker P E. Three-dimensional magnetotelluric modeling using difference equations-theory and comparisons to integral equation solutions. Geophysics , 1993, 58(2): 215-226. DOI:10.1190/1.1443407
[7] 蔡军涛, 陈小斌. 大地电磁资料精细处理和二维反演解释技术研究(二)-反演数据极化模式选择. 地球物理学报 , 2010, 53(11): 2703–2714. Cai J T, Chen X B. Refined techniques for data processing and two-dimensional inversion in magnetotelluric (Ⅱ):which data polarization mode should be used in 2D inversion. Chinese J. Geophys. (in Chinese) , 2010, 53(11): 2703-2714.
[8] Bostick F X Jr. A simple almost exact method of magnetotelluric analysis.//Ward S ed. Proceedings Workshop on electrical methods in geothermal exploration. U. S. Geol. Surv., 1977: 174-183.
[9] 徐世浙, 刘斌. 大地电磁一维连续介质反演的曲线对比法. 地球物理学报 , 1995, 38(5): 676–682. Xu S Z, Liu B. The curve comparsion method of MT inversion for one-dimensional continuous medium. Chinese J. Geophys. (in Chinese) , 1995, 38(5): 676-682.
[10] Constable S C, Parker R L, Constable C G. Occam's inversion: A practical algorithm for generating smooth models from electromagnetic sounding data. Geophysics , 1987, 52(3): 289-300. DOI:10.1190/1.1442303
[11] 陈小斌, 赵国泽, 汤吉, 等. 大地电磁自适应正则化反演算法. 地球物理学报 , 2005, 48(4): 937–946. Chen X B, Zhao G Z, Tang J, et al. An adaptive regularized inversion algorithm for magnetotelluric data. Chinese J. Geophys. (in Chinese) , 2005, 48(4): 937-946.
[12] Parker R L, Booker J R. Optimal one-dimensional inversion and bounding of magnetotelluric apparent resistivity and phase measurements. Phys. Earth Planet. Int. , 1996, 98(3-4): 269-282. DOI:10.1016/S0031-9201(96)03191-3
[13] Rodi W, Mackie R L. Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion. Geophysics , 2001, 66(1): 174-187. DOI:10.1190/1.1444893
[14] Smith J T, Booker J R. Rapid inversion of two-and three-dimensional magnetotelluric data. J Geophys. Res , 1991, 96(B3): 3905-3922. DOI:10.1029/90JB02416
[15] DeGroot-Hedlin C, Constable S C. Occam's inversion to generate smooth, two-dimensional models from electromagnetic sounding data. Geophysics , 1990, 55: 1613-1624. DOI:10.1190/1.1442813
[16] Siripunvaraporn W, Egbert G. An efficient data-subspace inversion method for 2-D magnetotelluric data. Geophysics , 2000, 65(3): 791-803. DOI:10.1190/1.1444778
[17] 李墩柱, 黄清华, 陈小斌. 误差对大地电磁测深反演的影响. 地球物理学报 , 2009, 52(1): 268–274. Li D Z, Huang Q H, Chen X B. Error effects on magnetotelluric inversion. Chinese J. Geophys. (in Chinese) , 2009, 52(1): 268-274.
[18] 陈小斌, 赵国泽. 自动构建大地电磁二维反演的测点中心网格. 地球物理学报 , 2009, 52(6): 1564–1572. Chen X B, Zhao G Z. Automatic construction of Site-centered Grid (SCG) for 2D magnetotelluric inversion. Chinese J. Geophys. (in Chinese) , 2009, 52(6): 1564-1572.