簡介一份面向地球物理全波形反演與優化算法研究者的示例工程演示如何在SEISCOPE優化工具箱中調用截斷牛頓Truncated Newton算法進行波形反演計算。程序實現對應Metivier等人2013年發表在SIAM Journal on Scientific Computing上的經典方法適合需要復現該算法或開展反演實驗的科研與工程人員參考。包內共22個文件核心為15個Fortran 90源文件涵蓋算法主程序與子模塊另含兩個數據文件、工程配置文件及Visual Studio相關文件整體壓縮包僅31KB結構精簡便于快速閱讀與調試。該資源已有239人學習瀏覽雖是小型示例但完整呈現了截斷牛頓法在全波形反演中的落地流程有助于理解工具箱調用方式、數據組織與迭代求解細節可作為二次開發或算法對比的起點。1. 為什么全波形反演需要截斷牛頓法而不是多跑幾輪 L-BFGSFWI 的迭代優化器選擇往往比正演精度更容易被低估。很多人在拿到一個波動方程正演代碼后第一反應是用 L-BFGS 往上堆迭代次數結果在強散射或深部目標上梯度方向被幾何擴散和多重散射拖慢幾十次迭代后成本函數仍在平臺期。截斷牛頓法Truncated Newton把一個完整的牛頓步拆成外迭代和內迭代外層做線搜索、更新模型內層用共軛梯度近似求解牛頓方程全程不需要顯式組裝海森矩陣只額外要求你提供一次海森向量積。這正是 SEISCOPE OPTIMIZATION TOOLBOX 里 TRN 示例的價值。如果你想在自己的 2D/3D 全波形反演流程里換一個收斂更扎實的優化器這個壓縮包里的代碼和運行日志值得從頭到尾拆一遍。適合的讀者是已經有正演和伴隨代碼、但對優化器原理還停留在“梯度下降加個擬牛頓”的人。2. TRN.ZIP 項目解剖SEISCOPE 的算法骨架與兩個迭代記錄文件2.1 壓縮包里的工程骨架先區分可編譯與可忽略TRN.ZIP不是傳統意義上的規范源碼包更像開發者直接打包的工作目錄。TRN.sln和TRN.vfproj是 Windows 上 Visual Studio Intel Visual Fortran 的工程入口src與trn_src保存核心源碼common放公共數據結構和常量根目錄的iterate_TRN.dat與iterate_TRN_CG.dat是上一次運行留下的外層與內層迭代日志。TRN.u2d、TRN.suo都是 IDE 生成的會話文件不參與編譯可以忽略tets從命名上看是 tests 的筆誤多半是作者做回歸驗證的腳本目錄。這種布局的好處是你直接能看到程序的運行產物壞處是構建順序和依賴項必須自己判斷。我拿到這類項目后第一件事是按“可編譯、可運行、可忽略”三類重新歸檔否則很容易在無關的調試文件上浪費時間。2.2 optim_type.h優化算法和用戶代碼之間的協議optim_type.h是理解 SEISCOPE 優化框架的鑰匙它把優化算法和具體反演問題解耦算法層只認結構體不關心你的正演是用時域有限差分還是頻率域 Born 近似。簡化開來看結構體長這樣/* 從 SEISCOPE 工具箱接口風格簡化的定義 */ typedef struct optim_type { int verbose; /* 是否打印每輪日志 */ int nouter_max; /* 外層截斷牛頓最大迭代次數 */ int ncg_max; /* 內層 CG 最大迭代次數 */ double cg_tol; /* 內層 CG 相對殘差容差 */ double c1, c2; /* Wolfe 線搜索參數 */ void (*cost)(const double *m, double *fcost, void *ctx); void (*grad)(const double *m, double *g, void *ctx); void (*hess_vec)(const double *m, const double *v, double *Hv, void *ctx); } optim_type;在 SEISCOPE 實際版本里字段名和組織方式不一定完全一樣但你要準備的三樣東西不會變計算目標函數、計算梯度、給一個向量 v 返回 Hv。這個設計是截斷牛頓法能落地的關鍵。傳統牛頓法需要顯式組裝 N×N 的海森矩陣FWI 里 N 是模型參數數量2D 網格到 10^5 已經令密矩陣存儲壓力很大3D 更不可能而用二階伴隨或散射正演只需要一到兩次等效正演就能得到 Hv。TRN 外層通過線搜索保證下降內層通過 CG 控制牛頓方向精度兩者都只依賴 Hv所以這套接口能撐起大參數規模的反演。2.3 外層 iterate_TRN.dat 與內層 iterate_TRN_CG.dat 分別記錄什么這兩個文件不是下載占位符而是算法主動寫出的診斷數據比終端輸出可靠。iterate_TRN.dat記錄外層截斷牛頓迭代典型列包括迭代號、目標函數值、梯度 L2 范數、線搜索步長、本輪 CG 次數。iterate_TRN_CG.dat記錄某輪外迭代中內層 CG 的每次迭代信息常見列是 CG 迭代號、殘差范數或方向修正量。我一般先看這兩個文件不看終端因為終端會滾動且通常只打印摘要。快速驗證是否正常可以用# 查看外層目標函數與梯度范數的變化趨勢 head -5 iterate_TRN.dat tail -5 iterate_TRN.dat # 查看內層 CG 記錄有多少條是否有異常長迭代 wc -l iterate_TRN_CG.dat這段命令里head/tail是 Linux/macOS 命令Windows 上對應 PowerShell 的Get-Content -Head。目標函數應大致單調下降梯度范數的數量級應從初始值下降幾個量級wc -l統計行數用于判斷內層 CG 是否頻繁突破預期迭代上限。2.4 從 Debug 目錄和 tets 目錄反推運行習慣Debug目錄說明作者使用了 Visual Studio 默認的 Debug 配置構建根目錄出現iterate_TRN_CG.dat則說明程序的工作目錄就是工程根目錄輸入輸出都用相對路徑。復現時盡量不要修改這個約束否則程序很容易找不到模型文件。tets目錄提醒我們TRN 這類與 CG 殘差、Wolfe 條件糾纏的算法必須有幾個可重復的小模型做回歸用例。你改一行容差可能讓一個在 A 模型上收斂良好的測試在 B 模型上悄悄退化沒有回歸用例很難意識到。2.5 TRN 與 L-BFGS 的本質差異從近似逆到近似解L-BFGS 的本質是用過去的梯度差和模型差去逼近海森矩陣的逆好處是每次迭代只需要額外的向量運算壞處是逼近質量依賴歷史迭代的多樣性和問題非線性程度。TRN 則不同它在每個外層步用 CG 去直接求解牛頓方程CG 的每次迭代都調用精確的 Hv 算子因此只要問題本身海森是正定的它獲得的搜索方向比 L-BFGS 的“替換品”更接近真實牛頓方向。這也是為什么在強散射介質、長偏移距數據下TRN 能把 FWI 從幾十次迭代壓到十幾次而 L-BFGS 常常卡在梯度平臺區。代價是每次外層迭代的正演次數顯著增加所以它不是替代 L-BFGS而是在預算允許時更高階的選擇。3. 從 TRN.sln 到反演結果Windows 構建、跨平臺編譯與迭代日志判讀3.1 在 Windows 上把 .sln 跑通的關鍵不是 VS而是 Intel FortranTRN.vfproj是 Intel Visual Fortran 工程文件打開TRN.sln后如果項目顯示為空通常是因為本機沒有安裝 Intel Visual Fortran 的 VS 集成。安裝 Intel oneAPI 時選擇 Fortran 編譯器與 Visual Studio 集成組件即可。構建時 Debug/Release 的選擇也影響性能Release 下編譯優化不會被調試符號拖慢FWI 正演熱區在 Debug 下可能慢一個量級。若出現無法解析的外部符號先檢查項目屬性的庫目錄是否指向 MKL 的對應架構路徑。命令行構建可以這樣call %ONEAPI_ROOT%\compiler\latest\env\vars.bat intel64 msbuild TRN.sln /p:ConfigurationRelease /p:Platformx64%ONEAPI_ROOT%是你的 oneAPI 安裝根目錄/p:Configuration選擇構建配置/p:Platform要與本機架構及 MKL 庫一致。注意vars.bat后需要加intel64或ia32這取決于終端架構如果混用編譯鏈接階段往往會出現環境變量沖突。3.2 跨平臺手動編譯兜底方案如果只有 Linux 或者不想依賴 Visual StudioFortran 源碼本身可以跨平臺只差一份構建規則。常見做法是把源碼列表拉出來然后直接編譯find common src trn_src \( -name *.f90 -o -name *.F90 -o -name *.f -o -name *.F \) sources.txt ifort -O3 -I common -I src -mkl -o trn_run $(cat sources.txt)find命令收集常見 Fortran 擴展名文件到sources.txt括號需要轉義防止被 shell 解釋-I common -I src把模塊目錄和公共頭文件目錄交給編譯器-mkl是 Intel 編譯器連接 MKL 的快捷開關gfortran 下不存在需要改成-lblas -llapack$(cat sources.txt)把文件列表展開成多個源文件參數。如果編譯報錯說找不到某個模塊多半是-I路徑順序不對Fortran 模塊依賴要求編譯順序按照 use 關系從底層往上排。最穩妥的做法是讓編譯器自動加載.mod文件而不是手動調源文件順序。3.3 主程序流程TRN 掉進優化器之前要準備好什么從src/trn_src的代碼組織看主程序并不是把優化器函數一股腦塞進去而是先準備好四個輸入初始模型、觀測數據、目標函數接口、梯度接口。簡化流程如下! 主程序流程示意 call read_model(m0) ! 初始速度/參數模型 call load_observed_data(dobs) ! 觀測波形數據 call misfit_initial(m0, fcost, grad) ! 初始誤差和伴隨梯度 optim_type%cost cost_wrapper optim_type%grad grad_wrapper optim_type%hess_vec hess_vec_wrapper call TRN_run(m0, fcost, grad, optim_type)這里cost_wrapper、grad_wrapper和hess_vec_wrapper是你自己的 Fortran 函數TRN 求解器通過函數指針回調。注意hess_vec_wrapper接收的擾動向量 v 來自 CG 迭代內容沒有任何物理含義可能包含棋盤格噪聲如果你的正演代碼對 v 做了平滑或截斷就會返回一個不一致的 Hv 給 CG。這是很多 TRN 實現跑起來比理論慢的原因——Hv 沒按“未處理的原向量”計算。3.4 跑完以后先看數據而不是看終端程序退出碼正常不代表反演正確。我通常立刻對比兩個迭代文件# 外層是否單調下降 awk {print $1, $2, $3} iterate_TRN.dat | head -30 # 內層 CG 殘差是否降到設定容差以下 awk $NF 1e-6 {print NR, $NF} iterate_TRN_CG.dat | head第一條awk命令取前三列分別對應迭代號、目標函數值、梯度范數如果第三列的數量級在下降說明外層搜索方向是有效的。第二條命令把內層 CG 記錄中最后一列小于 1e-6 的行打印出來用來確認內層迭代確實在收斂而不是被ncg_max硬截斷。如果打印為空說明 CG 基本沒達到默認容差需要去檢查預條件和 Hessian 向量積。注意不同版本的程序寫出文件的列數和排列可能不同先用head -1看表頭再決定 awk 的列號不要照搬這里的列索引。4. 調參截斷牛頓線搜索、CG 容差與 Hv 的一致性4.1 外層線搜索Wolfe 條件怎么給才算不過分保守TRN 的線搜索目標不是精確找到極小點而是找到一個滿足充分下降和曲率條件的步長即 Wolfe 條件。SEISCOPE 示例里的默認參數通常取 c11e-4c20.7 到 0.9。c2 太小會讓線搜索把每一步試滿成本飆升c2 太大則容易接受過大的步長導致下一輪 CG 起點變差。判斷線搜索狀態的指標是目標函數曲線中是否出現“連續小步長密集下降”。如果步長序列大量落在 0.01 以下但目標函數還在緩慢下降這通常不是 c1、c2 需要調而是梯度方向包含嚴重的數值噪聲這時要檢查正演是否用了足夠的模板精度或吸收邊界而不是繼續壓線搜索參數。4.2 內層 CG 容差跟著外層梯度走不要設死值截斷牛頓的“截斷”就是指內層 CG 沒有真正解到機器精度。在 L. Metivier、R. Brossier、J. Virieux、S. Operto 于 2013 年發表在 SIAM Journal on Scientific Computing 上的論文中核心討論之一就是如何控制內層 CG 的殘差。常見做法是采用 Eisenstat-Walker 策略初始給一個較大容差隨著外層梯度范數減小逐步收緊。用偽代碼表達就是double eta 0.5; // 初始容差 if (outer_iter 0) { eta fmin(0.5, sqrt(grad_norm / grad_norm_initial)); eta fmax(eta, 0.1); // 下限保護 } cg_solve(..., eta);這里的eta是內層 CG 的相對殘差目標。開始迭代時模型距離最優解很遠方向稍微粗糙點也能接受所以容差可以放到 0.1~0.3等到外層梯度縮小到初始的百分之一如果還保持 0.3 的容差CG 給出的方向可能離牛頓方向偏移很大所以要讓容差按梯度下降比例收縮。fmax下限 0.1 是為了避免在數值噪聲主導時把 CG 逼到無意義的精細求解。我一般在這個框架基礎上保留ncg_max作為硬上限因為即使容差很小正演算子的舍入誤差也會限制 CG 的實際可達精度。4.3 預條件與模型維度TRN 不是免預條件的很多用戶把 TRN 直接套用到時間域 FWI完全不加預條件然后發現內層 CG 永遠在極限迭代算法表現比 L-BFGS 還差。問題不在截斷牛頓而在于海森條件數太差。常見做法是給內層 CG 加一個對角預條件子用模型空間振幅歸一化或平滑算子近似海森對角項。要注意模型網格超過 10^7 量級時每次 Hv 至少需要一次正演一次伴隨3D 時間域代價仍然很高。因此 TRN 最適合的區間是“一次正演成本適中、但梯度方向需要更準確曲率”的中間規模問題。在超大網格上我寧可把 CG 內迭代次數壓到 5 次以下讓 TRN 退化成一種含自適應正則的梯度法也比硬上完整牛頓方向劃算。4.4 常見失敗模式與排查順序現象可能原因優先排查外層 cost 上升且步長趨近 0梯度或 Hv 不一致梯度有限差分校驗內層 CG 殘差長期不變Hv 實現與梯度不同源單獨測試 Hvcost 下降但梯度范數不降正則項壓掉了梯度信息檢查目標函數歸一化單輪時間幾乎全在 CGcg_tol 過嚴或沒有預條件放寬容差、加對角預條件換網格尺寸后收斂性突變數據殘差未按網格體積縮放檢查觀測誤差歸一化這張表是我在實際 FWI 中總結出來的排查順序。第二行要特別強調Hv 必須和梯度使用同一套正演與伴隨實現任何網格差異、邊界條件差異或近似差分階數不一致都會讓 CG 殘差像噪聲一樣停在某個水平線上。這類問題不會在梯度檢驗中發現需要手動構造一個小擾動方向來自測。注意調參順序應當是先驗證梯度與 Hv、再調線搜索參數、最后動 CG 容差反過來會浪費大量計算。5. 在自定義 FWI 代碼里復用 TRN 思想梯度與 Hv 的一致性命門5.1 梯度有限差分模板在把 TRN 接入自己的 FWI 流程前先做一次梯度驗證。中心差分是常規做法def check_gradient(m, misfit, grad, eps1e-6): mplus m.copy(); mplus[0] eps mminus m.copy(); mminus[0] - eps num_g (misfit(mplus) - misfit(mminus)) / (2 * eps) rel_err abs(num_g - grad[0]) / (abs(num_g) 1e-30) print(frel_err{rel_err:.3e})這個函數只檢驗第一個參數實際應當隨機抽一批分量并保證擾動不超過正演離散誤差允許的范圍。eps選太大會引入截斷誤差選太小則浮點對消開始支配結果。我一般先在 1e-7、1e-6、1e-5 三檔測試選取相對誤差穩定且不隨eps劇烈變化的一檔作為正式測試配置。5.2 Hv 的 Taylor 展開檢驗梯度驗證通過后繼續驗證 Hessian 向量積。思路是對擾動向量 v定義 φ(t)g(mt·v)那么 φ(0)Hv可以用差分近似側的數值導數來對比你的hess_vec。示例代碼如下import numpy as np t 1e-6 g_plus compute_gradient(m t * v) g_minus compute_gradient(m - t * v) num_Hv (g_plus - g_minus) / (2 * t) ana_Hv hess_vec(m, v) rel np.linalg.norm(num_Hv - ana_Hv) / np.linalg.norm(ana_Hv) print(rel)這里的v必須取隨機方向不能取網格基向量隨機方向的覆蓋更廣泛能抓出類似“只有某些分量有符號錯誤”的 bug。另一個細節是如果正演代碼中有邊界吸收條件擾動向量如果碰到邊界區域數值導數會出現來源不明的異常那不代表 Hv 錯而是邊界類實現沒有對擾動向量做相同處理。測試時把擾動限制在計算域內部縮小到遠處。5.3 量綱與歸一化正式接入時最值得注意的一件事當梯度與 Hv 檢驗都通過后TRN 才能真正被信任。實際接入時最容易引發問題的不是優化器內部而是模型物理量綱。以速度模型為例單位取 km/s 時速度值在 1.5 到 6.0 之間數據殘差如果是地震振幅數量級可能在 1e-3 到 1e2 之間。兩者相乘產生的梯度和 Hessian 條件數會讓內層 CG 的殘差在真正收斂之前就觸碰浮點極限。建議在主循環開始前把模型向量的最大值歸一化到 1反演迭代完成后把增量乘回原始量綱這一步往往比任何 cg_tol 和 c2 都更能影響收斂速度。本文還有配套的精品資源點擊獲取