地球物理学报  2012, Vol. 55 Issue (2): 671-682   PDF    
基于遗传算法的大地电磁阻抗张量分解方法研究
尹曜田1 , 魏文博1,2 , 叶高峰1,2 , 金胜1,2 , 董浩1,2     
1. 中国地质大学地下信息探测技术与仪器教育部重点实验室, 北京 100083;
2. 中国地质大学地球物理与信息技术学院, 北京 100083;
3. 中国地质大学地质过程与矿产资源国家重点实验室, 北京 100083
摘要: 本文针对Groom-Bailey分解法在消除地表电性不均匀体对MT数据的畸变效应时存在的问题引入遗传算法,对基于传统线性最优化方法的GB分解算法进行改进,提出基于遗传算法的大地电磁阻抗张量分解方法.通过对理论合成数据以及三维/二维模型正演数据的分解试验,并且对青藏高原北缘阿尔金断裂地区实际MT数据进行分解处理,证明了基于遗传算法的GB分解能够更加有效地校正三维/二维情况下近地表三维电性不均匀体所造成的畸变影响.最后,在已有算法基础上,研究了基于遗传算法的多频率GB分解算法和MT数据静校正算法,并通过实际MT数据的处理证明了这些算法的有效性.
关键词: 大地电磁      阻抗张量      地质构造走向      GB分解      遗传算法      多点多频      静校正     
An improved GB decomposition method based on genetic algorithm
YIN Yao-Tian1, WEI Wen-Bo1,2, YE Gao-Feng1,2, JIN Sheng1,2, DONG Hao1,2     
1. Geodetection Laboratory of the Ministry of Education, China University of Geosciences, Beijing 100083, China;
2. School of Geophysics and Information Technology, China University of Geosciences, Beijing 100083, China;
3. State Key Laboratory of Geological Processes and Mineral Resources, China University of Geosciences, Beijing 100083, China
Abstract: The galvanic distortion analysis approach advocated by Groom and Bailey has become the most adopted method. As proposed by Groom and Bailey,one must determine the appropriate frequency-independent telluric distortion parameters and geoelectrical strike direction by fitting the seven-parameter model on a frequency-by-frequency and site-by-site basis independently. Although this approach has the attraction that one gains a more intimate understanding of the data set,it is rather easy to be tapped into a local minimum before reaching the global minimum by using a traditional optimization methods (e.g. Least square method or Gauss Newton method) when conducting the repetitive application. In this paper,we improved GB decomposition method by Genetic Algorithm firstly,then tested the reliability of this new GB decomposition method by decomposing the synthetic data and data generated by forward modeling. Lastly,we applied it into the studying of the structure strike of an area in Northern Qinghai Tibetan Plateau and obtained results consistent with the results of previous geological and geophysical work. What's more, multifrequency GB decomposition and MT static shift correction method have also been studied and their validity has been verified by practical MT data test..
Key words: Magnetotelluric sounding      Impedance tensor      Geoelectrical strike direction      GB decomposition method      Genetic algorithm      Multi-frequency      Static shift     
1 引言

在大地电磁(MT)资料分析中出现的相当一部分问题,都是因为近地表电性水平非均匀性使二维或三维的区域构造中水平均匀的正常电磁场发生畸变,从而对MT 阻抗张量产生影响[1-2].故研究MT阻抗张量畸变理论及分解方法对于反演解释顺利进行具有重要的意义.

传统的一种MT 阻抗张量分解方法叫做Swift分解或者Swift旋转[3],对具有明确走向的近二维构造模型,适应性最好.但当有近地表三维电性非均匀体产生的畸变效应存在时,用该方法得到的电性主轴方向往往与真实二维构造情况相差甚远[4].而这种三维/二维构造在实际MT 资料处理中是普遍存在的[5].

Groom 和Bailey假设表层局部三维异常体覆盖在二维区域构造之上,提出的处理近地表局部非均匀体引起的电场电流畸变的张量分解法(GB 分解)[6-7]是目前最广泛应用的一种阻抗张量分解方法.该方法实际上是求解由分解定义式推导出的非线性超定方程组.传统的GB 分解求解非线性方程组方法的最优化策略是非线性问题的线性化,容易导致迭代过程收敛过快陷入局部最小而不能得到准确的分解结果.

本文针对这一问题,本文改进GB 分解方法,采用遗传算法求解非线性方程组,通过对理论合成数据和三维/二维正演模型数据的分解试验,以及对青藏高原北缘阿尔金断裂地区实际MT 数据进行分解处理,验证了该方法的有效性.并且针对当前MT数据处理中的热点问题,对基于遗传算法的多频率GB分解以及用上述方法消除MT 数据静位移进行了研究.

