2014, Vol. 57
1 引言
利用大地电磁阻抗张量分解技术不仅可以压制近地表的三维小异常体的影响、为二维反演提供更为可靠的区域阻抗张量数据，还可以获得随测点、频率变化的区域电性主轴方位，为分析研究电性构造的几何特征分布提供了基本的数据支持.这些特征使得阻抗张量分解技术成为大地电磁资料处理的核心技术，在国内外得到了广泛的应用(Fischer and Masero， 1994；赵国泽等，1996；晋光文等，2003；Toh等，2006；晋光文和孔祥儒，2006；Heise等，2007，2008；Weckmann等，2007；Stagpoole等，2009；詹艳等，2011；Ledo等，2011；Anahnah等，2011).
迄今，已出现了多种阻抗张量分解技术(Swift，1967；Bahr， 1988，1991；Groom and Bailey， 1989，1991；Caldwell等，2004；Bibby等，2005；蔡军涛等，2010).然而，要对阻抗张量分解结果进行有效分析，并不容易.长期以来，阻抗张量分解主要基于单频点数据进行，并一般采用电性主轴-频率散点图进行单测点分析.由于张量分解所能提供的参数较多，数据结构复杂，对单测点、单频点的阻抗张量分解结果采用确定性的方式进行分析，无法充分利用张量分解所能提供的数据信息，从而影响了这一技术的有效应用.由于受环境电磁噪声、构造维性等因素的影响，单频点的阻抗张量分解结果——包括电性主轴、畸变因子、维性参数、区域阻抗张量等数据，具有较大的不稳定性，经过张量分解获得的区域阻抗张量，数据质量上往往比分解前更差，最终难以在反演中使用.
通过多测点-多频点的统计技术进行分析，能够在纷繁复杂的数据集中提取出主要特征，因而是进行阻抗张量分解结果分析的有效工具.此外，也可以在多测点-多频点的统计分析的基础上，选择特征较为一致的测点和频点，利用最优化技术确定一致的主轴方位，从而获得同一主轴方位的畸变因子和区域阻抗.Jones和他的学生基于GB分解技术曾经提出了一个多测点-多频点分解算法(McNeice and Jones， 2001).
由于方法的差异、分解结果的不确定性、算法程序使用功能方面的限制等原因，采用传统的方法进行阻抗张量分解时，即便是经验丰富的大地电磁专家之间也会产生很大的分歧.最近，Jones课题组和Anahnah课题组就北非地区阿特拉斯山脉的岩石圈结构展开了一场非常激烈的争论.两个课题组在同一地区反演得到的岩石圈结构差异很大(Ledo等，2011；Anahnah等，2011).Jones课题组认为这种差异是由于二者所采用的电性主轴方位的不同造成的(Jones等，2012)：在地壳部分，二者的电性主轴相差30°，而岩石圈地幔中，相差达60°.Jones等认为Anahnah等的反演结果是建立在一个不正确的主轴方位上的.Farida等人在对数据进行重新处理和反演以后，但依然坚持自己原来的主轴方位的判定(Anahnah等，2012).这一争论表明两点：其一、电性主轴方位对于大地电磁二维反演结果影响很大；其二、利用阻抗张量分解技术确定电性主轴方位，不是一件容易的事情.
造成阻抗张量分解技术使用困难的主要原因是，我们将一个不确定的结果施以确定性的分析.无论是电性主轴-频率散点图，还是近年来风行的相位张量椭圆图，给出的都是单测点-单频点的信息.而在统计玫瑰图中，一般给出的是所有测点、所有频点的统计结果，忽视了数据集中内部结构的差异：电性主轴在不同剖面位置、不同频段可能存在的差异.此外，数据质量、构造维性的差异也可能对统计结果造成影响，这些因素目前尚未得到有效的考虑.针对上述存在的问题，本文提出了基于阻抗张量分解技术的多测点-多频点的统计成像分析方法.
2 统计成像技术的算法基础与软件实现 2.1 大地电磁阻抗张量分解的基本原理大地电磁阻抗张量分解结果中，除区域阻抗张量外，还涉及5个实参数：4元素的畸变矩阵和电性主轴.式(1)中，Z 为观测阻抗张量，Z r为主轴方位上的区域阻抗张量，R 为旋转矩阵，C 为局部畸变矩阵.业已证明(Groom and Bailey， 1989；Mcneice and Jones， 2001；Bibby等，2005)，式(1)中，静位移因子gx和gy为非确定性参数，是张量分解技术所不能确定的；畸变因子a、b，电性主轴方位角θ为确定性参数，可以利用张量分解技术求解.相应地，通过阻抗张量分解技术可以很好地确定区域阻抗张量的相位，而视电阻率曲线则无法准确获得，这是导致张量分解后视电阻率曲线质量下降的主要原因.
在当前大多数大地电磁阻抗张量分解算法中，有5个确定性参数，分别是畸变因子a、b，电性主轴方位角θ、2个阻抗元素的相位(分别对应于TE、TM极化模式)；另外存在4个不确定性参数：静位移因子gx、gy再加2个极化模式的视电阻率值.不过，静位移因子和视电阻率之间是不独立的，只要确定了静位移因子，视电阻率值也就确定了，因而，真正的不确定性参数是两个静位移因子gx、gy.大地电磁阻抗张量分解结果的统计分析主要针对5个确定性参数进行，同时也考察在一定先验约束条件下获得的不确定性参数的频率、剖面分布情况.此外，基于张量分解技术还可以计算一维偏离度S1d、二维偏离度S2d等维性参数，可以对测区构造进行维性分析.
2.2 电性主轴的统计成像新技术：频率分布云图和测点分布云图电性主轴方位角θ是阻抗张量分解所确定的最为关键的参数.关于电性主轴的分析是张量分解结果分析中最为重要的一步.为展示电性主轴的分布特征，存在多种图示方式，其中，统计玫瑰图的应用最为广泛.尤其是在进行多测点-多频点电性主轴数据分析时，统计玫瑰图是一种强有力的统计成像手段.
一般说来，在利用统计玫瑰图进行多测点-多频点分析时，为了了解数据集中的内部结构，将测点和频点分别分成多个组进行统计.但是，如何确定这些分组中各组测点范围、频率范围的合理取值呢？已有的方法都是按照一定的步长，等分测点和频率成组(詹艳等，2011).然而，无论是频率分组，还是测点分组，由于无法有效确定分组的合理性，预先确定的步长范围很可能横跨在两个或多个不同主轴方位的数据段上，因而很难获得好的统计结果.因此，这种做法容易破坏数据内部结构的一致性.
为了克服分组不合理的问题，作者提出了新的统计成像思路.将包含一定数目测点、所有频率范围的统计玫瑰图沿着频率方向展开，在一个较小步长的主轴方位区间和频率区间所组成的方格中，统计数据点频次，可以获得一种新的剖面式统计云图.这个统计云图精细地展示了一定数目的测点下，电性主轴沿着频率轴的统计分布特征，为我们分析了解电性主轴随频率的变化情况、合理确定频率分组范围提供了直观而有力的工具.为叙述方便，我们将这种统计成像称为“频率分布云图”.
同样，将包含所有测点、一定频率范围的统计玫瑰图沿着测点方向展开，在一个较小步长的测点序列区间和主轴方位区间所组成的方格中，统计数据点频次，可以获得类似的剖面式统计云图.该统计云图精细地展示了一定频率范围内，电性主轴沿着测点序列方向的统计分布特征，为我们分析了解电性主轴随测点序列的变化情况、合理确定测点分组范围提供了直观而有力的工具.我们将这种统计成像称为“测点分布云图”.
2.3 数据质量和构造维性的影响大量的大地电磁阻抗张量分解结果分析表明，数据质量和电性结构的维性分布对电性主轴的统计结果具有很大的影响.一般情况下，噪声过大的数据，其张量分解的结果是不正确的，在统计中应予以剔除.此外，一维性强、三维性强的数据，通过张量分解得到的电性主轴方位不稳定，并且易受噪声的影响.电性主轴应主要由“纯二维性”强的数据来确定.
对于噪声大的数据，在我们所编制的软件模块中，类似于反演数据的选择那样，通过可视化的方式，手动设定删除标记，使其不参与统计分析.为克服构造维性所带来的影响，本文根据传统的一维偏离度S1d、二维偏离度S2d，提出了二维有效因子e2d的定义：
从式(2)可以看出，当S1d小、即一维性强的时候，e2d也小(一般情况下，S1d≥S2d)；当S1d大、即一维性弱的时候，如果S1d-S2d小，则意味着S2d也比较大，即三维性比较强，此时e2d也小，如果此时S1d-S2d大，则意味着S2d比较小，是一维性和三维性都比较弱，即纯二维的情况，此时，由式(2)计算的e2d会比较大.亦即，二维有效因子描述了结构的“纯二维性”，其值高，表明该数据对应的结构二维性好，但一维性和三维性都不好；其值小，则表明纯“二维性”弱，一维性或者三维性好.
由于可以描述大地构造的纯二维性，二维有效因子e2d用于电性主轴的统计加权中，可以提高二维性好的数据的权重，而降低其他两种数据的权重，从而增加统计结果的可靠性.另一方面，断裂带、岩性突变带附近往往呈现为线性构造，“纯二维性”强，因此有可能利用e2d成像进行追踪分辨.
需要注意的是，一般情况下，e2d远小于1，故加权统计结果不仅不再是一个整数，而且其值较原始统计频次(整数)看起来要小很多.
2.4 关于电性主轴的90°模糊性众所周知，利用阻抗张量确定的电性主轴方位具有90°的模糊性，既可以是电性构造的走向，也可以是电性构造的倾向.由于张量分解所给出的只有一个数值结果，而这个值在另外三个象限里也是可能存在的，故常常使得随频率变化的最佳主轴方位曲线变得非常散乱，导致传统的逐测点、逐频点确定性分析无法获得规律性的结果.但对于统计分析而言，这种困难不再存在.既然这个数值有可能存在于4个象限里，那就对其在4个象限里都进行频次统计.由于所有的电性主轴都具有4个统计量，机会均等，因此分析结果在统计上没有差异.这一技术下的统计玫瑰图具有4象限分布特征，完整地描述了构造的走向、倾向相互垂直的几何状态.
2.5 共主轴方位的多测点-多频点阻抗张量分解由于二维模型假设的需要，一条剖面应只有一个构造走向.根据单频点自由分解的统计结果(统计玫瑰图、统计云图等)可以基本确定反演剖面的优势方位，可以将此优势方位作为电性剖面的构造主轴.这一结果是单纯从主轴方位分布的特征来看的，其是否符合阻抗数据分布的需要呢？利用共主轴的多测点-多频点阻抗张量分解结果进行对比分析，可以解决这个问题.
阻抗张量分解可视为一个简单的反演问题.对于多测点-多频点技术而言，由于数据量、参数量的增大，再加上固定参数的组合情况，要具体去实现类似McNeice和Jones(2001)的最优化算法具有一定的难度.在传统的Swift方法中，除了采用Swift 1967年给出的最优化公式以外，还可以采用枚举方位角度的方式来获得电性主轴.这一方法虽然笨拙，但容易实现.
与Swift方法有所不同，在本文所提出的利用枚举法实现多测点-多频点的阻抗张量分解时，需要事先进行单频点分解，得到每个测点的区域相位数据；然后，在0~90°范围内，按一定角度步长逐渐递增，进行固定主轴(主轴方位已知)的张量分解，将其得到的相位与先前自由分解得到的相位进行拟合分析，在计算残差均方根的过程中需要考虑相位90°模糊性、数据质量、构造维性等因素的影响，选择拟 合最好的那个方位角作为多频点-多测点的分解结果.
2.6 其他参数的剖面成像除电性主轴和区域阻抗张量以外，张量分解还可以获得畸变因子a、b，一维偏离度S1d、二维偏离度S2d、二维有效因子e2d以及估算静位移因子gx、gy等参数.这些参量由于不存在周期性反复等问题，其值自小到大，物理意义不存在二义性，因而可以直接在频率域沿剖面绘制剖面分布图.这些图件可以进一步展示沿剖面不同深度、不同位置的构造复杂程度、畸变情况，可为反演结果的地质解释提供独立的依据.由于篇幅所限，本文仅对S1d、S2d、e2d三个维性因子进行了讨论.
2.7 统计成像技术的可视化实现Jones和他的学生在其多测点-多频点分解算法的基础上，研发了基于GB分解(Groom and Bailey， 1989，1991)著名的Strike程序(McNeice and Jones， 2001). 但是，由FORTRAN语言编制、依赖DOS命令行进行数据交换的Strike程序使用起来极不方便.极为低效的工作效率限制了Strike程序的推广应用.
基于可视化程序设计技术，本文第一作者在上述算法的基础上，研发了一个完整的大地电磁阻抗张量分解模块，集成于MT-Pioneer软件平台(陈小斌等，2004)中.在该程序模块中，实现了基于Swift旋转(Swift，1967)、Bahr分解(Bahr， 1988，1991)、GB分解(Groom and Bailey， 1989，1991)、相位张量(CBB)(Caldwell等，2004)、共轭阻抗法(蔡军涛-陈小斌-赵国泽法，CCZ)(蔡军涛等，2010)、五参数优化1)等多种方法的多测点-多频点的阻抗张量分解及统计成像技术.利用这个可视化模块，可以非常方便地实现多种方法的张量分解及其结果的相互比较.本文所有的研究工作都是在这个模块上进行的.
1)陈小斌提出的最新张量分解算法，尚未发表. 3 算例研究：理论模型分析下面以两个理论模型为基础，讨论本文提出的阻抗张量分解的多测点-多频点统计成像技术的分析过程和效果.正如前述，张量分解技术有多种，限于篇幅问题，本文主要讨论我们在本研究系列(一)(蔡军涛等，2010)中所提出的共轭阻抗法(CCZ)的分解结果.
3.1 模型1：3D/2D叠加模型阻抗张量分解技术要克服的是近地表的局部小异常体的影响.其最初得以提出的物理基础是建立在一种被称为3D/2D叠加模型的基础之上的.所谓3D/2D叠加模型，指的是区域构造为二维、近地表的浅部叠加了三维小异常体的这样一种模型.得以广泛应用的Bahr分解、GB分解都是基于这样一种模型.后来发展的相位张量法(CBB)、共轭阻抗法(CCZ)等则不需要区域构造为二维的假设，因而具有更加广泛的适用性.
本文所采用的3D/2D叠加模型如图 1所示.一个简单的区域二维构造上叠加了4个三维小异常体模型.利用Mackie等(1993)研发的三维正演程序进行了三维正演计算，获得了该模型的三维响应，为我们研究近地表小异常体所产生的局部畸变效应、静位移效应及其校正技术提供了基本的数据支持.
 | 图 1 二维区域构造背景上叠加有4个小异常体模型(3D/2D叠加模型)：左为剖面图；右为平面图(部分) Fig. 1 A theoretical 3D/2D model with four small abnormities superimposed on the 2D regional structure |
