目前, 电磁勘探理论技术发展已经比较成熟, 数据处理解释已经形成比较全面、系统的方法技术系统, 一些新的处理技术和反演方法在实际应用中发挥了关键作用。随着电磁勘探应用领域的不断扩大, 面临的地质构造等复杂程度增加, 解释精度要求高, 解释难度加大, 因而对处理技术和解释方法提出了更高的要求, 对处理解释软件的要求也相应提高。

现有的电磁处理解释软件有国外的WinGlink和GeoTools, 软件功能比较全面, 相对比较成熟, 但是价格昂贵, 且与国内系统存在兼容性问题, 使用不便。其中WinGlink已经不再对外发布更新, 现有的版本存在两大短板, 一是缺少资料分析方法, 二是反演方法单一, 很难满足现阶段高精度解释的要求。国内许多科研院所也在软件开发方面做了大量的研究, 取得了一定的成果。最早由魏胜等^[1]推出的电磁解释工作站, 在一定程度上推动了国内电磁资料处理解释软件的发展。现在比较常用的有王绪本等^[2]开发的MtSoft和陈小斌等^[3-4]推出的MTPoineer, 都得到了很好的推广。经过十几年的沉淀和积累, 国内的软件在功能和体验上取得了很大的进步, 获得了从业者的肯定。但是, 传统的软件多是利用现有的可执行程序打包集成, 更多的是一个软件包, 因而软件执行效率较差, 而且仅开发了32位软件版本, 在二维大数据量计算方面效率较低, 很难满足大数据量反演和精细解释的需求。

当前, 电磁法已经迈入精细勘探和精细解释阶段, 对数据处理和反演算法提出了更高的要求, 需要我们开发能够满足当前处理解释需要的一体化软件, 提供不同特点的算法以满足不同数据的需要和各种类型的地质解释。实践证明, 不同的分析算法适应不同的数据类型, 不同的反演方法对于各类地质构造、异常体等反演能力也有一定的差距。基于此, 本文以C++和Fortran混编的方式, 基于动态链接库(DLL)集成, 开发了基于Windows系统(Windows XP/7/10)的大地电磁处理解释系统(MTI)。各类算法均基于DLL集成, 尤其是各种二维反演算法的集成, 相比于单独运行子程序提高了执行效率, 对同一数据实现多种方法的快速对比解释, 保证了反演的可靠性。

1 软件关键技术

MTI软件基于C++和Fortran混合编程, 充分利用了微软基础类库(MFC)的可视性和Fortran强大的计算能力。界面基于MFC开发, 所有的算法都采用Fortran编写或者利用现有的算法代码改编。软件的两大关键技术是人机交互的可视化实现和动态链接库的使用, 人机交互通过MFC来实现, 主要包括软件的系统设计和可视化操作两方面, 而动态链接库涉及相关算法和代码的优化以及通信等问题。

1.1 系统设计

参照Windows程序的设计模式, 结合电磁处理解释的需要, 采取工程化的软件设计理念。基于数据量方面的考虑, 所有的操作都在工程中进行, 最大程度保证操作的流畅性和软件运行的稳定性。为了达到更好的人机交互效果, 充分利用MFC的可视性, 软件整体采用分层设计, 无论是主界面还是各子系统界面, 操作按键与图形显示一目了然。各个子系统独立分布, 便于后期维护和更新, 子系统间通过公用数据类进行数据交换, 在降低内存需求的同时保证数据流的安全。

1.2 可视化操作

在人机交互操作中, 最核心的技术是解决图像的实时动态显示问题, 因为每一次鼠标操作都将带来图像的改变和重绘, 造成图像的闪烁。因此采取了基于双缓冲的成图机制。所谓双缓冲就是创建一个内存设备, 然后在这个内存设备上成图, 完成后再将它整体拷贝到屏幕上, 这样能提高成图的速度, 同时解决闪屏问题。由于MTI支持Windows 10操作系统, 可以安装在平板上操作, 希望借此便携的设备以及更快的速度, 将数据处理解释由室内带到野外, 带来新的进步。

1.3 动态链接库

目前, 混合编程一般有两种实现方式^[5]。一是以子程序的形式单独进行计算, 与主程序之间通过文件或者其它方式实现数据传递。这种方式的优点是简单易实现, 但是效率较低, 目前多数电磁处理解释软件都是通过这种方式实现的。另一种方式是以DLL的形式集成到主程序中, 需要改写现有的Fortran程序, 软件开发工作量较大, 但是软件运行效果最优。MTI系统全部采用DLL方式集成, 只有一个主程序, 功能模块或者计算模块全部由DLL完成, 没有利用任何一个现有的可执行程序或者软件包, 高度集成, 从根本上提高了软件的执行效率。