2 基于遗传算法的GB分解 2.1 GB分解基本原理[7-8]

Groom 和Bailey假设表层局部三维异常体覆盖在二维区域构造之上,将实测阻抗张量Zm分解为:

(1)

式中C为与频率无关的局部畸变矩阵,C=gTSAg为常数,R为旋转矩阵,RTR的转置矩阵,T为扭变矩阵,S为剪切矩阵,A为分裂张量,Z2D为区域真实的二维构造阻抗张量.将它们展开分别为:

其中θ 为旋转角(走向角),t为扭转(twist)因子,e为剪切(shear)因子,s为分裂比例因子.

需要注意的是,gAZ2D 在数学上有非惟一确定的困难,gA的物理含义相当于发生静位移,它们与Z2D的乘积gAZ2DZ2D在形式上完全一样,无法区分,因此又叫静位移因子,故并入Z2D中,而改用Z2表示gAZ2D.这样Z就可写成:

(2)

2.2 遗传算法[8]

遗传算法(Genetic Algorithm, 简称GA),是由Holland依据“适者生存"原理而创立的具有“生存+检测"的迭代过程的搜索算法.遗传算法从一组随机产生的称为“种群(Population)"的初始解开始搜索过程.其操作对象是一群二进制串,称为“个体"或“染色体(Chromosome)",每个个体都对应问题的一个解,利用转移概率规则帮助指导搜索.以适应度(Fitness)为依据,采用适当的选择策略在当前群体中选择个体,进行“繁殖(Reproduction)",保留拟合好的成员而去掉拟合差的成员,通过“交换(Crossover)"和“变异(Mutation)"进行随机部分交换和改变串中的某位来产生“后代(Offspring)",它将被评价、选择和繁殖.重复上述过程直至模型群体的拟合函数值变得很小.由于适值高的染色体被选中的概率较高,这样经过若干代之后,算法收敛于最好的染色体,它很可能就是问题的最优解或次优解.

在MT 阻抗张量数据存在噪声和区域结构受到近地表三维电性非均匀体影响导致数据畸变严重时,GB分解的非线性方程组的病态得到了放大[9].传统的最优化方法的策略是非线性问题的线性化,主要利用目标函数的梯度信息,通过反复迭代,寻找最优解,当初始模型远离真实模型时,解容易陷入局部极值[10-11].在实际大地电磁法工作中实测数据的阻抗张量具有很大的不确定性,因此不可能利用准确的先验信息为其提供接近真实模型的初始模型.而遗传算法由于是平行随机的评价模型空间各部分的拟合函数并对它们进行比较,所以只要模型群体的成员数大小、交换和变异的概率选择合适,这种算法一般不会使解陷入拟合函数的局部最优解中[10],从而能有效地求解非线性最优化问题.

2.3 基于遗传算法的GB分解

为了更好的确定非线性最优化问题的目标函数,首先对式(2)按实部和虚部展开得到的非线性方程组,用类似地球物理反演的方法将各已知量和未知量表示为已知数据向量和待求模型参数向量.

已知数据空间中的8个分量分别为实测阻抗张量四个分量的实部与虚部,用向量表示为:

(3)

待求模型参数向量中的7个分量,表示为模型参数向量为

(4)

其中,arbr 分别为所恢复的较大主阻抗a的实部与较小主阻抗b的实部;ai, bi 分别为ab的虚部,θ为较大的视电阻率的方位角(或走向角),e为剪切度,t为扭转度.

根据式(2)式可以确定数据空间和模型空间之间的函数关系,将该函数关系表示为:

(5)

则根据计算值与观测值之差的平方和最小为原则,求模型参数mi,使目标函数

(6)

为最小.

用遗传算法进行迭代时需要对模型参数变化范围加以限定[12].根据实测数据可知,MT 阻抗张量元素值变化范围相当大,不同元素间可能相差几个数量级,故无法对每一模型参数加以固定的搜索范围.为解决这一问题,本文中作者首先通过Swift分解粗略的确定模型参数主阻抗各元素以及θ 的较可靠初始值,确定模型参数数量级,并对该初值乘以不同的权值以确定各模型参数变化范围上下限.

图 1是用GA 算法拟合GB 分解参数的主要步骤.值得指出的是,由于该非线性问题存在固有的多解性.如果没有充分的利用先验信息,迭代效果必然不好.因此,做好拟合前的准备工作十分重要.最好能够利用各种已有地质或地球物理资料和方法来确定参数的搜索范围然后进行迭代.这样在一定程度上减小迭代的多解性,更好的保证精度和准确性.