图 2给出了Swift方法和共轭阻抗法(CCZ)得到的L17测线(位置见图 1)上所有测点、所有频点的电性主轴统计成像.从图中可以看出，由于地表小异常体的影响，Swift方法给出的统计结果较为混乱.从测点分布云图中可以看出，在小异常体及其附近区域，Swift方法给出的电性主轴很乱，无法确定正确的区域主轴方位，而在偏离小异常体一定距离后，可以给出正确的区域阻抗主轴方位.这表明，这种尺度的小异常体的影响范围是有限的.频率分布云图则显示Swift方法的统计结果不够集中.相对而言，共轭阻抗法的统计结果显著性更强，近地表三维小异常体的影响得到了有效压制.
 | 图 2 测线L17上所有测点、频点电性主轴统计成像，未施加维性加权 上：Swift法；下：共轭阻抗法(CCZ).由左至右依次为统计玫瑰图、测点分布云图、频率分布云图，虚线框内为三维小异常体影响区域. Fig. 2 Statistical images of geo-electrical strikes for L17 line with all sites-all frequencies data calculated without dimensional weighting The upper images are for Swift method， and the lowers are for CCZ method. From left to right，the graphs are respectively rose diagram，site-based cloud diagram and frequency-based cloud diagram. |
图 3给出了利用二维有效因子进行加权后的电性主轴统计成像结果.由于Swift方法给出的一维偏离度S1d、二维偏离度S2d受近地表的三维小异常体的影响，因此加权没有起到应有的作用，反而使统计结果变得更加凌乱.共轭阻抗法(CCZ)分解在进行加权统计以后，主轴方位更加集中一些，说明在共轭阻抗法中，二维有效因子e2d确实能够突出模型的纯二维性.从加权处理后的统计云图中还可以看出一些维性信息：在经过加权处理后，剖面左侧、高频部分(浅部)、低频部分(深部)的主轴方位的显著性明显降低，表明这些部位的纯二维性较弱.这和模型的实际情况是一致的，图 4给出的维性分布结果也印证了这一点.由于区域构造是纯二维的，故而二维有效因子的形状与一维偏离度基本一致，即在此模型背景下，一维性弱的地方，就意味着二维性好.
 | 图 3 测线L17上所有测点、频点的电性主轴统计成像，施加了维性加权 上：Swift法；下：共轭阻抗法(CCZ).由左至右依次为统计玫瑰图、测点分布云图、频率分布云图，虚线框内为三维小异常体影响区域. Fig. 3 Statistical images of geo-electrical strikes for L17 line with all sites-all frequencies data calculated with dimensional weighting The upper images are for Swift method， and the lowers are for CCZ method. From left to right，the graphs are respectively rose diagram，site-based cloud diagram and frequency-based cloud diagram. |
 | 图 4 测线L17的构造维性分布云图 由左至右依次为一维偏离度(S1d)、二维偏离度S2d和二维有效因子(e2d). Fig. 4 The cloud maps of dimensionality From left to right，the maps are 1D skew，2D skew and 2D effective factor，respectively. |
