本研究的通讯作者为吉林大学地球探测科学与技术学院的Yin Chang-Chun教授（邮箱：yinchangchun@jlu.edu.cn），第一作者为Cao Xiao-Yue博士研究生。合作单位包括加拿大纽芬兰纪念大学地球科学系。研究成果发表于《Applied Geophysics》期刊2018年第15卷第3-4期（9-12月刊），页码556-565，共包含10幅插图，DOI号为10.1007/s11770-018-0703-8。研究获得国家自然科学基金（41774125）、国家自然科学基金重点项目（41530320）、国家重点研发计划（2016YFC0303100和2017YFC0601900）以及中科院战略先导专项（XDA14020102）的资助。
大地电磁法（Magnetotelluric, MT）作为地球物理勘探的重要方法，在油气资源勘查和深部构造研究中具有广泛应用。传统三维MT反演多基于结构化网格，在模拟起伏地形和复杂构造时存在精度限制。随着计算机技术的发展，基于非结构化网格的有限元方法（Finite-Element Method, FEM）因其能精确拟合复杂几何形态的优势而备受关注。然而，大规模三维反演面临内存需求大、计算效率低的挑战。本研究旨在开发结合非结构化有限元网格和L-BFGS（Limited-memory Broyden-Fletcher-Goldfarb-Shanno）优化算法的三维MT反演方法，解决复杂地形条件下的高精度反演问题。
研究采用基于四面体单元的非结构化矢量有限元法求解电磁场控制方程。控制方程采用忽略位移电流的矢量亥姆霍兹方程：
∇×(∇×E) + iωμ₀σE = 0
其中μ₀为真空磁导率，E为电场，σ表示电导率。通过Galerkin有限元离散化，将计算域Ω划分为四面体单元，利用矢量基函数φ对控制方程进行弱形式处理，最终形成复线性方程组：
AE = B
通过求解该方程组获得计算域内的电场分布，再根据法拉第定律计算磁场，最终通过不同极化方向的场值计算阻抗张量。
反演目标函数采用Tikhonov正则化形式：
Φ(m) = Φ_d(m) + λΦ_m(m)
其中Φ_d(m)为数据拟合项，Φ_m(m)为模型光滑约束项，λ为正则化因子。数据拟合项定义为观测数据与模拟数据阻抗张量分量的差异：
Φ_d(m) = 1⁄2 Σ|(Z_n^obs - Z_n)/ε_n|²
模型光滑约束项通过模型参数与先验模型的差异构建：
Φ_m(m) = (m-m₀)^T W^T W(m-m₀)
采用L-BFGS优化算法求解该反演问题，其核心在于通过存储有限组向量对（s_k, y_k）来近似Hessian矩阵的逆，避免显式计算和存储完整的Hessian矩阵。算法流程包括： 1. 设置误差阈值ε和目标均方根误差r₀ 2. 初始化对角矩阵H_k⁰，计算缩放因子γ_k 3. 计算梯度∇Φ_k和RMS误差，满足收敛条件则停止 4. 通过递归公式更新近似Hessian逆矩阵H_k 5. 计算搜索方向d_k = -H_k ∇Φ_k 6. 采用Wolfe条件进行线性搜索确定步长αk 7. 更新模型参数m{k+1} = m_k + α_k d_k
特别地，除第一次迭代外，本研究将牛顿步长固定为1，满足充分下降条件，显著提高了计算效率。
设计两个相邻导电体（10 Ω·m）埋置于100 Ω·m半空间中的模型。反演结果显示（图2），在36次迭代后目标函数从19.60降至2.33，RMS从4.42降至1.03。恢复的异常体电阻率为10-20 Ω·m，位置与真实模型吻合良好。值得注意的是，除第一次迭代需要2次正演计算外，其余迭代仅需1次正演计算，在2.6 GHz Xeon处理器工作站上耗时5小时12分钟。
在梯形山体模型（高3 km）下测试算法性能。反演结果（图5）显示，36次迭代后目标函数从10.80降至1.93，RMS从3.28降至1.02。尽管存在显著地形影响，算法仍能准确恢复地下异常体位置和电阻率特征（10 Ω·m），计算耗时7小时58分钟。
针对更接近实际的复杂模型——倾斜30°的椭圆形矿体（10 Ω·m）位于高阻层（500 Ω·m）与基底（100 Ω·m）之间，并考虑地形影响。反演结果（图8,10）显示： - 上层高阻层反演结果为200-700 Ω·m - 基底电阻率恢复约100 Ω·m - 异常体走向和位置与真实模型高度一致 - 界面形态恢复良好 在31次迭代后，目标函数从15.40降至1.06，RMS从3.87降至1.03，耗时16小时31分钟。
本研究成功实现了基于非结构化有限元和L-BFGS算法的三维MT反演方法，得出以下重要结论：
方法学创新：非结构化网格有限元法能有效处理地形起伏和复杂构造的反演问题，无需进行地形校正即可获得可靠结果。L-BFGS算法通过避免显式计算Hessian矩阵，将内存需求降低至仅存储2×m×m个向量（m通常取3-7），显著提升了大规模反演的可行性。
计算效率突破：与传统方法相比，本算法在第一次迭代后采用固定牛顿步长（=1），每次迭代仅需1次正演计算，使计算效率提升约50%。三个模型测试表明，该方法在保证精度的同时具有稳定的收敛性。
应用价值：该方法特别适用于具有地形影响和复杂地下结构的实际勘探场景，为矿产资源勘查、深部构造研究等提供了新的技术手段。椭圆矿体模型的成功反演证明了该方法处理实际复杂地质问题的能力。
算法创新性：首次将非结构化有限元与L-BFGS算法系统结合应用于三维MT反演，解决了传统结构化网格方法处理复杂几何的局限性。
计算效率优化：通过固定牛顿步长和L-BFGS内存优化，实现了”一次迭代一次正演”的高效计算模式，为大规模三维反演提供了实用解决方案。
验证系统性：通过简单模型→地形模型→复杂地质模型的递进验证，全面证明了方法的可靠性和适应性，特别是对倾斜异常体和层状结构的恢复能力。
工程实用价值：所有测试均在普通工作站（2.6 GHz Xeon处理器）完成，最长计算时间控制在17小时内，具有较好的工程应用前景。
未来研究方向包括将该方法应用于实际野外数据反演，以及结合并行计算技术进一步提升计算效率。本研究为复杂条件下的大地电磁数据解释提供了新的方法论基础，对推动三维电磁勘探技术的发展具有重要意义。