CN114895360A - 一种任意测点旋转角的大地电磁三维反演方法 - Google Patents

Abstract

本发明公开了一种任意测点旋转角的大地电磁三维反演方法，包括以下步骤：S1、对反演测点逐一旋转，通过数据曲线形状是否连续光滑、视电阻率和相位是否匹配等，找到各测点数据质量最好的方位，在该方位挑选频点数据得到反演观测数据；S2、依据反演测点空间展布方位构建模型网格，得到初始模型和对应的模型方位角，构成正演计算坐标系；S3、通过观测数据、初始模型和相应反演参数进行任意测点旋转角的大地电磁三维反演计算；S4、反演中，每次迭代将每个测点、每个频点的正演响应从正演计算坐标系旋转到该测点该频点数据方位角确定的观测数据坐标系进行数据拟合，得到模型修正量，更新模型。循环迭代，直到反演满足迭代终止条件，结束反演。

Description

一种任意测点旋转角的大地电磁三维反演方法

技术领域

本发明涉及固体地球物理和地质勘探领域，提出一种任意测点旋转角的大地电磁三维反演方法。

背景技术

大地电磁测深法(Magnetotellurics,MT)是以天然电磁场为场源来研究地球内部电性结构的一种重要的地球物理手段。其基本原理是：依据不同频率的电磁波在导电媒质中具有不同趋肤深度的原理，在地表测量由高频至低频的地球电磁响应序列，经过相关的资料处理来获得大地由浅至深的电性结构。它有不受高阻层屏蔽、对高导层分辨能力强，横向分辨能力较强等优点，也有纵向分辨能力随着深度增加而迅速减弱等缺点。

实践中，每一个大地电磁测点配置5道电磁场数据，包括2道电场组成的水平电场矢量和3道磁场组成完整的磁场矢量，均可在地表观测得到。由于观测装置的限制，常规大地电磁无法观测地表垂直电场分量。大地电磁响应函数阻抗张量Z表示为水平电场矢量E和水平磁场矢量H之间的转换关系：

E＝ZH

而另一响应函数倾子矢量T表示垂直磁场分量H_z与水平磁场矢量H_l之间的转换关系：

H_z＝TH_l

阻抗张量通过旋转公式

可以进行旋转，得到任意方位的阻抗张量。式中Z′为旋转后的阻抗张量，Z为旋转前的阻抗张量，R为旋转矩阵，θ为旋转角。

同理，倾子矢量也可以通过旋转公式

T′＝RT

得到任意方位的倾子矢量。式中T′为旋转后的倾子矢量，T为旋转前的倾子矢量。

实际观测中，一般在在观测方位进行的数据编辑，大部分情况下在观测数据方位的质量是比较好的，对观测数据的旋转，一般来说会导致数据质量下降，数据连续性减低(图1)。

三维正演求取全阻抗张量和倾子矢量，它们可以进行旋转从而获得任意方向坐标系的响应。由于绝大多数情况下，正演计算的误差很小，且没有方向偏好，因此正演响应满足阻抗张量和倾子矢量旋转的一切条件，可以任意旋转而不造成数据质量的下降(图2)。

反演解释中，选用阻抗张量、倾子矢量、视电阻率、相位等数据，利用反演方法软件进行计算得到研究区域的地下电性结构，然后基于反演结果，结合现有的地质和其他地球物理资料，对研究区域进行地球物理和地质解释。目前最先进效果最好的反演是MT三维反演。

目前使用广泛的三维反演程序主要有WSINV3DMT和ModEM。WSINV3DMT是基于数据空间的Occam反演，虽然计算量和内存占用小于传统的Occam反演，但是因为还是需要显式计算并储存完整的雅可比矩阵，所以计算量和内存占用还是非常大，计算耗时较长，不太适合大规模数据集的三维反演。ModEM采用共轭梯度方式直接对反演目标函数最小化，避免了计算和储存雅可比矩阵，极大地减小了内存占用和计算量。不像Occam反演在单次迭代中需要多次正演，ModEM反演一次迭代仅需两次附加的正演和一次数据拟合的正演，也就是总共三次迭代，计算速度要快很多。再加上很好的模块化结构和相对较小的内存需求，ModEM开源后迅速成为了标准化的MT三维反演方法程序，基于该程序已经取得了许多成果。