第二个模型是一个立交桥模型，如图 5所示.模型浅部(0.5 km以上)有一个三维斜柱体，走向为45°；模型深部(0.7 km以下)的二维棱柱体的走向为0°.本文利用此模型来进一步演示阻抗张量分解多测点-多频点统计成像新技术的应用.
 | 图 5 立交的斜棱柱模型：左为立体视图，右为平面视图 Fig. 5 An intersection model with two prisms |
图 6给出了立交斜棱柱模型的最佳主轴的统计成像结果.共轭阻抗法(CCZ)很好地确定了该模型深浅不同的构造走向特征.经过二维加权处理后的统计结果更加突出了模型深部的主构造信息，而压制了一维性或三维性强的结果统计.由图可知，若单纯根据统计玫瑰图的结果，极有可能将测区构造走向方位确定为0°角，而忽略了浅部的45°角走向的构造.整条剖面的多测点-多频点共主轴张量分解的结果印证了这一点.两种统计云图为我们提供了更加丰富的信息，从中可以清楚地发现，在剖面中部的高频部分存在一个电性构造走向为45°角的浅部构造，这与真实模型是完全吻合的.
 | 图 6 测线L44上所有测点和频点的电性主轴统计成像(CCZ方法) 上：未施加维性加权，下：施加了维性加权.由左至右依次为统计玫瑰图、测点分布云图、频率分布云图、共主轴张量分解结果. Fig. 6 Statistical maps of geo-electrical strikes for L44 line with all sites and all frequencies data calculated without dimensional weighting(the upper) and with dimensional weighting(the lower)by using of CCZ method From left to right，the graphs are respectively rose diagram，site-based cloud diagram，frequency-based cloud diagram and rose diagram which results from impedance decomposition to a fixed strike azimuth based on all sites and all frequencies data. |
