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