尽管现有三维反演技术已经较为成熟，但实测数据的反演实践中发现，现有方法尚存在一些可以改进提升的地方。在进行野外观测时，为了方便数据管理，大部分情况下会沿正南北-东西方位正交布极，得到统一以正北为0°角的观测数据。而在测点布设时，为增强数据对地下结构的约束，提高数据的利用效率，需要尽量垂直于构造走向布设，大部分情况下多个构造不会正好沿着正南北东西走向，这就导致测点的展布方位与测点的布极方位不一致(如图3A所示)。在MT三维反演中，以上情形可以归纳为两个坐标系：测点布极坐标系和测点展布坐标系。两个坐标系在大部分情况下是不一致的，二者存在一定的夹角，反演时需要基于其中一个坐标系构建网格模型(如图3B、图3C所示)。测点布极坐标系和测点展布坐标系的不一致将可能导致：1、反演中模型的正演响应的方位与野外测点布极方位不一致；2、反演区域无谓增大，从而反演网格数增多、计算量增大、非唯一性增强等问题。对于这两个问题，采用现有的三维反演技术只能解决其一，而不能同时克服。

另外，在进行实际数据的反演应用中，我们还发现大多测点在非布极方位视电阻率和阻抗相位曲线形态更好、更加连续光滑，数据质量更好(如图4、图5所示)，所有参与反演的测点进行旋转挑选数据后，每个测点的观测布极方位可能都不一致，而现有三维反演技术对此无能为力。

综上，现今三维反演存在以下问题：1、反演网格坐标系与测点展布坐标系一致，网格适配测点空间分布，但大多情况下正演响应不适配数据布极方位，需要将所有测点数据旋转一个固定的角度，数据质量会变差；2、反演网格坐标系与观测布极坐标系一致，正演响应适配数据布极方位，但大多情况下网格不适配测点空间分布，需要在大量的没有测点覆盖的非探测目标区剖分网格，网格数变多；3、我们发现，每个测点数据质量最好的方向不一定是其观测布极方位，且大多不一致，这就要求反演中每个测点旋转角不一样。

发明内容

本发明目的是针对上述问题，提供一种操作简单、提高反演效果的任意测点旋转角的大地电磁三维反演方法。

为了实现上述目的，本发明的技术方案是：

一种任意测点旋转角的大地电磁三维反演方法，所述方法包括以下步骤：

S1、对参与反演的测点进行逐一旋转，通过数据曲线形状是否合理、是否连续光滑、视电阻率和阻抗相位是否相匹配等，找到每个测点数据质量最好的方位，在该方位上进行频点数据挑选，得到各测点参与三维反演的观测数据；

S2、依据参与反演的测点空间展布方位构建模型参数网格，得到反演的初始网格模型和对应的模型网格方位角，构成正演计算坐标系；

S3、通过参与三维反演的观测数据、反演初始网格模型和相应的反演参数进行任意测点旋转角的大地电磁三维反演计算；

S4、在反演计算过程中，每一步迭代均需将每个测点、每个频点的正演响应数据从由模型网格方位角确定的正演计算坐标系旋转到该测点对应频点的由数据方位角确定的观测数据坐标系求取雅可比矩阵并进行数据拟合，得到模型参数修正量，更新模型参数；依次循环迭代，直到反演满足迭代终止条件，结束反演计算。

进一步的，所述步骤S3中，进行任意测点旋转角的大地电磁三维反演计算需要构建目标函数；任意测点旋转角的大地电磁三维反演的目标函数为：

其中，φ(m)为目标函数，m为模型参数向量，d为观测数据向量，F_R(m)为模型方位计算时旋转到测点数据方位的正演响应向量，

