1. 項目概述從二維到三維的跨越搞地球物理勘探的同行尤其是做電磁法的對“瞬變電磁法”這個名字肯定不陌生。它就像給大地做CT通過向地下發射一個短暫的脈沖電流然后“聽”大地回應的電磁信號來推斷地下幾米到上千米深度的電性結構。過去十幾年二維正反演技術已經相當成熟成了礦產勘查、水文地質調查的常規武器。但現實世界是三維的當地下構造稍微復雜一點比如遇到傾斜的礦體、不規則的采空區或者城市里縱橫交錯的管線二維模型那“切片式”的簡化就開始捉襟見肘解釋結果往往失真甚至漏掉關鍵信息。我這次要聊的就是我們從二維舒適區跳出來啃下“瞬變電磁三維正演”這塊硬骨頭的研發歷程。正演是整個反演系統的基石。簡單說就是給定一個地下三維電性模型計算出地面觀測點會接收到什么樣的瞬變電磁響應。這聽起來像一道數學物理題但真正做起來你會發現它是一場對計算資源、算法穩定性和工程實現能力的極限挑戰。為什么非得做三維因為精度需求倒逼。現在項目甲方要的不再是“大概有個異常”而是“異常體的頂底板埋深、傾向、傾角、規模到底是多少”二維解釋給不了這個答案。三維正演算得準后續的反演才能靠譜最終的報告才有說服力。2. 核心思路與技術選型在精度與效率間走鋼絲開發一個實用的三維正演引擎本質上是在求解麥克斯韋方程組在時間域下的數值解。這中間有無數個岔路口每個選擇都直接關系到最終程序是“玩具”還是“工具”。2.1 控制方程與離散化方法有限差分法的務實之選時間域電磁場的擴散過程通常用電場擴散方程來描述。我們放棄了雖然精度高但網格適應性差、內存消耗巨大的有限元法選擇了在矩形網格上更易實現的時域有限差分法。FDTD的核心思想很簡單把連續的空間和時間用網格劃分開用差分近似微分。但魔鬼在細節里。對于瞬變電磁這種源突然關閉后場隨時間衰減的問題直接顯式時間推進雖然簡單但穩定性條件苛刻時間步長被網格最小尺寸死死限制計算會慢到無法忍受。我們采用了Du Fort-Frankel格式的一種改進方案來處理時間導數。這是一種顯式格式但通過巧妙的中心差分構造獲得了更好的穩定性。當然它也不是無條件穩定但相比經典顯式格式允許我們采用更大的時間步長。在離散化時我們采用Yee氏網格來交錯放置電場和磁場分量這樣天然滿足法拉第定律和安培定律的離散形式保證了算法的物理基礎扎實。注意網格剖分是第一個大坑。為了用有限的網格去模擬無限的地下空間必須在模型區域外圍添加足夠厚的吸收邊界層。我們試過PML但發現在晚期衰減信號很弱時PML邊界容易產生數值反射污染計算結果。后來換成了擴展網格加指數拉伸坐標的方法簡單粗暴但有效通過逐漸拉大外圍網格的尺寸讓場在邊界處“自然衰減”到近乎為零雖然多算了一些網格點但換來了晚期道數據的純凈。2.2 激發源的處理從“理想”到“接地”很多教科書和開源代碼喜歡用水平電偶極子作為發射源因為它公式簡潔計算方便。但在實際勘探中尤其是深部找礦或工程探測我們大量使用的是接地長導線源比如長達上千米的導線兩端接地或者大定回線源。這兩種源的場分布與電偶極子源有顯著差異尤其在近區。如果正演程序只用偶極子源那么你辛辛苦苦開發出來的引擎根本無法直接擬合野外實測數據因為“源”不對。我們必須實現接地線源的模擬。這里的技巧在于將長導線離散成一系列首尾相連的電偶極子然后進行疊加。但這會帶來巨大的計算量因為每個偶極子都要計算一次場然后累加。我們通過離散波數變換技術將空間域的疊加轉換到波數域進行利用FFT加速最終將計算復雜度從O(N2)降到了O(N log N)這才讓長導線源的正演在普通工作站上變得可行。2.3 并行計算策略擁抱多核與GPU一個中等規模的三維模型比如100x100x50個網格一次正演計算可能就需要幾個小時甚至幾天。不做并行化這程序就沒有實用價值。我們的并行化是分層級的頻率并行對于頻率域轉換法后面會提到不同頻率點的計算是完全獨立的可以完美并行。我們直接用OpenMP在CPU多核上并行循環這是“免費的午餐”效率提升立竿見影。網格區域分解對于單個頻率點下的大型線性方程組求解這是最耗時的部分我們采用了基于MPI的區域分解。將整個三維網格在空間上切割成若干個子區域分給不同的CPU進程計算進程間通過邊界交換數據。這里最大的挑戰是負載均衡如果地下模型電性反差大有的區域收斂快有的慢會導致一些進程早早就閑著等別人。GPU加速嘗試線性方程組求解的核心是大型稀疏矩陣的運算這正是GPU的強項。我們嘗試將最耗時的預條件共軛梯度法求解器移植到CUDA上。初期效果并不理想因為數據在CPU和GPU之間來回搬運的開銷抵消了GPU的計算優勢。后來我們重新設計了數據流盡可能讓整個求解流程待在GPU上只在一頭一尾進行數據傳輸終于獲得了3-5倍的加速比。但這部分代碼對硬件依賴強我們將其作為可選模塊供擁有高性能GPU顯卡的用戶使用。3. 核心算法實現細節穿越“數值不穩定”的雷區有了理論框架真正編碼實現時才是問題集中爆發的階段。三維正演不是把公式翻譯成代碼那么簡單它需要大量的“工程調優”。3.1 時間域與頻率域方法的抉擇直接時間步進求解時域方程直觀但就像前面說的穩定性是噩夢。我們走的是更主流的路線頻率域轉換法。即先在頻率域計算多個頻點的諧變場響應然后通過正弦或余弦變換合成出時間域的瞬變響應。這個方法優勢明顯頻率域方程是Helmholtz方程形式更簡單且每個頻率點獨立計算易于并行。但難點在于頻點選擇選多少頻點選哪些頻點選少了變換回時間域時精度不夠特別是早期和晚期信號失真選多了計算量劇增。我們采用了一種自適應頻點選取算法根據觀測時間窗口和大地電導率范圍動態確定需要計算的頻率范圍和采樣密度在保證精度的前提下將頻點數減少了約30%。數值積分從頻率域變換到時間域需要進行一個從零到無窮的積分。我們用的是數字濾波法特別是Guptasarma和Singh的濾波系數。這里的關鍵是濾波系數的長度和采樣間隔。系數太短精度差太長計算慢且可能引入震蕩。我們通過大量模型測試固化了一套適用于一般地電條件的濾波參數。3.2 大型稀疏線性方程組的求解這是整個正演計算的心臟也是最耗時的部分。頻率域下每個頻點的計算最終都歸結為求解一個形如Ax b的大型稀疏復線性方程組。其中A是系統矩陣規模可達幾十萬甚至上百萬階。迭代法 vs 直接法直接法如LU分解求解穩定一次分解后可快速求解多個右端項對應多個發射源位置。但對于三維問題矩陣A的填充元會爆炸式增長內存根本吃不消。我們毫無懸念地選擇了迭代法特別是Krylov子空間迭代法。預條件器的藝術迭代法收斂快慢幾乎完全取決于預條件器的好壞。一個糟糕的預條件器迭代可能幾千步都不收斂。我們嘗試了多種雅可比對角預條件最簡單但效果一般尤其當模型電性反差大時如圍巖和礦體對角線元素差異巨大效果很差。不完全LU分解效果很好能顯著加速收斂但構造ILU本身也有計算和存儲開銷。多重網格預條件這是我們最終采用的方案。它的思想非常巧妙在細網格上難以平滑掉的誤差轉移到粗網格上會變得很“光滑”容易消除然后再將修正量傳回細網格。我們實現了一個幾何多重網格預條件器與穩定雙共軛梯度法迭代器搭配。實測表明對于我們的問題它比ILU預條件快2-4倍而且內存占用更可控。收斂準則設置迭代什么時候停止不能光看殘差范數是否小于某個閾值。因為晚期道的信號幅值可能比早期道小6-8個數量級如果統一用絕對殘差會導致早期道算得不夠準而晚期道又過度計算。我們采用了相對殘差和絕對殘差相結合的自適應準則確保不同時間道的計算都達到合理的精度水平。3.3 場值計算與觀測系統模擬解出網格節點上的電場或磁場值只是第一步。我們實際觀測的是特定位置、特定方向上的磁場隨時間的變化率dB/dt或者感應電動勢。這需要從網格值進行插值。磁場計算在FDTD的Yee網格上磁場天然定義在網格棱邊中心。但如果接收點不在網格棱邊上就需要插值。我們采用線性插值從最近的幾個棱邊磁場值加權平均得到接收點處的磁場。對于回線源內部的點還需要對磁場進行面積分這又涉及到數值積分精度的控制。時間導數計算我們觀測的是dB/dt。直接從頻率域變換得到的就是時間域場值B(t)然后數值求導。數值求導是個噪聲放大器特別是對晚期低幅值信號。我們采用了三點中心差分結合平滑濾波的方法來求導在犧牲一點點時間分辨率的前提下換取了信號的信噪比。多分量與多裝置一個實用的正演程序必須能模擬各種裝置中心回線、重疊回線、偶極-偶極、動源等。同時要能計算多分量響應Hz, Hx, Hy。我們在程序架構設計時就將“發射源”和“接收器”抽象成獨立的模塊通過配置文件驅動可以靈活組合各種觀測系統。這為后續反演中同時擬合多種裝置類型數據打下了基礎。4. 精度驗證與性能調優用已知答案檢驗未知算法程序寫出來了算得飛快但結果對嗎這是最讓人忐忑的階段。我們建立了一套多層次的驗證體系。4.1 解析解對比測試這是最硬核的檢驗。我們尋找一切有解析解或半解析解的場景來對比。均勻半空間這是基礎測試。將我們的三維程序設置成一個均勻半空間模型對比其響應與經典一維解析解如Wait公式或Kaufman公式的差異。我們要求相對誤差在早期道小于1%晚期道小于5%考慮到數值誤差累積。第一次測試失敗了晚期誤差高達20%。排查后發現是吸收邊界在晚期不夠“吸收”調整了邊界層厚度和拉伸系數后通過。層狀大地我們實現了基于漢克爾變換的一維正演作為基準。將三維程序模擬一個水平層狀模型與一維結果對比。這里主要檢驗垂直方向網格剖分是否足夠精細特別是薄層的模擬能力。簡單三維體比如地下一個立方體異常體。這類模型很少有解析解但我們找到了國外某研究機構公開發表的高精度數值解用完全不同的方法計算如積分方程法作為參照。對比結果令人鼓舞在異常體上方我們的結果與參考解吻合得非常好但在異常體邊界附近由于網格離散的階梯效應存在一些局部差異這在預期之內。4.2 收斂性分析數值方法必須證明其收斂性。我們設計實驗對一個固定模型逐步加密網格比如從粗網格20x20x10加密到40x40x20再到80x80x40觀察計算結果的變化。當網格加密到一定程度后計算結果的變化應小于某個閾值比如1%此時可以認為解已經收斂。我們繪制了誤差隨網格尺寸變化的對數圖其斜率應該與理論收斂階一致比如二階方法斜率應接近-2。這個測試不僅驗證了程序正確性也幫助我們確定了針對不同勘探深度和精度要求應該使用多大的網格密度避免了盲目加密網格帶來的計算浪費。4.3 實際數據擬合“實戰”是最終的試金石。我們選取了幾處地質情況相對清楚、已有鉆探驗證的礦區老資料。用我們的三維正演程序根據已知地質信息構建初始三維電性模型計算理論響應然后與當年的實測瞬變電磁曲線進行對比。 這個過程極其痛苦因為野外數據包含各種噪聲人文干擾、天然電磁場噪聲等而我們的正演是“純凈”的。我們需要在正演程序中加入一個簡單的噪聲模型比如早期道加百分比噪聲晚期道加固定電平噪聲來模擬實際情況。通過反復調整模型細節如礦體邊界微調、圍巖電導率微調使得理論曲線與實測曲線在形態、幅值、衰減趨勢上達到最佳擬合。當看到那些起伏的野外曲線能被我們的程序生成的平滑理論曲線緊緊“跟隨”時那種成就感是無與倫比的。這證明我們的程序不僅數學上正確物理上也是合理的。5. 開發中的典型問題與實戰排坑記錄寫代碼、調參數的過程就是不斷踩坑、爬坑的過程。下面記錄幾個讓我印象深刻的“坑”。5.1 晚期數據震蕩與溢出問題現象在計算某些高阻模型時時間域晚期道的響應曲線會出現非物理的高頻震蕩甚至數值溢出變成NaN。排查過程首先懷疑是時間步長太大導致顯式格式不穩定。減小時間步長后問題依舊。檢查頻率域到時間域的變換過程。發現出問題的模型其頻率域響應在低頻部分非常小且變化劇烈。數字濾波法在處理這種“陡峭”的低頻譜時容易產生吉布斯現象導致變換后時間域信號震蕩。進一步分析根本原因在于頻率域求解時對于高阻模型低頻下的波數很大導致離散化后的系統矩陣條件數變得極差迭代求解本身就不準確輸出了含有較大數值誤差的頻率域解。解決方案改進預條件器針對高阻模型強化預條件器的效果確保即使在低頻也能獲得相對準確的頻率域解。我們為多重網格預條件器增加了針對高阻區域的特殊平滑算子。增加低頻采樣點在自適應頻點選取中強制在高阻模型對應的頻段內增加采樣密度用更多的點去刻畫劇烈變化的頻譜。后處理平滑在時間域結果上對晚期道施加一個輕量的滑動平均濾波作為最后一道保險。同時在程序中加入自動檢測機制當發現晚期道數據出現劇烈震蕩時給出警告并建議用戶檢查模型電阻率設置或加密網格。5.2 并行計算中的“幽靈”數據不同步問題現象在使用MPI進行區域分解并行時程序偶爾非每次會計算出錯誤結果且錯誤每次出現在不同的網格區域像幽靈一樣隨機。排查過程這是最棘手的并發bug。首先排除了算法錯誤因為單進程運行結果始終正確。懷疑是MPI通信不同步。在每一個MPI發送/接收操作后添加了同步路障問題頻率降低但未根除。使用調試工具檢查內存。發現某個數組在迭代求解過程中一個進程本該讀取鄰居進程發來的邊界數據但偶爾讀到的卻是過時的舊值。根本原因在于我們為了減少通信開銷使用了MPI的非阻塞通信MPI_Isend,MPI_Irecv。在調用非阻塞發送后立即對發送緩沖區進行了修改準備下一輪計算而接收方進程可能尚未完成接收操作導致數據污染。或者接收方在確保數據到達前未調用MPI_Wait就使用了接收緩沖區。解決方案嚴格同步在非阻塞通信后必須在對相關緩沖區進行任何操作之前使用MPI_Wait或MPI_Test確保通信完成。我們重新梳理了所有通信邏輯繪制了通信依賴圖。引入通信緩沖區副本對于需要重復使用的發送數據不再直接操作原數組而是先拷貝到一個專用的通信緩沖區再發送這個緩沖區。雖然增加了一次內存拷貝但徹底杜絕了數據競爭。壓力測試編寫了專門的測試用例用數百個進程對復雜模型進行反復計算確保并發bug完全消失。5.3 復雜地形模擬失真問題現象當模型包含劇烈起伏的地形如山谷、山脊時在地形突變處附近的計算點響應曲線出現明顯的畸變與物理直覺不符。排查過程我們的初始網格是規則的笛卡爾網格地形是通過一種“階梯”狀網格來逼近的。即將地表以上的空氣層網格電阻率設為極高的值如1e8 Ω·m地表以下的網格按真實電阻率賦值。當地形起伏時這個“空氣-大地”界面在網格中就是鋸齒狀的。問題根源這種階梯近似在平坦地形時還好但當地形陡峭時會產生兩個問題1) 網格的階梯狀邊界會人為引入虛假的電荷積累影響電場計算2) 在地形頂點或谷底網格尺寸可能無法精細刻畫曲率導致場值計算不準確。解決方案地形網格拉伸我們引入了非結構化網格的前處理思路但仍在FDTD框架內實現。具體做法是在保持網格邏輯結構仍是矩形的前提下對靠近地表的網格層進行坐標變換將物理坐標系中起伏的地形映射到計算坐標系中的一個水平面。這樣在計算坐標系中網格是規則的所有差分公式照常使用只需在雅可比矩陣中考慮坐標變換的系數。這個方法被稱為共形網格或地形擬合網格技術。空氣層處理優化對于空氣層我們不再簡單地設為高阻而是精確賦予其真實電阻率~1e14 Ω·m和介電常數。同時確保空氣層足夠厚以避免邊界反射影響地表附近的場。對于地形劇烈區域局部加密水平方向的網格。效果驗證我們用一個傾斜山坡下的低阻體模型測試。改進后山坡上下的響應曲線過渡自然異常形態清晰與基于有限元法的商業軟件結果一致性很好。6. 工程化封裝與用戶體驗優化一個強大的計算內核還需要一個好用的外殼才能交付給解釋人員使用。我們花了很大力氣在工程化上。6.1 輸入輸出接口設計輸入文件我們采用了JSON格式而不再是傳統的、難以閱讀和修改的卡片式文本。JSON結構清晰支持嵌套可以方便地描述復雜的三維模型、觀測系統參數。我們還為JSON配置文件編寫了詳細的模式定義用戶編輯時能有語法提示和錯誤檢查。 輸出方面除了標準的文本格式每個測點每條曲線一個文件我們主要支持NetCDF和HDF5這兩種科學數據格式。它們支持并行讀寫可以高效地存儲大規模的三維場值數據、模型參數和元數據。用戶可以用PythonnetCDF4, h5py庫或MATLAB輕松讀取和可視化。6.2 可視化與調試工具“黑盒”程序沒人敢用。我們開發了一系列配套工具模型可視化器可以三維顯示網格剖分、電阻率模型分布、發射源和接收點位置。讓用戶一目了然地檢查模型設置是否正確。正演結果瀏覽器可以按測點、按分量、按時間道快速繪制理論曲線并支持與實測曲線的疊合對比。可以方便地切換線性坐標和對數坐標。運行時監控程序運行時會輸出關鍵的迭代信息如殘差下降曲線、計算耗時等。對于大型任務我們還提供了一個簡單的Web監控頁面可以遠程查看計算進度和資源占用情況。6.3 性能剖析與調優指南我們使用gprof和Intel VTune等工具對程序進行了深度性能剖析。發現超過70%的時間花在了線性方程組求解的矩陣-向量乘法和預條件器應用上。基于此我們給出了針對不同硬件平臺的編譯和運行建議CPU平臺建議使用Intel編譯器搭配MKL數學庫并設置合適的OpenMP線程數通常等于物理核心數。混合平臺如果使用GPU加速建議將最耗時的單個頻率點計算特別是大型模型分配給GPU而多個頻率點的并行任務由CPU承擔形成流水線。內存優化我們提供了“內存優化模式”選項該模式下會使用單精度浮點數進行計算而非默認的雙精度并將一些不頻繁訪問的中間變量交換到磁盤。這可以將內存占用量減少近一半計算速度也有提升代價是損失一些精度適用于對精度要求不是極端高的快速模擬場景。開發三維正演的過程是一個不斷在理論理想與工程現實之間尋找平衡點的過程。沒有一種算法是完美的關鍵是要深刻理解每種方法背后的假設和局限然后針對我們要解決的具體問題瞬變電磁法做出最務實的選擇和優化。它不僅僅是一個數學物理問題的求解器更是一個融合了高性能計算、軟件工程和地球物理專業知識的復雜系統。當看到它成功模擬出復雜地質模型產生的電磁響應并與野外數據吻合時你會覺得所有在深夜調試bug、在集群前等待結果的時間都是值得的。這個正演引擎成為了我們后續攻克更艱難的三維反演問題的堅實跳板。