李展辉、黄清华（北京大学地球与空间科学学院地球物理学系，造山带与地壳演化教育部重点实验室）与王彦宾合作的研究论文《三维错格时域伪谱法在频散介质井中雷达模拟中的应用》发表于《地球物理学报》。该研究针对井中雷达（Borehole Radar, BHR）在低频段（如10 MHz以下）探测时面临的频散介质模拟难题，提出了一种基于交错网格时域伪谱法（Staggered-Grid Pseudo-Spectral Time-Domain, SG-PSTD）的高效数值模拟方法，显著提升了计算效率与精度。
学术背景
井中雷达是探测地下异常体（如裂隙、矿藏）的重要工具，但低频电磁波在频散介质（如含水电导地层）中传播时会产生显著的波包扩散效应。传统时域有限差分法（Finite-Difference Time-Domain, FDTD）在模拟大范围模型时需要极细的网格划分，导致计算成本高昂。伪谱法（Pseudo-Spectral Method, PSM）因每个波长所需网格数更少而具有效率优势，但常规伪谱法存在Gibbs现象（高频振荡）和界面误差问题。本研究旨在开发一种能兼顾精度与效率的三维井中雷达模拟方法，并解决频散介质中的波场计算稳定性问题。
研究方法与流程
算法设计
- 扩展坐标下的Maxwell方程组：引入完全匹配层（Perfectly Matched Layer, PML）吸收边界条件，通过复数坐标拉伸算子处理频散介质的电磁场方程。
- 频散介质建模：采用Debye模型描述介电常数的频散特性，通过分段线性递归卷积（Piecewise Linear Recursive Convolution, PLRC）方法离散化频散方程，替代传统的隐式中心差分，简化计算流程。
- 交错网格伪谱法：在空间导数计算中采用实数傅里叶变换，并通过交错网格布局（电场与磁场分量位于不同网格点）抑制Gibbs现象。时间步进采用二阶显式Adams-Bashforth（AB2）方法，其稳定性在低电导率介质中已得到验证。
数值实现
- 激励源模拟：以半波偶极子天线为发射源，电流脉冲采用截断三正弦波形，避免数值模拟中的电荷积累问题。
- 稳定性控制：通过分析一维非频散模型的CFL条件（Courant-Friedrichs-Lewy），确定时间步长上限（如Δ𝑡 ≤ 0.9Δ𝑥/√3𝑐，𝑐为波速）。
- 吸收边界优化：PML参数采用三次方渐变形式（如𝜎(𝑥) = 𝜎₀(𝑥/𝑑)³），减少界面反射。
主要结果
精度对比
- 在球形异常体模型中（半径3 m与1 m，中心频率50 MHz），SG-PSTD的模拟结果比FDTD更清晰地识别反射信号，且无尾波干扰（图3）。常规伪谱法的Gibbs现象在交错网格下完全消除（图4）。
- 裂隙模型（宽度放大至网格尺度）的反射相位反转现象与理论吻合：含水裂隙（ε𝑟=80）反射波相位与入射波相反，而干燥裂隙（ε𝑟=1）相位不变（图9），验证了算法对介电特性差异的敏感度。
效率优势
- 相同精度下，SG-PSTD的网格数比FDTD少4倍，计算内存需求降低约60%。例如，模拟50 MHz信号时，FDTD需Δ𝑥 ≤ 0.1 m，而SG-PSTD仅需Δ𝑥 ≤ 0.4 m。
结论与意义
本研究提出的SG-PSTD方法通过结合交错网格、PLRC频散处理与AB2时间迭代，实现了井中雷达三维模拟的高效性与鲁棒性。其科学价值在于：
1. 理论创新：首次将显式AB2方法应用于频散介质电磁场求解，突破了传统隐式算法的复杂性瓶颈。
2. 应用价值：为大尺度地下探测（如油气储层、地热资源）提供了可靠的数值工具，尤其适用于低频井中雷达数据的解释。
亮点与创新
- 方法学突破：通过交错网格与PLRC的结合，解决了伪谱法在非均匀介质中的精度问题。
- 工程实用性：AB2显式格式简化了代码实现，适合实际勘探场景的快速建模。
- 跨学科意义：算法框架可扩展至地震波模拟等其他波动问题。
补充发现：研究还指出，PML参数的三次方渐变形式可将边界反射衰减至入射幅度的1‰以下（图1），这一优化对低频模拟尤为重要。