图 1 基于遗传算法的GB分解程序流程图 Fig. 1 Program flow chart of GB decomposition method based on GA
3 理论模型实验算例 3.1 理论公式合成数据的GB分解

该实验中所用的二维模型数据为AlanJones和Gary McNeice在北美中央平原所测得的一组典型的二维阻抗张量数据[13].通过GB 分解理论公式Zm=RTSZ2RT 人为的在该数据中加入畸变干扰,然后分别用不同方法对该畸变数据进行恢复,以验证它们之间对比效果(表 1).

表 1 理论公式合成数据的GB分解效果表 Table 1 Results of different decomposition methods for the data synthesized by theoretical formula

表 1表明,基于马奎特法和遗传算法的GB 分解均能较准确的恢复出区域二维阻抗分量、走向角度以及畸变参数.为了对比这两种GB 分解的稳定性,分别对同一合成数据分解1000次,观察重复试验的统计效果.图 2 为分解出的区域走向角的统计玫瑰图(仅绘制第一象限部分),可以看出,虽然马奎特法在宏观统计上也能反映出较真实的区域构造走向,但是其结果分布空间较分散,远没有遗传算法的结果稳定,甚至于相当一部分点分解失败.分析原因是因为收敛过早,迭代停止时目标函数已经达到精度要求,但陷入局部最优解而不能达到全局最优.大量的实验表明,遗传算法能够稳定并准确的恢复出各模型参数.并且,因为采取了遗传基因多点交换、逐步缩小参数搜索范围以及引入S形函数的加速收敛的措施,不仅很好的克服了传统非线性最优化方法收敛速度慢、效率低的缺点,而且能够使迭代过程对模型空间的搜索更快也更彻底,也提高了搜索结果的精度(表 2).

图 2 理论公式合成数据1000次重复GB分解试验恢复区域走向角效果对比统计玫瑰图 Fig. 2 Rose diagram of the regional strike results of1000 times repeated test of GB decomposition methods based on GA and Marquette approach for the data synthesizedby theoretical formula
表 2 理论公式合成数据GB分解重复试验的平均效率及精度对比 Table 2 Comparisons of the efficiency and precision between repeated test of GB decomposition methods based on GA and Marquette approach for the data synthesizedby theoretical formula
3.2 正演模型数据的GB分解实验

为了进一步验证基于遗传算法的GB 分解消除近地表三维电性非均匀体引起的畸变以恢复区域二维构造特征的效果,本文使用泰国Mahidol大学Weerachai Siripunvaraporn教授编写的大地电磁三维正反演程序WSINV3DMT(Siripunvaraporn et al.,2005)进行正演数值模拟试验.

具体做法为:首先,构造一个二维地质模型,利用三维正演程序得到一组相应的响应;然后,在此二维地质模型的基础上,在地表处增加若干个三维异常体,这样就构造了一个三维/二维地质模型,计算出相应的响应,接着用笔者设计的基于遗传算法的阻抗张量分解校正程序对正演结果进行校正以对比效果.本文所用正演模型具体参数如图 3所示,区域构造走向角为N30°W.

图 3 包含地表三维畸变体的二维构造模型示意图 (a)三维/二维模型示意图;(b)测点及地表三维电性不均匀体分布位置图. Fig. 3 2-D structural model for forward modeling with local near-surface 3-D electrical inhomogeneities (a)Forward model of 3-D/2-D structure; (b)Locations of sites and near-surface 3-D inhomogeneities.

图 4是正演模型中5 号测点分别用不包含(图 4a)和包含(图 4b)了地表三维不均匀体的二维构造模型进行正演所得到的响应曲线以及各频点阻抗张量元素极化图.图 4a中的电阻率曲线和阻抗张量元素极化都表现出明显的二维性,并且极化方向能很好的指示出二维构造走向,为-30°(设顺时针方向为正角度).图 4b中能看出,地表三维不均匀体无论是对视电阻率曲线还是对阻抗张量元素极化都产生了比较严重的畸变影响,表现出比较强烈的三维极化特征.极化方向与真实的构造走向相去甚远.如果此时用传统的Swift分解做构造走向分析必将得出错误的结果,所以必须用阻抗张量分解的方法对其进行校正.

