簡介pyGIMLi 例子是一套開源地球物理建模與反演軟件的教程示例合集面向地球物理、地質工程及相關領域的初學者與進階使用者。壓縮包共11個文件全部為Python腳本大小僅19KB便于快速下載與本地運行。內容按 _tutorials_auto_python 目錄組織涵蓋基礎操作、數(shù)值建模、反演與可視化等主題。基礎部分演示 gimli 核心對象使用、網(wǎng)格插值與質量檢查建模部分包含有限元法正演如直流電阻率、二維/一維熱傳導反演部分提供直流電阻率一維平滑與塊狀反演示例。此外還包括指數(shù)擬合、多項式擬合等入門腳本幫助用戶理解從數(shù)據(jù)到模型的全流程。資源已有1069人學習下載適合希望快速上手 pyGIMLi、掌握地質建模與反演實戰(zhàn)技能的讀者。示例腳本均提供完整可運行代碼配合注釋便于逐行理解通過運行這些腳本可以直觀學習代碼寫法與參數(shù)設置為后續(xù)開展科研或工程應用提供可直接改造的模板。1. 拆解 pyGIMLi 例子從正演到反演的一套完整地球物理工作流很多做電法、電阻率層析成像的從業(yè)者手里缺的不是數(shù)據(jù)而是一套能直接跑通的代碼。pyGIMLi 這個庫恰好補上了這塊短板它把網(wǎng)格生成、正演計算、反演迭代、結果可視化和模型評估全串在一條流水線上。我最初接觸它是因為同事丟來一個高密度電法數(shù)據(jù)集用傳統(tǒng)腳本處理完還要手動轉格式、手動挑反演參數(shù)折騰一下午才出一個像樣的剖面。換成 pyGIMLi 之后從加載數(shù)據(jù)到拿到反演結果整個過程被壓縮到一個 Python 腳本里還能順手把靈敏度、誤差分布這些常規(guī)工具不好出的東西一起畫出來。適合誰用做地質勘察、水文探測、工程檢測的工程師以及需要批量處理電阻率數(shù)據(jù)的研究生和研究員。這份資源的核心價值在于它不止給了能跑的 examples而是把地球物理反演里最容易被忽略的邊界條件、網(wǎng)格密度、正則化參數(shù)這些坑直接在代碼注釋里標了出來。2. 為什么選 pyGIMLi反演框架不是黑匣子是半透明的工具箱2.1 模塊化架構正演、反演、網(wǎng)格各自獨立又互相咬合pyGIMLi 的設計邏輯和很多商業(yè)軟件不一樣。商業(yè)軟件往往把正演和反演焊死在一個圖形界面后面你點幾個按鈕就出圖但中間發(fā)生了什么完全看不到。pyGIMLi 是 Python 庫底層用 C 實現(xiàn)核心計算上層暴露統(tǒng)一 API正演、反演、網(wǎng)格、可視化各是各的模塊串起來用的時候又能互相傳遞數(shù)據(jù)。從我的使用經驗來看這種設計最實際的好處是排查方便。正算出問題就在正演模塊里查反演不收斂就在反演參數(shù)里調網(wǎng)格質量差就直接看網(wǎng)格對象。舉一個具體場景做二維電阻率剖面時如果地表地形起伏很大商業(yè)軟件處理起來往往要額外買模塊pyGIMLi 里直接在網(wǎng)格生成階段把地形高程數(shù)組傳給createGrid或者通過Mesh類的節(jié)點修改實現(xiàn)。代碼層面大概長這樣import pygimli as pg import pygimli.meshtools as mt # 定義測線長度和地形高程數(shù)組 x pg.utils.grange(start0.0, stop100.0, n101) topo pg.Vector(101, 0.0) # 先用平地表占位 for i in range(len(x)): topo[i] 5.0 * pg.math.sin(x[i] / 20.0) 2.0 # 構造一個起伏地形 # 創(chuàng)建包含地形的網(wǎng)格 geom mt.createWorld( x, topo, worldBoundary[-20, -10, 120, 20], # 四個值左、下、右、上邊界 boundaryMarkers[1, 2, 3, 4] # 給四周邊界打標記 )這段代碼的邏輯是先用pg.utils.grange生成從 0 到 100 米的等間距測線坐標然后基于公式構造一個起伏地形數(shù)組最后createWorld用這些地形點生成網(wǎng)格幾何模型。worldBoundary四個參數(shù)控制網(wǎng)格外邊界范圍boundaryMarkers給邊界編號后續(xù)加邊界條件時直接引用標記號就行。參數(shù)層面有個細節(jié)topo數(shù)組的長度必須和x一致否則createWorld會直接拋異常。如果你的實際測線坐標不是從 0 開始記得先做平移不然地形和測線錯位反演結果會非常難看。2.2 反演引擎Gauss-Newton 迭代與正則化參數(shù)怎么配合pyGIMLi 反演核心用的是 Gauss-Newton 方法配合不同類型的一階或二階平滑約束。這套東西聽起來很數(shù)學實際用起來需要關心的其實就是兩個量正則化參數(shù)lam和最大迭代次數(shù)maxIter。lam控制模型平滑度和數(shù)據(jù)擬合度之間的權衡調大了結果更平滑調小了更貼數(shù)據(jù)但容易出現(xiàn)虛假異常體。我先給一個最常規(guī)的反演調用框架# 以電阻率反演為例 from pygimli.physics import ert # data 是之前加載并處理好的 ERT 測量數(shù)據(jù) mgr ert.ERTManager(data) model mgr.invert( lam20.0, maxIter10, verboseTrue, robustDataTrue, zWeight0.5 )lam20.0是很多標準例子的起始值如果反演出來的模型太碎、異常體邊緣鋸齒嚴重就調大lam如果數(shù)據(jù)擬合差、殘差始終降不下來就適當減小。robustDataTrue表示用 L1 范數(shù)處理數(shù)據(jù)中的離群點實測數(shù)據(jù)有異常跳變時強烈建議開。zWeight0.5是垂直方向平滑權重當?shù)刭|體傾向水平層狀時經常設成 0.5 以下。還有個被很多人忽略的點ERTManager會自動做數(shù)據(jù)誤差估計。默認按 3% 相對誤差算如果實際儀器精度達不到反演結果的信噪比會虛高。我會在數(shù)據(jù)加載后手動重置誤差見后面避坑章節(jié)。2.3 網(wǎng)格質量對反演結果的決定性影響網(wǎng)格是反演的容器。容器不對再好的算法也白搭。pyGIMLi 里網(wǎng)格生成用的是無結構三角網(wǎng)格相比規(guī)則矩形網(wǎng)格它能把地形和異常體邊界貼合得更自然。但無結構網(wǎng)格也意味著你要對網(wǎng)格密度分布有概念電極附近的網(wǎng)格必須加密遠離測線的區(qū)域可以放粗。常見做法是先用createWorld生成一個大范圍背景網(wǎng)格再在電極位置用createPolygon或者createCircle嵌入局部加密區(qū)。下面這段代碼演示了如何把加密區(qū)域和背景網(wǎng)格合并# 創(chuàng)建一個圓形的局部加密區(qū)域 circle mt.createCircle( pos[30.0, -5.0], # 圓心坐標 radius3.0, # 半徑 marker5, # 區(qū)域標記 area0.2 # 該區(qū)域最大單元格面積 ) # 合并背景網(wǎng)格和加密區(qū)域 merged mt.merge(geom, circle) # 生成最終的計算網(wǎng)格 mesh mt.createMesh(merged, quality34.0, area2.0)area參數(shù)控制加密區(qū)的網(wǎng)格單元面積上限值越小網(wǎng)格越密。這里circle的area0.2意味著加密區(qū)內每個三角單元面積不超過 0.2 平方米而背景網(wǎng)格area2.0就疏松得多。quality34.0是最小三角形質量角低于這個角度的單元會被重劃分通常 30-34 之間比較安全。一個經驗值電極間距是 1 米時電極附近網(wǎng)格單元面積設到 0.1-0.3 比較合理間距 5 米時加密區(qū)面積可以放到 1.0-2.0。網(wǎng)格過密會讓計算時間成倍增加反演一次從幾十秒變成十幾分鐘但精度提升很有限這個平衡點需要自己根據(jù)數(shù)據(jù)量試。3. 手把手跑通一個二維電阻率成像例子完整流程與參數(shù)說明3.1 數(shù)據(jù)準備從文本文件到 pyGIMLi 數(shù)據(jù)對象拿到任何實測數(shù)據(jù)第一步永遠是統(tǒng)一格式。pyGIMLi 的 ERT 模塊支持幾種常見格式最常用的是 Protocol 格式每行一條測量值包含電流電極 A、B電位電極 M、N視電阻率值或者電流電壓對。真實數(shù)據(jù)往往沒這么干凈可能電極編號從 1 開始坐標原點不在測線起點甚至電極排列不是等間距的。我先展示怎么手動構造一個數(shù)據(jù)對象import pygimli as pg from pygimli.physics import ert # 手工構造電極編號從1開始4個電極組成的測量序列 elecs [(0.0, 0.0), (5.0, 0.0), (10.0, 0.0), (15.0, 0.0)] # 每個測量四元組: (A, B, M, N) 對應電極數(shù)組下標從0開始 measurements [ (0, 1, 2, 3), (0, 2, 1, 3), (1, 2, 0, 3) ] # 對應的視電阻率值單位: 歐姆米 rho_values [120.5, 135.2, 118.7] # 轉換為 pyGIMLi 數(shù)據(jù)結構 data ert.createData(elecs, measurements, rho_values)這段構造代碼的關鍵點在于measurements里的四元組索引是從 0 開始而實際儀器記錄電極號通常從 1 開始轉換時記得減一這是個一不留神就翻車的細節(jié)。createData內部會把電極坐標和測量序列綁定成DataContainer對象后面所有正演反演都復用這個對象。如果你的數(shù)據(jù)量很大比如高密度電法幾百條測線、上千個電極手工構造不現(xiàn)實。這時直接用ert.ERTManager(datafilename)加載文件更方便它會自動識別常見格式并完成解析。加載這塊我在避坑章節(jié)會再補充一個關于單位的問題。3.2 正演計算驗證網(wǎng)格與測量方案是否匹配反演之前先做正演這個習慣幫我擋掉很多數(shù)據(jù)問題。正演的核心作用是給一個已知電阻率模型計算理論上的測量響應然后和實測數(shù)據(jù)對比。如果差得離譜說明網(wǎng)格或電極方案有問題早發(fā)現(xiàn)早處理。正演調用非常直接# 創(chuàng)建一個均勻半空間模型電阻率固定為100歐姆米 rhomap [[1.0, 100.0]] # 網(wǎng)格標記1對應電阻率100 forward_scheme mgr.createForwardScheme(mesh) # 執(zhí)行正演返回模擬數(shù)據(jù)對象 sim_data mgr.forward(mesh, rhomap, schemeforward_scheme) # 對比模擬數(shù)據(jù)和實測數(shù)據(jù)的視電阻率范圍 print(實測數(shù)據(jù)范圍:, min(rho_values), max(rho_values)) print(模擬數(shù)據(jù)范圍:, min(sim_data[rhoa]), max(sim_data[rhoa]))這里rhomap的寫法是[網(wǎng)格標記, 電阻率值]的列表網(wǎng)格標記在createWorld和createCircle階段定義。forward方法返回的sim_data[rhoa]是模擬視電阻率數(shù)組。如果模擬值和實測值差一個數(shù)量級以上先檢查單位pyGIMLi 默認長度單位是米電阻率單位是歐姆米如果你的電極坐標或電阻率值用的是厘米或毫西門子每米結果必然對不上。還有一個常見坑電極位置必須落在網(wǎng)格節(jié)點上。用createWorld生成背景網(wǎng)格后如果電極坐標超出了網(wǎng)格覆蓋的 x 范圍正演計算會跳過這些電極導致模擬數(shù)據(jù)不完整。檢查方法很簡單打印mesh.nodes()的最大最小值跟電極坐標對照一下。3.3 反演執(zhí)行參數(shù)矩陣與收斂判據(jù)正演驗證通過后進入反演階段。反演的本質是反復正演、計算殘差、更新模型直到殘差不再下降。pyGIMLi 的ERTManager.invert把這一整套迭代過程封裝了起來但理解內部邏輯對調參非常關鍵。一個完整的反演執(zhí)行片段mgr ert.ERTManager(data) model mgr.invert( lam20.0, maxIter15, startModelpg.Vector(1, 100.0), # 初始模型均勻100歐姆米 lamTolerance0.1, # 正則化參數(shù)下降閾值 fitTolerance0.01, # 數(shù)據(jù)擬合相對變化閾值 verboseTrue ) # 查看反演結果摘要 print(mgr.inv.invModel) print(均方根殘差:, mgr.inv.absRMS)參數(shù)含義拆開看startModelpg.Vector(1, 100.0)構造一個只有 1 個元素的向量初始值 100 歐姆米pyGIMLi 會把它自動擴展到所有網(wǎng)格單元。lamTolerance0.1表示當正則化參數(shù)的變化率小于 10% 時認為達到收斂fitTolerance0.01表示數(shù)據(jù)擬合度的相對改善小于 1% 時停止迭代。這兩個容差是常見的收斂判據(jù)如果反演一直跑不完可以把fitTolerance放寬到 0.03。注意startModel的選值最好接近測區(qū)背景電阻率。如果測區(qū)是含水地層背景可能只有幾十歐姆米你卻給了 1000反演的迭代次數(shù)會明顯增加甚至陷入局部極小值。沒有先驗信息時用所有視電阻率的對數(shù)平均值作為初始模型是最穩(wěn)妥的。4. 避坑指南反演跑不通、結果跳變、網(wǎng)格錯位這幾個真實踩坑記錄4.1 電極坐標單位不匹配正演數(shù)據(jù)全是 NaN現(xiàn)象加載實測數(shù)據(jù)后正演計算出的模擬數(shù)據(jù)全是nan或者反演進行到一半直接報錯退出。原因這是我一個同事踩過的真坑。他手里的電極坐標以公里為單位而 pyGIMLi 默認按米處理測線長 1.2 公里網(wǎng)格范圍卻只有幾十米電極全部跑到網(wǎng)格外面去了。正演在電極位置找不到對應的網(wǎng)格單元計算自然失敗。解決加載數(shù)據(jù)后先檢查坐標范圍。如果坐標數(shù)值很大或很小主動縮放單位。data.setSensorPositions([(p[0] * 1000.0, p[1] * 1000.0) for p in data.sensorPositions()])把坐標統(tǒng)一到米以后重新打印范圍確認一下再進入正演環(huán)節(jié)。4.2 反演結果出現(xiàn)條紋狀虛假異?，F(xiàn)象反演完成的剖面圖上地表淺部出現(xiàn)交替排列的高阻-低阻條紋看起來像斑馬線。這些條紋在實測數(shù)據(jù)里沒有明確對應。原因網(wǎng)格加密區(qū)域和背景網(wǎng)格疏密差異過大導致反演時淺部單元權重失衡。更常見的原因是lam值設得太小模型過度擬合了數(shù)據(jù)中的噪聲。解決遇到條紋先別急著改網(wǎng)格優(yōu)先把lam調大試試。我會按 20 → 50 → 100 的順序遞進測試同時觀察均方根殘差變化。如果殘差變化不大而條紋明顯消失說明是正則化約束不夠。如果殘差顯著上升才考慮網(wǎng)格問題。4.3 反演不收斂殘差停在 10% 上下反復橫跳現(xiàn)象迭代到第 5 到第 8 次時均方根殘差不再單調下降而是在 10% 附近小幅震蕩verbose 輸出里能看到數(shù)值反復跳。原因數(shù)據(jù)中存在離群點單個測量值的誤差遠大于其他值Gauss-Newton 迭代被這個異常值持續(xù)拉扯。另一個可能原因是初始模型遠離真實模型迭代在局部極小值附近抖動。解決打開robustDataTrue這個參數(shù)會切換為 L1 范數(shù)約束顯著降低離群點的影響。同時把maxIter調大到 20 以上讓迭代有足夠空間跳出抖動區(qū)。4.4 網(wǎng)格邊界條件錯誤導致深部模型失真現(xiàn)象反演結果的深部區(qū)域出現(xiàn)明顯的高阻異常即使測線下方實際是低阻含水層也會被反演成高阻。原因createWorld的worldBoundary設置得太窄底部邊界離測線深度太近邊界條件對深部區(qū)域產生了強烈約束。解決把底部邊界放到目標探測深度的 3 倍以上。例如測線長度 100 米、最大探測深度按 30 米算worldBoundary的底部至少放到 -100 米。geom mt.createWorld( x, topo, worldBoundary[-20, -100, 120, 20], boundaryMarkers[1, 2, 3, 4] )這個修改會略微增加網(wǎng)格單元總數(shù)但深部可靠性顯著提升。5. 進階技巧用 pyGIMLi 做模型質量評估與可視化輸出5.1 靈敏度圖識別反演結果哪些區(qū)域可信哪些是陪跑反演結果不能拿到手就畫圖交差。我建議每次反演完都額外輸出一張靈敏度圖。靈敏度表示模型參數(shù)變化對觀測數(shù)據(jù)的響應程度靈敏度高的區(qū)域反演結果可信度高靈敏度低的區(qū)域基本屬于數(shù)值填充沒有實際地質意義。pyGIMLi 里計算靈敏度并不復雜sens mgr.inv.sensitivity() # 歸一化到 0-1 范圍 sens_norm (sens - sens.min()) / (sens.max() - sens.min()) # 和反演結果一起可視化 ax1 pg.show(mesh, model, coveragesens_norm, labelResistivity (Ohm m))coverage參數(shù)控制模型顯示時被靈敏度閾值遮罩的程度。coverage0.1表示靈敏度低于 10% 的區(qū)域將以半透明方式顯示這樣剖面上不可信的區(qū)域一目了然。實際交付報告時我通常把coverage設到 0.05 到 0.15 之間太低遮不住噪聲區(qū)太高會把淺部有效信息也糊掉。5.2 多測線批量反演腳本化處理的一個閉環(huán)范例當你有幾十條平行測線要處理時逐條手工配置和等待就是不現(xiàn)實的。pyGIMLi 的腳本化特性在這里能體現(xiàn)最大價值。我一般把一個批處理流程組織成三步加載配置、循環(huán)反演、匯總輸出。import os import glob # 假設測線數(shù)據(jù)文件按 std_01.dat, std_02.dat 排列 dat_files sorted(glob.glob(std_*.dat)) for idx, f in enumerate(dat_files): mgr ert.ERTManager(f) model mgr.invert( lam20.0, maxIter12, robustDataTrue, verboseFalse ) # 保存反演結果和靈敏度 mgr.saveResults(f.replace(.dat, _inverted.png)) mgr.inv.save(f.replace(.dat, _invModel.npy)) print(f第{idx1}/{len(dat_files)}條測線完成殘差: {mgr.inv.absRMS:.2f}%)saveResults會生成帶色標和地形顯示的反演剖面圖save把模型向量直接存成 numpy 格式方便后續(xù)批量讀取做統(tǒng)計對比。一個建議批處理時verboseFalse關掉逐次迭代輸出只在每個文件結尾打印殘差摘要日志會干凈很多。absRMS的單位是百分比表示數(shù)據(jù)擬合的均方根相對誤差這個值也是你判斷每條測線數(shù)據(jù)質量最直接的指標。5.3 與地質解釋銜接剖面圖疊加鉆孔信息反演剖面最終要能跟地質分層對上話。我經常在 pyGIMLi 生成的圖上疊加鉆孔巖性分界點用來校驗反演邊界和實際地質界線是否一致。import matplotlib.pyplot as plt # 假設一個鉆孔在 x45 米處巖性分界在 -12 米 ax pg.show(mesh, model, labelResistivity (Ohm m)) ax[0].plot([45.0, 45.0], [0.0, -25.0], k--, linewidth1.5, labelborehole) ax[0].plot([45.0], [-12.0], rv, markersize10, labelinterface) ax[0].legend() plt.show()坐標對應關系要注意pyGIMLi 的 y 軸向下為正所以深度 -12 米在垂直方向上正好顯示在下方。鉆孔位置如果落在測線范圍之外plot還是會把線畫出來但已經超出網(wǎng)格覆蓋區(qū)需要在圖上留意別誤解。疊加圖的意義在于快速驗證反演結果的可靠性為最終地質解釋提供一層獨立證據(jù)。從那以后我每次跑完測線數(shù)據(jù)都會強制走一遍「正演驗證 → 反演 → 靈敏度評估 → 鉆孔對比」這個流程少選條反演直接交圖的情況就很少發(fā)生了。希望幫到你。本文還有配套的精品資源點擊獲取