李學(xué)蘭董艷輝李帝銓③

（①中南大學(xué)，地球科學(xué)與信息物理學(xué)院 長沙 410083）

（②中國科學(xué)院地質(zhì)與地球物理研究所 北京 100029）

（③中南大學(xué)，有色金屬成礦預(yù)測教育部重點實驗室 長沙 410083）

基于地下水?dāng)?shù)值模擬的地球物理反演解釋優(yōu)化方法研究*

李學(xué)蘭①董艷輝②李帝銓①③

（①中南大學(xué)，地球科學(xué)與信息物理學(xué)院 長沙 410083）

（②中國科學(xué)院地質(zhì)與地球物理研究所 北京 100029）

（③中南大學(xué)，有色金屬成礦預(yù)測教育部重點實驗室 長沙 410083）

地球物理的多解性無處不在，面對這種情況我們往往結(jié)合地質(zhì)資料來對反演結(jié)果進(jìn)行合理的地質(zhì)解釋。然而在偏遠(yuǎn)地區(qū)地質(zhì)、深部鉆探資料非常缺乏，野外已有的井孔也往往難以獲取巖性編錄資料。當(dāng)這些資料比較少時，反演結(jié)果依然存在難以解釋或者不確定等問題。是否可以借助于地下水位與水質(zhì)信息包含的地質(zhì)體水力特性，來輔助地球物理反演過程？基于此思路探索了基于地下水模型的地球物理反演解釋優(yōu)化方法：以大地電磁法為實現(xiàn)手段，通過研究已知理想地質(zhì)模型和多個基于反演結(jié)果解釋得到的地質(zhì)模型，建立地下水?dāng)?shù)值模型，對比地下水水流和溶質(zhì)運移數(shù)值模擬結(jié)果，最終確定相對最佳地質(zhì)模型。并通過敏感性分析驗證了運用地下水?dāng)?shù)值模擬搜索相對最佳地質(zhì)模型的可行性。

地球物理反演 地質(zhì)模型 地下水流動數(shù)值模擬 溶質(zhì)運移

0 引 言

地球物理資料的反演和解釋（王家映，2002）是用來了解地球內(nèi)部物理性質(zhì)以及礦產(chǎn)資源分布的。在科學(xué)研究和工程技術(shù)領(lǐng)域中，凡是涉及根據(jù)觀測數(shù)據(jù)推求模型或者體系時，都無法回避反演問題（魏超等，2008）。目前，自然地質(zhì)活動與人文地質(zhì)活動不斷增加，使得目標(biāo)體所處的地質(zhì)環(huán)境越來越復(fù)雜，使用傳統(tǒng)單一的地球物理探測方法并不能很好地獲得高質(zhì)量的野外采集數(shù)據(jù)，加上各類地球物理方法都具有一定的局限性（陳思宇，2014），如：電法（電磁法）勘探分辨率的高低往往受探測深度的影響，有時獲得地質(zhì)體異常響應(yīng)可能會表現(xiàn)出信號微弱、輪廓不清晰等特點（戴前偉等，2013），屏蔽層的存在使得地震勘探也無能為力，空間分布特征的探測對地球物理測井也是一種挑戰(zhàn)，非唯一性對重力及磁法勘探的影響更不容小覷。單一地球物理勘探方法的局限性及場的等效性是造成反演問題多解性的一部分，它是無法避免的。除此之外，地球物理問題反演多解性與野外數(shù)據(jù)采集過程中噪聲強弱、儀器精度高低、數(shù)據(jù)量多少等因素也息息相關(guān)，這些因素均影響了地球物理勘探成果的可靠度。