图 4 三维/二维模型正演结果(号测点) (a)不包含地表三维不均匀体的二维构造模型;(b)受到地表 3D不均匀体畸变作用的2D构造模型. Fig. 4 Forward modeling results for the 3D/2D structural model (a) 2D structural model without local near-surface 3D electrical inhomogeneities; (b) 2D model with 3D electrical inhomogeneities.

从多个频点数据分解结果分析中可以看出(图 5),遗传算法GB分解能够更加稳定的分解出真实的构造走向及畸变参数,而基于普通最优化方法的GB分解却很不稳定,具有较大的误差.

图 5 不同最优化方法对正演数据GB分解结果对比(号点) Fig. 5 Comparison of the recovery effects for GB decomposition methods based on different optimization methods
4 实际资料分解

以青藏高原北部阿尔金断裂区域的实际MT测深资料为例,对基于遗传算法的GB 分解方法在MT 数据分析中的应用效果进行验证.

阿尔金构造系位于青藏块体与新疆块体间,为两个块体的分界构造,总体走向为北东向60°~70°[14-16].该地区4000 号测线横跨新疆维吾尔自治区若羌县和青海省芒崖镇之间的几条主要断裂上(图 6).查阅已有的地质以及地球物理资料[14-19]可知,该区域内各断裂带两侧岩性及电性构造区分明,在深部应该表现为较为明显的二维电性构造.但该区岩石圈纵向呈层状,具断块结构,水平特别是浅层水平电性非均匀质强烈[18].4410~4450号测点阻抗张量的维性参数和元素极化图表现出强烈的三维性,推断受到近地表复杂的电性三维不均匀性的影响,故用基于遗传算法的GB分解对其进行校正.分解结果如图 7所示,为了反映分解效果,我们把相邻测点的分解结果绘制于同一坐标系下.不难看出,相邻测点间恢复出的区域构造走向角度曲线形态和数值在某一频段上可以表现出很强的相关性(用图 7中灰色点线表示该趋势),这是因为相邻测点虽然在地理位置上可能相差几十甚至上百公里,但是却反映了相同区域的背景构造的维性及走向特征.

图 6 青藏高原北部阿尔金断裂系地区实际MT测点位置及工区内主要地质构造分布(崔军文等,《阿尔金断裂系》,1999)[14] Fig. 6 Sites locations of the practical MT data collected in Northern Qinghai Tibet Plateau andthe distribution of main regional structures in this area[14]
图 7 4000线各测点阻抗张量数据进行基于遗传算法GB分解的结果 (a)4410 点和 4420 点;(b) 4425、4430 和 4440 点;(c) 4445 和 4450 点. Fig. 7 The result of GB decomposition based on GA for the impedance tensors of 4000 line (a) Site 4410 and 4420; (b) Site 4425,4430 and 4440; (c) Site 4445 and 4450.

为了更好的验证基于遗传算法的GB 分解的实际效果,我们绘制了这些测点走向角度的统计玫瑰图(图 8b)直观的反映不同频段(即不同深度)处的区域构造走向,并与前人(鲁新便等,1995)在该区域所做MT 剖面测量结果(图 8a)[18]进行对比.将4000线测点投影到鲁新便等所做MT 剖面上,通过位置对比可知,4410 与4420 号测点主要位于Fv-4断裂附近,4425、4430 及4440 三个测点位于Fv-5断裂附近,4445 与4450 号测点位于Fv-6 断裂附近.根据鲁新便等的解释,Fv-5是一条断面南倾的岩石圈断裂.以此断裂为界,南、北两侧岩石圈电性结构明显不同,断裂以北地壳上部电阻率为n×102 Ωm, 而以南为n×103 Ωm.在1∶100万中国地质图上,该断裂呈北东走向,地表延伸数百公里,从其位置判断,相当于阿尔金北缘断裂.阿尔金北缘断裂是前人(崔军文等,2002)[17]提出的祁连地块与塔里木地块间5条北东东向断裂之一.该断裂总体走向N60°E,延伸较平直[17].而对4425、4430 和4440 三个测点的做统计玫瑰图所确定出的区域构造走向位于N45°E至N65°E 之间,基本上与前人推断出的构造走向相同.Fv-4断裂是鲁新便等推断出的一条较浅的隐伏断裂,其断面两侧电性结构差异较大,走向与Fv-5断裂基本一致,4410 与4420 点的构造走向统计玫瑰图也能较好的反映出这一信息.至于4445 与4450两点所反映出的区域构造走向角,却集中位于N10°E~N20°E 之间,分析其原因,可能是因为该断裂两侧岩石电性结构差别不大,区域电性构造表现出一定的一维性,这两点的维性参数与阻抗张量极化图也能验证这一假设.