图 7为依据图 6的频率统计云图所给出全剖面高、低频分段统计结果.可以看出，在高频段(50～1000 Hz)，整条剖面上电性主轴方位的一致性很好，为45°左右；而在低频段，电性主轴方位很明显就是0°.高频段基本上看不到其他主轴方位的构造.在低频段，从统计玫瑰图上无法分辨出偏离0°主轴方位的信息，但从测点分布云图上则可以看出，在剖面中部，也就是模型中的两个柱体的交叉部分，存在既非0°也非45°的电性主轴方位.这一主轴方位信息是虚假的，是由于模型的三维特性造成的.在分析二维反演结果时，需要考虑这些部位的三维性可能造成的影响.此外，从图 7中还可以看出，无论高频段还是低频，多测点-多频点的统计成像结果与多测点-多频点共主轴方位张量分解结果具有非常好的一致性，表明二者的结果可相互印证.
 | 图 7 测线L44上电性主轴多频段统计成像(CCZ方法) 上：高频段(50~1000 Hz)，下：低频段(0.001~50 Hz).由左至右依次为统计玫瑰图、测点分布云图、共主轴张量分解结果.施加了二维加权. Fig. 7 Statistical maps of geo-electrical strikes for L44 line with two different frequency b and s data(CCZ method) The upper maps are those for the high frequency b and (50~1000 Hz) and the lower maps are those for the low frequency b and (0.001~ 50 Hz). From left to right，the graphs are respectively rose diagram，site-based cloud diagram and rose diagram for the result from impedance decomposition to a fixed strike azimuth based on multi-frequencies data. |