針對反演問題的多解性，解決方案通常有兩種。其中之一是通過已有的可靠的地質(zhì)信息來約束地球物理反演過程，引導(dǎo)反演過程趨于合理、趨向真實的野外地質(zhì)情況（郭冬等，2014）。它是指對于先驗地質(zhì)信息確定的單元，迭代時保持其物性資料的不變，即在目標(biāo)函數(shù)中引入含有先驗地質(zhì)信息項，得到考慮數(shù)學(xué)和地質(zhì)合理性的反演結(jié)果（Yeh et al.，2002；祁光，2013；陳思宇，2014）。它能降低反演非唯一性，提高分辨率，但是先驗地質(zhì)信息往往比較匱乏，在這種情況下盡可能擴大數(shù)據(jù)量，從數(shù)學(xué)角度出發(fā)，讓反演方程不斷趨向完善與合理是行之有效的一種方法（郭冬等，2014）。因此眾多學(xué)者在某些經(jīng)典算法上進(jìn)行改進(jìn)（魏超等，2008；Yuan et al.，2009），如魏超等（2008）提出了基于量子蒙特卡羅的地球物理反演方法，這種方法適合于非線性、多極值的地球物理反演問題，在收斂速度和避免陷入局部極小等方面有著一定的優(yōu)勢，具有較強的通用型，但是基于量子蒙特卡羅的反演方法中試探波函數(shù)以及動能項的設(shè)置是根據(jù)量子系統(tǒng)來設(shè)置的，對于實際地球物理反演問題而言，這種設(shè)置可能并不那么合理。與此同時，地球物理聯(lián)合反演（敬榮中等，2003；于鵬等，2006；何委徽等，2009；Wu et al.，2012）被廣泛應(yīng)用于反演中，其實質(zhì)為數(shù)據(jù)融合。聯(lián)合反演發(fā)展主要分為以下兩個方向，其一為已知地質(zhì)信息和地球物理模型統(tǒng)一起來的聯(lián)合反演，或者是不同地球物理勘探方法所得觀測數(shù)據(jù)之間的聯(lián)合，如地震和電法，地震和重力，重力和大地電磁法等，但是不同物性之間存在著相關(guān)的內(nèi)在聯(lián)系是聯(lián)合的基礎(chǔ)（林達(dá)明等，2014），通過這種方法能達(dá)到相互約束、補充的目的，從而減小反演多解性。其二為數(shù)學(xué)、地質(zhì)、地球物理模型統(tǒng)一起來的聯(lián)合反演，它利用波動場和擴散場間的相互聯(lián)系，通過數(shù)學(xué)變換，建立多種地球物理勘探方法的數(shù)學(xué)物理模型，并將其統(tǒng)一起來，得到共同數(shù)學(xué)物理模型，從而進(jìn)行統(tǒng)一的反演成像（陳潔等，2007）。聯(lián)合反演使得反演結(jié)果有更高的分辨率，并能較好地解決非唯一性問題，但其要緊密結(jié)合巖石物理性質(zhì)的統(tǒng)計規(guī)律，才能取得兼容地質(zhì)與綜合方法的應(yīng)用效果（楊文采，2002）。盡管不同的地球物理反演方法均有一定的局限性，但是隨著反演約束條件的增加，反演所得的結(jié)果與真實情況的吻合度逐漸增加。源于電流場與滲流場的相似性（喻永祥等，2006），地球物理與地下水緊密結(jié)合在一起。目前已有學(xué)者嘗試用各種地球物理方法來間接求取滲透系數(shù)，如喻永祥等（2006）實現(xiàn)了運用高密度電阻率成像法獲得非均質(zhì)多孔介質(zhì)的滲透系數(shù)。Andersen et al.（2013）將地下水與地球物理模型結(jié)合起來研究山谷流域的捕獲區(qū)、補給區(qū)、地下水年齡，盡管文中指出不能太依賴地球物理反演結(jié)果，但是它對本文的研究有很好的啟發(fā)性。Singha et al.（2006）提出了一種建立電阻率與濃度物性關(guān)系的方法，通過物性關(guān)系可以將電磁法中所得反演剖面轉(zhuǎn)換成濃度剖面，很好地從反演結(jié)果中提取了地下水的相關(guān)信息。Yeh et al（2002）提出了將先驗地質(zhì)信息引入到反演中，并由電阻率得出美國新墨西哥州索科羅STVZ區(qū)域含水分布曲線。

當(dāng)?shù)刭|(zhì)資料特別是鉆探資料豐富時，地球物理反演結(jié)果可以得到較好的優(yōu)化，可以得到地層結(jié)構(gòu)較準(zhǔn)確的刻畫。但在偏遠(yuǎn)的、缺少地質(zhì)資料地區(qū)，深部鉆探資料非常缺乏，野外已有的井孔也往往難以獲取巖性編錄資料。由于地質(zhì)體的水力特性，如滲透系數(shù)，與地層巖性有著直接關(guān)聯(lián)，那么，是否可以通過地質(zhì)體中的地下水系統(tǒng)循環(huán)特征，如水位、水質(zhì)信息來反應(yīng)地層巖性并輔助地球物理反演解釋呢？基于此思路，本文以大地電磁法為實現(xiàn)手段，探索地質(zhì)數(shù)據(jù)缺乏時，基于地下水模型優(yōu)化地球物理反演解釋的方法，通過已知理想地質(zhì)模型及敏感性分析來驗證運用地下水?dāng)?shù)值模擬搜索相對最佳地質(zhì)模型的可行性。

1 研究方法

1.1 地球物理正反演

地球物理正演是給定場源的分布，求解場值的大小。反演是正演的逆過程，即由場源觀測數(shù)據(jù)求出場源的物性分布（徐世浙，1994）。如果把模型空間中的一個點定義為m，把數(shù)據(jù)空間中的一個點定義為d，按照物理定律，可以把相關(guān)關(guān)系寫成：

其中，G為模型空間M到數(shù)據(jù)空間D的一個映射。從式（1）、式（2）可知正演是從模型空間到數(shù)據(jù)空間的過程，反演是從數(shù)據(jù)空間到模型空間的過程。用流程圖描述如下（圖1）：

圖1 正反演流程圖Fig.1 The flowchart of forward and inversion

1.2 地下水?dāng)?shù)學(xué)模型與數(shù)值求解

1.2.1 地下水流動數(shù)學(xué)模型及數(shù)值方法

依據(jù)三維地下水流動的水文地質(zhì)概念模型，地下水三維非穩(wěn)定流方程可以描述為（Panday et al.，2013）

式中，K為滲透系數(shù)；h為水頭高度；Ss為儲水系數(shù)；t為時間；W為單位體積的源匯項，上述數(shù)學(xué)模型采用的數(shù)值解法為有限差分法（陳崇希等，2009）。

1.2.2 溶質(zhì)運移數(shù)學(xué)模型及數(shù)值方法

飽和帶溶質(zhì)運移的對流－彌散方程為（薛禹群等，2007）

