地～井方位激电（IP）人机交互解释软件开发及应用

吕玉增，，彭苏萍

（1.中国矿业大学 煤炭资源与安全开采国家重点实验室，北京 100083；2.桂林理工大学 地球科学学院，桂林 541004）

地～井方位激电（IP）人机交互解释软件开发及应用

吕玉增1，2，，彭苏萍1

（1.中国矿业大学 煤炭资源与安全开采国家重点实验室，北京 100083；2.桂林理工大学 地球科学学院，桂林 541004）

基于MSR（Modified Sparse Row）非零元素存储和SSOR－PCG方程求解技术，编制了三维IP的有限元法快速正演程序，在此基础上，利用Delphi等开发了可视化地～井方位IP人机交互解释软件。在地～井方位IP的数据解释中，遵循“先井口，后其他方位”的分析方法，先通过井口方位观测曲线初步断定异常体的电性特征、大致埋深和倾向；然后对比分析其它方位曲线的异常特征，推断出异常的赋存方位。对于异常复杂情况，该软件可利用测井数据、钻孔柱状图等，形成分层背景，通过差值法和比值法数据预处理技术以突出井旁异常，提高异常分辨能力。实测数据的解释结果表明，针对地～井方位IP观测数据量少，充分利用测井、钻孔等信息，用人机交互正演拟合反演方式解释是可行的。

地～井方位观测；有限单元法；IP；人机交互解释；软件

地～井方位激电观测是井中激电的工作方式之一，指在地面供电，井中测量的井中激发极化法，主要用来发现井旁盲矿，确定其埋藏深度及方位，追踪和圈定矿化带等。由于地～井方位观测是在井中观测，能更好地接近深部的目标异常体，因此，在深部隐伏矿产的勘查以及老矿山的“探边摸底”等方面有着广泛的应用。

近年来，国内、外对井地、井井电法的研究非常重视，取得了丰富的研究成果［1～13］，但对地～井电法的研究成果偏少［14、15］。一方面由于地～井方位观测野外施工难度大；另一方面地～井方位观测正反演所用的模型单元数多，而可利用的实测数据少，正演模拟和数据解释难度相对较大［16］。针对地～井方位IP的实际情况，作者提出人机交互解释思路，并开发相应软件［17～19］。

1 地～井方位IP的快速正演

地～井方位IP的人机交互解释需建立在快速、准确的三维正演基础之上。有限单元法是研究复杂三维地电模型正演的最有效数值模拟方法之一，但它需要在区域内进行网格剖分。对于地～井方位观测而言，由于沿井方向的网格剖分数量大导致整个网格剖分数量巨大，数值计算上存在存储、计算速度和精度等多个矛盾。

下面介绍解决上述矛盾，实现三维地－井方位IP有限元快速正演的技术。

1.1 MSR系数矩阵非零元素存储

有限元法求解地～井方位IP的正演问题，通过网格剖分、单元分析、矩阵合成等步骤，最终形成系线性方程组：

其中 K是系数矩阵；u是全部节点的电位ui组成的列向量；p是与点电源项有关。

有限元法正演形成的系数矩阵K是对称、稀疏矩阵，其中绝大部份元素值为零，但每行零元素的分布和个数不同，与网格剖分和节点编号有关例如对于直角坐标系中网格节点（x×y×z）＝（5 ×5×5）的三维四面体剖分，网格节点个数125（5 ×5×5），系数矩阵K的大小125×125＝15 625，常规变带宽存储这个下三角阵需要3 225个元素，而非零元素个数仅有810个，可见非零元素存储能明显减少内存存储。

MSR（Modified Sparse Row）是一种改进的行压缩索引存储技术，是以行为单位，顺序存储系数矩阵中的非零元素。考虑到在方程组的求解过程中，系数矩阵对角线元素均不为零，且对角线元素的访问和运算最频繁，把对角线元素提出来单独存储（见图1）。该存储方式仅需要两个数组：一个实型数组AA和一个整型数据ICOL，具体数组说明见表1。

1.2 SSOR－PCG迭代法求解方程组

SSOR－PCG迭代法求解方程组（1）的过程为［19］：对称系数矩阵K可写为对角阵D和严格下三角阵E的和形式：K＝D－E－ET，而M为K的 SSOR预条件矩阵，则：

其中 ω∈［0，2］为松弛因子，通过试算确定SSOR－PCG求解的整个迭代过程如下：

由于预条件矩阵M与系数矩阵K具有相同的稀疏性，同样也采用MSR压缩存储，整个迭代过程只是系数矩阵非零元素的运算，计算量小，计算速度快。

图1 四面体网格剖分和系数矩阵下三角非零元素分布Fig.1 Tetrahedral mesh and lower triangular distribution of nonzero elements in the coefficient matrix

表1 MSR存储说明Tab.1 MSR storage instruction

2 地～井方位IP人机交互解释软件设计

2.1 软件功能和执行流程