图 8展示了从剖面分段方面进一步统计成像的结果.根据图 6的测点序列统计云图，将整条测线分为三个部分.由图 8可以看出，剖面左段(-1.5～0.4 km)和右段(0.4～1.5 km)部分，电性主轴分布一致性很好，在频率域均可分为高频段和低频段两个部分，高频段表现了浅部45°走向的斜柱体的影响，低频段体现了深部0°角走向的柱体的影响.由于深部0°走向构造为模型的主构造，故左、右两个剖面段的全频段共主轴方位的张量分解结果均体现了深部0°走向的电性构造，而浅部45°走向的电性构造基本没有体现出来.
 | 图 8 测线L44上多测点-全频段的电性主轴统计成像(CCZ方法) 左、中、右分别为不同剖面范围(见图上所注)的统计结果；上、中、下依次为多测点-多频点统计玫瑰图、 频率分布云图和多测点-多频点共主轴张量分解结果. Fig. 8 Statistical maps of geo-electrical strikes for the three profile sections of L44 line(CCZ method) From left to right，the graphs are statistical results for different profile section(showing in the upper of each graph). the upper，the middle and the lower graphs are rose diagram for all the frequencies data，frequency-based cloud diagram，rose diagram rose diagram for the result from impedance decomposition to a fixed strike azimuth based on multi-sites data，respectively. |
剖面中间部分(-0.4~0.4 km)的电性主轴较为复杂，出现了三种显著性差别不大的方位：高频段的45°(或-45°)、中频段的66°(或-24°)、低频段的5°(或95°)的走向方位.其中高频段45°走向方位的频率范围延伸得更低，到20 Hz左右.全频段共主轴张量分解结果为75°，体现了高、中、低三个不同频段的综合影响.从图 9可以看出，在这个剖面中部的中频段，一维偏离度S1d和二维偏离度S2d都较大，而二维有效因子e2d相对较小.浅部斜柱体和深部的直柱体在此交叉重叠，故剖面中部这一较为复杂的主轴方位统计特征主要是该处模型较强的三维性造成的.
 | 图 9 测线L44的维性分布云图(CCZ方法) 由左至右依次为一维偏离度(S1d)、二维偏离度S2d、二维有效因子(e2d). 图中白色块状是被删除掉的质量差的频点数据. Fig. 9 The cloud maps of dimensionality From left to right，the maps are 1D skew，2D skew and 2D effective factor，respectively. |