式中，c為溶液中某種組分的濃度；→u為實際平均流速矢量；→D為水動力彌散系數(shù)。大多數(shù)解決溶質(zhì)運移的數(shù)值方法為歐拉法、拉格朗日法、歐拉—拉格朗日混合法，詳細(xì)求解可參考（Zheng et al.，1999）。

1.3 基于地下水?dāng)?shù)值模擬的地球物理反演解釋優(yōu)

化方法

為了了解地球內(nèi)部物理性質(zhì)及其礦產(chǎn)資源分布（魏超等，2008），往往會通過地球物理資料的反演來對其進(jìn)行合理的解釋，電法、電磁法、地震、測井等地球物理方法為常用手段。基于地下水模型的地球反演地質(zhì)解釋優(yōu)化方法大致實現(xiàn)過程為（圖2）：

（1）采用合理的地球物理反演方法對野外實測數(shù)據(jù)進(jìn)行反演成像；

（2）結(jié)合已知地質(zhì)、鉆探資料，對反演成像的結(jié)果做出合理的地質(zhì)解釋，由于反演的多解性，往往能得到多個解釋結(jié)果；

（3）為了得到更合乎實際情況的地質(zhì)解釋，將（2）中所得的多個地質(zhì)解釋結(jié)果（包含巖性信息）轉(zhuǎn)換成地下水模型（包含水力特性參數(shù)信息）；

（4）通過水文地質(zhì)調(diào)查，確定邊界定解條件、源匯項等相關(guān)水文信息，對（3）中地下水模型進(jìn)行地下水?dāng)?shù)值模擬，得到不同地質(zhì)解釋結(jié)果對應(yīng)地下水模型的水位、濃度信息；

（5）將（4）中模擬的水位、濃度信息與水文地質(zhì)中水位、水質(zhì)監(jiān)測結(jié)果進(jìn)行對比，最終確定更符合野外實際情況的地質(zhì)解釋。

本文從已知理想地質(zhì)體出發(fā)，以大地電磁法為實現(xiàn)手段，通過數(shù)值模擬實現(xiàn)上述優(yōu)化方法。已知水位、濃度的獲取為已知理想地質(zhì)體地下水?dāng)?shù)值模擬的結(jié)果。首先，采用MTsoft2D軟件完成地球物理正演，得到反演數(shù)據(jù)輸入并進(jìn)行地球物理反演，反演過程中主要采用不同模式的選擇、道間距以及插值頻點的改變得到反演結(jié)果，進(jìn)一步將其進(jìn)行地質(zhì)解釋得到基于反演解釋的地質(zhì)模型。其次，將多個基于反演解釋的地質(zhì)模型轉(zhuǎn)換成地下水模型，通過地下水?dāng)?shù)值模擬得到地下水模型的模擬水位、濃度分布。最后，統(tǒng)計多個基于反演解釋的地質(zhì)模型對應(yīng)的模擬水位、濃度與已知理想地質(zhì)體數(shù)值模擬的水位、濃度的標(biāo)準(zhǔn)偏差、均方根誤差、歸一化誤差、平均歸一化均方根誤差、相關(guān)系數(shù)的大小，并比較各項誤差大小，最終確定相對最佳地質(zhì)模型，達(dá)到優(yōu)化地球物理反演解釋的目的。

圖2 基于地下水模型的地球物理反演解釋優(yōu)化方法實現(xiàn)圖Fig.2 The chart of an optimization method for interpreting ERT data based on groundwater numericalmodeling

2 已知理想地質(zhì)模型的建立

針對1.3節(jié)所述的基于地下水模型的地球反演解釋優(yōu)化方法，采用MTsoft2D建立給定地質(zhì)模型Model i（圖3）。模型長2000m，深300m，表層為25m厚的黏土覆蓋層，中間為充填中砂的U型谷，U型谷中間為一條厚度為50m的黏土層，基巖為花崗巖。參照幾種巖石電阻率值的分布范圍曲線（李金銘，2005），取黏土電阻率為10Ω·m，砂巖為500Ω·m，花崗巖為10000Ω·m，得到構(gòu)造地電模型。

3 大地電磁的正反演計算及反演解釋優(yōu)化方法

圖3 已知理想地質(zhì)模型Model iFig.3 The known geologicalmodel Model i

以上述已知理想地質(zhì)模型為例，通過大地電磁法對其進(jìn)行研究。電磁法常用的正演模擬軟件有MTsoft2D，反演軟件有MTsoft2D、SCS2D、GME_3DI（V4.1）等，本次MT正演模擬采用成都理工大學(xué)開發(fā)的MTsoft2D，數(shù)值解法為有限單元法。反演計算采用中南大學(xué)自主開發(fā)的重磁電三維反演成像解釋一體化系統(tǒng)GME_3DI（V4.1），反演成像方法為二維連續(xù)介質(zhì)廣義逆快速反演（戴世坤等，1997）。

3.1 大地電磁正演計算

通過設(shè)置不同的網(wǎng)格剖分方式，根據(jù)模擬敏感性分析結(jié)果，采用垂直、水平剖線0.05、0.01km／格，并添加自動邊界剖線最合適。水平方向共布置了41個測點，間距50m，頻率為1.1250～7680.589Hz（共40個頻點）。