动态链接库是微软提供的一种共享函数的方式, 在Windows系统中得到了广泛使用。其本质可以理解为一个可被重复调用的函数, 最大的特点是可被多种语言共享调用且安全性极高。传统电磁处理解释软件核心算法都交由子程序来独立执行, 而DLL调用需要主程序与子程序共同完成, 涉及到主程序与DLL的即时通信问题, 需要给出最合适的技术解决方案。

1.4 多线程工作模式

MTI系统采用了多线程的工作方式, 尤其是在二维反演计算中, 开发了工作线程、用户界面线程分别负责冗长的后台计算和处理结果反馈。工作线程没有消息循环, 只负责后台的数据计算; 界面线程可以在计算过程中保证程序持续运行, 继续接受消息队列, 处理用户反馈以及系统响应等。多线程的工作方式保证了复杂的计算独立于主程序之外, 对于提升软件性能和用户体验效果具有很大的帮助。

1.5 内存优化

计算量是制约电磁处理解释软件发展的关键因素, 随着精细化勘探解释的需要, 数据量会越来越大, 而又很难从算法角度进行改进, 所以必须从内存方面对代码进行优化。MTI采取了两种方式, 一是对现有的Fortran代码进行优化, 将大数组改为动态存储或者是采用数据文件进行缓存; 二是充分利用计算机内存, 由于32位的Windows系统最大支持4G内存, 而64位的系统对内存没有限制, 于是对原32位的软件和动态链接库进行了改写, 编译了64位版本, 最大限度地使用计算机内存, 对大数据二维反演是很好的支撑。内存优化带来的优势不仅是处理量的提升, 更是执行效率的提高。

2 软件功能开发

基于大地电磁数据处理、解释的需要, MTI软件主要分为4个功能模块, 包括数据管理、数据分析、数据处理和数据解释, 各个功能模块又下设多个子系统, 相应的功能由各子系统来完成, 详细功能如图 1所示。

2.1 数据管理

数据管理主要实现从工程建立到测线测点管理以及坐标、剖面等数据的输入、输出等。MTI系统采取工程、测线、测点三级管理, 结构明了, 层次清晰, 所有操作均在工程中进行, 简单安全。MTI系统主要支持两种电磁数据格式, 分别为功率谱数据文件和阻抗倾子数据文件, 从而支持了GMS系列、凤凰V5V8以及EH4等主流电磁仪器采集的数据, 基本支持所有的EDI数据文件。

2.2 数据分析

在实际资料中, 由于局部畸变等影响, 观测到的数据往往不能最大限度地反映地下真实信息, 在进行二维反演的时候需要将数据旋转到构造主轴方向。MTI系统中数据分析模块主要分析阻抗张量和倾子矢量, 主要手段包括张量分解和极化模式识别, 核心目的是求得最能反映地下真实信息的数据用于反演, 两大难点是电性主轴的选取和极化模式的判别, 二者直接影响最终的反演结果。

张量分解的目的是从观测的大地电磁响应中消除局部畸变的影响, 获得区域构造阻抗、电性主轴和走向等参数^[6], 常用的方法有Swift分解^[7]、Bahr分解^[8-10]、相位张量^[11-12]和GB分解^[13-14], 不同的方法都有其优缺点。MTI系统提供了3种张量分解方法(Swift分解、Bahr分解、相位张量(CBB)), 基于统计玫瑰图进行快速定性分析, 如图 2所示。利用阻抗张量极坐标图可以分析构造走向、维性及极化模式(图 3)。

2.3 数据处理

MTI系统中, 数据处理和反演为两大核心模块, 人机交互操作的可视化程度将直接影响软件体验效果和处理水平。MTI将成图模块与处理模块统一放在视图类中, 采取双缓存绘图的形式解决人机交互过程中的闪屏问题。

坐标测线管理：可以实现地理坐标系与平面坐标系之间的转换, MTI实现了WGS84, 北京54和西安80三种坐标系统之间的任意转换, 实现了测点自动拟合、手绘测线自动投影和按照设计方向投影3种测点投影方式。

数据处理：进行解释时, 需要将数据尽可能地平滑以满足反演拟合的需要, 设计了三点、五点、七点线性平滑以及二次、三次等7种圆滑算法, 满足不同质量数据的需要。对于质量较差的数据, 设计了基于交互式操作的人工挑点、单频点删除、多频点删除等, 基于反演需要, 可以完成测线整体截频。结合数据分析, 设计了主轴旋转、极化模式选择、数据交换、静校正(手动、自动)等整体化数据处理功能。详细功能如图 4和图 5所示。

2.4 数据反演

