CN104375195A - 时频电磁的多源多分量三维联合反演方法 - Google Patents

Abstract

本发明是一种地球物理数值模拟和反演的可控源时频电磁的多源多分量三维联合反演方法。依据已知电阻率测井数据和地震勘探资料确定时频电磁三维反演初始电阻率模型，计算出水平方向x、y方向的最大和最小坐标，确定出三维反演在水平方向的反演范围，选择时频电磁多源多分量三维反演的垂直方向网格尺寸，计算时频电磁发射源在初始模型中的一次场，计算地下每个六面体之间的格林张量；计算每个时频电磁场源的导数，采用共轭梯度的迭代算法使目标函数最小化，迭代次后完成多场源多分量时频电磁三维反演。本发明实测数据处理后得到了工区地下介质的电阻率三维分布，满足了构造、断层、圈闭的解释需要。

Description

时频电磁的多源多分量三维联合反演方法

技术领域

本发明属于地球物理勘探的电磁法勘探数据处理领域，是一种地球物理数值模拟和反演的可控源时频电磁的多源多分量三维联合反演方法。

背景技术

上世纪90年的建场测深法用于油气藏探测并取得了非常好的应用效果。在此方法基础上，结合时间域瞬变电磁测深和频率域电磁测深的特点，2005年采用了时频电磁法。这种方法具有时间域电磁法和频率域电磁法的特点。

时频电磁法测量Ex和Hz两个分量，具有时间域瞬变电磁和频率域测深特点。近年，随着勘探目标体越来越复杂、勘探的程度越来越精细，勘探由二维（2D）转变成三维（3D），测线由一条测线变成网状的布设的多条测线。传统的处理方法以时间域和频率域一维（1D）反演、定性的振幅和相位参数求取为主，这些处理方法都不能进行三维数据处理，网状布设的时频电磁数据只有进行了三维的反演，才能对目标存储、圈闭的空间展布和范围了解的更清楚，这对提高时频电磁法在油气勘探中的应用效果非常重要，因此，研究时频电磁三维反演就非常有必要。

时频电磁法在进行网状测线布设时候，并不是所测线共用一个发射，实际上是多个发射。因为在发射源位置固定后，测线距离发射源的最大偏移距是基本固定的，偏移距超过最大偏移距，测线上的接收器是无法采集到发射的信号。因此网状的测线必然存在多个发射源，当对网状测线上所有的测点的数据做三维反演时候，实际上是一个多场源多分量的三维电磁反演问题。

对于时频电磁法的多源多分量的三维反演问题而言，目前国内外都没有开展过相关的研究，相关的文献也没有。时频电磁法属于可控源电磁法，可控源电磁法中处理技术发展最快的就数海洋可控源电磁法，目主要用的是一维（1D）、二维（2D）和三维（3D）反演处理技术，也涉及到多个发射源的问题，但是主要处理电分量，磁分量是用来计算大地电磁的视电阻率。而且海洋电磁的野外观测方式和陆地上的时频电磁法是不相同的。

目前时频电磁主要用一维处理方法。随着勘探目标越来越复杂，用一维反演去刻画三维复杂目标体是不能满足理解释的精度要求。一个场源，几条测线的三维反演也存在同样的问题，反演结果只能是工区的一部分，不能一次反演出整个工区地下电阻率的分布。而且我们实际观察方式就是多个发射源，多个测线，多个分量的，所以只有采用多源、多测线、多分量的三维反演方法，才能将整个工区地下电阻的分布一次全部反演出来。

由于存在多个发射场源，发射场源的长度不同，所以不同发射场源的观测数据相差很大，在反演过程中如何使不同场源的数据起到的作用一样就显的非常关键。时频电磁法的三维反演中需要大量计算雅克比导数矩阵，如果采用传统的差分法，也就是两次正演相差除以电阻率的变化量，这种反演算法计算速度非常慢，很难实现快速成像。快速计算三维反演的雅克比矩是时频电磁法多源多分量三维反演的关键之一。如何解决这些技术问题，在已经发表的论文中没有发现有关的线索。

发明内容

