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