数据反演是MTI的核心模块, 目前主要的解释方法是二维反演。一维反演主要用于电性结构分析, 便于对地下结构有整体的认识, 而三维反演方法^[15-20]虽然已经进入实用化阶段, 但是由于计算效率低等不足, 全面推广尚需一定的时间。目前常用的解释软件中一般只提供一种二维反演方法, 并以非线性共轭梯度反演(NLCG)居多。考虑到反演的多解性, MTI系统实现了3种一维反演和3种二维反演。由于全部采用了动态链接库的集成方式, 各种反演方法所需的数据、模型等参数格式不再受代码限制, 使用同一个模型可以完成所有的反演解释, 极大地提高了反演的效率, 并且可以在相同条件下进行快速对比解释, 以尽可能地保证解释的准确性。

2.4.1 一维反演

MTI可以完成Bostick反演、奥可姆(OCCAM)反演^[21]和广义逆反演。Bostick反演简单快速, 但精度较差; OCCAM反演稳定性和精度较好; 广义逆反演分层性较好, 适合定量解释, 但是高度依赖初始模型。3种反演方法具有不同的适应性, 可以各取所长, 互补所短, 针对不同的地质情况提供了多种选择。基于初始模型对于反演的重要性, MTI提供了4种初始模型设计方案, 分别为自动模型、基于初始量的自动模型、分层模型和人机交互模型, 能够满足不同的解释需要。

MTI系统也对算法本身进行了改进, 在OCCAM算法中增加了一个反演因子, 该反演因子直接参与到反演过程中, 能够控制最终的模型粗糙度。经过大量数据实测, 总结出了模型优化算法, 当检测到所设置的初始模型容易引起假异常或者影响收敛时, 会利用相应的优化算法自动对初始模型进行优化, 以使其更利于反演。

2.4.2 二维反演

网格剖分是二维反演的基础, 尤其是对初始条件要求高的反演方法, 会直接影响反演结果。人机交互的可视化操作对于剖分必不可少, MTI可以基于高程数据自动生成带地形的模型, 通过交互式操作实现对模型的修改。同时, 模型对于所有反演方法都通用, 每个方法既可以使用当前模型, 也可以使用其它方法已经反演过的模型继续反演, 实现多种方法的联动反演。

NLCG和OCCAM作为目前主要的两种反演方法, 代表了当今电磁反演的精髓。MTI中可以完成OCCAM^[22], 简化基奥可姆(REBOCC)^[23], NLCG^[24-25]3种反演, 其中NLCG包含了带地形的反演。这3种反演方法中, OCCAM反演比较稳定, 初始条件要求低; REBOCC反演速度较快, 但是反演精度不如OCCAM; 而NLCG反演精度较高, 反演速度较快, 尤其是带地形反演对于地形起伏较大的数据解释优势明显, 但是对初始条件依赖程度较高。现在计算机硬件水平已经完全能够满足二维快速反演的需要, 因此多种方法的联合解释也变成了一项基本要求, 这既降低了反演的多解性, 增加了解释的可信度, 又能够对反演结果进行快速评价。

3 算例分析

为了检验MTI中各个子系统的有效性和动态链接库的准确性, 设计了多个模型验证各反演方法对各种地质模型的适应性。同时基于相同的初始条件, 与经典的电磁处理解释软件WinGlink进行对比研究, 保证了对比的公正性和结果的有效性。本文以断层模型和层状模型为例进行说明。

3.1 模型一

模型一为一个断层模型, 如图 6所示, 主要用来检验反演方法对于层位的识别效果和对断层参数的反演精度。图 7到图 9给出了采用MTI软件中3种方法和WinGlink软件的NLCG方法对断层模型在不同模式下的反演结果。从图 7到图 9可以看出, MTI中的3种反演方法均能较准确地反演出模型的相关参数(边界、电阻率), 能较好地识别断层, 尤其是OCCAM和NLCG反演精度较高, 不逊于WinGlink软件的反演结果, 相比而言, REBOCC的反演精度稍差。从反演细节看, OCCAM的反演结果优于两种NLCG反演结果, 对各项构造参数的识别均比较接近实际, 尤其是层位识别, 优势较为明显, 说明OCCAM反演方法对于断层的反演效果更优。就3种模式而言, TE模式和联合模式的反演结果更接近于设置的模型。

3.2 模型二

模型二为一个层状模型, 总共分为5层, 第1层和第3层较薄, 第4层和第5层为厚度递增的高阻层, 如图 10所示。首先进行了一维反演, 如图 11所示, 3种一维反演结果对于层状划分效果均不错, 各个反演结果的拟合程度也较高。OCCAM反演对于每一层均有较好的识别效果, 拟合精度很高; 而广义逆反演分层性更加明确, 层位参数更加准确, 在反演过程中其厚度会参与反演, 基本与我们设置的模型参数一致, 说明其对于层状模型的地电结构划分效果更好。

