井震约束的广域电磁数据反演及其在保靖页岩气储层评价中的应用

王延忠,沈金松，周正武,符　超，路肖尧

(1.神华地质勘查有限责任公司，北京 100083；2.中国石油大学 地球物理与信息工程学院，北京 102249；3.中国石油集团公司 物探重点实验室，北京 102249；4.长江大学 地球物理与石油资源学院,湖北 武汉 434023)



井震约束的广域电磁数据反演及其在保靖页岩气储层评价中的应用

王延忠1,沈金松2,3，周正武1,符超1，路肖尧4

(1.神华地质勘查有限责任公司，北京 100083；2.中国石油大学 地球物理与信息工程学院，北京 102249；3.中国石油集团公司 物探重点实验室，北京 102249；4.长江大学 地球物理与石油资源学院,湖北 武汉 434023)

目的是研究一种井震约束的广域电磁数据反演方法，用以处理湖南保靖地区采集的广域电磁数据,获取地层电阻率反演结果，然后用以得到研究区的地层分布格架、断层展布特征和页岩储层的纵横向分布。文中所用的方法包括三个方面。首先，用测井数据标定地震资料，对二维地震剖面进行处理解释，得到地震层位和断层分布；然后，用地震层位结果建立初始电阻率分布模型, 用测井分层电导率作为上下界约束，对广域电磁数据进行井震约束反演，得到研究区地下地层的电阻率分布；最后，对电阻率反演结果进行综合地质解释，建立页岩气储层的纵横向分布剖面，给出页岩气富集的有利分布位置。通过与3口测井及试气结果对比，初步证实了井震约束的广域电磁数据反演方法在保靖页岩气储层评价中的应用效果。

井震约束；广域电磁数据；电磁场灵敏度计算；正则化反演；页岩气储层

1　引　言

页岩气是典型的非常规天然气资源，产自极低孔渗以富有机质页岩为主的储集岩系中，以自生自储方式或以游离气，或以吸附气赋存于地下储层中。我国的页岩气勘探尚处于起步阶段，与美国等西方国家的页岩气储层相比，我国页岩气赋存区受构造破坏更严重,地表和地下地质条件更为复杂，致使勘探开发成本高，风险更大，急需综合地球物理勘探手段优选和精细评价有利分布区。鉴于页岩储层与非页岩区存在的电性差异大，除了必备的地震勘探外，许多页岩气勘探区还将电磁勘探作为补充手段。在本研究的目标区保靖区块选用了何继善[1-3]提出的广域电磁法。

广域电磁法是何继善从场的统一性出发，将“近区”、“过渡区”和“远区”有机地结合而提出的,他在深入分析可控源音频大地电磁(CSAMT)测量及其视电阻率定义的基础上，提出广域电磁测量方法和全区视电阻率的定义[2,3]。该方法改善了非远区的畸变效应，使得测深能在广大的、不局限远区的区域进行，在同等收发距上可勘探的深度增大。同时，用适合于全域且无需简化的公式计算的视电阻率，拓展了人工源电磁法的观测范围，提高了观测速度、精度和野外工作效率。在广域电磁测量响应的数值模拟和反演方面，国内学者也开展了多方面的研究。余云春[4]实现了各种水平电流源广域电磁法一维正演，并从理论响应分析了各种水平电流源和场的不同分量对不同地电模型的响应能力，比较了广域电磁法与其他电磁测深法的探测性能。李帝诠等[5]采用积分方程法实现了广域电磁响应的三维正演模拟，给出4种不同装置的测量结果，对比广域视电阻率与Cagniard视电阻率对典型三维地质体的分辨能力。

本文在对湖南保靖页岩气勘探区块的地质和地球物理特征分析的基础上，采用测井标定地震数据，对二维地震剖面进行处理解释，得到地震层位和断层构造分析结果。接着，用地震层位结果建立初始电导率分布模型，测井分层电导率结果作为上下界约束，对广域电磁数据进行反演，得到研究区地下地层的电阻率分布。最后，参考测井、地震的处理解释结果，对电阻率反演结果进行综合地质解释，刻画页岩储层的纵横向分布。实现了井震约束的广域电磁数据反演在湖南保靖地区电磁测量的反演处理，得到了研究区的地层分布格架、断层展布特征和页岩储层的分布剖面。