圖4分別為TE、TM模式正演模擬結(jié)果。從TE模式結(jié)果看（圖4a），整個剖面的電阻率值在0～5000Ω·m之間，與模型所設(shè)計的電阻率相比有一定的差異，最高電阻率沒有達(dá)到預(yù)期值10000Ω·m。從整個剖面上來看它分為3層，整體電阻率為低－中－高的趨勢，并且能夠看到U型谷的存在，但是并不能刻畫U型谷中間低阻層的存在。而TM模式正演結(jié)果剖面（圖4b）與TE模式相比整體趨勢兩者結(jié)果相似，但是基底的電阻率遠(yuǎn)高于TE模式，U型谷刻畫得比較明顯。

3.2 大地電磁反演

不同模式測量、插值頻點數(shù)目以及道間距可以獲得不同的反演結(jié)果，本文中采用不同模式測量、插值的道間距以及頻點數(shù)描述反演解釋多解性問題，并非描述地球物理反演算法的優(yōu)劣。本次反演中采用了3種極化模式（TE、TM、TEM）以及3組插值后頻率和間距（40，25）、（40，50）、（50，25）（原本頻組為40頻點，間距50m），反演深度300m，得到9個二維反演結(jié)果（圖5）。

圖5的Model 1-3、4-6、7-9分別為TE、TM、TEM模式頻間組合為（40，25）、（40，50）、（50，25）的反演結(jié)果。從圖5可知TE模式下不同參數(shù)反演斷面圖Model 1-3均能大致刻畫U型谷中低阻帶的存在，但詳盡不一，3個反演結(jié)果與已知理想地質(zhì)模型相比，存在低阻帶與基巖不完全相連。Model 1刻畫較完整，Model 2在低阻帶的上部出現(xiàn)了帶狀高阻，Model 3的低阻帶主要分布在U型谷的右半部分。TM模式Model 4、Model 5 U型谷中低阻隱約能顯示，但Model 4頂部的低阻覆蓋層和基巖部分已經(jīng)不完整，尤其是Model 6、Model 9中間的低阻帶已經(jīng)缺失，失去原本模型的大致結(jié)構(gòu)。TEM模式中的低阻帶的Model 7、Model 8低阻也與基巖部分不相連，結(jié)果和Model 3相似。

3.3 大地電磁反演結(jié)果解釋－地質(zhì)模型

圖6是基于大地電磁反演結(jié)果解釋所得的地質(zhì)模型，主要依據(jù)電阻率（李金銘，2005）的分區(qū)劃分。由于黏土電阻率一般比較低，當(dāng)電阻率低于50Ω·m劃分為黏土。砂巖的電阻率變化范圍比較大，50～1000Ω·m劃分為砂巖。當(dāng)電阻率大于1000Ω·m劃分為花崗巖。

3.4 基于地下水?dāng)?shù)值模擬的優(yōu)化方法

圖4 TE模式（a）和TM模式（b）正演結(jié)果Fig.4 The results of TE（left）and TM（right）mode forward modeling

由上文可知，通過改變地球物理反演參數(shù)（插值的頻點數(shù)、道間距）或者采用TE、TM、TEM不同測量模式，可以得到多個可供解釋的大地電磁反演結(jié)果，但是具體哪個反演結(jié)果更符合已知理想地質(zhì)模型，野外實測數(shù)據(jù)的反演到底哪一個更吻合地質(zhì)解釋，這就需要通過地下水模型來進(jìn)行下一步的研究。

圖5 大地電磁反演結(jié)果Fig.5 The results of MT inversion

3.4.1 地下水?dāng)?shù)值模型建立

本次數(shù)值模型的建立采用數(shù)值模擬軟件Visual MODFLOW 2011.1，該軟件由于其優(yōu)良的可視化界面以及豐富的功能被廣泛用于地下水?dāng)?shù)值模擬研究中（魏亞強等，2015；鄭燦政等，2015），它集成了MODFLOW，MODPATH，MT3DMS、ZoneBudget、PEST模塊，模塊間相互獨立，是三維地下水流動和污染物運移模擬實際應(yīng)用中功能最為完整、操作比較簡單、圖形界面直觀的專業(yè)軟件。

地下水?dāng)?shù)值模型的建立主要包括以下幾個步驟：

（1）空間剖分：研究模型長2000m，寬300m，剖面面積為2000×300m2，采用有限差分的方法進(jìn)行矩形剖分，單元長10m，寬5m，共剖分了1行200列60層，單元總數(shù)12000個。

（2）邊界、定解條件處理：所有地下水模型底部處理為隔水邊界，左右邊界設(shè)置為距離頂板厚度為10m定水頭，水位分別為290m、315m，其他單元格的初始水位為290m。

（3）水文地質(zhì)參數(shù)處理：參照不同巖性滲透系數(shù)經(jīng)驗值（中國地質(zhì)調(diào)查局，2012），將不同巖性的區(qū)域賦相應(yīng)的值，黏土、中砂、花崗巖依次為2.9e－6m·d－1、1.2e－4m·d－1、5e－13m·d－1，其他地質(zhì)參數(shù)均采用系統(tǒng)默認(rèn)值。

（4）源匯項來源：定水頭邊界的補給；降雨入滲（100mm·y－1）；在X方向500m，Z方向230～250m處設(shè)置注水井（注水量150m3·d－1），X方向1500m，Z方向90～110m處設(shè)置抽水井（抽水量200m3·d－1），每口注水井加入2個點污染源，濃度分別為400mg·L－1、500mg·L－1。