由于可视化的同步实现，本文所发展的多测点-多频点的统计成像技术已多次应用于大地电磁(MT)、音频大地电磁(AMT)、张量可控源音频大地电磁(TCSMT)等实测资料的处理分析之中.这里以一个张量可控源音频大地电磁法(TCSMT)实例来简单演示这一技术在实测资料处理分析中的应用情况.
在南北地震带中南段的地震活断层探查中，我们联合北京欧华联公司开展了一次张量可控源技术的活断层探测试验.由于测区干扰大，加之对张量可控源新技术缺乏野外施工经验，观测的数据质量一般.此外，观测数据中还存在近场的影响.直接对原始数据进行张量分解分析，结果会很不可靠.图 10显示了采用本文提出的去噪、二维性加权等手段处理前后的主轴方位的统计结果对比.可以看到，原始数据的统计结果主轴方位的显著性很低，而且其主 要方位与最终得到的方位不同.去除质量不好的数据以后，统计结果的轮廓显得较为清晰一些；进一步采用二维性加权以后，得到了非常理想的统计玫瑰图，北东东向(结合地质构造排除90°的模糊性)的电性构造非常清楚.统计结果与多测点-多频点共主轴张量分解结果(NEE61°)非常一致.
 | 图 10 所有测点、频点的电性主轴统计玫瑰图(CCZ方法) 由左至右依次为原始数据统计结果、去除低质量频点数据统计结果、去除低质量数据 并进行二维性加权统计结果、多测点-多频点共主轴张量分解结果. Fig. 10 The statistical rose diagrams(CCZ method) From left to right are rose diagram for the raw data with all sites and all frequencies，for the good quality data，for the good quality data with structure dimensional weighting， and for the result from impedance decomposition to a fixed strike azimuth based on multi-sites，multi-frequencies data，respectively. |
从图 11可以看出，在未施加二维性加权以前，依据测点分布云图可将剖面分为4个构造分区，其 中最东面(测点序号0—10，右为东)构造区没有明显的主轴方位，中间偏东构造区(测点序号11—28)主轴方位约为30°，中间偏西构造区(测点序号29—39)主轴方位约为70°，最西部构造区(测点序号40—52)主轴方位约为55°；依据频率分布云图可将测区分为高、低两个频段.高频段(>1500 Hz)在75°附近主轴方位有一定的集中，但显著性不强；低频段(<1500 Hz)的主轴方位主要集中在30°～70°之间，也较为分散.在施加二维性加权以后，两种统计云图都明显集中.沿剖面可分为东西两个构造分区，西部(测点序号31—52)主轴方位更加集中于65°左右，而东部(测点序号0—30)主轴方位不明显.沿频率的分区和位置均未改变，但高频段(>1500 Hz)显著性更弱，低频段(<1500 Hz)则更集中于65°附近.
 | 图 11 电性主轴统计成像(CCZ方法) 由左至右依次为未施加维性加权的测点分布云图、频率分布云图和施加了维性加权的 测点分布云图、频率分布云图.其中，低质量数据已被删除. Fig. 11 The statistical cloud maps of geo-electrical strike(CCZ method) From left to right，the maps are site-based cloud image(no dimensional weighting)，frequency-based cloud image(no dimensional weighting)，site-based cloud image(using dimensional weighting)，frequency-based cloud image(using dimensional weighting). Data pointes with poor quality are deleted in all of these works. |
主轴方位的空间分布与构造维性具有较好的对应性.主轴方位显著的地方，意味着线性构造发育，二维性较强.因此，依据上述的分析，我们可推知测区东部地区二维性弱，一维性或者三维性较强，而西部地区的深部二维性较好.这一结果得到了图 12所示的维性分析结果的验证.一维偏离度S1d分布图显示，东部及整个剖面有效数据的高频部分，一维性很好.二维偏离度S2d分布图则显示，剖面最西端具有较强的三维特性.二维有效因子则综合了二者数据，突出了二维性强的区域.可以看出二维有效因子最大的区域正好对应于图 11中主轴方位的统计显著性最强的剖面和频率范围.
 | 图 12 沿剖面的维性分布云图 由左至右依次为一维偏离度(S1d)、二维偏离度S2d、二维有效因子(e2d). 图中白色块状是被删除掉的质量差的频点数据. Fig. 12 The cloud maps of dimensionality From left to right，the maps are 1D skew，2D skew and 2D effective factor，respectively. |