2　广域电磁视电阻率数据的预处理

2.1广域电磁视电阻率简介

为便于理解后续的电磁响应正演模拟中源的处理，文中采用的广域电磁数据采集时使用的电流源为一对接地电极A和B，且当正、负电流源A和B之间的距离dL(m)与观测点到场源的距离r(m)相比很小时，将之看成一对电偶极子。水平电偶极源的场分布用解析关系表达[1-3]。实际工作中使用的长导线作为发射源，AB距离较大，不能简化成点偶极源，需要将长导线剖分成若干小段，将每个小段用点偶极源近似，再通过点偶极源激发的电磁场响应的积分求解。由于在均匀大地表面布设激发电流I(A)，长为dL(m)的电偶源的电场关系前人做过详细的研究和分析[7-9]，长导线电磁场的积分计算也有何继善等[1-3]开展过系统研究，在此不再赘述。

需要说明的是，广域视电阻率的计算较之经典的卡尼亚(Cagniard)视电阻率计算更为复杂。一方面，它也需要将观测到的电位差、发送电流以及有关的几何参数代入视电阻率的定义式计算[1]，另一方面，由于电磁效应函数中也含有未知电阻率ρ(m/S)[3],需要采用迭代法逐次逼近求解[3,4,6]。同时，由于该视电阻率定义对观察点到发送源的距离没有任何限制，它适合于广大测量区域，而不是像可控源音频大地电磁(CSAMT)那样只适用于远区，因此，把该视电阻率称为广域视电阻率[2,3,10]。

2.2视电阻率数据的预处理

广域电磁测量数据中不可避免含有部分噪声，以及受到近地表非均匀体的静态位移影响。因此需要对野外测量的原始数据进行预处理。在广域电磁法施工中由于各种电磁噪音的存在，不可避免地给实测资料带来一定的误差，严重时使曲线的形态发生变化，图1给出了4个测点的实测原始电磁响应计算的视电阻率曲线，图中看到，不同频率上存在干扰造成的测量误差。所以在资料解释前必须进行去噪处理，一般需要观察原始记录曲线形态，剔除畸变和压制噪声。

由于工区地表起伏不平及近地表介质电阻率局部横向变化较大，需要对数据实施静态效应校正，以免后续的反演解释得出错误的结果[11,12]。本文是通过给出频率—视电阻率等值线图，通过对频率—视电阻率等值线图的分析，确定测线数据是否存在静态效应的影响，并对受到影响的数据进行空间滤波校正[13,14]。图2给出了6个测点的电磁响应计算的视电阻率曲线经过静态校正的结果，从图2中看到，经过空间滤波基本消除了随机干扰和静态效应影响。

图1　四个测点的实测原始电磁响应计算的视电阻率曲线(注:ρ/Ω·m=ρ/m·S-1)Fig.1　The apparent resisitivity curves field EM measurements of 4 measuring points (note:ρ/Ω·m=ρ/m·S-1)

图2　六个测点的电磁场响应计算的视电阻率经过静态校正的结果(注:ρ/Ω·m=ρ/m·S-1)Fig.2　The static correction results of the apparent resistivity derived from the EM responses of 6 stations (note:ρ/Ω·m=ρ/m·S-1)

3　广域电磁数据井震约束反演方法

文中广域电磁测量数据的井震约束反演是将在地表测得的视电阻率数据通过正演模拟和地震构造和层位信息约束的迭代反演方法，获得各测点地下不同深度的电阻率值，给出勘探剖面对应的地下地层的电阻率分布剖面。由于在地面的实测电磁场数据只是电流源激发的空间电磁场分布的极小部分信息，与地下未知的电阻率分布相比，仍是一个数据信息少得多的欠定问题，必须通过其他已知信息的约束及适当的正则化反演求解。