图 12到图 14给出了采用MTI软件中3种方法和WinGlink软件的NLCG方法在不同模式下对层状模型的反演结果。可以看出各个方法的分层性均不错。OCCAM的反演精度最高, WinGlink的NLCG反演效果整体表现不如MTI。而仔细对比两种NLCG反演结果可以发现, 两种软件的NLCG反演结果中均有一定的边界效应, 说明NLCG反演对初始条件要求较高, 其稳定性较OCCAM低, 建议与OCCAM反演结果进行联合解释。

4 实测数据对比分析

选取2016年在青海某地采集的大地电磁数据进行处理解释, 旨在查明10km以上的地层结构、构造特征及花岗岩的电性结构等。该测线包含66个测点, 分别用MTI和WinGlink进行了带地形的反演, 用imagem的Bostick进行了简单反演, 反演结果如图 15所示。

3种软件不同方法的反演结果整体较为一致, 基本查明了第四系、新近系和古近系的厚度, 下伏花岗岩受F1, F2断层控制逆冲至浅地表; 测线中段受断裂控制形成一局部地堑构造, 新近系和古近系沉积厚度约6km。两种NLCG的反演结果无论是在构造参数还是电阻率方面都具有很强的一致性, 而OCCAM与NLCG对F2断裂的反演有细微差别, 与算法本身和其对模型的依赖程度有一定的关系。相比而言, Bostick由于是简单反演, 其精度受到一定的限制。但3种软件的反演结果都反映了地下电性结构特征, 也证明了MTI软件的可靠性和稳定性。

MTI和WinGlink软件3种反演方法的反演参数如表 1所示。可以看出OCCAM的单次反演速度较慢, 但是其达到设计精度所需的迭代次数较少, 反演稳定性较高, 能够得到较好的反演结果。NLCG算法计算优势明显, 但是对于初始模型要求较高, 迭代收敛性和反演粗糙度较差。综合对比两种软件, NLCG算法基本一致, 但是MTI的OCCAM算法更加稳定, 在拟合差和粗糙度方面的寻优能力更强。

基于反演结果, 设计了能够反映地下电性结构的地质模型, 模型由两个主要断裂控制, 如图 16a所示。基于此地质模型进行正演, 再分别用MTI和WinGlink进行了反演解释, 一方面为了验证实测数据反演结果的准确性, 另一方面可以再次对比反演软件的可靠性。

从图 16中可以看出, 3种反演结果均较好地反映了地下电性结构特征, 包括古近系、新近系和第四系的厚度分布、断裂参数等, 电阻率结果也均与模型接近。基于地质模型的反演结果, 能够辅证实测数据解释结果的可靠性。整体看, 3种反演方法对于复杂模型的反演能力均较好, 反演精度符合预期要求。从细节看, OCCAM的反演稳定性较好, 反演得到的构造参数与模型设置更加接近, 相比而言两种软件的NLCG反演结果对于断裂的识别稍有偏差, 但是两种软件的反演结果具有很强的一致性。表 2给出了MTI和WinGlink软件3种反演方法反演的参数。由表 2可见, OCCAM在收敛性和模型稳定性方面优势明显, MTI和WinGlink的NLCG算法在迭代收敛性和计算速度上基本一致。

综合对比两种软件, 对于较为复杂的地质解释, OCCAM算法更加稳定, 对初始模型要求较低, 拟合寻优和控制模型粗糙度的能力更强; NLCG算法计算速度较快, 但是模型粗糙度较大, 对初始模型要求较高。对于复杂构造的地质解释, 建议采用对比解释或者OCCAM和NLCG联合反演, 充分利用二者的优点, 取长补短, 相辅相成。两种软件的NLCG算法表现基本一致, 对不同的地质模型其反演能力存在一定的细微差别, 主要是由于初始模型的设计模式不同引起的。由于MTI设计了64位的软件模式, 能够承载更多测点的计算, 满足较大数据量的反演要求, 较WinGlink存在一定的内存优势。

5 结论

基于DLL集成的MTI软件集成度高, 对内存进行了一定的优化, 执行效率高。基于MFC开发的人机交互系统, 能够完成电磁资料(MT/AMT/EH4)的数据处理解释一体化操作。系统包含了完整的数据分析方法和目前主流的经典反演算法, 功能比较全面, 可以实现快速的对比分析和解释, 尤其是OCCAM和NLCG两大经典算法的加入, 保证了反演结果的可靠性, 从技术层面降低了解释结果的多解性。开发了64位软件版本, 在大数据量的处理解释方面具有很大的优势。

目前, 软件在功能方面依旧需要进一步完善, 随着电磁法和计算机的快速发展, 我们将对MTI继续改进更新, 加入新的电磁方法和新的计算机技术, 以使MTI满足电磁勘探技术发展的需要。