该软件主要有两大功能：①地～井方位IP三维有限元正演模拟；②地～井方位IP数据的人机交互反演。软件是基于三维有限元正演基础上开发的，适合于地～井多方位观测方式下，任意的三维地电体的直流激电正演模拟和正演拟合交互反演解释。三维有限元正演软件包含测井基本信息输入、三维异常体模型构建、有限元正演、显示结果曲线、保存结果数据等功能，软件执行流程见图2。野外数据的反演拟合步骤大致分为实测数据输入、图示数据、测井数据分层、建立层状模型、三维模型构建、曲线对比拟合等步骤，详见图3流程图。

2.2 软件菜单及功能项介绍

地～井方位IP人机交互解释软件是在Delphi 7.0和Fortran 4.0开发平台下开发完成的可视化软件，软件设置了文件、数据输入、数据切换、模型构建等菜单项，各菜单项功能和对应的快捷方式见图4和下页表2。

3 软件几个主要功能的实现

3.1 测井工程文件的建立

建立测井工程文件是地～井方位IP正演和交互解释的第一步，用户根据地～井方位的模型和实际工作情况，通过“工作参数输入”对话框输入（编辑）方位测井的工作参数（见后面图5），包括地面方位点的供电位置、井孔深度、测点距等信息。当完成工作参数输入后，软件系统会根据工作参数信息，绘出地～井方位测井示意图。这时用户可以进行三维建模等操作，还可以通过功能项保存当前的模型到工程模型文件（.GP）中，便于下一次的直接打开和编辑。

表2 菜单选项、工具条及其功能介绍Tab.2 Menu options，toolbars and features

3.2 三维建模和修改

测井工作参数输入完成后，系统自动对三维空间进行水平分层。三维建模的思路是：将每个深度层看作是一定厚度的板状体，通过对各不同的深度层建模，建立复杂的三维模型。具体编辑模型操作步骤如下：

（1）选择要编辑的深度层。

（2）用不同的颜色区分不同的异常体，通过鼠的颜色示意在主界面右侧深度层选择区域上。

3.4 数据预处理和交互解释

数据预处理的目的是为下一步的分析、解释提供帮助。根据方位测井的实际经验，软件提供了比值法、差值法二种预处理方法，以突出旁侧异常。比值法（差值法）预处理就是把所有方位观测数据统一除以（减去）井孔方位观测数据值，以压制孔壁局部影响，突出旁侧异常信息。

根据预处理后的曲线结果，结合原始数据曲线特征和钻探等地质信息，遵循“先井口、后其他方位”的分析方法，推断异常体的大致深度、方位和异常体电性特征。并初步建立预测初始模型，进行正演试算。将模型正演结果与实测数据对比，不断修正模型并再次正演，直到结果满足要求为止。

4 应用实例

广西某矿山钻孔深度493m，孔径约100mm，采用了地～井方位激电观测，目的是了解已知低阻高极化矿带在该钻孔附近的位置和隐伏情况。勘测单位进行了钻孔四周四个方位探测（A1，A2，A3和A4），各方位距离相等，井中观测电极距MN＝10m，相邻测点距10m，井中实际观测位置从深度为100m～480m。各方位点相对坐标、高程和方位距离见表3，视电阻率、视极化率等实测数据保存在数据文件中。为了便于对比分析，我们用钻孔四周四个方位探测（A1、A2、A3和A4）数据的平均值当作是井口方位（A0）的探测结果。

表3 各方位点相对坐标、高程和方位距离Tab.3 Relative coordinates and azimuth distances of all sites

人机交互正演拟合反演解释过程如下。

（1）工作参数和实测数据导入。进入“工作参数输入”对话框，导入供电点坐标、输入井深、直径和观测参数（见图8（a））。当基本参数输入完成后，导入视电阻率、视极化率等实测数据到主界面左侧相应数据表格中，并绘制出四个方位的观测曲线（见图8（b））。

图8 工作参数与实测数据输入和观测曲线显示Fig.8 Working parameters（a），the measured data input an result curves show（b）

（2）根据实测曲线异常特征建立初始模型。从图8方位测井视电阻率和视极化率曲线上，有二处较明显的异常段：①深度300m－350m区段，四个方位都有明显低阻高极化异常，且视电阻率和视极化率异常对应较好，推断矿体大致埋深320m另外，四个方位激电异常不对称，对比发现A1和A2方位视电阻率和视极化率异常相近，而A3和A4方位观测的视电阻率和视极化率异常相近且幅值较大，初步断定矿化体靠近A3、A4方位，结合地质资料，可初步推断该异常体为近似水平板状体②在钻孔的底部450m附近，低阻高极化异常，初步推断引起该激电异常的是位于钻孔底部附近另一个矿化体。

通过上述分析，建立初始模型，背景电阻率和极化率近似取观测数据平均值，分别为410Ω·m和8%；依据该地区的矿石标本物性参数，异常体电阻率和极化率取值分别为50Ω·m和20%。