通过上述多测点-多频点张量分解结果的统计成像分析，我们不仅获得了主轴方位和区域响应数据，同时对于测区电性构造的空间分布也有了相当多的了解.这些由电磁观测数据独立给出的构造几何信息，与已知的地质结果具有非常好的对应性.这说明在二维反演以前，通过本文给出的多测点-多频点统计成像技术，可以分析和掌握测区构造的几何分布特征，为二维反演结果的合理解释提供了新的依据.此外，根据图 11，还可以进行更细致的多测点-多频点统计成像分析，但限于篇幅，这些工作不再赘述.
5 讨论与结论在大地电磁法(包括音频大电磁法和张量可控源音频大地电磁法)中，电性主轴的确定是二维反演解释以前数据处理与分析的最为关键的一个步骤.从Swift开始，围绕阻抗张量所发展的一系列数据处理技术——多种阻抗张量分解算法，已经成为大地电磁数据处理的核心技术.然而，由于其理论上较为复杂，涉及到参数较多，分析不方便等，使得这一技术所具备的大量信息难以得到充分的利用.本文提出了两种新的统计分布云图，定义了二维有效因子e2d，考虑了统计过程中数据观测误差的影响，实现了新的多测点-多频点阻抗张量共主轴分解算法，并完成了整个分析过程的可视化，最终建立了一套大地电磁阻抗张量分解多测点-多频点的统计成像分析技术.
利用本文所提出的主轴方位频率分布云图，在去除低质量频点数据，并施加二维性影响以后，可获得沿频率分布的电性主轴的统计成像结果；进而，分频段执行共主轴多测点-多频点分解计算，获得高、低频各自统一的主轴方位，从而可区分高、低频主轴方位是否一致.因此，采用本文所提出的统计成像新技术，为解决Jones课题组与Anahnah课题组有关北非阿特拉斯山脉岩石圈电性主轴的争论提供了良好的途径.
从本文所给出的研究结果来看，数据质量对电性主轴的统计成像结果具有很大的影响.本文通过可视化手动删除来解决这一问题，与反演前的数据编选采用了类似的手段.实践表明这样做是很有效的.不过，本文所给出的算例统计中，还没有考虑数据观测误差值大小的影响.采用类似反演的处理方式，使用误差棒做分母进行加权，效果可能会更好.
本文引入的新的维性描述参数——二维有效因子e2d带来了“纯二维性”的思想.这一参数的定义主要是基于，电性结构的二维性越强，其计算的电性主轴就越稳定，而二维性强的数据是二维反演中所主要依赖的数据，应该在统计过程中给予更高的权值.通过本文的算例可以看出，二维有效因子e2d能很好地压制一维结构和三维结构，不仅在主轴方位统计成像过程中加权使用，起到滤波的作用，还可能用来直接分析地下电性结构的二维性情况，发现线性构造的分布，为反演结果的地质解释提供辅助性信息.
通过方位角枚举法拟合单频点张量分解得到的相位，本文提出并实现了新的多测点-多频点共主轴阻抗张量分解算法.这一算法构成了多测点-多频点统计成像技术的最后一个步骤.实践表明，在主轴方位一致的前提下，采用张量分解技术易获取可用的区域视电阻率曲线数据，因此，最终用于二维反演的区域阻抗数据就是来自该算法的处理结果.本文算例结果已经验证了这一技术的有效性.
综上所述，本文所提出的大地电磁阻抗张量分解的多测点-多频点统计成像技术，能够为我们提供更多可靠的电性结构的几何信息和维性分布情况，为大地电磁资料的精细处理提供新的技术支撑力量.
致谢 本研究所用的张量可控源音频大地电磁数据(TCSMT)张量可控源数据由中国地震局地质研究所和欧华联公司联合观测提供，朱自串、朱学会、马云等人参与了野外数据采集工作.| [1] | Anahnah F, Galindo-Zaldivar J, Chalouan A, et al.2011.Deep resistivity cross section of the intraplate Atlas Mountains (NW Africa): New evidence of anomalous mantle and related Quaternary volcanism.Tectonics, 30(5), doi: 10.1029/2010TC002869. |
| [2] | Anahnah F, Galindo-Zaldivar J, Chalouan A, et al.2012.Reply to the comment by A.G.Jones et al.on "Deep resistivity cross section of the intraplate Atlas Mountains (NW Africa): New evidence of anomalous mantle and related Quaternary volcanism".Tectonics, 31(5), doi: 10.1029/2010TC003116. |
| [3] | Bahr K.1988.Interpretation of the magnetotelluric impedance tensor: regional induction and local telluric distortion.J.Geophys., 62: 119-127. |
| [4] | Bahr K.1991.Geological noise in magnetotelluric data: a classification of distortion types.Phys.Earth Planet.In., 66(1-2): 24-38, doi: 10.1016/0031-9201(91)90101-M. |
| [5] | Bibby H M, Caldwell T G, Brown C.2005.Determinable and non-determinable parameters of galvanic distortion in magnetotellurics.Geophys. J.Int., 163(3): 915-930, doi: 10.1111/j.1365-246X.2005.02779.x. |
| [6] | Cai J T, Chen X B, Zhao G Z.2010.Refined techniques for data processing and two-dimensional inversion in magnetotelluric I: Tensor decomposition and dimensionality analysis. Chinese Journal of Geophysics (in Chinese), 53(10): 2516-2526, doi: 10.3969/j.issn.0001-5733.2010.10.025. |
| [7] | Caldwell T G, Bibby H M, Brown C.2004.The magnetotelluric phase tensor. Geophys.J.Int., 158(2): 457-469, doi: 10.1111/j.1365-246X.2004.02281.x. |
| [8] | Chen X B, Zhao G Z, Zhan Y.2004.A visual integrated windows system for MT data process and interpretation.Oil Geophysical Prospecting (in Chinese), 39(Suppl.): 11-16. |
| [9] | Fischer G, Masero W.1994.Rotational properties of the magnetotelluric impedance tensor: the example of the Araguainha impact crater, Brazil. Geophys.J.Int., 119(2): 548-560, doi: 10.1111/j.1365-246X.1994.tb00141.x. |
| [10] | Groom R W, Bailey R C.1989.Decomposition of magnetotelluric impedance tensors in the presence of local three-dimensional galvanic distortion. J.Geophys.Res., 94(B2): 1913-1925, doi: 10.1029/JB094iB02p01913. |
| [11] | Groom R W, Bailey R C.1991.Analytic investigations of the effects of near-surface three-dimensional galvanic scatterers on MT tensor decompositions.Geophysics, 56(4): 496-518, doi: 10.1190/1.1443066. |
| [12] | Heise W, Bibby H M, Caldwell T G, et al.2007.Melt distribution beneath a young continental rift: The Taupo Volcanic Zone, New Zealand.Geophys.Res.Lett., 34(14), doi: doi: 10.1029/2007GL029629. |
| [13] | Heise W, Caldwell T G, Bibby H M, et al.2008.Three-dimensional modelling of magnetotelluric data from the Rotokawa geothermal field, Taupo Volcanic Zone, New Zealand. Geophys.J.Int., 173(2): 740-750, doi: 10.1111/j.1365-246X.2008.03737.x. |
| [14] | Jin G W, Sun J, Bai D H, et al.2003.Impedance tensor distortion decomposition of magnetotelluric data from West Sichuan-East Xizang. Chinese Journal of Geophysics (in Chinese), 46(4): 547-552. |
| [15] | Jin G W, Kong X R.2006.The Distortion and Decomposition of Magnetotelluric Impedance Tensor (in Chinese).Beijing: Seismological Press. |
| [16] | Jones A G, Kiyan D, Fullea J, et al.2012.Comment on "Deep resistivity cross section of the intraplate Atlas Mountains (NW Africa): New evidence of anomalous mantle and related Quaternary volcanism".Tectonics, 31(5), doi: 10.1029/2011TC003051. |
| [17] | Ledo J, Jones A G, Siniscalchi A, et al.2011.Electrical signature of modern and ancient tectonic processes in the crust of the Atlas mountains of Morocco.Phys. Earth Planet.In., 185(3-4): 82-88, doi: 10.1016/j.pepi.2011.01.008. |
| [18] | Mackie R L, Madden T R, Wannamaker P E.1993.Three-dimensional magnetotelluric modeling using difference equations;theory and comparisons to integral equation solutions.Geophysics, 58(2): 215-226, doi: 10.1190/1.1443407. |
| [19] | McNeice G W, Jones A G.2001.Multisite, multifrequency tensor decomposition of magnetotelluric data. Geophysics, 66(1): 158-173, doi: 10.1190/1.1826619. |
| [20] | Stagpoole V M, Bennie S L, Bibby H M, et al.2009.Deep structure of a major subduction back thrust: Magnetotelluric investigations of the Taranaki Fault, New Zealand. Tectonophysics, 463(1-4): 77-85, doi:10.1016/j.tecto.2008.09.035. |
| [21] | Swift C M.1967.A magnetotelluric investigation of an electrical conductivity anomaly in the southwestern United States[Ph.D.thesis].Massachusetts: Massachusetts Institute of Technology. |
| [22] | Toh H, Baba K, Ichiki M, et al.2006.Two-dimensional electrical section beneath the eastern margin of Japan Sea.Geophys.Res.Lett., 33(22), doi: 10.1029/2006GL027435. |
| [23] | Weckmann U, Ritter O, Jung A, et al.2007.Magnetotelluric measurements across the Beattie magnetic anomaly and the Southern Cape Conductive Belt, South Africa. J.Geophys.Res., 112(B5), doi:10.1029/2005JB003975. |
| [24] | Zhan Y, Zhao G Z, Wang L F, et al.2011.Deep structure in Shijiazhuang and the vicinity by magnetotellurics. Seismology and Geology (in Chinese), 33(4): 913-927, doi:10.3969/j.issn.0253-4967.2011.04.015. |
| [25] | Zhao G Z, Tang J, Liu T S, et al.1996.Preliminary interpretation of MT data along profile from Yanggao to Rongcheng: application of decomposition of MT impedance tensor. Seismology and Geology (in Chinese), 18(1): 66-74. |
| [26] | 蔡军涛, 陈小斌, 赵国泽.2010.大地电磁资料精细处理和二维反演解释技术研究(一)——阻抗张量分解与构造维性分析.地球物理学报, 53(10): 2516-2526, doi:10.3969/j.issn.0001-5733.2010.10.025. |
| [27] | 陈小斌, 赵国泽, 詹艳.2004.MT资料处理与解释的Windows可视化集成系统. 石油地球物理勘探, 39(增): 11-16. |
| [28] | 晋光文, 孙洁, 白登海等.2003.川西一藏东大地电磁资料的阻抗张量畸变分解. 地球物理学报, 46(4): 547-552. |
| [29] | 晋光文, 孔祥儒.2006.大地电磁阻抗张量的畸变与分解.北京: 地震出版社. |
| [30] | 詹艳, 赵国泽, 王立凤等.2011.河北石家庄地区深部结构大地电磁探测. 地震地质, 33(4): 913-927,doi:10.3969/j.issn.0253-4967.2011.04.015. |
| [31] | 赵国泽, 汤吉, 刘铁胜等.1996.山西阳高～河北容城剖面大地电磁资料的初步解释-阻抗张量分解技术及其应用. 地震地质, 18(1): 66-74. |