本发明目的是提供一种以满足解释断层、构造、基地起伏形态、圈闭的含油气性评价需要的时频电磁的多源多分量三维联合反演方法。

本发明通过以下步骤实现：依据实测时频电磁发射频率的分布范围和分量类型，选择参加反演的数据；

所述的参加反演中时频电磁多源多分量的发射频率的范围为0.01-100Hz。

所述的参加反演中时频电磁多源多分量的分量为与场源平行的电场分量Ex和垂直磁场分量Hz。

2）依据已知电阻率测井数据和地震勘探资料确定时频电磁三维反演的初始1D电阻率模型，层厚度和电阻率值；

所述的电阻率测井数据采用深侧向的电阻率数据，依据地震剖面的反射同相轴、电测井数据的高低关系，确定1D背景模型的层厚度和层数，1D背景电阻率值大小通过电测井数据确定，同时保证测井数据的总纵电导曲线和1D背景电阻率模型的总纵电导曲线重合。

3）依据地面所有测点的水平方向x、y方向的坐标、计算出水平方向x、y方向的最大和最小坐标，确定出三维反演在水平方向的反演范围，依据目标体的深度再加上2000m得到三维反演的深度；用x，y，z三个方向的反演范围除以各自剖分网格尺寸，得到三维反演x，y，z三个方向的剖分网格数；

所述三维反演的水平方向x，y的反演范围包含所有的测点，网格尺寸是100m。

所述三维反演的垂直方向z的最大反演深度大于目标体的深度2000m，垂直方向z方向的网格尺寸由逐渐增加，第一个网格的尺寸是50m，最后一个网格的尺寸是500m。

4）选择正确的时频电磁多源多分量三维反演的垂直方向z网格的尺寸；

所述的选择是依据一个网格中只有一个电阻率值的原理，三维反演的垂直z方向的任何一个网格不能同时位于1D背景层两个地层中；

5）计算时频电磁发射源在1D电阻率初始模型中的一次场；

6）计算地下每个六面体对地面接收点的格林张量，计算地下每个六面体之间的格林张量；

7）计算每个时频电磁场源的Frechet导数，所述的计算时频电磁Frechet导数矩阵通过解析公式计算：

式中，E^a是水平电场Ex，H^a是垂直磁场，σ是电阻率的导数，和时电场和磁场的格林张量，E^b背景电场。r_j是到测点的距离，r是到剖分网格的距离，是与格林张量有关的系数矩阵，g是与背景场E^b和格林张量有关的系数矩阵。

8)依据正则化反演原理计算目标函数，所述的计算如下：