为数据误差相关矩阵，λ为正则化因子，m₀为参考模型向量，

为模型协方差矩阵。

进一步的，所述步骤S4中雅可比矩阵的计算公式为：

式中：

其中，Z为模型方位的正演阻抗张量，Z_R为观测数据方位的正演阻抗张量，R为旋转矩阵。

与现有技术相比，本发明具有的优点和积极效果是：

本发明针对当前三维反演技术所遇到的问题，提出了一种任意测点旋转角的大地电磁三维反演方法。该方法以旋转正演响应为核心，每个测点每个频点观测数据在不同的旋转角下进行三维反演，在反演中通过把模型响应旋转到观测数据方位进行数据拟合，实现了数据与网格的分离。在测点展布方位构建反演模型，在每个测点观测数据曲线最光滑方位挑选反演数据反演，既保证了参与反演的测点数据质量，又减少了无效的网格单元，减小了反演操作的整体计算量，提升了数据拟合水平，从而得到更可靠精准的反演结果，短了计算时间，提高了计算效率和结果准确性。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为观测数据旋转前后视电阻率曲线示意图；其中，上排为三个测点在观测布极方位的视电阻率曲线，下排为对应测点在旋转-45°到测点展布方位后的视电阻率曲线；

图2为正演响应数据旋转前后视电阻率曲线示意图；其中，上排为三个测点在观测布极方位的视电阻率曲线，下排为对应测点在旋转-45°到测点展布方位后的视电阻率曲线；

图3为三维反演的数据选择示意图；其中，A为某研究区域测点分布图，测点观测布极方位为0°，测点展布方位为-45°；B为基于测点布极坐标系构建的反演模型核心区域(51x56)；C为基于测点展布坐标系构建的反演模型核心区域(31x58)，实心箭头表示测点布极方位，虚线箭头表示测点展布方位；

图4为观测方位相位与其它方位的相位曲线对比示意图；

图5为观测方位的视电阻率与其它方位的视电阻率曲线对比示意图；

图6为正演计算合成数据的理论模型示意图；其中，a为测点展布方位0°的浅部两异常体顶部0km的水平切片；b为测点展布方位0°的深部异常体顶部5km的水平切片；c为黑线处垂直切片；d为测点展布方位45°的浅部两异常体顶部0km的水平切片；E为测点展布方位45°的深部异常体顶部5km的水平切片；黑点表示测点，黑线表示垂直切片在水平切片位置；

图7为不同坐标系构建模型网格核心区示意图；其中，A为在测点展布方位45°构建的网格核心区，水平网格单元大小1767m，网格数为81x81x82；B为观测布极方位0°构建的网格核心区，水平网格单元大小2500m，网格数为51x51x82；

图8为视电阻率相位数据反演RMS对比示意图；其中，NR表示基于观测布极方位构建网格，不旋转观测数据和正演响应，用传统三维反演方法计算；RD表示基于测点展布方位构建网格，旋转观测数据，不旋转正演响应，用传统三维反演方法计算；RDR表示基于测点展布方位构建网格，旋转观测数据和正演响应，用新方法计算；

图9为合成数据结果模型显示图；其中，黑点表示测点，黑线表示垂直切片在水平切片位置；NR为基于观测布极方位构建网格，不旋转观测数据和正演响应，用传统三维反演算法反对角线视电阻率相位数据反演的结果模型；RDZ为基于测点展布方位构建网格，旋转观测数据，不旋转正演响应，用传统三维反演算法反对角线视电阻率相位数据反演的结果模型；RDR为基于测点展布方位构建网格，旋转观测数据和正演响应，用新算法计算反对角线视电阻率相位反演的结果模型；

图10为RDR反演两相邻测点拟合曲线图；其中，左侧为测点L4_7拟合曲线，数据旋转角为77°，右侧为测点L4_8拟合曲线，数据旋转角为-14°；