圖6 基于大地電磁反演的地質(zhì)解釋Fig.6 The geological interpretation of MTmethod inversion

3.4.2 地下水流和溶質(zhì)運移模擬結(jié)果

地下水流、溶質(zhì)運移模擬分別使用的是Visual MODFLOW2011.1中MODFLOW、MT3DMS程序包。模擬非穩(wěn)定流水位、濃度變化情況，時間為2000d。

在整個剖面上設(shè)立22口監(jiān)測井監(jiān)測激勵源產(chǎn)生的響應(yīng)。將已知理想地質(zhì)模型對應(yīng)的地下水流監(jiān)測井的水位、濃度作為實測水位、濃度（Observation），地球物理解釋對應(yīng)的地下水模型的模擬水位、濃度作為計算水位、濃度（Calculation），統(tǒng)計觀測值與計算值之間的各項誤差，分別見表1、2。其中平均歸一化均方根誤差為歸一化均方根誤差與參與模擬計算模型歸一化均方根誤差平均值的比值，當(dāng)模型間的差值大于5％，模型間具有很好的區(qū)分度。

從表1可知Model 4、Model 6、Model 7、Model 9的各項誤差相比于其他模型均比較大，平均歸一化均方根誤差依次為169.7％、184.8％、163.1％、110.7％，遠(yuǎn)遠(yuǎn)高于平均值，根據(jù)相關(guān)系數(shù)的值可知這幾個模型與實測值的相關(guān)性均比較低。Model 1、Model 3的各項誤差都比較小，平均歸一化均方根誤差分別為43.7％、48.6％，但是平均歸一化均方根誤差差值小于5％，因此Model，1、Model，3并不具有很好的區(qū)分度。Model，2、Model，5的各項誤差相當(dāng)接近且Model，8的結(jié)果優(yōu)于Model，2、Model，5。根據(jù)表1初步確定與Model，i的擬合優(yōu)劣順序為Model，1（3）、8、2（5）、9、7、4、6。

表1 不同地下水模型模擬水位與實測水位誤差統(tǒng)計Table 1 The errors between observation and calculation water level ofmultiple geologicalmodels

從表2中可知Model 6、7、9的平均歸一化均方根誤差為131.2％，與已知理想地質(zhì)模型的相關(guān)系數(shù)值接近于0。Model 3、5、8的誤差的統(tǒng)計結(jié)果也相當(dāng)接近，平均歸一化均方根誤差分別為79.5％、79.5％、77.8％（誤差差值小于5％），在溶質(zhì)運移中并不具有很好的區(qū)分度，Model 4各項誤差遠(yuǎn)大于3、5、8。Model 1各項誤差均為最小，與實測濃度擬合的最佳，與其他模型的平均歸一化均方根誤差最小差值也為16.2％。根據(jù)表2確定與Model i的擬合優(yōu)劣順序為Model 1、3（5、8）、2、4、6（7、9）。

從不同模型模擬水位、濃度與實測結(jié)果相比可知，不管是水流模擬還是溶質(zhì)運移的結(jié)果，Model 6、7、9均被證明與已知理想地質(zhì)模型差異性比較大，相關(guān)性幾乎為0。Model 4在溶質(zhì)運移中誤差值也明顯高于Model 1、3、5、8，在模擬水位誤差中更容易體現(xiàn)。Model 1與2、3、5雖然在模擬水位上的各項誤差值極為接近，但是在溶質(zhì)運移中具有明顯的區(qū)分度，類似有Model 3、5、8在溶質(zhì)運移的結(jié)果相似，但是在模擬水位上存在明顯的差異。以水流模擬為基礎(chǔ)，當(dāng)模型間的平均歸一化均方根誤差差值低于5％時，進(jìn)一步利用溶質(zhì)運移結(jié)果，從而得出與已知理想地質(zhì)模型的擬合優(yōu)劣順序為：Model 1、3、8、5、4、9、7、6。

表2 不同地下水模型濃度觀測值與計算值誤差統(tǒng)計Table 2 The errors between observation and calculation concentration ofmultiple geologicalmodels

3.4.3 數(shù)值模擬中抽注水井靈敏度分析

由于3.4.2只分析一種抽注水情景，如何確定所得最佳擬合結(jié)果是偶然因素造成還是實際上就存在這么一種響應(yīng)？本節(jié)主要來探討這個問題。對于上述問題的研究，通過抽取上一小節(jié)中的4個地下水模型Model 1、3、5、8，敏感性分析通過更換抽注水井的位置來實現(xiàn)（其他條件保持一致），分別設(shè)置以下兩種情景：

（1）情景一：X方向650m，Z方向80～100m處設(shè)置注水井（注水量150m3·d－1），X方向1500m，Z方向230～250m處設(shè)置抽水井（抽水量200m3·d－1）。

（2）情景二：X方向500m，Z方向230～250m處設(shè)置注水井（注水量150m3·d－1），X方向1000m，Z方向90～110m處設(shè)置抽水井（抽水量300m3·d－1），X方向1500m，Z方向230～250m處設(shè)置注水井（注水量100m3·d－1）。

類似于3.4.2，統(tǒng)計了情景一、二下不同地下水模型模擬水位（濃度）與實測水位（濃度）值誤差（表3、表4）。