在本研究的广域电磁数据反演中，首先根据现有的地质认识和地震解释层位建立一个初始地层电阻率模型；其次，采用2.5维有限元方法计算地层模型的理论广域电磁响应，估算实际测量值与正演模拟理论值之间的差异；再次，采用激发源产生的电磁场与接收点假想源电磁场的互换原理计算理论电磁响应对地下电阻率参数的灵敏度，结合实测电磁场与模拟电磁场的差异估算地下电阻率模型参数的修正量；最后，重复正演模拟响应与实测电磁场误差计算与电磁场对电阻率参数的灵敏度估算过程，修改地电模型的电阻率值，使实测电磁场响应与正演模拟的理论值之差的平方和达到最小，得到最终的地下电阻率反演结果。为了文章的严谨和完整性，下文就电磁场响应正演模拟、电磁场对地下电阻率的灵敏度和井震约束正则化反演算法加以简单介绍。

3.1电磁场响应正演模拟

对于电流源的广域电磁法，频率域麦克斯韦方程为[15]：

(1)

针对保靖的二维广域电磁测量数据，用图3所示的二维介质模型模拟。设y方向与走向平行，在该方向上，电导率σ，介电常数ε和磁导率μ不变，在x-z平面内这些参数可以变化。将方程(1)中电场E和磁场H表示成3个分量，并消去Ex、Ez及Hx、Hz，得到关于Ey、Ey分量耦合方程，再对变换后的方程和源项作y方向的Fourier变换，即

(2)

(3)

(4)

(5)

(6)

图3　二维介质模型(y方向与走向平行)Fig.3　Two dimensional model of the conductivity (y direction is along with the strike)

KF=S

(7)

3.2灵敏度计算

电磁场响应对电导率的灵敏度定义为单位电导率变化引起的电磁响应变化，用(7)式的电磁场离散关系来考察电磁响应的灵敏度时，对于电流源的广域电磁法，由于源项与电导率参数无关，电磁响应的灵敏度可以写成[14,16-18]：

(8)

式中，p是电导率参数；∂F/∂p是三个电场或磁场分量对电导率的灵敏度；∂K/∂p是刚度矩阵对电导率的变化，只有反演单元上的单元刚度矩阵的元素不为零。

方程(8)看到，右端等效源项是由原电磁场与反演单元上的单元刚度矩阵变化的内积得到。因此,反演计算中,需要计算相当于反演单元数的正演问题,这在实际电磁数据的反演中是难于实现的。事实上，在实际反演计算中，只需要接收器位置的电磁场对电导率的灵敏度，因此可以采用电磁场的互易原理简化接收器位置电磁场灵敏度的计算[18]。

3.3井震约束的正则化最小二乘反演方法

Φ=‖WD(Dc-Do)-WDJAΔm‖2·

‖Wm(m-mp)‖2

(9)

其中，式中，D0是观测数据；Dc是计算的广域电磁波理论响应；JA是广域电磁理论响应对电阻率的偏导数组成的雅克比(Jacobian)矩阵；Δm是迭代反演中电阻率(电导率的倒数)的变化量；WD和Wm是测量数据和模型参数的加权矩阵，通常取各自的标准方差的倒数组成的矩阵；mp为先验信息确定的电阻率参数。

在实际资料的计算中，模型参数是未知的，其标准差无法得到，考虑到介质电性参数的渐变特性，文中取Wm为等间隔离散的拉普拉斯二维差分算子。对于测量数据加权矩阵WD，文中采用视电阻率的倒数组成的对角阵。式(9)的极小化目标泛函对应的正则化方程为：

(10)

式中，Δm迭代反演中电阻率(电导率的倒数)的变化量。

考虑到地震层位和测井电阻率的约束，方程(7)加上两个约束的解即为下边矩阵的最小二乘解[20]：

(11)

m=m(0)+Δm

(12)

其中，m(0)是上一次迭代的更新值。

考虑到式(11)的刚度矩阵条件数较大，在实现反演过程中采用Gauss-Newton最速下降方法和带先验信息约束的加权L2模正则化法[20]求解，以保证迭代反演的每一步具有一定的误差减小。

4　保靖区块广域电磁数据反演结果的页岩气储层评价应用

4.1工区位置及地质地球物理概况

