2018, Vol. 61
2. 防灾科技学院电子科学与控制工程学院, 河北三河 065201;
3. 中国地质大学(北京)地球物理与信息技术学院, 北京 100083
2. School of Electronic Science and Control Engineering, Institute Disaster of Prevention Science and Technology, Hebei Sanhe 065201, China;
3. School of Geophysics and Information Technology, China University of Geosciences, Beijing 100083, China
电阻率法是一种应用于矿产、能源、环境、灾害、水文和地质构造等领域的经典地球物理方法(Yuval and Oldenburg, 1996；Hauck，2002；Diaferia et al., 2006；Krautblatter and Hauck, 2007；汤井田等，2007；Jones et al., 2010；Revil et al., 2010；Barde-Cabusson et al., 2013；李术才等，2015；范国栋，2017).电阻率法灵活性高，拥有多种观测方式，如二极、三极、偶极-偶极、对称四极等数据采集装置，而且其供电电极和接收电极还可以与钻井搭配使用，如地-井、井-地、井-井采集装置.但是由于在地表获取的数据中往往含有浅地表局部异常体的干扰信息，导致更深部的目标体反演结果不理想(Li and Oldenburg, 1992，2000；Ellis and Oldenburg, 1994；Oldenburg and Li, 1994).因此，减小浅地表局部异常体对目标体反演结果精度的影响显得十分重要.
为了精准的恢复目标体电阻率结构，学者们研究了各种反演算法.混合非线性最小二乘反演对初始模型要求比较严格，比较重构法反演虽然成像质量较高，但是要求高精度的观测数据(白登海和于晟，1995)；共轭梯度反演算法避免求解灵敏度矩阵，为三维反演节约大量时间(Zhang et al., 1995；吴小平和徐果明，2000)；非线性多目标遗传反演算法无需计算灵敏度矩阵，且对初始模型的依赖性较小(Schwarzbach et al., 2005)；二维数据空间Occam反演能够有效的降低反演矩阵维数，节省存储空间(Boonchaisuk et al., 2008)；二维电阻率快速优化的最小复杂神经网络反演算法，减小人工神经网络模型复杂程度(Singh et al., 2010).除此之外，还有岭回归的最小二乘法、α中心法、近似法、非线性共轭梯度法等反演算法(Pelton et al., 1978；Petrick et al., 1981；Park and Van, 1991；Ma et al., 2015).虽然以上算法都精准的恢复出目标体电阻率，但是没有考虑浅地表局部异常体对目标体反演结果精度的影响.
在地表实施电阻率法勘探过程中，浅地表局部电性异常体会影响电流场分布.当浅地表存在局部低阻异常体时，低阻体吸引电流，影响电场分布，进而使得实测视电阻率数据中大多是浅地表局部低阻异常体信息，而仅包含少量地下勘探目标信息，其预期勘探深度大大减小.为此，不得不增加极距来获得更大勘探深度的信息，从而导致需要更大功率的发电机增加电流强度，耗费更多人力物力.相对地表电阻率法而言，浅地表局部低阻异常体对井中电阻率法数据信息影响较小，但是由于钻井位置的局限性，井中装置采集的数据信息覆盖范围不广.为此，在不增加极距和电流的前提下，本文提出将地表装置和井中装置采集的视电阻率数据联合反演的策略，以达到减小浅地表局部异常体对电场分布影响的目的，而且该方法可以获得更多勘探目标数据信息，从而减小反演多解性，更准确地恢复勘探目标的电阻率结构.
另外，这种联合反演的策略虽然减小了反演多解性，但是随着数据量增加，会使得反演计算效率降低.近年来，计算机并行计算在地球物理数值模拟领域发展迅速，很多国内外地球物理学者通过利用MPI并行计算提高算法执行效率：在大地电磁正演模拟算法中，可以利用MPI将多个频率的计算任务分给多个进程同时进行计算，或者将MPI并行用于QMR解方程中系数矩阵和向量的乘积运算(Tan et al., 2006；李焱等，2010；Weiss，2013)；也可以使用区域分解和MPI并行技术，实现地震预测中的LURR逼近，加速比可达到20~40倍(Feng et al., 2008)；在可控源电磁法(CSEM)正反算法中，采用MPI并行算法加速计算灵敏度矩阵的转置和向量的乘积(Kelbert et al., 2014)；利用粗细粒度结合的MPI多通讯域并行算法实现地震波动方程正演模拟(李祎昕，2014).所以，针对上述联合反演数据量大的特点，可以利用MPI将电阻率法多个点电源正演模拟任务和灵敏度矩阵的相关计算任务分给多个进程同时进行计算，从而达到提高程序执行效率的目的.
为了减小浅地表局部异常体对目标体反演结果精度的影响，本文开发了一套联合地表、地-井、井-地和井-井装置数据体的三维正反演算法.由于这种方法使得数据量大幅增加，导致本来就比较复杂的灵敏度矩阵变得更为庞大，求解更困难，耗时更长.为此，选择NLCG反演算法避免计算灵敏度矩阵，同时在NLCG反演算法的基础上利用MPI并行算法加速正演算法和灵敏度矩阵的相关计算，最后通过理论合成数据和实测数据论证该正反演算法的有效性，并讨论MPI并行算法的执行效率.
1 理论部分 1.1 正演算法点电源电阻率法满足的偏微分方程及其边界条件为
|
(1) |
其中，σ是岩矿石电导率(S·m-1)；φ是电位向量(V)；I是供电电流强度(A·m-3)；δ是狄拉克函数；n是边界点的向外法向量；r是点电源与边界点的距离；θ是边界点向外法向量n和场源到边界点距离r之间的夹角.
该偏微分方程的边界条件分为：1)在地表空气边界采用第二类边界条件(Newman边界条件)；2)在无穷远边界处采用第三类边界条件(混合边界条件).本文采用七点格式有限差分法离散方程(1)，形成大型线性方程组，其矩阵形式表示为
|
(2) |
其中，A是系数矩阵；b是含有场源信息的向量.本文采用基于不完全LU分解的稳定双共轭梯度法(BICGSTAB)解方程获得电位(vander Vorst，1992；熊杰等，2012)，求得电位φ后即可得到视电阻率ρa，公式为
|
(3) |
其中，α是在所有网格节点电位φ中选择接收电极位置处电位的向量；U是二极装置M点采集的电位；ΔU是三极或四极装置MN的电位差；K为装置系数(m).
由于本文介绍的数据类型涉及井中数据，那么其装置系数也会相应的变化，通过镜像法(施国良和张国雄，2008)推导的点电源装置系数通用公式为
|
(4) |
其中，AM、AN、BM、BN分别是场源A、B至接收点M、N的距离；AM′、AN′、BM′、BN′分别是虚场源A′、B′至接收点M、N的距离.下面分别说明该装置系数公式的三种情况：
1) 当使用二极或三极装置采集时，B放置在无穷远，则BM=BM′=BN=BN′= ∞；2)当场源在地表时，虚场源和场源重合，所以AM=AM′，BM=BM′，AN=AN′，BN=BN′；3)当场源在井中时，场源和虚场源不重合，直接将各个量代入公式(4)即可.所以，无论是地表装置还是井中装置(地-井、地-井和井-井)视电阻率数据，都可以直接通过公式(4)计算的装置系数代入公式(3)中获得.
1.2 反演算法本文采用NLCG算法反演三维电阻率结构，其目标函数ϕ为
|
(5) |
其中，dobs是观测视电阻率向量；d是正演视电阻率向量；Cd是随机误差矩阵；λ是拉格朗日乘子，在反演迭代过程中用于平衡数据误差和模型正则化的影响；mref是电导率先验模型向量；m是每次迭代所更新的电导率模型向量，在算法中取模型向量的自然对数ln(σ)，即电导率的自然对数；L是拉普拉算子，对空间X、Y、Z三个方向的二阶偏导数(王家映，2002)；上标T表示矩阵的转置.
在电阻率法反演算法中，灵敏度矩阵是电位对电导率模型求偏导数.方程(2)两端对电导率求偏导数，再利用方程(3)即可导出灵敏度矩阵J为
|
(6) |
其中，i和j分别代表第i个数据和第j个网格单元.非线性共轭梯度法的优势是避免求解灵敏度矩阵，直接计算灵敏度矩阵J或者其转置和一个向量的乘积(Jp或JTe)，能够减少占用的计算机内存.可以表示为(Zhang et al., 1995；吴小平和徐果明，2000；吴小平和汪彤彤，2001)：
|
(7) |
|
(8) |
其中，p是反演算法中的搜索方向，e是数据残差向量.若要求解公式(7)和(8)，首先需要求解两个公式中的V向量，这个过程相当于求解正演线性方程组的过程，由于其右端项不同于正演数值模拟，所以称之为“拟正演”.如公式(9)所示，通过BICGSTAB法解线性方程组，求得V向量.再将V代入(7)和(8)中，得到相应结果.公式(9)为
|
(9) |
MPI主要用于计算机进程之间传递消息，适用于多个独立计算任务.由于电阻率法NLCG反演算法中，各场源正演和“拟正演”数值模拟互相独立且耗时，所以本文利用MPI主从模式(主进程分配任务，从进程执行任务并返还结果给主进程)开发一套基于MPI并行的电阻率法三维联合反演算法.下面介绍正演和“拟正演”MPI并行算法.
公式(5)中数据d的正演过程与公式(7)中Jp的“拟正演”计算过程类似，其MPI并行算法如图 1所示.假设有4个场源，每个场源供电时观测10个数据，利用MPI开辟两个进程进行计算，每个进程都分配两个场源的计算任务(进程一计算d1至d20；进程二计算d21至d40)，计算完毕后，主进程使用MPI_AllGatherv函数收集各进程计算结果，整合到向量d中.
|
图 1 MPI_AllGatherv并行算法原理图 Fig. 1 Principle of parallel algorithm MPI_AllGatherv |
公式(8)计算灵敏度矩阵的转置和数据残差向量的乘积，其结果是一个一维向量，向量长度和剖分网格个数相同.然而，当存在多个场源时，其结果需将每个场源的结果相加.假设有4个场源，利用MPI开辟两个进程计算，每个进程都分配两个场源的计算任务(进程一计算前两个场源的JTe；进程二计算计算后两个场源的JTe)，计算完毕后，主进程调用MPI_Reduce函数将两个进程计算结果进行累加，得到最终JTe，算法流程如图 2.最后，将上述MPI并行算法嵌入NLCG反演算法中，形成一套适用于电阻率法多种装置组合数据的三维联合反演MPI并行算法，其流程如图 3所示.
|
图 2 MPI_Reduce并行算法原理图 Fig. 2 Principle of parallel algorithm MPI_Reduce |
|
图 3 电阻率法三维NLCG反演并行算法流程图 Fig. 3 Flow chart of parallel algorithm for 3D NLCG inversion in resistivity method |
为了测试正演算法的精度，将正演算法的电位结果与均匀半空间解析解、预处理共轭梯度(CGPC)正演算法(Spitzer，1995)结果进行对比.假设有一个100 Ωm均匀半空间，数据观测装置分别采用地表二极装置、地-井二极装置、井-地二极装置和井-井二极装置.如图 4所示，圆点是场源，倒三角是接收电极，五角星是井口位置，(-80, 0, 0)是接收井，(80, 0, 0)是场源井.如图 5所示，解析解和正演模拟结果吻合较好，两者的最大相对误差为0.027%(图 6).
|
图 4 模型和观测系统图 (a) XOZ剖面；(b)地表投影. Fig. 4 Model and observation system (a) XOZ profile; (b) Surface projection. |
|
图 5 正演结果与解析解对比图 (a)地表装置电位曲线；(b)地-井装置电位曲线；(c)井-地装置电位曲线；(d)井-井装置电位曲线. Fig. 5 Comparison of forwardmodeling results and analytical solutions (a) Potential curve of surface array; (b) Potential curve of surface-borehole array; (c) Potential curve of borehole-surface array; (d) Potential curve of borehole-borehole array. |
|
图 6 正演结果与解析解相对误差图 (a)地表装置电位相对误差图；(b)地-井装置电位相对误差图；(c)井-地装置电位相对误差图；(d)井-井装置电位相对误差图. Fig. 6 Relative errors between forward modeling results and analytical solutions (a) Relative error of surface potential; (b) Relative error of surface-borehole potential; (c) Relative error of borehole-surface potential; (d) Relative error of borehole-borehole potential. |
随后，在上述均匀半空间中设置一个几何尺寸为100 m×100 m×70 m，埋深为70 m，电阻率为10 Ωm的低阻棱柱体模型，数据观测系统保持不变(图 4).本文的正演模拟结果与Spitzer(1995)的CGPC正演模拟结果一致性较好(图 7)，最大相对误差为0.0727%(图 8)，通过以上结果对比，证明了本文三维正演算法的有效性.
|
图 7 正演结果与CGPC结果对比图 (a)地表装置电位曲线；(b)地-井装置电位曲线；(c)井-地装置电位曲线；(d)井-井装置电位曲线. Fig. 7 Comparison of forward modeling and CGPC results (a) Potential curve of surface array; (b) Potential curve of surface-borehole array; (c) Potential curve of borehole-surface array; (d) Potential curve of borehole-borehole array. |
|
图 8 正演结果与CGPC结果相对误差图 (a)地表装置电位相对误差图；(b)地-井装置电位相对误差图；(c)井-地装置电位相对误差图；(d)井-井装置电位相对误差图. Fig. 8 Relative error between forward modeling and CGPC results (a) Relative error of surface potential; (b) Relative error of surface-borehole potential; (c) Relative error of borehole-surface potential; (d) Relative error of borehole-borehole potential. |
由于在浅地表常常存在电阻率异常干扰，所以在背景为1000 Ωm的均匀半空间中放置五个电阻率异常体(Ellis and Oldenburg, 1994)，其中三个顶界面位于地表的棱柱体(S1，S2，S3)用于模拟浅地表的电阻率异常干扰，另外两个顶界面埋深分别为50 m和95 m的棱柱体(B1，B2)模拟勘探目标体，具体参数见表 1.为了进行效果对比，以下三个算例都采用相同的反演初始参数：整个研究区域的网格剖分为17496(27×27×24)个三维网格块；X、Y、Z三个方向的剖分范围分别是(-1500，1500)、(-1500，1500)和(0，475)；视电阻率数据加入2%的高斯误差参与反演；初始模型选用1000 Ωm的均匀半空间；拉格朗日乘子λ=10-5.
|
表 1 棱柱体物性参数表 Table 1 Physics parameters of prism |
数据采集装置为地表二极装置，点电源供电时，所有接收电极采集视电阻率数据(图 9)，数据量为10690个.图 10a和图 10c是X=±100 m真实模型剖面图，图 10b和图 10d是X=±100 m反演结果剖面图.反演结果中，S2和S3浅地表电阻率异常恢复较好，但是高阻目标体B1和低阻目标体B2的反演结果偏离真实模型中心，边界不清晰，尤其是X= 100 m剖面B1目标体的反演结果.而且X=-100 m剖面反演结果中，B2目标体的电阻率值只恢复至300 Ωm左右.该算例经118次反演迭代，RMS误差为5.98(图 11a)，从三维结果(本文合成数据算例中的三维结果都保留100 Ωm至300 Ωm以及大于1700 Ωm部分)中也能看出(图 11b)，仅使用地表装置数据进行反演，受到浅地表局部异常体干扰，该算例低阻目标体反演结果不好，浅地表存在假异常.
|
图 9 异常模型位置和数据观测系统图 (a)三维异常模型和钻井位置图；(b)地表数据观测系统图. Fig. 9 Location of anomaly model and observation system (a) Location of 3D anomaly model and boreholes; (b) Surface observation system. |
|
图 10 地表二极装置数据YOZ剖面反演结果 (a) X=-100 m剖面模型图；(b) X=-100 m剖面反演结果；(c) X=100 m剖面模型图；(d) X=100 m剖面反演结果. Fig. 10 Inversion results with pole-pole data of surface in YOZ profile (a) model of X=-100 m profile; (b) Inversion result of X=-100 m; (c) model of X=100 m profile; (d) Inversion result of X=100 m. |
|
图 11 地表二极装置数据反演拟合差及三维结果 (a)拟合差曲线；(b)三维结果. Fig. 11 RMS and 3-D results of inversion with pole-pole data of surface (a) RMS curve; (b) 3-D result of inversion. |
添加接收电极于两口接收钻井(图 9黑色虚线)中，利用地-井二极装置采集的视电阻率数据和上一算例的地表视电阻率数据联合反演.在两口井中分别放置11个接收电极，深度分布范围是10~380 m.反演结果如图 12所示，与上述算例的反演结果比较，X=±100 m(图 12b和图 12d)的目标体B1反演结果比算例一有所改善，与真实模型位置一致，但是目标体B2仍然没有恢复至真实电阻率值.该算例经119次反演迭代，RMS误差为10.22(图 13a)，虽然也受到浅地表局部异常体干扰，浅地表存在假异常，但是该算例低阻目标体反演结果优于上一算例的结果(图 13b).
|
图 12 地表和地-井二极装置数据YOZ剖面反演结果 (a) X=-100 m剖面模型图；(b) X=-100 m剖面反演结果；(c) X=100 m剖面模型图；(d) X=100 m剖面反演结果. Fig. 12 Inversion results with pole-pole data of surface and surface-borehole in YOZ profile (a) model of X=-100 m profile; (b) Inversion result of X=-100 m; (c) model of X=100 m profile; (d) Inversion result of X=100 m. |
|
图 13 地表和地-井二极装置数据反演拟合差及三维结果 (a)拟合差曲线；(b)三维结果. Fig. 13 RMS and 3-D results of inversion with surface and surface-borehole data (a) RMS curve; (b) 3-D result of inversion. |
在算例二观测系统的基础上，添加发射源于两口钻井(图 9粉色虚线)中，利用地-井、井-地、井-井二极装置采集的视电阻率数据体和算例一视电阻率数据体联合反演.在两口井中分别放置发射源11个，其对应深度相同，分布范围是10~380 m.从X=±100 m剖面反演结果中可知(图 14b和图 14d)，无论是异常范围边界，还是电阻率值，该算例反演结果都明显优于算例一和算例二，尤其是低阻目标体B2的电阻率值已恢复至接近100 Ωm.该算例经115次反演迭代，RMS误差为7.23(图 15a).如图 15b所示，使用地表、地-井、井-地和井-井装置数据进行联合反演，浅地表局部异常体、高阻和低阻目标体都恢复较好，而且地表不存在假异常.
|
图 14 地表、地-井、井-地和井-井二极装置数据YOZ剖面反演结果 (a) X=-100 m剖面模型图；(b) X=-100 m剖面反演结果；(c) X=100 m剖面模型图；(d) X=100 m剖面反演结果. Fig. 14 Inversion results with the pole-pole data of surface, surface-borehole, borehole-surface and borehole-borehole in YOZ profile (a) model of X=-100 m profile; (b) Inversion result of X=-100 m; (c) Model of X=100 m profile; (d) Inversion result of X=100 m. |
|
图 15 地表、地-井、井-地和井-井二极装置数据反演拟合差及三维结果 (a)拟合差曲线；(b)三维结果. Fig. 15 RMS and 3-D results of inversion with surface, surface-borehole, borehole-surface and borehole-borehole data (a) RMS curve; (b) 3-D result of inversion. |
为了测试三维电阻率法正反演MPI并行算法的执行效率，以上述合成数据体的算例二为例，确保网格剖分、反演初始模型、观测数据误差和拉格朗日乘子不变的情况下，分别用2个进程、6个进程、12个进程和16个进程执行并行程序，并与串行程序执行时间对比.表 2展示了串行程序和并行程序的运行时间和相对加速比，当开辟12个进程执行程序时，获得最大加速比4.51.
|
|
表 2 合成数据算例二反演执行时间表 Table 2 Runtimes for inversion of synthetic data for example 2 |
从表 2中还可以发现以下两个特点：首先，与串行程序的相对加速比不是整数倍关系.这是由两个原因造成的：1)在并行程序运行过程中，进程间数据传递会占用时间；2)场源总数与进程数不是整数倍关系，在计算任务少的进程完成任务后，会等待任务多的进程，从而导致时间消耗.其次，开辟16个进程比开辟12个进程的程序执行时间相差不多，而且16个进程的执行时间长，加速比相对较小.造成这种情况主要有两个原因：1)12进程和16进程的最大任务量(两个场源计算任务)和最小任务量(一个场源计算任务)相同，这就导致了计算任务少的进程有近似的等待时间；2)开辟16个进程进行计算，其分配的等待进程数为10个，而开辟12个进程的等待进程数为2个，相对而言，开辟16个进程进行计算，不仅增加了进程间的通讯时间，而且增加了更多的等待进程数量.
3.3 实测数据算例测区位于某沙漠边缘区域，为了勘查污染物进入地下水情况，在二维观测剖面上使用地表三极装置进行勘查，将其视电阻率资料进行三维反演.三极装置的无穷远极(B极)放在测区中心5 km以外，工区共273个点电源(A极)，MN为10 m，当一个场源供电，则该场源所在的测线进行数据采集.如图 16所示，星号代表点电源，圆圈代表接收电极，东西方向5条测线，南北方向1条测线.
|
图 16 观测系统图 Fig. 16 Observation system |
将所有观测数据代入本文反演算法中，正反演网格剖分大小为97×97×44，拉格朗日因子选取100，初始模型为40 Ωm的均匀半空间，由于工区地形起伏较小(不到2 m)，故没有考虑地形影响.图 17为测线Y=160 m，场源位置为(X=95 m，Y=160 m)的数据拟合曲线，正演模拟结果和实测数据拟合较好，证明了反演结果的可靠性.拟合差曲线表明反演过程稳定下降，共迭代89次，拟合差RMS下降至2.66(图 18a).在三维结果(图 18b)中发现了几处局部浅地表低阻异常，这可能是由于数据噪声引起的假异常.根据水文钻井资料，证实低电阻率区域为地下水.从X=160 m、Y=160 m和Y=210 m三条剖面的反演结果可知(图 19)，该工区浅层为不连续的高阻覆盖，电阻率在260 Ωm至1000 Ωm之间，不存在从地表入侵地下水的低阻污染物，说明后期治理效果较好，在20 m至40 m深度位置为水平低阻层，电阻率在0.1 Ωm至1 Ωm之间，这是由于前期污染，地下水中含有较多金属离子，使得反演电阻率低至0.1 Ωm.
|
图 17 数据拟合曲线 Fig. 17 Data fitting curve |
|
图 18 实测数据反演拟合差及三维结果 (a)拟合差曲线；(b)三维结果. Fig. 18 RMS and 3-D results of inversion with measured data (a) RMS curve; (b) 3-D result of inversion. |
|
图 19 实测数据反演结果图 (a) X=160 m剖面反演结果；(b) Y=160 m剖面反演结果；(c) Y=210 m剖面反演结果. Fig. 19 Inversion results of measured data (a) Inversion result of X=160 m profile; (b) Inversion result of Y=160 m profile; (c) Inversion result of Y=210 m profile. |
本文以有限差分法正演和NLCG反演为理论基础，利用MPI并行技术，开发了电阻率法多种装置数据的三维联合反演算法.合成数据反演算例说明了多种装置数据体联合反演，不仅有利于恢复异常目标体的电阻率值、位置和边界范围，而且还可以减弱浅地表局部电阻率异常体的影响，实测数据算例更进一步证实了本文反演算法的有效性.MPI并行算法测试表明：并不是开辟进程数越多，执行时间越短，还要根据电阻率法场源数量具体分析合适的进程数量，才能发挥MPI并行计算效果.由此可见，像这样MPI开辟进程数量的问题并不只存在于电阻率法中，在电磁法和人工地震中也常常遇见，比如频率域电磁法中可以按照频率数量分配进程，人工地震中可以按照震源数量分配进程.本文的创新点有两点：1)在正演算法中，整理出了地表、地-井、井-地和井-井装置的装置系数通用计算公式，反演算法不仅可以反演单一装置电阻率法的视电阻率资料，还可以反演多种装置组合的视电阻率资料；2)利用MPI并行算法加速正演模拟和灵敏度矩阵或其转置与一个向量乘积的运算，提高反演效率.当开辟12个进程时，获得最大加速比4.51.在今后的研究过程中，希望可以获取三维地-井、井-地和井-井装置的实测资料，进一步论证本文算法的有效性，同时研究其他高效的并行算法.
Bai D H, Yu S. 1995. Theory and methods of resistivity tomography. Progress in Geophysics (in Chinese) (in Chinese), 10(1): 56-75. |
Barde-Cabusson S, Bolós X, Pedrazzi D, et al. 2013. Electrical resistivity tomography revealing the internal structure of monogenetic volcanoes. Geophysical Research Letters, 40(11): 2544-2549. DOI:10.1002/grl.50538 |
Boonchaisuk S, Vachiratienchai C, Siripunvaraporn W. 2008. Two-dimensional direct current (DC) resistivity inversion:Data space Occam's approach. Physics of the Earth and Planetary Interiors, 168(3-4): 204-211. DOI:10.1016/j.pepi.2008.06.022 |
Diaferia I, Barchi M, Loddo M, et al. 2006. Detailed imaging of tectonic structures by multiscale earth resistivity tomographies:the Colfiorito normal faults (central Italy). Geophysical Research Letters, 33(9): L09305. DOI:10.1029/2006GL025828 |
Ellis R G, Oldenburg D W. 1994. The pole-pole 3-D DC-resistivity inverse problem:a conjugategradient approach. Geophysical Journal International, 119(1): 187-194. DOI:10.1111/j.1365-246X.1994.tb00921.x |
Fan G D. 2017. Application of resistivity method in Peiling geological conditions of shale gas drilling platform.//China Petroleum Institute Geophysical Technology Workshop (in Chinese). Tianjin, 880-883.
|
Feng Y D, Chi X B, Wang W, et al. 2008. Fastcomputing for LURR of earthquake prediction. Pure and Applied Geophysics, 165(3-4): 749-759. DOI:10.1007/s00024-008-0313-0 |
Hauck C. 2002. Frozen ground monitoring using DC resistivity tomography. Geophysical Research Letters, 29(21): 12-1. |
Jones K A, Ingham M, Pringle D J, et al. 2010. Temporal variations in sea ice resistivity:Resolving anisotropic microstructure through cross-borehole DC resistivity tomography. Journal of Geophysical Research:Oceans, 115: C11023. DOI:10.1029/2009JC006049 |
Kelbert A, Meqbel N, Egbert G D, et al. 2014. ModEM:A modular system for inversion of electromagnetic geophysical data. Computers & Geosciences, 66: 40-53. DOI:10.1016/j.cageo.2014.01.010 |
Krautblatter M, Hauck C. 2007. Electrical resistivity tomography monitoring of permafrost in solid rock walls. Journal of Geophysical Research:Earth Surface, 112: F02S20. DOI:10.1029/2006JF000546 |
Li S C, Nie L C, Liu B, et al. 2015. Advanced detection and physical model test based on multi-electrode sources array resistivity method in tunnel. Chinese Journal of Geophysics (in Chinese) (in Chinese), 58(4): 1434-1446. DOI:10.6038/cjg20150429 |
Li Y, Hu X Y, Kim K, et al. 2010. Research of 1-D magnetotelluric parallel computation based on MPI. Progress in Geophysics (in Chinese) (in Chinese), 25(5): 1612-1616. DOI:10.3969/j.issn.1004-2903.2010.05.012 |
Li Y G, Oldenburg D W. 1992. Approximate inverse mappings in DC resistivity problems. Geophysical Journal International, 109(2): 343-362. DOI:10.1111/gji.1992.109.issue-2 |
Li Y G, Oldenburg D W. 2000. 3-D inversion of induced polarization data. Geophysics, 65(6): 1931-1945. DOI:10.1190/1.1444877 |
Li Y X. 2014. Seismic exploration in 2D frequency-domain forward based on parallel algorithms[Master's thesis] (in Chinese). Beijing: China University of Geosciences (Beijing).
|
Ma H, Tan H D, Guo Y. 2015. Three-dimensional induced polarization parallel inversion using nonlinear conjugate gradients method. Mathematical Problems in Engineering, 2015(Article ID 464793). DOI:10.1155/2015/464793 |
Oldenburg D W, Li Y G. 1994. Inversion of induced polarization data. Geophysics, 59(9): 1327-1341. DOI:10.1190/1.1443692 |
Park S K, Van G P. 1991. Inversion of pole-pole data for 3-D resistivity structure beneath arrays of electrodes. Geophysics, 56(7): 951-960. DOI:10.1190/1.1443128 |
Pelton W H, Rijo L, Swift C M Jr. 1978. Inversion of two-dimensional resistivity and induced-polarization data. Geophysics, 43(4): 788-803. DOI:10.1190/1.1440854 |
Petrick W R Jr, Sill W R, Ward S H. 1981. Three-dimensional resistivity inversion using alpha centers. Geophysics, 46(8): 1148-1162. DOI:10.1190/1.1441255 |
Revil A, Johnson T C, Finizola A. 2010. Three-dimensional resistivity tomography of Vulcan's forge, Vulcano Island, southern Italy. Geophysical Research Letters, 37: L15308. DOI:10.1029/2010GL043983 |
Schwarzbach C, Börner R U, Spitzer K. 2005. Two-dimensional inversion of direct current resistivity data using a parallel, multi-objective genetic algorithm. Geophysical Journal International, 162(3): 685-695. DOI:10.1111/gji.2005.162.issue-3 |
Shi G L, Zhang G X. 2008. Macroscopic Field Theory (in Chinese). Wuhan: China University of Geosciences Publishing House.
|
Singh U K, Tiwari R K, Singh S B. 2010. Inversion of 2-D DC resistivity data using rapid optimization and minimal complexity neural network. Nonlinear Processes in Geophysics, 17(1): 65-76. DOI:10.5194/npg-17-65-2010 |
Spitzer K. 1995. A 3-D finite-difference algorithm for DC resistivity modeling using conjugate gradient methods. Geophysical Journal International, 123(3): 903-914. DOI:10.1111/gji.1995.123.issue-3 |
Tan H D, Tong T, Lin C H. 2006. The parallel 3D magnetotelluric forward modeling algorithm. Applied Geophysics, 3(4): 197-202. DOI:10.1007/s11770-006-4001-5 |
Tang J T, Zhang J F, Feng B, et al. 2007. Determination of borders for resistive oil and gas reservoirs by deviation rate using thehole-to-surface resistivity method. Chinese Journal of Geophysics (in Chinese) (in Chinese), 50(3): 926-931. |
van der Vorst H A. 1992. Bi-CGSTAB:A fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems. SIAM Journal on Scientific and Statistical Computing, 13(2): 631-644. DOI:10.1137/0913035 |
Wang J Y. 2002. Inverse Theory in Geophysics (in Chinese). 2nd ed. Beijing: Higher Education Press.
|
Weiss C J. 2013. Project APhiD:A Lorenz-gauged A-φ decomposition for parallelized computation of ultra-broadband electromagnetic induction in a fully heterogeneous Earth. Computers & Geosciences, 58: 40-52. DOI:10.1016/j.cageo.2013.05.002 |
Wu X P, Wang T T. 2001. Study on some problems for 3D resistivity inversion using conjugate gradient. Seismology and Geology (in Chinese) (in Chinese), 23(2): 321-327. |
Wu X P, Xu G M. 2000. Study on 3-D resistivity inversion using conjugate gradient method. Chinese Journal of Geophysics (in Chinese) (in Chinese), 43(3): 420-427. |
Xiong J, Zou C C, Meng X H. 2012. Using the BICGSTABalgorithm with the incomplete LU factorization precondictioning to implement 2D FDFD induction logging fast forward modeling. Geoscience (in Chinese) (in Chinese), 26(6): 1283-1288. |
Yuval D, Oldenburg W. 1996. DC resistivity and IP methods in acid mine drainage problems:results from the Copper Cliff mine tailings impoundments. Journal of Applied Geophysics, 34(3): 187-198. DOI:10.1016/0926-9851(95)00020-8 |
Zhang J, Mackie R L, Madden T R. 1995. 3-D resistivity forward modeling and inversion using conjugate gradients. Geophysics, 60(5): 1313-1325. DOI:10.1190/1.1443868 |
白登海, 于晟. 1995. 电阻率层析成象理论和方法. 地球物理学进展, 10(1): 56-75. |
范国栋. 2017.电阻率法在涪陵页岩气钻井平台地质条件勘查中的应用.//中国石油学会2017年物探技术研讨会.天津, 880-883.
|
李术才, 聂利超, 刘斌, 等. 2015. 多同性源阵列电阻率法隧道超前探测方法与物理模拟试验研究. 地球物理学报, 58(4): 1434-1446. DOI:10.6038/cig20150429 |
李焱, 胡祥云, 金钢燮, 等. 2010. 基于MPI的一维大地电磁并行计算研究. 地球物理学进展, 25(5): 1612-1616. DOI:10.3969/j.issn.1004-2903.2010.05.012 |
李祎昕. 2014.地震勘探频率域二维正演并行算法研究[硕士论文].北京: 中国地质大学(北京).
|
施国良, 张国雄. 2008. 宏观场论. 武汉: 中国地质大学出版社.
|
汤井田, 张继锋, 冯兵, 等. 2007. 井地电阻率法歧离率确定高阻油气藏边界. 地球物理学报, 50(3): 926-931. DOI:10.3321/j.issn:0001-5733.2007.03.035 |
王家映. 2002. 地球物理反演理论. 2版. 北京: 高等教育出版社.
|
吴小平, 汪彤彤. 2001. 利用共轭梯度方法的电阻率三维反演若干问题研究. 地震地质, 23(2): 321-327. DOI:10.3969/j.issn.0253-4967.2001.02.028 |
吴小平, 徐果明. 2000. 利用共轭梯度法的电阻率三维反演研究. 地球物理学报, 43(3): 420-427. DOI:10.3321/j.issn:0001-5733.2000.03.016 |
熊杰, 邹长春, 孟小红. 2012. 不完全LU分解预条件BICGSTAB算法实现感应测井二维FDFD快速正演模拟. 现代地质, 26(6): 1283-1288. DOI:10.3969/j.issn.1000-8527.2012.06.022 |