图 8 北阿尔金地区4000线实采MT数据地质解释效果图 (a)前人(鲁新便等,1995年)[18]在该区域所做MT剖面电性及地质解释结果;(b)用基于遗传算法的 GB分解所恢复区域电性主轴统计玫瑰图及与实际区域构造对比. Fig. 8 Geological and geoelectrical interpretation of practical MT data of 4000 line in North Altyn Tagh region (a) Geological and geoelectrical rnterpretation of previous practical MT profile; (b) The rose diagrams of regional strike angles recovered by GB decomposition based on GA.

图 9为对4425 号测点周期为1~2000s频段数据进行基于遗传算法GB 分解结果.从中可以看出,不论是恢复的区域构造走向角度,还是三维畸变参数-剪切角和扭转角均保持平稳.该频段电性主轴方向角分布于40°~60°之间,剪切角分布于20°~30°之间,扭转角大多分布于-30°~-20°之间.这说明基于遗传算法的GB 分解在该频段作用效果稳定.相比之下,在实际处理过程中基于普通最优化方法的GB分解发挥很不稳定,分解结果具有很大的均方根误差.并且通过以上与前人在该区域工作的地质解释结果对比可以看出,基于遗传算法的GB分解恢复得到的电性主轴方位角与当地区域构造走向角度基本吻合.故认定基于遗传算法的GB 分解较好的恢复了真实的区域构造走向角度,具有良好的实际应用效果.

图 9 4425号测点1〜2000 s频段数据基于遗传算法GB分解结果 Fig. 9 The result of GB decomposition based onGA for site 4425 of 4000 line(1 〜2000 s)
5 基于遗传算法的单点多频率GB 分解及在MT 静校正中的应用 5.1 基于遗传算法的单点多频率GB分解方法原理及实现

Gary McNeice和AlanJones(2001)等的研究表明,多点多频的GB 分解无论是在精度、效率,还是在对静位移的校正方面相比单点单频的方法都具有较明显的优势[13].在Gary McNeice 和AlanJones多点多频的GB 分解算法中,对于S个测点,N个频率值的MT 阻抗张量数据,数据空间中一共有S×N×8个已知的阻抗张量分量,根据GB 分解计算公式,需要用该数据空间去拟合S×(N×4+2)+1 个未知量.待求未知量包括S×N个区域复阻抗分量(AB),S×2个畸变参数值(te)以及二维区域构造走向角(θ)[13].该算法适用的前提假设是,这些测点的区域构造背景相对单一,拥有共同的走向角度,并且对于同一测点,畸变参数te与频率相关.很显然,这种假设并不符合中国地区实际MT 工作实际情况.通过对阿尔金地区以及其他工区实际MT 数据的单点单频GB 分解结果可知,大部分地区构造背景复杂,不同测点间的构造情况及维性特征可能具有较大差异.但对于同一测点,不同频率间特别是相邻频率间的畸变参数(te)却相差不是很大.故笔者根据实际工作需求对上述多点多频算法进行改造.在实际计算过程中笔者发现,对单站点5个频率的数据进行GB 分解效果最佳,既能保证效率和精度,又能避免在最优化计算过程中模型空间过大而出现的迭代结果不稳定的现象.

与单点单频的GB 分解相似,多点多频GB 分解的计算过程也是令特定的目标函数最小化的过程.对于单点5频率GB分解,可以认为在这5个频率之间,畸变参数te及走向角度θ 与频率无关.这样待求的模型空间中未知量一共为23 个(AB各分量及teθ),而此时的数据空间中一共有40个已知阻抗张量分量,即用40 个已知数据去拟合23个模型参数.故如果实际地质条件满足该方法的前提假设,可以大大的提高GB分解的效率.此时的目标函数可写为

(7)

不同于单点单频的GB 分解,(7)式中的mi为一个长度为7的一维向量.通过遗传算法求解上述最优化问题即可得出最优解.

对上述阿尔金地区4425号测点进行单点多频的GB分解与单点单频算法的效果对比如图 10 所示.显然,基于遗传算法的单点多频GB 分解效果较好,与单点单频的算法相比,基本上能够得出准确的分解结果,其精度完全能够达到实际工作的要求.不仅如此,因为该算法拥有相对更少的未知量,故拟合过程效率更高.对4425号点所有80个频率阻抗张量进行GB 分解,用单点单频的方法所需时间为6.7469s, 而用单点多频的方法仅需2.4175s, 效率提高一倍以上.同时,基于普通最优化算法的单点多频GB分解却因为模型空间中拥有更多的未知量而更容易陷入局部最小得到了更加混乱的分解结果(图 10c).因此可得出结论,遗传算法是适合多点多频GB分解的一种最优化方法.特别是对于数据量巨大的实际工作,使用单点多频算法将节省大量的运算时间.