图11为两个模型网格核心区示意图；其中，A为在观测布极方位0°构建的网格核心区，模型网格数为53x49x124；B为在测点展布方位42°构建的网格核心区，模型网格数为46x49x120；

图12为不同反演结果切片示意图；其中，NR为在观测布极方位构建模型不旋转观测数据和正演响应的传统方式反演，10km、15km和20km深度的结果切片；RD为在测点展布方位构建模型旋转观测数据不旋转正演响应的传统方式反演，10km、15km和20km深度的结果切片；RDR为各测点观测数据旋转到最佳方位，在测点展布方位构建模型旋转正演响应的任意测点旋转角反演，10km、15km和20km深度的结果切片。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

本发明提出了一种任意测点旋转角的大地电磁三维反演方法，在反演过程中，正演计算完整的大地电磁传输函数响应(阻抗张量、倾子矢量)，通过旋转传输函数去拟合观测数据。在正演响应旋转过程中，每个测点、每个频点的旋转角可以是任意的，因此，采用该算法进行反演，每个测点、每个频点的观测数据可以不同，可以在任意的旋转角下进行三维反演；同时，正演和反演各自自由的旋转，使得数据布极方向与网格展布方位分离，可以同时保证反演模型空间极小化和反演数据质量极优化，提高了反演计算的效率和反演结果的可靠性。

本发明所述的任意测点旋转角的大地电磁三维反演算法，通过旋转正演响应去拟合观测数据，而数据可以旋转到任意的方位，只要这个方位能使得数据质量最好。

任意测点旋转角的大地电磁三维反演目标函数如下

其中，φ(m)为目标函数，m为模型参数向量，d为观测数据向量，F_R(m)为正演响应向量，

为数据误差相关矩阵，λ为正则化因子，m₀为参考模型向量，

为模型协方差矩阵。

基于链式求导法则，正演向量F_R(m)的雅可比矩阵为：

式中：

其中，Z为模型方位的正演阻抗向量，Z_R为观测数据方位的正演阻抗向量，R为旋转矩阵。对倾子T的旋转也可以做类似的操作。

由以上公式发现，阻抗对模型偏导数的旋转与阻抗本身的旋转在形式上完全一致，因此，我们可以在已有的成熟的反演算法程序的基础上很容易地实现任意测点旋转角的反演算法。

本发明的任意测点旋转角的大地电磁三维反演流程包括以下步骤：

步骤一：对参与反演的测点逐一旋转找到数据质量最好的方位，在该方位上进行数据挑选，保存测点观测数据的方位角。

步骤二：在测点展布方位构建反演模型，保存模型方位角。

步骤三：准备好数据文件、模型文件和参数文件，用任意测点旋转角的大地电磁三维反演算法程序进行反演。在反演迭代中，将每个测点的响应数据从模型方位旋转到观测数据方位进行拟合，求出模型修正量，直到最后得到结果模型。

下面通过具体的实例并结合附图6～图12对本发明做进一步的详细描述。

首先构建一个包含局部异常体的理论模型，如图6A、图6B和图6C所示。理论模型背景电阻率为100Ωm，包含了三个局部异常体，浅部异常体电阻率值为10Ωm(20kmx40kmx5km)和100Ωm(20kmx40kmx5km)，深部异常体电阻率为1Ωm(35kmx40kmx5km)。测点以5km测点据按照16x16的网状分布，共256个测点，设置频率40个从0.0001Hz～1000Hz，用传统三维反演算法进行全阻抗正演计算得到0°方位的响应数据，然后将测点展布旋转45°(图6D和图6E)，同时响应数据继续保持在0°方位作为观测数据。

然后设计三组反演实验：一组为传统的在观测布极方位构建模型的反演；二组为传统的在测点展布方位构建模型的反演；三组为在测点展布方位构建模型并任意旋转观测数据模拟不同测点有不同最佳旋转角的反演。三组实验用于测试在测点展布方位与观测布极方位不一致的情况下，按照传统方式进行反演和进行任意测点旋转角的反演的过程与结果区别。每组测试都采用视电阻率相位数据进行反演。