$P\left(m\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{W}_d^i\·{\mathrm{\φ}}^i\left(m\right)+\mathrm{\αS}\left(m\right)=\min ---\left(3\right)$

式中是第i个时频电磁发射源的数据的拟合误差；是第i个时频电磁发射源的权函数，S(m)是稳定化函数；α是正则化参数；Aⁱ正演算子；dⁱ是i个时频电磁发射源时，时频电磁观测点的电场E_x和磁场H_z分量。

9）采用共轭梯度的迭代算法使目标函数最小化；

所述的目标函数最小化计算时，在最小化第一次迭代过程中正则化参数α值为0.5，α的取值范围为0.05-5。在迭代过程中自动调节α的大小，如拟和误差增加，增加α；如果拟合误差减小，α减小。

所述的迭代算法中采用如下公式计算时频电磁数据加权矩阵

$W_d^i=\max (\mathrm{sum}\left({\mathrm{\φ}}_0^i\right)\·\·\·\mathrm{sum}\left({\mathrm{\φ}}_n^{\mathrm{TFEM}}\right))/\mathrm{sum}\left({\mathrm{\φ}}_0^i\right)---\left(4\right)$

式中：是第i个时频电磁发射源的初始模型的拟合误差；是第n个时频电磁发射源的初始模型的拟合误差；max是计算最大值函数；sum是求和函数。

10）重复步骤9），当多个场源的时频电磁反演的目标函数的拟合误差达到设定的误差标准或者迭代次数超过了设定的最大迭代次数后，完成多场源多分量时频电磁三维反演。

本发明对模型合成和实测数据做了实验处理，模型合成的数据反演结果与理论模型基本一样，说明本发明有效；实测数据处理后得到了工区地下介质的电阻率三维分布，满足了构造、断层、圈闭的解释需要。

附图说明

图1时频电磁电场和磁场振幅曲线，（left）Ex，（right）Hz。

图2测井曲线和总纵电导曲线，（left）测井曲线，（right）纵电导曲线。

图3一维背景电阻率模型。

图4某个工区测线和测点位置图，（left）旋转和平移前，（right）旋转和平移后。

图5某工区实测数据的时频电磁多源多分量3D联合反演结果的3D显示。

图6深度在6000m的水平方向电阻率平面分布图。

图7测线1的电阻率断面图。

具体实施方式

以下结合附图详细说明本发明。

1）依据实测时频电磁发射频率的分布范围和分量类型，选择参加反演的数据；时频电磁多源多分量的发射频率的范围为0.01-100Hz。图1是时频电磁电场和磁场振幅曲线，（left）Ex，（right）Hz。

2）依据已知电阻率测井数据和地震勘探资料确定时频电磁三维反演的初始1D电阻率模型，层厚度和电阻率值。图2测井曲线和总纵电导曲线，（left）测井曲线，（right）纵电导曲线，图2（left）是依据电阻率测井曲线变化规律，建立1D背景的电阻率值和层厚度。图2（right）是电测井数据的纵电导曲线与建立的1D背景模型的总纵电导曲线对比，保证使二者之间拟合。图3是一维背景电阻率模型。

3）依据地面所有测点的水平方向x、y方向的坐标、计算出水平方向x、y方向的最大和最小坐标，确定出三维反演在水平方向的反演范围，依据目标体的深度再加上2000m得到三维反演的深度；用x，y，z三个方向的反演范围除以各自剖分网格尺寸，得到三维反演x，y，z三个方向的剖分网格数。

图4是测线和测点位置图，（left）旋转和平移前，（right）旋转和平移后。经过旋转和平移后，最终的x方向的反演范围为-5000m—12600m，y方向的反演范围为0m—12400m，z方向的反演范围为0m—10000m。x、y和z方向的网格单元的尺寸都为200m，所以三维反演的网格数位89×63×51。

4）选择正确的时频电磁多源多分量三维反演的垂直方向z网格的尺寸（参数）。

5）计算时频电磁发射源在1D电阻率初始模型中的一次场。

6）计算地下每个六面体对地面接收点的格林张量，计算地下每个六面体之间的格林张量。

7）计算每个时频电磁场源的Frechet导数。Frechet导数矩阵的计算公式见公式（1）、（2）。

8）依据正则化反演原理计算目标函数，所述的计算公式见公式（3）。

9）采用共轭梯度的迭代算法使目标函数最小化，在最小化第一次迭代过程中正则化参数α值为0.5，α的取值范围为0.05-5。在迭代过程中自动调节α的大小，如拟和误差增加，增加α；如果拟合误差减小，α减小。

10）重复步骤9），当多个场源的时频电磁反演的目标函数的拟合误差达到设定的误差标准或者迭代次数超过了设定的最大迭代次数后，完成多场源多分量时频电磁三维反演。图5是某工区实测数据的时频电磁多源多分量3D联合反演结果的3D显示，图6是经过20次迭代反演后，得到了工区地下介质的三维电阻率分布。图6是深度在6000m的水平方向电阻率平面分布图，图7是测线1的电阻率断面图。从上面的三维反演结果可以看出，基本的电性层已经反演出，与测井数据吻合。三维反演清楚的将基地起伏形态反演出，

为精确划定潜山顶面埋深，提高地震的偏移成像和解释精度。

Claims (10)

1.一种时频电磁的多源多分量三维联合反演方法，特点是通过以下步骤实现：

1）依据实测时频电磁发射频率的分布范围和分量类型，选择参加反演的数据；

2）依据已知电阻率测井数据和地震勘探资料确定时频电磁三维反演的初始1D电阻率模型，层厚度和电阻率值；

3）依据地面所有测点的水平方向x、y方向的坐标、计算出水平方向x、y方向的最大和最小坐标，确定出三维反演在水平方向的反演范围，依据目标体的深度再加上2000m得到三维反演的深度；用x，y，z三个方向的反演范围除以各自剖分网格尺寸，得到三维反演x，y，z三个方向的剖分网格数；

4）选择正确的时频电磁多源多分量三维反演的垂直方向z网格的尺寸；

5）计算时频电磁发射源在1D电阻率初始模型中的一次场；

6）计算地下每个六面体对地面接收点的格林张量，计算地下每个六面体之间的格林张量；

7）计算每个时频电磁场源的Frechet导数；

8)依据正则化反演原理计算目标函数，所述的计算如下：

