一种快速计算电磁场扰动的网格规划方法

1.本发明涉及电磁法勘探技术领域，具体涉及一种快速计算电磁场扰动的网格规划方法。


背景技术：

2.页岩水力压裂、煤层气压裂、油气运移等过程中，储集层会发生电导率变化。电磁法通过观测地下储集层的电导率变化导致的电磁场扰动，有利于准确描述地下储集层的裂隙与压裂液空间分布，为后续施工与井网设计提供依据，有效减少开采成本。
3.传统电磁法数值模拟需要用单元网格对整个地质勘探区域进行空间离散，分别计算地下电导率变化前后的电磁场，通过计算地下电导率变化后的电磁场与地下电导率变化前的电磁场的差值，即可得到地下电导率变化导致的电磁扰动场。但是在对大面积的地质勘探区域进行离散时，往往形成的离散单元数量大、待求解方程组规模大，带来较大的计算负担。
4.经研究发现，水力压裂等过程往往只改变地下电导率变化的储集层周边区域的电导率，这意味着地下其他区域的电导率并不会变化；同时，地下的电磁场变化的集中区域也是处于储集层附近，远离储集层的区域几乎不受到电磁场扰动的影响。
5.利用水力压裂等过程的上述特点，提出一种快速计算电磁场扰动的网格规划方法以解决现有技术中存在计算负担大的问题。


技术实现要素：