一组为传统的观测布极方位(0°)构建模型的反演，构建初始模型网格，设置水平网格单元大小为1767m(2500xsin(45°))构建原始模型，其核心区如图7a所示，设置水平延展网格10个，增长步长为1.5，模型网格数为81x81x82，设置模型电阻率为100Ωm。反演参数方面，设置正则化因子为1000，变化步长为5，目标RMS为1。反演数据方面，所有测点都参与反演，采用0.001Hz～437Hz共31个频率进行计算。设置视电阻率门槛误差为2％的相对误差，相位门槛误差为0.573°。

二组为传统的测点展布方位(45°)构建模型的反演，将观测数据从布极方位(0°)旋转到展布方位(45°)，在测点展布方位(45°)构建初始模型网格，设置水平网格单元大小为2500m构建原始模型，其核心区如图7b所示，设置水平延展网格10个，增长步长为1.5，模型网格数为51x51x82，设置模型电阻率为100Ωm。反演参数和反演数据设置与组一一致。

三组为测试任意测点旋转角的反演，将观测数据进行随机旋转来模拟数据在曲线最光滑的方位时的情况。在测点展布方位(45°)构建初始模型网格，网格参数与组二一致。反演参数和反演数据设置与组一一致。

从表1可以看到，因为测点的展布方位为45°，NR(表示基于观测布极方位构建网格，不旋转观测数据和正演响应，用传统三维反演算法计算)反演整体网格数过多，计算量加大，计算耗时增加，反演的非唯一性增加，造成结果的不稳定。而RD(基于测点展布方位构建网格，旋转观测数据，不旋转正演响应，用传统三维反演算法计算)和RDR(基于测点展布方位构建网格，旋转观测数据和正演响应，用新算法计算)整体网格数不到传统方法的网格数的一半，计算量减少，计算耗时急剧地减小。最后的RMS都小于1，达到收敛条件。

表1反演收敛参数比较(反演数据：反对角视电阻率相位)

反演形式	网格数	计算时长	迭代次数	RMS
					NR	538002	45:21:05	76	0.980
RD	213282	21:57:28	93	0.997
					RDR	213282	09:33:11	71	0.995

表1中：NR表示基于观测布极方位构建网格，不旋转观测数据和正演响应，用传统三维反演算法计算；RD表示基于测点展布方位构建网格，旋转观测数据，不旋转正演响应，用传统三维反演算法计算；RDR表示基于测点展布方位构建网格，旋转观测数据和正演响应，用新算法计算；

图8为反对角线视电阻率相位数据反演并计算各个测点RMS的分布图，NR、RD和RDR的反演都达到了RMS≤1的拟合要求。总体来看，任意测点旋转角的反演很好的恢复了理论模型。

合成数据不同反演的结果如图9所示，可以看出，所有反演都基本还原了原始模型，都达到了预期，而且，因为网格构建更加合理，RD和RDR的反演对浅部高导体的还原要更好。针对RDR的反演，选取同一条测线上两相邻测点查看拟合情况(如图10所示)。这两个测点的相对旋转角很大，接近90°，可以看出，不同旋转角的观测数据都拟合很好。

以某研究区域三维反演为例子来验证新反演算法结果可靠性。

测点数为84，数据频率范围从320Hz～6.9E-5Hz共46个频率，测点观测布极方位为0°，测点展布方位为42°。设计三个反演计算，分别为NR(基于观测布极坐标系构建反演模型，观测数据与模型网格保持一致，用传统三维反演算法进行计算)、RD(基于测点展布坐标系构建反演模型，观测数据与模型网格保持一致，用传统三维反演算法进行计算，作为对照组)和RDR(在测点展布方位构建模型，每个测点观测数据都旋转到最合适的角度，用新算法进行计算)。基于测点展布和观测布极坐标系构建的模型核心区域如图11所示，网格单元都为2000m，基于测点展布方位构建的模型网格数为46x49x120，基于观测布极方位构建的模型网格数为53x49x124，都采用100Ωm的均匀半空间模型作为初始模型。数据方面，在对数据进行挑选后，采用160Hz～2.75E-4Hz共40个频率的视电阻率相位数据进行反演，设定门槛误差为5％。