$P\left(m\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{W}_d^i\·{\mathrm{\φ}}^i\left(m\right)+\mathrm{\αS}\left(m\right)=\min ---\left(3\right)$

式中是第i个时频电磁发射源的数据的拟合误差；是第i个时频电磁发射源的权函数，S(m)是稳定化函数；α是正则化参数；Aⁱ正演算子；dⁱ是i个时频电磁发射源时，时频电磁观测点的电场E_x和磁场H_z分量；

9）采用共轭梯度的迭代算法使目标函数最小化；

10）重复步骤9），当多个场源的时频电磁反演的目标函数的拟合误差达到设定的误差标准或者迭代次数超过了设定的最大迭代次数后，完成多场源多分量时频电磁三维反演。

2.根据权利要求1的方法，特点是步骤1）所述的参加反演中时频电磁多源多分量的发射频率的范围为0.01-100Hz。

3.根据权利要求1的方法，特点是步骤1）所述的参加反演中时频电磁多源多分量的分量为与场源平行的电场分量Ex和垂直磁场分量Hz。

4.根据权利要求1的方法，特点是步骤2）所述的电阻率测井数据采用深侧向的电阻率数据，依据地震剖面的反射同相轴、电测井数据的高低关系，确定1D背景模型的层厚度和层数，1D背景电阻率值大小通过电测井数据确定，同时保证测井数据的总纵电导曲线和1D背景电阻率模型的总纵电导曲线重合。

5.根据权利要求1的方法，特点是步骤3）所述三维反演的水平方向x，y的反演范围包含所有的测点，网格尺寸是100m。

6.根据权利要求1的方法，特点是步骤3）所述三维反演的垂直方向z的最大反演深度大于目标体的深度2000m，垂直方向z方向的网格尺寸由逐渐增加，第一个网格的尺寸是50m，最后一个网格的尺寸是500m。

7.根据权利要求1的方法，特点是步骤4）所述的选择是依据一个网格中只有一个电阻率值的原理，三维反演的垂直z方向的任何一个网格不能同时位于1D背景层两个地层中。

8.根据权利要求1的方法，特点是步骤7）所述的计算时频电磁Frechet导数矩阵通过以下解析公式计算：

式中，E^a是水平电场Ex，H^a是垂直磁场，σ是电阻率的导数，和时电场和磁场的格林张量，E^b背景电场。r_j是到测点的距离，r是到剖分网格的距离，是与格林张量有关的系数矩阵，g是与背景场E^b和格林张量有关的系数矩阵。

9.根据权利要求1的方法，特点是步骤9）所述的目标函数最小化计算时，在最小化第一次迭代过程中正则化参数α值为0.5，α的取值范围为0.05-5。在迭代过程中自动调节α的大小，如拟和误差增加，增加α；如果拟合误差减小，α减小。

10.根据权利要求1的方法，特点是步骤9）所述的迭代算法中采用如下公式计算时频电磁数据加权矩阵

$W_d^i=\max (\mathrm{sum}\left({\mathrm{\φ}}_0^i\right)\·\·\·\mathrm{sum}\left({\mathrm{\φ}}_n^{\mathrm{TFEM}}\right))/\mathrm{sum}\left({\mathrm{\φ}}_0^i\right)---\left(4\right)$

式中：是第i个时频电磁发射源的初始模型的拟合误差；是第n个时频电磁发射源的初始模型的拟合误差；max是计算最大值函数；sum是求和函数。