图4给出了工区的测线位置和构造纲要, 纵横坐标为相对坐标，所有与地名相关的信息均采用了覆盖处理。该页岩气区块位于云贵构造活化区南部边缘，后期构造破坏强烈，主要构造线呈北北东向展布，向斜倾角平缓，背斜倾角较大，同时伴生断裂发育，主要断裂走向与褶皱轴线基本一致。所布置的测线基本与构造走向垂直。除北北东向构造外,还有北东向构造及其伴生断裂，少部分断裂走向为北西向。研究区页岩分布层位为古生界下志留统龙马溪组(S1ln)以及下寒武统牛蹄塘组(∈1n)。两套目标层埋深大约为0～2 500 m和0～5 000 m范围内。工区近地表普遍灰岩区覆盖，地下高陡地层发育。

图4　保靖区块广域电磁勘探测线分布Fig.4　Configuration of exploration lines of the wide field electromagnetic measurements in Baojing area

4.2广域电磁测量数据的定性分析

资料的定性分析是针对频率—视电阻率剖面进行的，依据不同地质构造、电性分布特征的大地电磁响应规律，分析提取原始资料中的地质信息，定性地把握地下电性层分布特征、地层起伏变化情况、局部构造、构造单元划分等，为进一步的定量解释提供依据，同时评价、检验、落实定量解释成果的可靠性。在广域电磁法中，对实测曲线类型的分析、比较是资料定性认识解释更准确地获得测区地质结论的重要组成部分。曲线类型定性地反映地下电性层的分布特征，如电性层数，相对埋深和各电性层间电阻率相对变化情况。特别是通过对测区内曲线类型的分析比较，可对测区的地质构造单元进行划分与归类，给出测区地质构造的定性概念。

本文利用各测点相应频率上的视电阻率值勾绘等值线图，得到图5所示的L5测线的频率—视电阻率拟断面图。分析该频率—视电阻率拟断面图，可以定性了解测线上的电性分布、基底的起伏、断层的分布、电性层的划分等断面特征。一般而言，在深部(低频)高视电阻率等值线的起伏形态与基底起伏相对应，而视电阻率等值线密集、扭曲和畸变的地方又往往与断层有关，断层越浅，这种特征越明显。

4.3井震约束反演电阻率剖面地质解释

图6给出了L5线在测井曲线标定下地震解释剖面，图中看到基本的分层和对应的断层较为清晰确定，为广域电磁数据反演电阻率分布奠定了基础。图7给出了L5线广域电磁视电阻率数据在地震解释剖面约束下得到的反演电阻率剖面，图中看到反演电阻率也较好地反映了三个断层的存在。

图5　L5线频率—视电阻率拟断面Fig.5　The pseudo section of the apparent resistivity as a function of frequency at line L5

图6　L5线在测井曲线标定的地震解释剖面Fig.6　The seismic interpretation result scaled with well logs of the line 5

图7　L5线广域电磁响在地震解释剖面约束下得到的反演电阻率剖面Fig.7　The inverted resisitivity profile of the wide field electromagnetic responses constrained with well logs and seismic stratum

4.4电阻率反演结果的储层评价应用

图8给出了L3测线广域电磁数据反演电阻率剖面解释断层分布，从图8中可以看到，剖面中分布有2条主要断层分别为F27、F29，推断为逆断层，从地层地表出露情况及倾角判断L3线垂直于构造走向，穿越向斜构造带，整体上呈两侧地层为老地层，自西向东经过向斜核部、向斜右翼。受构造作用影响，地层倾向倾角出现变化，地层之间出现断层并产生错动，地层分布出现变化。

图9 L5测线广域电磁数据反演电阻率剖面解释断层分布，从图9中可以看到，电阻率分布表现出3条断层分别为F17、F13、F19，推断F17为逆断层，F13为逆断层，F19为逆断层。从地层地表出露情况及倾角判断L5线垂直于构造走向，穿越向斜构造带，整体上呈两侧地层为老地层，中间为向斜核部新地层。

经过测井标定的目标层反演结果如图10所示，剖面反演图中可见龙马溪组及牛蹄塘组。龙马溪组埋藏较浅，在测线中完整出现，两侧埋藏浅，中部埋藏深，中部埋藏深度不大于2 000 m；牛蹄塘组埋藏深度从3 000 m到5 000 m。L5龙马溪组较完整，且埋藏不深，牛蹄塘组埋藏相对较深。

图8　L3测线广域电磁数据反演电阻率剖面解释断层分布Fig.8　Interpretation results of the faults obtained from inverted resistivity of the wide field electromagnetic measurements at line 3