從表3中可知模型間平均歸一化均方根誤差差值均大于5％，因此Model 1、3、5、8具有很好的區(qū)分度，并不需要再運用地下水中的溶質(zhì)運移模擬即可確定與已知理想地質(zhì)模型擬合優(yōu)劣順序為Model 1、3、8、5。

從表4可以看出Model 1的各項誤差均為最小，Model 5為最大，Model 1與Model 8兩者平均歸一化均方根誤差差值8.5％，具有很好的區(qū)分度。而Model 3、8的平均歸一化均方根誤差差值2.3％，因此進(jìn)一步對其進(jìn)行溶質(zhì)運移模擬，在濃度各項誤差統(tǒng)計結(jié)果中模型間平均歸一化均方根誤差相差23.2％，所以通過溶質(zhì)運移能夠很好地將兩者區(qū)分出來。在情景二下得到擬合順序仍為Model 1、3、8、5。

表3 情景一：不同地下水模型模擬水位與實測水位值誤差統(tǒng)計Table 3 The scenario one：the errors between observation and calculation water level ofmultiple geologicalmodels

表4 情景二：不同地下水模型模擬水位與實測水位值誤差統(tǒng)計Table 4 The scenario two

4 模擬結(jié)果與分析

通過大地電磁反演得到多個反演結(jié)果（圖5），將其解釋為地質(zhì)模型（圖6）并轉(zhuǎn)換成地下水?dāng)?shù)值模型。當(dāng)運用水位進(jìn)行比較時，若存在兩個、甚至多個模型不具有良好區(qū)分度、可以通過溶質(zhì)運移進(jìn)行進(jìn)一步分析，最終確定與已知理想地質(zhì)模型的擬合順序，從而篩選出相對最優(yōu)結(jié)果。

從模擬結(jié)果可知Model 2、4、6、7在地下水模擬中與Model i相差比較大，通過水位能很好地將其區(qū)分；對比地質(zhì)模型Model 2、4、6、7與Model i的差異，從圖6也能看出這4個模型確實與Model i的相似性比較差。對于Model 1、3、5、8盡管從地質(zhì)模型的視覺效果上看，差異性可能并不是很大，但是通過地下水?dāng)?shù)值模擬可以很好地將其區(qū)分出來，并確定了Model 1為相對最佳地質(zhì)模型。

從地下水?dāng)?shù)值模擬結(jié)果可知不管抽注水井的位置如何變化，均能確定與已知理想地質(zhì)模型吻合最好的模型為Model 1，達(dá)到了優(yōu)化反演地質(zhì)解釋結(jié)果的目的。

5 結(jié)論與展望

本文通過學(xué)科交叉與融合，針對地球物理反演結(jié)果解釋的不確定性，探索了一種結(jié)合地下水模擬來優(yōu)化地球物理反演解釋思路。對于鉆探和地質(zhì)資料缺乏但具有較多可監(jiān)測的地下水井孔的地區(qū)，該方法可為地球物理反演結(jié)果解釋進(jìn)一步限定，提供一種新的思路和方法。主要得到以下兩點結(jié)論：

（1）提出了一種基于地下水模型的地球物理反演地質(zhì)解釋優(yōu)化方法。

首先，地球物理勘探數(shù)據(jù)的前期處理，并采用合理的地球物理反演方法對數(shù)據(jù)進(jìn)行反演成像。其次，結(jié)合已知地質(zhì)、鉆探資料，對反演成像的結(jié)果做出合理的地質(zhì)解釋；為了得到更合乎實際情況的地質(zhì)解釋，將反演成像所得的多個地質(zhì)解釋結(jié)果轉(zhuǎn)換成地下水模型，通過水文地質(zhì)調(diào)查，確定邊界定解條件、源匯項等相關(guān)水文信息，對地下水模型進(jìn)行地下水?dāng)?shù)值模擬，得到基于不同地質(zhì)解釋結(jié)果對應(yīng)地下水模型的模擬水位、濃度信息。最后，將模擬水位、濃度信息與水文地質(zhì)中水位、水質(zhì)監(jiān)測結(jié)果進(jìn)行對比，從而確定更符合野外實際情況的地質(zhì)解釋。

（2）通過地下水流、溶質(zhì)運移模擬結(jié)合的方式證明了方法的可行性。

文中以已知理想地質(zhì)模型和大地電磁法為例。首先，建立已知理想地質(zhì)模型，并對給定地質(zhì)模型進(jìn)行大地電磁正反演、得到多個基于大地電磁反演解釋的地質(zhì)模型。為了進(jìn)一步優(yōu)化反演解釋結(jié)果，將已知理想地質(zhì)模型以及基于反演解釋的地質(zhì)模型進(jìn)行水流和溶質(zhì)運移模擬，統(tǒng)計各項誤差，最終確定了相對最佳地質(zhì)模型。通過敏感性分析，不管抽注水井位置如何，均得到相同的結(jié)果，達(dá)到了優(yōu)化反演解釋的目的，很好地驗證了運用地下水流、溶質(zhì)運移數(shù)值模擬搜索相對最佳地質(zhì)模型的可行性。