图 10 4425号测点单点多频GB分解结果 (a)单点单频遗传算法GB分解结果;(b)单点多频遗传算法 GB分解结果;(c)单点多频普通最优化方法GB分解结果. Fig. 10 The results of single-site, multifrequency GB decomposition for site 4425 (a) Single-site, single-frequency GB decomposition; (b) Singlesite»multi-frequency GB decomposition bashed on GA; (c) Single-site, multi-frequency GB decomposition ba^ed on Marquette approach.

尽管如此,我们也不能忽视一个问题,那就是如果实际地质情况不符合单点多频算法的前提假设,即相邻5个频点间构造参数差异特别大,再使用该算法将会人为的加入系统误差.因此,我们需要认识到,多点多频的GB 分解更多反映的是区域构造的一种综合统计化的效果.对于大地电磁资料的精细处理,采用单点单频的GB 分解将是一种更加稳妥的方法.故在本文中多采用单点单频的GB分解方法.

5.2 基于遗传算法的GB分解在MT静校正中的应用

上文中提到,常规GB分解公式中的Z2 因为吸收了增益常数g和各向异性矩阵A,而并不是真实的二维构造阻抗张量矩阵Z2D .一般假设gA与频率无关,即它们对阻抗张量元素的振幅和走向主轴会起到一定拉伸和压缩的作用,但并不影响走向角度;表现在视电阻率曲线上,就是GB 分解能够正确的恢复区域视电阻率曲线形态,但是视电阻率曲线由于这两个与频率无关的参数作用而会发生一定程度的静位移.故单纯用常规的阻抗张量分解无法完全消除MT 数据的静位移[6].

AlanG.Jones曾经提出过一种使用电性主轴旋转消除MT 静校正的方法[20],笔者使用基于遗传算法的GB分解对其进行改进.

众所周知,在MT 勘探中由于趋肤效应,高频的勘探范围较小,其阻抗张量理论上满足一维地电模型条件下的特征,即对于阻抗张量元素有:Zxx=Zyy=0,Zxy+Zyx=0,反映在视电阻率曲线上即TE 和TM 模式曲线的首枝重合.实际上这也是我们判断资料中有无静态位移干扰的一个标准.

但在实际情况下,由于地表三维电性不均匀体的影响,高频MT 阻抗张量往往表现出二维或者三维地电模型特征.故首先对其进行GB 分解,求出Z2.这里Z2=gAZ2D,吸收了增益常数g和各向异性矩阵A,而Z2D是真实的二维构造阻抗张量.对于高频的阻抗张量,Z2D事实上是Z1D,即静态校正的预期结果.设静校正实系数张量Cs=(gA)-1,则有:

(8)

其中,

于是有

(9)

将其展开,Z1D中的各要素应满足

(10)

将上式中的复数ZTEZTM 展成虚实分量,可获得6个等式,再加上条件|Cs|=1,可采用最优化方法求出Cs 的4个主对角元素.理论上讲,Cxy= Cyx=0,Cs 为主对角矩阵.假设gA与频率无关,即可用Cs对同一测点的所有频率进行校正以消除这两个畸变参数的影响,从而也就校正了MT 数据的静位移.

运用上述方法对河南省开封市兰考县的11535S10测点实测MT 数据进行校正.该测点正好位于聊城-兰考隐伏断裂上,为标准的二维电性结构.但是从该测点数据的维性参数和阻抗张量元素极化图的分析中可知,该点数据三维特性明显.特别是其高频数据可以表现出强烈的静态位移特征(图 11a).根据前人(向宏发、王学潮等,2000年)在该地区所做浅层地震勘探以及钻探工作的解释成果可知,该区域距地表 1000 m 以内的第四纪地层中断错层大量发育[21-22],这些断错层特别是垂向的浅层构造表现出强烈的近地表三维电性非均匀,使MT测量数据不仅受到电场型畸变的影响,而且会因为垂直界面上的电荷积累而发生磁场型畸变.因此,需要对其进行静态位移的校正.校正效果如图 11b所示,TE 和TM 模式曲线在首枝重合,符合MT 勘探的规律.此时再对校正后的阻抗张量进行电性主轴分析,得到了相当准确的构造走向角度,与当地已有地质资料高度吻合,从而证明上述静校正方法效果显著.