图9　L5测线广域电磁数据反演电阻率剖面解释断层分布Fig.9　Interpretation results of the faults obtained from inverted resistivity of the wide field electromagnetic measurements at line 5.

图10　L5测线经过测井标定的广域电磁数据目标层反演结果Fig.10　Inversion resistivity of the target formation scaled with well logs and obtained from wide field electromagnetic measurements at line 5

5　结论与认识

文中采用广域电磁法在湖南保靖地区的测量数据开展了页岩气储层电阻率的约束反演和评价研究，初步证实了文中约束反演算法的有效性，明确了保靖地区地下储层的断层构造和碳质泥页岩的电性分布规律。采用井震约束的广域电磁法数据反演方法，得到了地下地质体的电阻率分布特征，较好地反映了研究区的构造和优质页岩地层分布。从广域电磁法测量数据的电阻率反演结果的评价分析，可以得到如下认识：

1)广域电磁法对地下地层电性参数的分辨率较高，通过井震约束的广域电磁数据电阻率反演结果可以较好地反映研究区的地层分布格架、断层展布特征和页岩储层的纵横向分布；

2)从L5线的反演结果还可以看到，研究区采用的E-Ex广域电磁法对地下高低阻异常体具有较好的分辨能力，为井震电阻率反演奠定了物质基础；

3)用地震层位结果建立初始电阻率分布模型，用测井分层电导率作为上下界约束，对广域电磁数据进行井震约束反演，得到研究区页岩气储层的纵横向分布剖面，可以清晰地给出页岩气富集的有利分布位置。

4) 初步应用表明，由于广域电磁法对高阻和低阻体的良好响应和对复杂地表良好的适应性，将在中浅层页岩气勘探中具有良好的推广应用前景。

致谢本文全部野外原始数据由中南大学继善高科有限公司采集，部分成果由国家自然科学基金(编号：41374141)和国家重点基础研究发展计划(973计划)(编号：2013CB228605)联合资助。

[1]何继善.广域电磁法和伪随机信号电法[M].北京:高等教育出版社,2010.

[2]何继善,李帝铨,戴世坤.广域电磁法在湘西北页岩气探测中的应用[J].石油地球物理勘探,2014,49(5):1 006-1 012.

[3]何继善.广域电磁测深法研究[J].中南大学学报(自然科学版),2010,41(3):1 065-1 072.

[4]余云春.广域电磁法一维正反演[D].长沙:中南大学地球科学与信息物理学院,2010.

[5]李帝铨,谢维,程党性.E-Ex广域电磁法三维数值模拟[J].中国有色金属学报,2013,23(9):2 459-2 470.

[6]汤井田,何继善.水平电偶极频率测深中区全视电阻率定义的新方法[J].地球物理学报,1994,37(4):543-552.

[7]底青云,王若.可控源音频大地电磁数据正反演及方法应用[M].北京:科学出版社,2008.

[8]Goldstein M A, Strangway D W. Audio-frequency magnetotellurics with a grounded-electric dipole source[J].Geophysics,1975,40(4):669-683.

[9]汤井田,何继善.可控源音频大地电磁法及其应用[M].长沙:中南大学出版社,2005.

[10]毛先进,鲍光淑.水平电偶源频率域电磁测深全区视电阻率的直接算法[J].中南工业大学学报(自然科学版),1996,27(3):253-256.

[11]Sasaki Y, Meju M A. Useful characteristics of shallow and deep marine CSEM responses inferred from 3D finite-difference modeling[J]. Geophysics,2009,74(5):F67-F76.

[12]Sasaki Y. Bathymetric effects and corrections in marine CSEM data[J].Geophysics,2011,76(3):139-146.

[13]Sternberg B K, Washburne J C, Pellerint L. Correction for the static shift in mag-netotellurics using transient electromagnetic soundings[J]. Geophysics,1988,53(11):1 459-1 468.

[14]卢文华,杨亚新,吴信民,等.叠加次数对瞬变电磁探测效果的试验研究[J].工程地球物理学报,2016,13(1):105-108.

[15]Smith J T. Conservative modeling of 3-D electromagnetic fields, PartI: Properties and error analysis[J]. Geophysics,1996,61(5)：1 308-1 318.