本研究中反演結(jié)果解釋部分仍然存在一定的主觀性、局限性，在未來研究中需要做進(jìn)一步的完善，如通過大量的數(shù)值模擬實驗、實驗室測量建立電阻率與滲透系數(shù)的映射關(guān)系，或者將地下水約束加入地球物理反演算法中，從算法上優(yōu)化反演結(jié)果。

Andersen T R，Poulsen SE，Christensen S，et al.2013.A synthetic study of geophysics-based modelling of groundwater flow in catchments with a buried valley［J］.Hydrogeology Journal，21（2）：491～503.

Chen C X，Lin M.2009.Groundwater Hydraulics［M］.Wuhan：China University of Geosciences Press：1～300.

China Geological Survey.2012.Hydrogeological manual［M］.Beijing：Geological Publishing House：1－874.

Chen J，Wen N，Chen B Y.2007.Joint inversion of gravity-magneticelectrical－seismic combination survey：progress and prospect［J］. Progress in Geophysics，22（5）：1427～1438.

Chen S Y.2014.Research and application of integrated geophysical interpretation method［D］.Chengdu：Chengdu University of Technology.

Dai QW，Chen D P，Chen Y X，etal.2013.The enhanced algorithmsand its implementation for the abnormal response characteristics in electrical exploration［J］.Coal Geology＆Exploration，41（3）：75～78.

Dai SK，Xu S Z.1997.Rapid inversion ofmagnetotelluric data for 2-D and 3-D continuousmedia［J］.Oil Geophysical Prospecting，32（3）：305～317.

Guo D，Yan JY，Lu Q T，et al.2014.3D Densitymapping constrained by geological information：model study and application［J］.Acta Geologica Sinica，88（4）：763～776.

He W H，Wang JL，Yu P.2009.Overview of the status and prospect of geophysical joint inversion［J］.Progress in Geophysics，24（2）：530～540.

Jing R Z，Bao G S，Chen S Q.2003.A review of the researches for geophysical combinative inversion［J］.Progress in Geophysics，18（3）：535～540.

Li J M.2005.Electric field and electric prospecting［M］.Beijing：Geological Publishing House：1～473.

Ling D M，Sang Y J，Jin W J，et al.2014.Preliminary study of rockphysics and physical properties of rockmass structure［J］.Journal of Engineering Geology，22（6）：1154～1158.

Panday S，Langevin C D，Niswonger R G，et al.2013.MODFLOW-USG version 1：An unstructured grid version of MODFLOW for simulating groundwater flow and tightly coupled processes using a control volume finite-difference formulation［R］.Washington：US Geological Survey.

Qi G.2013.3D gravity and magnetic inversion modeling constrained by geological information［D］.Jilin：Jilin University.

Singha K，Gorelick S M.2006.Hydrogeophysical tracking of threedimensional tracer migration：The concept and application of apparent petrophysical relations［J］.Water Resources Research，42（6）：1～14.

Wang J Y.2002.Inverse theory in geophysics［M］.Beijing：Higher Education Press：1～187.

Wei C，Li X F，Zhang M G.2008.The geophysical inversemethod based on quantum Monte Carlo［J］.Chinese Journal of Geophysics Chinese，51（5）：1494～1502.

Wei Y Q，Qiao X J，Li G M.2015.Influence on calculation accuracy by differentalgorithms and parameters in MODFLOW［J］.Hydrogeology＆Engineering Geology，42（1）：14～21.

Wu G，Hu X Y，Huo G P，et al.2012.Geophysical exploration for geothermal resources：An application of MT and CSAMT in Jiangxia，Wuhan，China［J］.Journal of Earth Science，23（5）：757～767.

Xu S Z.1994.Finite element method in Geophysics［M］.Beijing：Science Press：1～308.

Xue Y Q，Xie C H.2007.Numerical simulation of groundwater［M］. Beijing：Science Press：1～451.

Yang W C.2002.A perspective to development of geophysical inversions［J］.Earth Science Frontiers，9（4）：389～396.

Yeh T C J，Liu S，Glass R J，et al.2002.A geostatistically based inverse model for electrical resistivity surveys and its applications to vadose zone hydrology［J］.Water Resources Research，38（12）：14－1－14－13.

Yu P，Wang J L，Wu J S，et al.2006.Review and discussions on geophysical joint inversion.［J］.Progress in Exploration Geophysics，29（2）：87～93.

Yu Y X，Wu JC.2006.Application of electrical resistivity tomography data to estimate hydraulic conductivity in porous medium［J］. Hydrogeology＆Engineering Geology，33（2）：41～44，49.

Yuan SY，Wang SX，Tian N.2009.Swarm intelligence optimization and its application in geophysical data inversion［J］.Applied Geophysics，6（2）：166～174.

Zheng C X，Wang P P.1999.A modular three-dimensionalmulti-species transportmodel for simulation of advection，dispersion and chemical reactions of contaminants in groundwater systems；documentation and user's guide［R］.Vicksburg，Mississippi，USA：US Army Engineer Research and Development Center Contract Report SERDP-99－1.

Zheng C Z，Tong L Y，ZhangM F，etal.2015.Characteristics of effects on groundwater flow by underground structures［J］.Journal of Engineering Geology，23（S1）：553～557.

陳崇希，林敏.2009.地下水動力學(xué)［M］.武漢：中國地質(zhì)大學(xué)出版社：1～300.

陳潔，溫寧，陳邦彥.2007.重磁電震聯(lián)合反演研究進(jìn)展與展望［J］.地球物理學(xué)進(jìn)展，22（5）：1427～1438.

