本文由吉林大学地球探测科学与技术学院的徐凯军、李桐林、张辉和李建平共同完成，发表于《西北地震学报》（Northwestern Seismological Journal）2006年第28卷第2期。
大地电磁测深法（Magnetotelluric, MT）是研究地球电性结构的重要地球物理方法。随着油气勘探深度和难度的增加，该方法已成为解决复杂地质问题的重要工具。虽然大地电磁测深的二维正反演已趋于成熟，但三维模拟仍面临挑战。本研究针对均匀导电半空间中的三维大地电磁响应，采用积分方程法（Integral Equation Method）进行数值模拟，旨在为研究三维地电构造的大地电磁响应分布规律提供有效工具，并为三维反演算法奠定基础。
在均匀导电半空间含三维异常体的模型中，研究者建立了电磁响应的积分方程。通过Maxwell方程和电磁张量格林函数（Tensor Green’s Function），得到了描述均匀大地中三维异常体电磁响应的关键方程：
电场积分方程： e® = e_p® + ∫v δσ·g̃_e(r,r’)·e(r’)dv’
磁场积分方程： h® = h_p® + ∫v δσ·g̃_h(r,r’)·e(r’)dv’
其中，e®和h®为总场，e_p®和h_p®为一次场，δσ为电导率差异，g̃_e和g̃_h为并矢格林函数。
研究团队将积分方程转化为矩阵方程进行数值计算，具体步骤包括：
1) 模型剖分：将异常体区域剖分为n个小立方体单元，每个单元内的场值取其中心点值。为提高计算精度，采用二次剖分算法解决奇异值问题，每个小单元再细分为8个次单元。
2) 张量格林函数计算： - 电场张量格林函数分为一次部分（pγ̃_e）和二次部分（sγ̃_e） - 一次部分采用解析近似方法求解电流项，差分近似方法计算电荷项 - 二次部分采用结合连分式展开的高斯求积法处理含贝塞尔函数的积分项，替代传统的快速汉克尔变换
3) 磁场张量格林函数计算： - 基于麦克斯韦定理，通过电场张量格林函数的旋度求得 - 一次部分（pγ̃_h）采用三重积分直接计算 - 二次部分（sγ̃_h）先对电张量格林函数求旋度再积分
4) 方程组求解：采用全选主元高斯-约当消去法求解大型线性复系数方程组，得到各剖分单元的总电场值，进而计算总磁场值。
5) 视电阻率计算：通过公式ρ_ij=1/(ωμ)|E_i/H_j|^2计算视电阻率，其中i,j=x,y。
研究团队设计了一个二维低阻棱柱体模型（x方向4km，z方向2km，顶面埋深2km，电阻率5Ωm，围岩100Ωm），分别用二维有限元法和本文三维程序计算。结果显示在1.0Hz和0.1Hz频率下TE模式的两种方法计算结果高度吻合，验证了算法的正确性。
1) 单一低阻体模型： - 半空间电阻率100Ωm - 低阻体尺寸1km×4km×2km，电阻率5Ωm - 按200m线性尺寸剖分为1000个单元，二次剖分100m - 获得了不同频率下x轴方向地表中心剖面的视电阻率ρ_yx分布
2) 高低阻体共存模型： - 半空间电阻率100Ωm - 高阻体300Ωm和低阻体5Ωm共存，y方向延伸4km - 计算了不同频率下x轴方向地表中心剖面的视电阻率ρ_xy - 结果显示视电阻率能反映地下异常体分布，且异常幅度随频率升高而减小，逐渐接近背景场电阻率
方法优势：积分方程法只需对异常体剖分和求积，相比有限元和有限差分法具有计算速度快、物理意义明确、内存占用少等优点，特别适合有限大小三维体的电磁响应模拟。
技术创新：
应用价值：
算法创新：首次将二次剖分算法与结合连分式展开的高斯求积法相结合，解决了积分方程法中张量格林函数计算的精度和奇异值问题。
验证全面：通过二维模型对比和多种三维模型试算，系统验证了算法的正确性和适用性。
实用性强：针对实际勘探中常见的高低阻体共存情况进行了模拟，结果对野外数据解释具有指导意义。
计算效率：通过优化积分计算方法和方程组求解策略，在保证精度的同时提高了计算效率。
研究团队详细讨论了张量格林函数计算这一关键问题，比较了不同方法的优缺点。特别指出对于含有贝塞尔函数的积分项，传统快速汉克尔变换方法不适用，因为积分核函数不总是随贝塞尔函数参数增大而收敛。这一发现对后续研究具有重要参考价值。
此外，论文还提供了完整的数学模型推导过程和张量格林函数的具体表达式，为读者复现研究结果提供了充分依据。研究结果展示了频率变化对视电阻率响应的影响规律，这对实际勘探中选择合适的工作频率具有指导意义。