[16]Smith J T. Conservative modeling of 3-D electromagnetic fields, Part II: Biconjugate gradient solution and an accelerator[J]. Geophysics,1996,61(5):1 319-1 324.

[17]Newman G A, Alumbaugh D L. Three-dimensional massively parallel electromagnetic inversion: 1-Theory[J]. Geophysical Journal International,1997,62(12):345-354.

[18]Rodi W, Mackie R L. Nonlinear conjugate gradients algorithm for 2D magnetotelluric inversion[J]. Geophysics, 2001,66(1):174-187.

[19]Spies B R, Habashy T M. Sensitivity analysis of crosswell electromagnetics[J]. Geophysics,1995,60(3): 834-845.

[20]Hu Wenyi, Abubakar A, Habashy T M. Joint electromagnetic and seismic inversion using structural constraints[J]. Geophysics,2009,74(6):R99-R109.

[21]Abubakar A, Habashy T M, Druskin V L. 2.5D forward and inversion modeling for interpreting low-frequency electromagnetic measurements[J]. Geophysics,2008,73(4):F165-F177.

[22]唐新功,胡文宝,苏朱刘,等.我国南方地区页岩气电磁勘探技术初探[J].长江大学学报(自然科学版)，2011,8(11):33-36.

[23]许杰,何治亮,董宁,等.含气页岩有机碳含量地球物理预测[J].石油地球物理勘探,2013,48(增1):64-68.

[24]汪忠浩, 陈嗣, 李厚霖,等.泥页岩气层脆性特征的地球物理测井研究方法[J].工程地球物理学报,2016,13(1):7-13.

[25]刘振武,撒利明,杨晓,等.页岩气勘探开发对地球物理技术的需求[J].石油地球物理勘探,2011,46(5):810-818.

The Inversion of the Wide Field Electromagnetic Measurements with Well Logging and Seismic Constraints and Its Application to Evaluation of Shale Gas Reservoir in Baojing Area

Wang Yanzhong1，Shen Jinsong2,3，Zhou Zhengwu1，Fu Chao1，Lu Xiaoyao4

(1.ShenhuaGeologicalExplorationLimitedCoporation,Beijing100083,China;2.CollegeofGeophysicsandInformationEngineering,ChinaUniversityofPetroleum,Beijing102249,China;.CNPCKeyLabofGeophysicalExploration,CNPC,Beijing102249,China;4.SchoolofGeophysicsandResources,YangtzeUniversity,WuhanHubei434023,China)

The objective of the paper is to study the inversion of the resistivity distribution of subsurface formation from wide field electromagnetic measurements based on the well logging and seismic constraints, and the inversion results have been applied to the gas reservoir assessment in Baojing area, Hunan province. And the framework of the formation distribution, the extent of the fault, and the shale formation have been obtained from the interpretation of the resisitivity inversion result. And the method applied in the paper includes three aspects. First, the seismic data were scaled and depth matching through well logging data, and the seismic sequences and faults were characterized with geological interpretation tool of the seismic processing system. Next，the inertial resistivity model was established based on the seismic sequences interpretation results and the wide field electromagnetic measurements were inverted for the formation resisitivity distribution under the upper and lower bounds of well logs. Finally, the resistivity distribution was interpreted geologically and provided the vertical and lateral distribution of the shale reservoir. By comparing with the gas test results and well logging data of 3 boreholes, it is validated that inversion method of the wide field electromagnetic measurements with well logging and seismic constraints is effective and its application to the evaluation of the shale gas reservoir in Baojing area is sound.

constraints with well logs and seismic data; wide field electromagnetic measurement; electromagnetic sensitivity calculation; regularized inversion; shale gas reservoir

1672—7940(2016)03—0253—10

10.3969/j.issn.1672-7940.2016.03.001

国家自然科学基金(编号：41374141)；国家重点基础研究发展计划(973计划)(编号：2013CB228605)

王延忠(1969-)，男，博士，教授级高工，主要从事页岩气勘探开发研究工作。E-mail:geophyexplorer@163.com

沈金松(1964-)，男，博士，教授，主要从事电磁探测理论、方法和应用研究。E-mail：shenjinsongcup@163.com

P631.3

A

2015-12-30