在現今的地球物理勘探領域，瞬變電磁（TEM）方法因其對地下導電結構敏感而被廣泛應用于資源勘查、環境調查和工程探測。然而，無論是商業巨頭提供的專業軟件，還是科研社區用Python、MATLAB等現代語言開發的各類程序，都面臨著一些共同的“成長的煩惱”。一方面，商業軟件往往被費用、使用權和修改限制所“捆綁”；另一方面，雖然Python等語言憑借其強大的內置函數和高效的矩陣操作，非常適合編寫反演程序的主體，但其解釋執行的特性使得被反復調用的核心正演計算子程序成為整個程序執行速度的瓶頸。想象一下，在需要進行數十萬次正演計算的馬爾可夫鏈蒙特卡洛（McMC）反演中，或者在需要快速生成海量訓練集的人工智能網絡訓練中，一個“慢吞吞”的正演模塊將是多么令人頭疼。此外，隨著儀器數字化的發展，波形記錄方式日新月異，對儀器系統響應（即儀器自身特性對觀測數據的影響）進行精確建模的需求也日益迫切。現有的工具是否足以靈活、高效且精準地應對這些挑戰？為了克服計算速度與建模靈活性之間的矛盾，并為社區提供一個可靠、通用的計算核心，來自丹麥奧胡斯大學地球科學系的Niels B. Christensen開發并開源了一個FORTRAN子程序。

這項研究旨在提供一個開源、高效且功能全面的FORTRAN子程序，用于計算一維（1D）地球模型的TEM響應及其關于模型參數的導數。該子程序支持常見的儀器配置（如圓形和多邊形發射回線），能夠完整建模系統響應（包括接收機和放大器的帶寬限制、重復波形效應以及實際發射波形的卷積），并且可以選擇性地將基于Cole-Cole模型的激發極化（IP）效應納入正演計算。代碼設計為可被Python（通過f2py）、MATLAB等高級語言調用，或通過文件進行“傻瓜式”集成的外部編譯模塊，從而在保持高級語言易于進行反演算法開發的優勢的同時，極大提升核心正演計算的效率。該研究成果以論文形式發表在《Computers 》期刊上。

本研究主要應用了幾個關鍵的數值計算技術：1. 場方程求解：基于準靜態近似，推導了頻率域中垂直磁偶極子源在層狀大地模型上產生的磁場表達式，并通過漢克爾變換（Hankel Transform）將其轉換到空間域。研究中采用了快速漢克爾變換（FHT）方法，并集成了相應的開源濾波器系數計算例程。2. 時域轉換：采用基于Gaver-Stehfest算法的拉普拉斯逆變換將頻率域響應轉換到時域，該方法相比傅里葉變換在早期延遲時間更準確，且為純實數運算，速度更快。3. 系統響應建模：在拉普拉斯域中，將低通濾波器（一階或二階Butterworth濾波器）的傳遞函數與場方程的核函數相乘，來模擬儀器帶寬限制。通過時間偏移和疊加來模擬重復波形的影響。將實際發射波形視為分段線性，并通過其二次導數與階躍響應的卷積來高效計算波形效應。4. 導數計算：通過遞歸公式對核函數進行解析求導，從而得到響應關于層電導率和厚度的導數，確保了導數與正演響應具有相同的計算精度。5. IP效應集成：在計算中采用Cole-Cole模型來描述電阻率的頻散特性，通過修改層狀模型的波數遞歸關系，將IP參數（零頻電阻率ρ₀、充電率m、時間常數τ和頻率依賴系數c）集成到正演響應計算中。

研究明確了子程序的設計邊界。它支持對IP效應進行正演模擬（但暫不計算其對IP參數的導數），將發射波形視為分段線性處理，目前僅計算垂直磁場分量（因這是1D反演中最常用的），不包含近似響應選項，每次僅處理一個發射矩，并且不內置對時間門（gate）的積分，而是輸出密集采樣的響應由用戶自行處理，以保持靈活性和通用性。

子程序支持大多數常用配置：發射源（Tx）可以是圓形回線或多邊形分段線性回線，接收器（Rx）始終模擬為磁偶極子。為適應儀器發展，還允許包含兩個Tx和兩個Rx（由用戶選擇極性），盡管在當前版本中，所有Tx-Rx對的橫向距離必須相同。

論文詳細給出了垂直磁偶極子和水平電偶極子源在1D模型上產生的垂直磁場頻率域積分表達式。在準靜態近似下，通過FHT和Gaver-Stehfest拉普拉斯逆變換，將其數值計算為時域響應。程序可輸出三種類型的響應：基本階躍響應、基本脈沖響應以及包含全部系統響應元素的卷積響應。

研究完整實現了對儀器系統響應的建模，包括：1) 接收機和放大器/記錄系統的帶寬限制（用低通濾波器模擬）；2) 重復波形的影響（通過疊加當前及前幾個周期的響應實現）；3) 單個測量周期發射波形的效應（通過波形二次導數與校正后響應的卷積計算）。圖1通過一個三層模型的實例，直觀展示了逐步加入重復性、低通濾波和波形卷積等系統響應元素后，脈沖響應曲線的變化，強調了完整系統響應建模對于獲得可與實測數據直接對比的理論值的重要性。

子程序允許在正演響應中引入基于Cole-Cole模型的IP效應。電阻率表示為ρ = ρ₀· [1 - m· (1 - 1/(1 + (sτ)^c))]，其中s為拉普拉斯變量。圖2展示了在一個均勻半空間模型中，包含IP效應（充電率0.3、時間常數0.1、頻率系數0.5）與不包含IP效應時，基本階躍響應的顯著差異，特別是在晚期時間出現負值，這是IP效應的典型特征。

通過將子程序計算結果與均勻半空間模型的解析解進行對比，驗證了其精度。測試了中心回線（Tx半徑10 m）和偏移磁偶極子（Tx-Rx距15 m）兩種配置。如圖3所示，在從早期到晚期數個時間數量級的范圍內，數值解與解析解吻合得非常好，證明了子程序計算的準確性。

本研究成功開發并發布了一個開源、高效、功能全面的FORTRAN子程序，用于計算1D地球模型的TEM響應及其導數。該程序的核心優勢在于其速度（作為編譯代碼，特別適合被高階語言反演程序反復調用）、精度（采用FHT和Gaver-Stehfest算法確保數值穩定性）和完整性（全面支持系統響應和IP效應建模）。它有效填補了當前TEM建模工具在計算效率與建模靈活性之間的空白，為地球物理反演（特別是計算密集的McMC反演）、機器學習訓練以及開源軟件開發提供了強大的基礎計算引擎。通過將代碼以FORTRAN77源代碼形式在GitHub（https://github.com/hydrogeophysicsgroup/TEM1D）公開，并遵循開放科學的精神，本研究鼓勵全球用戶使用、修改并改進該代碼，旨在促進TEM數據解釋工具的普及與發展，特別是在資源有限的國家和地區。這項工作不僅為TEM社區貢獻了一個高質量的計算模塊，也體現了利用編譯語言優化高性能計算環節、同時保留高級語言快速開發優勢的現代科學計算理念。