图 11 兰考县11535S10测点实测MT数据用基于遗传算法GB分解做静校正效果图 (a)静校正前的视电阻率及相位曲线;(b)静校正之后的视电阻率及相位曲线. Fig. 11 The effect of static shift correction by GB decomposition based on GA for the practical MT data of site 11535S10 in Lankao, Henan province (a) Apparent resistivity and phase curves before static shitt correction; (b) Apparent resistivity and phase curves after static shitt correction.
6 结论和建议

本文通过对理论合成数据以及三维/二维正演模型数据的分解试验,并且对实际MT 数据进行分解处理,证明了基于遗传算法的GB 分解能够更加有效的校正三维/二维情况下近地表三维电性不均匀体对MT 阻抗张量造成的畸变影响,可以较好地确定区域构造特征,如走向等.并且,笔者针对单点单频阻抗张量分解中存在的问题,结合当前MT 数据处理研究的热点,开发了基于遗传算法的多频阻抗张量分解以及MT 静校正算法,并取得了良好的实际应用效果.

但是,由于方法本身的性质所决定,GB 分解都有一些不可确定的畸变参数,如gA,这就导致完成GB分解后的大地电磁响应函数中仍然包含静位移的影响.虽然采用本文中的基于阻抗张量分解的MT 静校正算法对这些畸变参数加以校正,并能够满足实际工作需求.但是,我们同时也应该注意到本文方法所确定出的静位移矩阵只是各种畸变效应的综合效果,我们依然不能确定出各种畸变的准确值.因此,我们需要结合其他的一些数学或者物理的方法,从阻抗张量数据中获得更多的信息,或者尽可能多的利用已有的地质资料和其他地球物理方法的资料,综合考虑以期达到较理想的效果.

致谢

本文是在张乐天博士和谢成良博士的热切关注和精心指导下完成的,在此表示衷心的感谢.

