本研究由张昆（中国地质科学院矿产资源研究所/中国地质大学（北京））、董浩（中国地质大学（北京））、严加永（中国地质科学院矿产资源研究所/中国地质大学（北京））等合作完成，发表于《地球物理学报》（Chinese Journal of Geophysics）2013年第56卷第11期，页码3922-3931，DOI：10.6038/cjg20131134。
本研究属于地球物理勘探领域，专注于大地电磁（Magnetotelluric，MT）三维反演方法的改进与优化。传统上，在三维地质环境中采集的MT数据常被近似解释为二维结果，但在复杂地质区域，这种简化方法往往无法满足精度要求。三维MT反演面临两大核心挑战：一是观测数据远少于网格节点导致的超欠定病态矩阵问题；二是由大规模网格节点带来的巨大计算量。
研究团队旨在开发一种高效、实用的三维MT反演算法，使其能够在普通计算机上运行，降低对硬件性能的要求。该研究基于Newman和Alumbaugh（2000）提出的三维非线性共轭梯度（Nonlinear Conjugate Gradient, NLCG）算法，以及Rod和Mackie（2001）的二维NLCG反演预处理方法，通过改进预处理方法和引入并行计算技术，实现了计算效率的显著提升。
研究首先建立了精确的三维正演模拟基础。采用交错网格有限差分法（Staggered Grid Finite Difference Method）求解Maxwell方程组。在笛卡尔坐标系下将地下空间划分为nx×ny×nz个小长方体网格单元，电场分量定义在网格边缘中点，磁场分量定义在网格面中心。通过离散化处理，得到了描述电磁场分布的线性方程组：
iωμ₀E = ∇×(1/σ∇×E)
其中ω为角频率，μ₀为真空磁导率，σ为电导率，E为电场强度。研究采用Smith（1996）提出的静态散度校正方法，有效提高了低频阻抗响应的计算精度。
反演问题采用Tikhonov正则化方法构建目标函数：
φ = Σ[(Z_obs - Z)/ε]² + λmᵀWᵀWm
其中Z_obs和Z分别表示观测阻抗和正演阻抗，ε为数据误差，W为正则化矩阵，λ为正则化因子。该目标函数同时考虑了数据拟合差和模型平滑度约束。
研究采用非线性共轭梯度法求解这一优化问题，其核心迭代过程包括： - 初始模型设定（m₀） - 线性搜索确定模型修改量（αₖ） - 更新搜索方向（pₖ）
本研究的关键创新之一是改进了预处理方法。传统NLCG算法中，每次线性搜索都需要计算Hessian矩阵，导致巨大的计算量。研究提出使用预处理矩阵C = (γI + λWᵀW)⁻¹代替Hessian矩阵计算，其中γ是与当前迭代模型电阻率和数据误差相关的参数。这一改进使得每次迭代每个频点仅需6次正演计算（1次响应计算+5次模型修改量计算），相比传统方法减少了4次正演计算。
为应对三维正演的巨大计算量，研究采用OpenMP并行计算平台，实现了分频点并行计算策略。具体实现方式为： 1. 对所有稀疏矩阵进行一维定位存储，减少内存消耗 2. 根据计算机线程数量自动分配各线程的工作量 3. 将不同频点数据分配给不同CPU线程进行并行计算
这种并行方式避免了模型分块导致的边界条件破坏问题，保证了计算精度。测试表明，使用8个线程时，计算效率可提升至串行计算的约8倍。
研究设计了一个包含4个异常体（2个低阻10Ωm，2个高阻1000Ωm）的理论模型进行算法验证。模型参数包括： - 频率范围：0.01-100Hz（8个频点） - 网格规模：30×30×38（含8层空气网格） - 测点数量：64个 - 数据误差：10%（附加1%噪声） - 正则化因子：0.01
使用日本Kayabe地区实测数据进行算法实用性验证： - 频率范围：64-2Hz（6个频点） - 网格规模：23×23×38 - 测点数量：150个 - 正则化因子：0.1
经过52次迭代后，反演达到终止条件（RMS=1）。结果显示： - 在1.5-5.5km深度范围内，4个异常体均能被清晰识别 - 低阻异常体中心电阻率约10Ωm，与理论值一致 - 高阻异常体中心电阻率约370Ωm，低于理论值 - 异常体上方出现轻微反向变化（低阻体上方略高，高阻体上方略低），差异<20%
计算效率方面，在2.6GHz的6核12线程计算机上，总计算时间为4858秒。并行效率测试显示： - 1线程：平均单次迭代758秒 - 3线程：平均单次迭代232秒 - 8线程：平均单次迭代93秒
Kayabe地区数据反演经过89次迭代后收敛（RMS=2.32），在2.6GHz的2核4线程计算机上耗时4小时。主要发现： - 地下500m范围内主要表现为低阻异常 - 50-200m深度：低阻异常由东向西收敛，呈现西南和东南两个低阻异常体 - 300m以下：西南部及中部低阻异常体逐渐变小，与东南部分离
将本研究的NLCG方法与以下方法进行对比： 1. REBOCC（Weerachai，2006）：三维数值域反演 2. CG（林昌洪，2010）：共轭梯度反演 3. RRI（谭捍东等，2003）：快速松弛反演
对比结果显示： - 四种方法在100-200m深度均表现出主体低阻和周边高阻特征 - NLCG与REBOCC结果相似度最高 - CG方法结果最平滑，分辨率最低 - RRI方法对中部高阻体范围和值估计最大
本研究成功开发并验证了一种高效并行的大地电磁三维反演算法，主要贡献包括：
算法创新：改进的预处理方法显著降低了计算量，使每次迭代每个频点仅需6次正演计算，同时降低了对初始模型的依赖性。
计算效率：通过OpenMP实现的频点级并行计算，使普通PC机也能完成大规模三维反演计算。测试显示8线程下计算速度提升约8倍。
实用价值：理论模型和实测数据验证表明，该方法能有效识别地下电性结构异常，为复杂地质条件下的资源勘探提供了可靠工具。
方法优势：相比传统方法，本算法在保证精度的同时，大幅降低了对计算机硬件的要求，有利于三维MT反演技术的推广应用。
创新的预处理方法：使用与当前迭代模型电阻率相关的参数代替Hessian矩阵计算，既保证了精度，又显著提高了计算效率。
高效的并行策略：采用频点级并行计算，避免了模型分块导致的边界问题，实现了真正的”无损失”并行加速。
实用的计算方案：通过稀疏矩阵存储和算法优化，使大规模三维反演在普通PC上成为可能，降低了技术应用门槛。
全面的验证体系：通过理论模型和实际数据双重验证，证实了方法的正确性和实用性，并与多种现有方法进行了对比分析。
本研究得到了多项国家级科研项目的支持，包括： - 国家科技专项”深部矿产资源立体探测技术与试验（SinoProbe-03）”（2008-2012） - 中国地质科学院基本科研业务费项目”大地电磁场静位移效应的机理分析和校正方法探索（K1318）”（2013-2014） - 国家自然科学基金重点项目”长江中下游成矿带深部动力学背景与成矿作用研究（40930418）”（2010-2013）
这些支持反映了该研究的学术价值和应用前景，也体现了其在深部资源勘探领域的重要意义。