6.本发明目的在于提供一种快速计算电磁场扰动的网格规划方法，旨在解决现有电磁法数值模拟计算电磁扰动场存在计算负担大的问题，具体技术方案如下：
7.一种快速计算电磁场扰动的网格规划方法，具体步骤如下：
8.步骤s1：圈定正常网格的计算边界；
9.步骤s2：根据电磁监测前获得的地质信息以及正常网格的计算边界，对计算边界内部的计算区域进行基于有限单元法的网格离散；
10.步骤s3：基于正常网格的框架，对包含储集层在内的区域生成局部网格；
11.步骤s4：设计沿水平方向压裂的储集层模型，赋予水平压裂前的储集层的电导率，以及水平压裂液的电导率；
12.步骤s5：基于电场双旋度方程，采用有限单元总场法在正常网格上计算未压裂前的初始场；
13.步骤s6：依据储集层的地理位置，将正常网格中储集层区域的初始场承继至局部网格中的储集层区域；
14.步骤s7：在局部网格上基于二次场法求解储集层区域的电磁扰动场。
15.以上技术方案中优选的，步骤s1中，正常网格的计算边界所圈定的计算区域包括场源、测点和地质勘探区域；并且根据第一类边界条件要求正常网格上所模拟的场在计算
边界处的幅值衰减为零。
16.以上技术方案中优选的，多个测点在地表上布设为二维的接受阵列；场源采用broadside发射形式。
17.以上技术方案中优选的，在步骤s3中：基于正常网格的框架，对于储集层区域以及储集层上方的接受阵列区域，在局部网格生成时直接承继正常网格的剖分单元；对于储集层的外部区域，通过稀疏化和计算边界的截断缩减处理生成局部网格。
18.以上技术方案中优选的，对于处在电磁扰动场满足第一类边界条件的边界范围之外的区域，直接进行截断。
19.以上技术方案中优选的，稀疏化过程中，合成的新的网格单元的电导率由相邻网格单元体积平均获取：
[0020][0021]
其中，为新网格单元的电导率，σn为电场采样点相邻网格单元的电导率，vn为电场采样点相邻网格单元的体积，n为稀疏化过程被合并为一个新网格单元的原正常网格单元的个数。
[0022]
以上技术方案中优选的，步骤s5中采用有限单元总场法在正常网格上计算初始场具体为：
[0023]
电场的双旋度方程为：
[0024][0025]
在(4)式中，e代表电场矢量，σ代表介质的电导率，w代表圆频率，μ0为真空中的磁导率，i为单位虚数，为旋度算子，jc为外加场源；
[0026]
利用加权余量法，将(4)式在每个网格单元中进行离散，得到：
[0027][0028]
在(5)式中，re是单个网格单元的单元加权余量，j、k为网格单元的棱边序号，和分别为第j和第k条边的矢量形函数，为总场法中获得的初始场，为单个网格单元的体积；
[0029]
对于单个网格单元，将(5)式写成矩阵的形式：
[0030]
re＝(ce iwme)ee iwbeꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)，
[0031]
对于(6)式的各矩阵中的元素分别通过(7)-(9)式得到：
[0032][0033]
[0034][0035]
为矩阵ce中的元素，为矩阵me中的元素，为矩阵be中的元素，为矩阵ee中的元素，re为矩阵re中的元素，σe为每个网格单元的电导率值，为与网格单元有关的场源项；
[0036]
令每个网格单元中(6)式的re为0，得到：
[0037]
(ce iwme)ee＝-iwbeꢀꢀꢀꢀꢀꢀꢀꢀ
(11)，
[0038]
依据有限元法中每个网格单元之间的连接关系，将每个网格单元的各部分单元矩阵组装整合成有关整个正常网格的整体矩阵，得到一个大型稀疏方程：
[0039]
(c iwm)e＝-iwb
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(12)，
[0040]
在(12)式中，c为单元矩阵ce进行组装后得到的整体矩阵，m为单元矩阵me进行组装后得到的整体矩阵，e为单元矩阵ee进行组装后得到的整体矩阵，b为单元矩阵be进行组装后得到的整体矩阵。
[0041]
最后，依据在储集层压裂前的地质模型，采用总场法根据(12)式求解每个单元网格的初始场。
[0042]
以上技术方案中优选的，场源在网格单元上包括两种加载方式，一种是场源直接落在网格单元的面上；另一种是场源落在网格单元的棱边上。
[0043]
以上技术方案中优选的，根据储集层的地理位置，分别索引其在正常网格和局部网格上的网格单元号，并将其一一对应，实现局部网格中储集层区域的初始场的直接承继。
[0044]
以上技术方案中优选的，步骤s7中在局部网格上基于二次场法求解储集层区域的电磁扰动场具体为：
[0045]
将电磁扰动场的双旋度方程与电场的双旋度方程相减得到：
[0046][0047]
其中，es为电磁扰动场，e
p
为初始场，σa为压裂后的地质模型与未压裂前的地质模型间的电导率差值，σ代表介质的电导率，w代表圆频率，μ0为真空中的磁导率，i为单位虚数，为旋度算子；
[0048]
利用加权余量法，将(15)式在每个网格单元中进行离散，得到：
[0049][0050]
其中，是单个网格单元的单元加权余量，j，k为网格单元的棱边序号，和分别为第j和第k条边的矢量形函数；为局部网格中的网格单元承继所获得的初始场；为二次场法获得的电磁扰动场，为单个网格单元的体积；
[0051]
对于单个网格单元，将(16)式写成矩阵的形式：
[0052]
[0053]
对于(16)式中各矩阵中的元素分别通过(18)-(20)式得到：
[0054][0055][0056][0057]
为矩阵中的元素，为矩阵中的元素，为矩阵中的元素，为矩阵中的元素，为矩阵中的元素，σe为每个网格单元的电导率值；
[0058]
令每个网格单元中(17)式的为0，得到：
[0059][0060]
依据有限元法中每个网格单元之间的连接关系，将每个网格单元的各部分单元矩阵组装整合成有关整个局部网格的整体矩阵，得到一个大型稀疏方程：
[0061]
(cs iwms)es＝-iwbsꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(22)，
[0062]
其中，cs为单元矩阵进行组装后得到的整体矩阵，ms为单元矩阵进行组装后得到的整体矩阵，es为单元矩阵进行组装后得到的整体矩阵，bs为单元矩阵进行组装后得到的整体矩阵；
[0063]
根据式(22)求解电导率变化导致的储集层区域的电磁扰动场。
[0064]
应用本发明的技术方案，具有以下有益效果：
[0065]
相对于正常网格，局部网格的数量得到了明显的减少。电磁监测往往需要连续性地多次对电磁扰动场进行计算，本发明使用局部网格来进行电磁扰动场的计算，避免了多次使用一套稠密的正常网格进行数值模拟所带来的庞大的压力，仅需在正常网格上进行一次初始场的计算，后续电磁扰动场的求取都在小规模的局部网格上进行。这种基于局部网格的数值模拟方式，可以大幅度提高电磁扰动场的计算效率，为水力压裂等电磁法监测数据实时解释提供依据。
[0066]
本发明中，初始场计算未采用常规的有限单元二次场法的处理方式(即导出所研究的地质勘探区域的层状背景模型，再基于汉克尔变换等方法求取层状背景模型的解析解作为初始场)，而是采用适应性更广泛的总场法，解决了有限单元二次场法对于含起伏地形的复杂地质模型难以寻求合适的层状背景模型的问题，对含起伏地形的复杂地质模型的初始场实现更为精确的计算。
[0067]
相对于初始场，电磁扰动场(即由储集层的电导率变化导致的二次场)具有衰减快，幅值小的特征。局部网格是基于电磁扰动场的特征来构建的，在采用第一类边界条件时，场的幅值在边界上需衰减至零，而电磁扰动场相对于总场(即由场源直接引起的一次场)的幅值要小得多，故局部网格的计算区域相对于依据总场特征所构建的正常网格可以实现一定程度的缩减；更重要的是，局部网格生成时对储集层外部的网格进行稀疏化，这些
被稀疏化的网格区域在监测过程中电导率不发生变化，具有对总场影响较大，而对电磁扰动场影响较小的特征。据此，局部网格在电磁监测的过程中可以较大程度地实现网格数量上的缩减，且无需像正常网格考虑除储集层外的其他地质体的剖分需求，实现更简便更高效的网格构建。
[0068]
除了上面所描述的目的、特征和优点之外，本发明还有其它的目的、特征和优点。下面将参照图，对本发明作进一步详细的说明。
附图说明
[0069]
构成本技术的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：
[0070]
图1是快速计算电磁场扰动的网格规划方法的流程图；
[0071]
图2是场源发射形式与测点组成的接受阵列的示意图；
[0072]
图3是受水平压裂后的储集层模型示意图；
[0073]
图4是有限单元二次场法对地质模型的背景模型进行层状近似的示意图。
具体实施方式
[0074]
为了便于理解本发明，下面将对本发明进行更全面的描述，并给出了本发明的较佳实施例。但是，本发明可以以许多不同的形式来实现，并不限于本文所描述的实施例。相反地，提供这些实施例的目的是使对本发明的公开内容的理解更加透彻全面。
[0075]
除非另有定义，本文所使用的所有的技术和科学术语与属于本发明的技术领域的技术人员通常理解的含义相同。本文中在本发明的说明书中所使用的术语只是为了描述具体的实施例的目的，不是旨在于限制本发明。
[0076]
实施例1：
[0077]
参见图1，一种快速计算电磁场扰动的网格规划方法，具体是一种支持快速计算地下电导率变化导致电磁场扰动的网格规划方法。本实施例中先对计算区域进行正常网格划分，其中计算区域包括储集层；并对储集层外部的区域进行稀疏化和计算边界的截断缩减处理，获得包含储集层在内的局部网格；然后通过有限单元总场法在正常网格上计算电导率变化前计算区域的初始场，并将正常网格中储集层区域的初始场承继至局部网格中的储集层区域；最后通过有限单元二次场法在局部网格上计算电导率变化导致的储集层区域的电磁扰动场。
[0078]
进一步的，所述方法具体包括以下步骤：
[0079]
步骤s1：圈定正常网格的计算边界；
[0080]
具体地，正常网格的计算边界所圈定的计算区域必须包括场源、测点和地质勘探区域；并且根据第一类边界条件要求正常网格上所模拟的场在计算边界处的幅值衰减为零，可表示为：其中，e表示正常网格上所模拟的场的幅值，表示计算边界。
[0081]
步骤s2：根据电磁监测前获得的地质信息以及正常网格的计算边界，对计算边界内部的计算区域进行基于有限单元法的网格离散；
[0082]
具体地，在电磁监测前，往往对地质勘探区域进行了大量的前期勘探工作，依据勘探所获得的地质信息(包括异常体的位置、分解面的位置、地下的背景电导率分布情况、储
集层的地质信息等)，为相应的地质体开辟相应的网格单元，并赋予相应的电导率，并使得网格单元的组合形态与真实的地质体的形态相拟合。
[0083]
进一步地，在网格剖分时，需考虑相邻网格单元之间的拉伸系数(即相邻网格单元的长或宽的比率)，使得拉伸系数的值介于1～1.5之间，可保证网格剖分的整体质量，从而提高数值模拟的精度。为进一步保证数值模拟的精度，我们需要对场源、测点、异常体、储集层、电导率反差较大的位置进行精细的剖分，一般来说，进行精细剖分时选用的网格间隔为20m。
[0084]
本实施例中，场源采用broadside发射形式，为更全面地对电磁扰动场进行描述，将测点布设为二维的接受阵列，接受阵列由水平间距为20m的测点组成，布设为地表上，如图2所示。
[0085]
步骤s3：基于正常网格的框架，对包含储集层在内的区域生成局部网格；
[0086]
具体地，基于正常网格的框架，对于储集层区域以及储集层上方的接受阵列区域，在局部网格生成时直接承继正常网格的剖分单元；对于储集层的外部区域，通过稀疏化和计算边界的截断缩减处理生成局部网格。
[0087]
对于储集层的外部区域，忽视外部区域的其他异常体的剖分需求，即假定异常体区域不存在异常体，依据此异常体与储集层的距离，逐次增加网格间隔，无需对异常体进行精细剖分，可实现对正常网格的稀疏化。
[0088]
在局部网格的稀疏化步骤中，会出现正常网格的多个单元因稀疏化而融合为一个单元的情况，在稀疏化过程中存在上下以及水平相邻网格单元合并为一个新的网格单元的过程，合成的新的网格单元的电导率由相邻网格单元体积平均获取，具体公式为：
[0089][0090]
其中，为新网格单元的电导率，σn为电场采样点相邻网格单元的电导率，vn为电场采样点相邻网格单元的体积，n为稀疏化过程被合并为一个新网格单元的原正常网格单元个数；例如，如果两个正常网格单元合并成一个新网格单元，则n为2。
[0091]
对于远离储集层的区域，这部分区域处在电磁扰动场满足第一类边界条件的边界范围外，则直接进行截断，即实现计算边界的缩减。
[0092]
步骤s3能够实现网格数量缩减的原因在于：(1)异常体对总场的影响大，故在正常网格中需要对异常体进行精细剖分来准确地传导异常体对总场的影响；而在电磁监测中，异常体的电导率未发生变化，电磁扰动场由发生电导率变化的储集层产生，异常体对电磁扰动场的影响很小，故在计算电磁扰动场的局部网格上，无需对异常体进行刻画；(2)有限单元法的网格剖分(即网格离散)中，网格质量是保证数值模拟精度的前提，而相邻网格单元的变化程度越大，网格质量就越差，故相同的长度，若多处需要精细剖分，则无需精细剖分的区域由于要与精细剖分的区域保持较小的网格变化程度，故其网格间隔也受到限制；局部网格无需考虑场源和储集层的剖分需求，局部网格的数量相对于正常网格的数量实现了大幅度的缩减。
[0093]
步骤s4：设计沿水平方向压裂的储集层模型，赋予水平压裂前的储集层的电导率，
以及水平压裂液的电导率。
[0094]
步骤s5：基于电场双旋度方程，采用有限单元总场法在正常网格上计算未压裂前的初始场。
[0095]
频率域中，在准静态近似条件下，选用时谐因子e-iwt
基于麦克斯韦方程组可以得到以下关系：
[0096][0097][0098]
e和h分别代表电场矢量和磁场矢量，σ代表介质的电导率，w代表圆频率，μ0为真空中的磁导率，i为单位虚数，为旋度算子，jc为外加场源。
[0099]
为了避免同时将电场和磁场代入有限元方程增加求解的复杂性，将(2)式两边求取旋度后代入(3)式得到电场的双旋度方程：
[0100][0101]
利用加权余量法，将(4)式在每个网格单元中进行离散，得到：
[0102][0103]
在(5)式中，re是单个网格单元的单元加权余量，j,k为网格单元的棱边序号，和分别为第j和第k条边的矢量形函数，为总场法中获得的初始场(即总场法需要求解的初始场)，为单个网格单元的体积。
[0104]
对于单个网格单元，将(5)式写成矩阵的形式：
[0105]
re＝(ce iwme)ee iwbeꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)，
[0106]
对于(6)式的各矩阵中的每一个元素(例如，表示为ce中第j行，第k列的元素)，分别通过(7)-(9)式得到：
[0107][0108][0109][0110]
对应为矩阵ce中的元素，对应为矩阵me中的元素，对应为矩阵be中的元素，对应为矩阵ee中的元素，re对应为矩阵re中的元素，σe为每个网格单元的电导率值，为与网格单元有关的场源项。
[0111]
进一步的，场源在网格单元上存在两种加载方式，(a)场源直接落在网格单元的面
上；(b)场源落在网格单元的棱边上。本实施例中采取场源落在网格单元棱边上的加载方式，可以将(9)式写成：
[0112][0113]
其中，rate是单元体积的加权系数项，idl为场源在单元上的电偶极矩。
[0114]
令每个网格单元中(6)式的re为0，可以得到：
[0115]
(ce iwme)ee＝-iwbeꢀꢀꢀꢀꢀꢀꢀꢀ
(11)，
[0116]
依据有限元法中每个网格单元之间的连接关系，将每个网格单元的各部分单元矩阵组装整合成有关整个正常网格的整体矩阵，得到一个大型稀疏方程：
[0117]
(c iwm)e＝-iwb
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(12)，
[0118]
在(12)式中，c为单元矩阵ce进行组装后得到的整体矩阵，m为单元矩阵me进行组装后得到的整体矩阵，e为单元矩阵ee进行组装后得到的整体矩阵，b为单元矩阵be进行组装后得到的整体矩阵。
[0119]
最后，依据在储集层压裂前的地质模型，采用总场法根据(12)式求解每个单元网格的初始场。
[0120]
步骤s6：依据储集层的地理位置，将正常网格中储集层区域的初始场承继至局部网格中的储集层区域。
[0121]
具体的，由于在电磁监测中仅存在储集层的电导率发生变化，除储集层区域，其他区域在压裂过程中没有发生电导率变化，故我们仅需对储集层的初始场进行考虑。由于正常网格与局部网格在储集层处的网格剖分情况是相同的，故我们可以根据储集层的地理位置，分别索引其在正常网格和局部网格上的网格单元号，并将其一一对应，即可实现局部网格中储集层区域的初始场的直接承继。
[0122]
步骤s7：在局部网格上基于二次场法求解储集层区域的电磁扰动场。
[0123]
有限单元二次场法是将电场分解为初始场e
p
和电磁扰动场es：
[0124]
e＝e
p
esꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(13)，
[0125]
传统的有限单元二次场法的初始场是由场源下的均匀半空间或水平层状介质模型导出，而本实施例的(13)式中的初始场e
p
是通过总场法在正常网格上计算得到的(即承继得到)。
[0126]
图4展示了起伏地形下的电测监测模型，椭圆形圈闭体为电导率发生变化的储集层区域，由于传统二次场法需要寻求一个层状背景模型，故需要对起伏的地层进行层状近似，图4中的虚线为近似后的地层分界线，这样会导致初始场的计算存在误差，而用总场法直接对场源进行加载，无需寻求层状背景模型，适应性更加广泛，解决了复杂地质模型缺少合适的层状背景模型的问题，能够导出更加精确的初始场。
[0127]
进一步的，电磁扰动场是由电导率变化的储集层所产生，建立场源与初始场的关系(即建立电磁扰动场的双旋度方程)：
[0128][0129]
式中，σ
p
为未压裂前的地质模型的电导率。
[0130]
将(14)式与(4)式相减得到：
[0131]
[0132]
式中，σa为压裂后的地质模型与未压裂前的地质模型间的电导率差值，称异常电导率。
[0133]
进一步的，利用加权余量法，将(15)式在每个网格单元中进行离散，得到：
[0134][0135]
式中，是单个网格单元的单元加权余量，j，k为单元的棱边序号，和分别为第j和第k条边的矢量形函数；为承继步骤s5所获得的初始场，即由步骤s5中所计算的通过步骤6承继而来；为二次场法获得的电磁扰动场(即二次场法需要求解的电磁扰动场)，为单个网格单元的体积。
[0136]
对于单个网格单元，将(16)式写成矩阵的形式：
[0137][0138]
其中：对于(16)式中各矩阵中的每一个元素(例如，表示为中第j行，第k列的元素)，可以通过(18)-(20)式得到：
[0139][0140][0141][0142]
对应为矩阵中的元素，对应为矩阵中的元素，对应为矩阵中的元素，对应为矩阵中的元素，对应为矩阵中的元素，σe为每个网格单元的电导率值。
[0143]
令每个网格单元中(17)式的为0，可以得到：
[0144][0145]
依据有限元法中每个网格单元之间的连接关系，将每个网格单元的各部分单元矩阵组装整合成有关整个局部网格的整体矩阵，得到一个大型稀疏方程：
[0146]
(cs iwms)es＝-iwbsꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(22)，
[0147]
在(22)式中，cs为单元矩阵进行组装后得到的整体矩阵，ms为单元矩阵进行组装后得到的整体矩阵，es为单元矩阵进行组装后得到的整体矩阵，bs为单元矩阵进行组装后得到的整体矩阵。
[0148]
根据式(22)求解电导率变化导致的储集层区域的电磁扰动场。如图3，依据所设置
的储集层模型在x方向上的压裂距离，对局部网格上的储集层电导率参数做出改变，即将已进行过水力压裂的压裂区的电导率设置为压裂液的电导率，将储集层未被压裂的区域设置为原储集层的电导率，依据(22)式，对储集层因受到压裂产生电导率变化而产生的电磁扰动场进行计算。
[0149]
本实施例中基于局部网格的电磁扰动场算法，仅需在正常网格上计算初始场，电磁扰动场则是在小规模的局部网格上计算，所形成的方程组规模将远小于只使用一套正常网格进行计算的传统数值模拟方法，极大提升地下电导率变化导致电磁场扰动的计算效率。
[0150]
以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。