（3）模型修改和交互拟合解释。在初始模型基础上，利用“共轭梯度法”进行模型快速正演计算并反复修改模型、正演计算，对比模型正演与实测曲线的拟合情况，正演拟合以曲线形态拟合为主图9为修改后最终模型有限元数值解与实测数据曲线的对比情况。

图9 模型正演（红）与实测数据（蓝）拟合对比Fig.9 Fitting comparison forward modeling（red）with the measured data（blue）

（4）解释结果。最终的拟合结果显示，在该钻孔330m～360m深度上，矿带比较靠近钻孔，且在A3、A4方位方向有一定延伸，延伸长度约40m；另外，推断在钻孔底部下方深度450m以下，赋存另一高阻高极化矿化体。

根据勘测单位钻探情况，该解释结果基本证实了地质、物探人员的判断，即该钻孔已经穿过浅部的矿脉，但还没有揭露深部的矿化带。同时，矿体埋深、赋存方位等推断解释结果也为勘测单位提供了帮助。

图10 正演拟合交互解释结果Fig.10 The result of forward interactive interpretation

5 结论

（1）利用MSR非零元素存储和SSOR－PCG方程求解技术实现了三维IP的快速正演，为实现地－井方位激电（IP）人机交互解释奠定了基础。

（2）利用Delphi等软件平台开发了可视化地～井方位IP人机交互解释软件，在适合于地～井多方位观测方式下，任意的三维地电体的直流激电正演模拟和正演拟合交互反演解释。

（3）在地～井五方位IP实测数据的推断解释中，可遵循“先井口，后其他方位”的分析方法，先通过井口方位观测曲线初步断定异常体的电性特征、大致埋深和倾向；然后对比分析其它方位曲线异常特征，推断异常的赋存方位。对于异常复杂情况，可利用测井数据、钻孔柱状图等，形成分层背景，通过差值法和比值法对比分析不同方位观测曲线特征，以突出井旁异常，提高异常分辨能力。

（4）实测数据的解释结果表明，针对地～井方位IP观测数据量少，充分利用测井、钻孔等信息，用人机交互正演拟合反演方式解释是可行的。

致谢：

感谢上海应用技术学院黄俊革教授给予的帮助。

［1］ ZHDANOV M S，YOSHIOKA K.Cross－well electromagnetic imaging in three dimensions［J］.Exploration Geophysics，2003（34）：34.

［2］ ZHOU B，GREENHALGH S A.Rapid 2－D/3－D crosshole resistivity imaging using the analytic sensitivity function.Geophysics［J］，2002，67（2）：755.

［3］ WILKINSON P B，CHAMBERS J E.Extreme sensitivity of crosshole electrical resistivity tomography measurements to geometric errors［J］.Geophys.J Int.，2008，173：49.

［4］ PARDO D，CALO V M，TORRES－VERDIN C，e al.Fourier series expansion in a non－orthogona system of coordinates for the simulation of 3D－DC borehole resistivity measurements［J］.Comput Methods Appl.Mech.Engrg.2008（197）：1906.

［5］NAM M J，PARDO D，TOORES－VERDIN C Simulation of DC dual－laterolog measurements in complex formations：A Fourier－series approach with nonorthogonal coordinates and self－adapting finit elements［J］.Geophysics，2009，74（1）：E31.

［6］ 安然，李桐林，徐凯军.井地三维电阻率反演研究［J］地球物理学进展，2007，22（1）：247.

［7］ 王志刚，何展翔，刘昱.井地直流电法三维数值模拟及异常规律研究［J］.工程地球物理学报，2006，3（2）87.

［8］ 柯敢攀，黄清华.井地电法的三维正反演研究［J］.北京大学学报：自然科学版，2009，45（2）：294.

［10］底青云，王妙月.积分法三维电阻率成像［J］.地球物理学报，2001，44（6）：844.

［11］沈平，强建科，李永军，等.井间视电阻率的几何成像方法［J］.中南大学学报：自然科学版，2010，41（3）1079.

［12］潘纪顺，葛为中，折京平.地面/井地/井间超高密度电阻率成像技术［J］.华北水利水电学院学报，2010，31（2）：74.

［14］郭文波，宋建平，李貅，等.层状介质井中电测数值计算及其应用研究［J］.地球物理学报，2006，49（5）1561.

［15］蔡柏林，黄智辉，谷守民.井中激发极化法［M］.北京地质出版社，1983.

［18］吕玉增.地－井、井－地IP三维快速正反演研究［D］.长沙：中南大学，2008.

book=7,ebook=7

1001—1749（2012）03—0358—07

P 631.3＋24

A

10.3969/j.issn.1001－1749.2012.03.22

吕玉增（1978－），男，河北邢台人，博士后副教授，主要从事电法数值模拟与反演成像研究。

国家自然科学基金（40774057）；全国危机矿山接替找矿项目新技术新方法示范项目（200799086）；广西自然科学基金（桂科自0832263）

2011－07－27改回日期：2011－10－31