经过反演后，反演收敛参数如表2所示。

表2反演收敛参数的比较(反演数据：反对角视电阻率相位)；

反演形式	网格数	数据数(％)	计算时长	迭代数	RMS
						NR	322028	9919(74.69％)	47:56:18	80	3.865
RD	270480	8991(66.90％)	28:04:48	83	3.874
						RDR	270480	10058(74.84％)	40:33:43	94	3.513

表2中:NR表示基于观测布极方位构建网格，不旋转观测数据和正演响应，用ModEM计算；RD表示基于测点展布方位构建网格，旋转观测数据，不旋转正演响应，用ModEM计算；RDR表示基于测点展布方位构建网格，旋转观测数据和正演响应，用新算法计算；

从表2中可看出，NR网格数大，反演数据数适中，计算时长最长；RD网格数小，反演数据数最小，计算时长最短；RDR网格数小，反演数据数最多，计算时长适中，RMS最小。

综上所述，反演的网格数最小，参与反演的数据数最多，最后RMS最小，我们认为RDR的结果是更加可信的。

本发明提出了一种任意测点旋转角的大地电磁三维反演方法。该方法以旋转正演响应为核心，每个测点每个频点观测数据在不同的旋转角下进行三维反演，在反演中通过把模型响应旋转到观测数据方位进行数据拟合，实现了数据与网格的分离。在测点展布方位构建反演模型，在每个测点观测数据曲线最光滑方位挑选反演数据反演，既保证了参与反演的测点数据质量，又减少了无效的网格单元，减小了反演操作的整体计算量，提升了数据拟合水平，从而得到更可靠精准的反演结果，短了计算时间，提高了计算效率和结果准确性。

Claims (3)

1.一种任意测点旋转角的大地电磁三维反演方法，其特征在于：所述方法包括以下步骤：

S1、对参与反演的测点进行逐一旋转，通过数据曲线形状是否合理、是否连续光滑、视电阻率和阻抗相位是否相匹配，找到每个测点数据质量最好的方位，在该方位上进行频点数据挑选，得到各测点参与三维反演的观测数据；

S2、依据参与反演的测点空间展布方位构建模型参数网格，得到反演的初始网格模型和对应的模型网格方位角，构成正演计算坐标系；

S3、通过参与三维反演的观测数据、反演初始网格模型和相应的反演参数进行任意测点旋转角的大地电磁三维反演计算；

S4、在反演计算过程中，每一步迭代均需将每个测点、每个频点的正演响应数据从由模型网格方位角确定的正演计算坐标系旋转到该测点对应频点的由数据方位角确定的观测数据坐标系求取雅可比矩阵并进行数据拟合，得到模型参数修正量，更新模型参数；依次循环迭代，直到反演满足迭代终止条件，结束反演计算。

2.如权利要求1所述的任意测点旋转角的大地电磁三维反演方法，其特征在于：所述步骤S3中，进行任意测点旋转角的大地电磁三维反演计算需要构建目标函数，其计算公式为：

其中，φ(m)为目标函数，m为模型参数向量，d为观测数据向量，F_R(m)为模型方位计算后时旋转到测点数据方位的正演响应向量，

为数据误差相关矩阵，λ为正则化因子，m₀为参考模型向量，

为模型协方差矩阵。

3.如权利要求2所述的任意测点旋转角的大地电磁三维反演方法，其特征在于：所述步骤S4中，雅可比矩阵的计算为：

式中

其中，Z为模型方位的正演阻抗张量，Z_R为观测数据方位的正演阻抗张量，R为旋转矩阵。

Images