陳思宇.2014.地質(zhì)－綜合地球物理聯(lián)合解釋方法研究及應(yīng)用［D］.成都：成都理工大學(xué).

戴前偉，陳德鵬，陳勇雄，等.2013.電法勘探中異常響應(yīng)特征的增強算法及其實現(xiàn)［J］.煤田地質(zhì)與勘探，41（3）：75～78.

戴世坤，徐世浙.1997.MT二維和三維連續(xù)介質(zhì)快速反演［J］.石油地球物理勘探，32（3）：305～317.

郭冬，嚴(yán)加永，呂慶田，等.2014.地質(zhì)信息約束下的三維密度填圖技術(shù)研究及應(yīng)用［J］.地質(zhì)學(xué)報，88（4）：763～776.

何委徽，王家林，于鵬.2009.地球物理聯(lián)合反演研究的現(xiàn)狀與趨勢分析［J］.地球物理學(xué)進(jìn)展，24（2）：530～540.

敬榮中，鮑光淑，陳紹裘.2003.地球物理聯(lián)合反演研究綜述［J］.地球物理學(xué)進(jìn)展，18（3）：535～540.

李金銘.2005.地電場與電法勘探［M］.北京：地質(zhì)出版社：1～473.

林達(dá)明，尚彥軍，金維浚，等.2014.巖石物性與巖體結(jié)構(gòu)物性研究進(jìn)展［J］.工程地質(zhì)學(xué)報，22（6）：1154～1158.

祁光.2013.地質(zhì)條件約束下重磁三維反演建模方法研究［D］.吉林：吉林大學(xué).

王家映.2002.地球物理反演理論［M］.北京：高等教育出版社：1－187.

魏超，李小凡，張美根.2008.基于量子蒙特卡羅的地球物理反演方法［J］.地球物理學(xué)報，51（5）：1494～1502.

魏亞強，喬小娟，李國敏.2015.MODFLOW不同算法及參數(shù)設(shè)定對計算精度的影響［J］.水文地質(zhì)工程地質(zhì)，42（1）：14～21.

徐世浙.1994.地球物理中的有限單元法［M］.北京：科學(xué)出版社：1～308.

薛禹群，謝春紅.2007.地下水?dāng)?shù)值模擬［M］.北京：科學(xué)出版社：1～451.

楊文采.2002.評地球物理反演的發(fā)展趨向［J］.地學(xué)前緣，9（4）：389～396.

于鵬，王家林，吳健生，等.2006.地球物理聯(lián)合反演的研究現(xiàn)狀和分析［J］.勘探地球物理進(jìn)展，29（2）：87～93.

喻永祥，吳吉春.2006.利用ERT數(shù)據(jù)推求非均質(zhì)多孔介質(zhì)滲透系數(shù)初探［J］.水文地質(zhì)工程地質(zhì)，33（2）：41～44，49.

鄭燦政，童立元，張明飛，等.2015.地下水結(jié)構(gòu)對地下水流場影響的特性規(guī)律［J］.工程地質(zhì)學(xué)報，23（S1）：553～557.

中國地質(zhì)調(diào)查局.2012.水文地質(zhì)手冊［M］.北京：地質(zhì)出版社：1～874.

AN OPTIM IZATION METHOD FOR INTERPRETING ERT DATA BASED ON GROUNDWATER NUMERICAL MODELING

LIXuelan①DONG Yanhui②LIDiquan①③
（①School of Geosciences and Info-Physics，Central South University，Changsha 410083）
（②Institute of Geology and Geophysics，Chinese Academy of Sciences，Beijing 100029）
（③Key Laboratory of Metallogenic Prediction of NonferrousMetals，Ministry of Education，Central South University，Changsha 410083）

Geophysical non-uniqueness are everywhere.Faced with this situation we usually tend to combine geological data to calibrate the interpretation results.However，the information on geology and deep drilling are deficiency in remote areas，and field existing wells are often difficult to obtain information on lithology catalog，which will cause that the inversion results are difficult to interpret.Whether hydraulic characteristics of the water level and water quality information contained in the geological body can be used to assist the geophysical inversion？Based on this idea，amethod for interpreting ERT data based on groundwater numericalmodeling is provided in this article.In this article，Magneto-telluric method is a way to accomplish numerical simulation of geophysical forward and inversion.Firstly，by studying the known geologicalmodel andmultiple geologicalmodels interpreted from ERTdata，a set of groundwater numericalmodels can be set up.Then numerical simulation of groundwater flow and solute transportwere conducted.Lastly，the optimal geologicalmodel can be obtained by comparing the result of numerical simulation of groundwater.Sensitivity analysiswere also conducted to verify themethod.

Geophysical inversion，Geologicalmodel，Numerical simulation of groundwater flow，Solute transport

P641.7

：A

10.13544／j.cnki.jeg.2016.06.037

2015－10－15；

2016－06－16.

國家自然科學(xué)基金（41202182）資助.

李學(xué)蘭（1993－），女，碩士生，從事地球探測信息與技術(shù)研究.Email：lixuelan＠csu.edu.cn

董艷輝（1980－），男，博士，副研究員，主要從事水文地質(zhì)方面研究.Email：lemondyh＠m(xù)ail.iggcas.ac.cn