参考文献
[1] 晋光文, 孔祥儒. 大地电磁阻抗张量的畸变与分解. 北京: 地震出版社, 2006 . Jin G W, Kong X R. The Distortion and Decomposition of Magnetotelluric Impedance Tensor (in Chinese). Beijing: Seismological Press, 2006 .
[2] 王书明. 表面局部三维大地电磁曲线畸变校正—MT畸变校正阻抗张量分解方法. 西北地震学报 , 1998, 20(4): 1–11. Wang S M. The correction of magnetotelluric curve distortion caused by surfical local three-dimension inhomogeneities: the impedance tensor decomposition technique for the correction of MT curves distortion. Northwestern Seismological Journal (in Chinese) , 1998, 20(4): 1-11.
[3] Swift C M. A magnetotelluric investigation of an electrical conductivity anomaly in the southwestern United States : Massachusetts Institute of Technology, 1967.
[4] 王立凤, 晋光文, 孙洁, 等. 一种简单的大地电磁阻抗张量畸变分解方法. 西北地震学报 , 2001, 23(2): 172–180. Wang L F, Jin G W, Sun J, et al. A simple decomposition method of distortion in magnetotelluric impedance tensor. Northwestern Seismological Journal (in Chinese) , 2001, 23(2): 172-180.
[5] 蔡军涛, 陈小斌, 赵国泽. 大地电磁资料精细处理和二维反演解释技术研究(一)—阻抗张量分解与构造维性分析. 地球物理学报 , 2010, 53(10): 2516–2526. Cai J T, Chen X B, Zhao G Z. Refined techniques for data processing and two dimensional inversion in magnetotelluricⅠ: Tensor decomposition and dimensionality analysis. Chinese J. Geophys. (in Chinese) , 2010, 53(10): 2516-2526. DOI:10.3969/j.issn.0001-5733.2010.10.025
[6] Groom R W, Bailey R C. Decomposition of magnetotelluric impedance tensors in the presence of local three-dimensional galvanic distortion. Journal of Geophysical Research , 1989, 94(B2): 1913-1925. DOI:10.1029/JB094iB02p01913
[7] Groom R W, Bailey R C. Analytic investigations of the effects of near-surface three-dimensional galvanic scatterers on MT tensor decompositions. Geophysics , 1991, 56(4): 496-518. DOI:10.1190/1.1443066
[8] Holland J H. Adaptation in Natural and Artificial Systems. Michigan: University of Michigan Press, 1975.
[9] 杨长福, 林长佑, 徐世浙. 大地电磁GB张量分解法的改进. 地球物理学报 , 2002, 45(增刊): 356–364. Yang C F, Lin C Y, Xu S Z. The improvement of the Groom-Baily's tensor decomposition in magnetotellurics. Chinese J. Geophys. (nn Chinese) (in Chinese) , 2002, 45(增刊): 356-364.
[10] 杨文采. 地球物理反演的遗传算法. 石油物探 , 1995, 34(1): 116–112. Yang W C. Genetic algorithm for the geophysical inversion. Geophysical Prospecting for Petroleum (in Chinese) , 1995, 34(1): 116-112.
[11] 师学明, 王家映. 地球物理资料非线性反演方法讲座四—遗传算法. 工程地球物理学报 , 2008, 5(2): 129–140. Shi X M, Wang J Y. Lecture on non-linear inverse methods in geophysics (4)—Genetic Algorithm Method. Chinese Journal of Engineering Geophysics (in Chinese) , 2008, 5(2): 129-140.
[12] 王小平, 曹立明. 遗传算法: 理论、应用及软件实现. 西安: 西安交通大学出版社, 2002 . Wang X P, Cao L M. Genetic Algorithms: Theory, Applications and Implementation of Software (in Chinese). Xi'an: Xi'an Jiaotong University Press, 2002 .
[13] McNeice G W, Jones A G. Multisite, multifrequency tensor decomposition of magnetotelluric data. Geophysics , 2001, 66(1): 158-173. DOI:10.1190/1.1444891
[14] 崔军文, 等. 阿尔金断裂系. 北京: 地质出版社, 1999 . Cui J W, et al. Altyn Tagh Fault System (in Chinese). Beijing: the Geological Publishing House, 1999 .
[15] 葛肖虹, 刘永江, 任收麦. 青藏高原隆升动力学与阿尔金断裂. 中国地质 , 2002, 29(4): 346–350. Ge X H, Liu Y J, Ren S M. Uplift dynamics of the Qinghai Tibet Plateau and Altun fault. Geology in China (in Chinese) , 2002, 29(4): 346-350.
[16] 许志琴, 杨经绥, 张建新, 等. 阿尔金断裂两侧构造单元的对比及岩石圈剪切机制. 地质学报 , 1999, 73(3): 193–205. Xu Z Q, Yang J S, Zhang J X, et al. A comparison between the tectonic units on the two sides of the Altun sinistral strike slip fault and the mechanism of lithospheric shearing. Acta Geologica Sinica (in Chinese) , 1999, 73(3): 193-205.
[17] 崔军文, 张晓卫, 李朋武. 阿尔金断裂: 几何学、性质和生长方式. 地球学报 , 2002, 23(6): 509–516. Cui J W, Zhang X W, Li P W. The Altun fault: its geometry, properties and growth pattern. Acta Geoscientia Sinica (in Chinese) , 2002, 23(6): 509-516.
[18] 鲁新便. 阿尔金山及邻近地区深部地电结构探讨. 新疆地质 , 1995, 13(3): 256–263. Lu X B. Deep-seated geoelectric structures in the Altun Mountains and its neighbouring area. Xinjiang Geology (in Chinese) , 1995, 13(3): 256-263.
[19] 鲁新便, 石彦, 田春来, 等. 塔里木盆地东南部及阿尔金地区深部地质构造分析. 新疆石油地质 , 1996, 17(4): 314–317. Lu X B, Shi Y, Tian C L, et al. Analysis of deep geological structures in southeastern Tarim basin and Altun area. Xinjiang Petroleum Geology (in Chinese) , 1996, 17(4): 314-317.
[20] Jones A G. Static shift of magnetotelluric data and its removal in a sedimentary basin environment. Geophysics , 1988, 53(7): 962-978.
[21] 向宏发, 王学潮, 郝书俭, 等. 聊城-兰考隐伏断裂的第四纪活动性——中国东部平原区一条重要的隐伏活动断裂. 中国地震 , 2000, 16(4): 307–315. Xiang H F, Wang X C, Hao S J, et al. Activity of Liaocheng-Lankao buried fault in Quaternary. Earthquake Research in China (in Chinese) , 2000, 16(4): 307-315.
[22] 向宏发, 王学潮, 虢顺民. 聊城一兰考隐伏断裂第四纪活动性的综合探测研究. 地震地质 , 2000, 22(4): 351–359. Xiang H F, Wang X C, Guo S M. Integrated survey and investigation on the Quaternary activity of the Liaocheng-Lankao buried fault. Seismology and Geology (in Chinese) , 2000, 22(4): 351-359.