崔益安 ，朱肖雄 ，陳志學 ，柳建新

(1. 中南大學 地球科學與信息物理學院, 長沙 410083； 2. 中南大學 有色資源與地質災害探查湖南省重點實驗室, 長沙 410083)

隨著我國經濟的持續(xù)高速發(fā)展，在工業(yè)化、城鎮(zhèn)化過程中的環(huán)境污染問題也日益嚴重。特別是越來越多的工礦廢水廢渣庫、生活和工業(yè)垃圾填埋場等，嚴重威脅周邊的土壤與地下水安全[1]，即時檢測和控制污染源擴散是保護土壤和地下水不受污染或少受污染的積極方法[2]。與鉆孔測試、定期采樣分析等傳統(tǒng)環(huán)境監(jiān)測手段相比，直流電阻率法等地球物理方法具有效率高、成本低、快速無損、無二次污染等優(yōu)點[3]， 在地下污染探測的應用中有較強的優(yōu)勢。王迪等[4]就采用電阻率成像法對北京某垃圾場進行了污染探測，取得不錯的效果；VAUDELET 等[5]也采用電阻率法成功地在城市地區(qū)實現(xiàn)對污染物的成像。類似的研究與應用還有很多[6-9]。隨著污染控制與治理的進一步需求，電阻率法在污染物探測的應用中正逐步由二維觀測向三維觀測發(fā)展，同時也由短時單次觀測向長期連續(xù)監(jiān)測發(fā)展，如尚浩等[10]和郭秀軍等[11]研制開發(fā)了三維分布式監(jiān)測系統(tǒng)用于垃圾填埋場滲漏的監(jiān)測；GASPERIKOVA 等[12]采用電阻率法對地下污染活動進行長期監(jiān)測。在這些電阻率法污染監(jiān)測的應用中，即便是采用了三維觀測，但其觀測數據的反演解釋還是基本以二維處理為主[4,10]。污染監(jiān)測數據的反演解譯至關重要，只有開展三維反演才能滿足污染監(jiān)測中的高精度數據處理要求。

因此，有必要針對直流電阻率法在地下污染監(jiān)測中的應用開展監(jiān)測數據的快速三維反演算法研究，以滿足實際應用中對數據解譯的精確性和時效性要求。正演模擬是反演計算的基礎，潘克家等[13]采用多重網格法以及黃俊革等[14]采用有限單元法進行了直流電阻率法正演模擬，獲得較高的模擬精度，但這類正演方法計算量巨大，比較耗時。ZHANG 等[15]采用模擬阻抗網絡開展直流電阻率法正演模擬研究，獲得的阻抗剛度矩陣除對角元素和次對角元素外全為0，求解簡便，便于提高正演計算速度。許多研究人員在直流電阻率反演方面開展了大量研究。戴前偉等[16]將智能優(yōu)化方法應用到電阻率層析成像的反演計算中。LI 等[17]和劉斌等[18]研究了通過約束來提高三維電阻率反演收斂速度的方法。但這些反演方法普遍耗時較長，不適合污染監(jiān)測數據實時快速解譯的需求。梯度法是基于求偏導數來尋找極值點，是一類迭代次數少、收斂快的反演方法。吳小平等[19]采用共軛梯度法開展電阻率三維反演研究，計算速度較快，取得不錯的反演效果。黃俊革等[20]更是針對坑道超前探測專門研究了直流電阻率的快速反演方法。污染監(jiān)測數據同樣需要實時快速解譯以便及時提供決策依據，因此，本文作者開展相關研究，在采用模擬阻抗網絡進行正演模擬提高計算速度的基礎上，通過共軛梯度快速迭代來實現(xiàn)對污染監(jiān)測電阻率數據的快速三維反演。

1 模擬阻抗網絡模型

正演模型是開展反演研究的基礎。一般情況下，大地介質中的直流電阻率法問題滿足式(1)

正演計算的目的是通過給定介質的三維電阻率分布ρ(x,y,z)以及電流密度J(x,y,z)來計算空間中的電勢分布V(x,y,z)。計算過程中需要先將模型離散化，然后進行求解。圖1 所示為模擬阻抗網絡模型。由圖1 可知，模擬阻抗網絡模型將地電模型分割成許多六面體電阻塊，每個塊在x、y、z 3 個方向的長度分別為 Δx 、Δz 。把空間電勢分布與電流密度分布均定義在每個塊的頂面中心點，稱之為網絡節(jié)點。除邊界外，每個節(jié)點在三維空間中都有6 個相鄰的節(jié)點通過阻抗臂與之連接，每條阻抗臂的阻抗由所在塊的電阻率以及節(jié)點之間的距離決定。

圖1 模擬阻抗網絡模型 Fig. 1 Model of analogy impedances network

如將圖1 中的中心立方塊在空間中的位置定義為(i, j, k)，電阻率為ρ(i, j, k)，其網絡節(jié)點電壓為V(i, j, k)，則節(jié)點(i, j, k)與節(jié)點(i-1, j, k)之間阻抗臂的阻抗Rx為

節(jié)點(i, j, k)與節(jié)點(i, j-1, k)之間阻抗臂的阻抗Ry為

節(jié)點(i, j, k)與節(jié)點(i, j, k-1)之間阻抗臂的阻抗Rz為

為滿足紐曼邊界條件(n·j=0)，處于地表的節(jié)點沒有與空氣相連的阻抗臂，地下邊界節(jié)點的設置則需要滿足第一類邊界條件。由于空間中各節(jié)點的電勢與電流均滿足基爾霍夫定律，將式(1)中電勢梯度用模擬電阻網絡中相鄰節(jié)點的電勢差分代替后有

式中：K 為表示阻抗分布結構的剛度矩陣。矩陣K 的對角元素為與其對應節(jié)點相連的所有電導之和，次對角元素為對應節(jié)點與相連節(jié)點之間的負電導，其他元素均為0。通過式(5)求解便可以獲得各節(jié)點的電勢 V(x, y, z)，即實現(xiàn)正演模擬。

2 共軛梯度反演算法

算數據f(m)之間的殘差； mΔ 為每次迭代的模型增量。

有許多方法可以用來求解迭代方程(9)，其中許多方法均需要完全計算靈敏度矩陣。而在采用共軛梯度的迭代法中，只需要計算靈敏度矩陣A 或其轉置矩陣與一個向量X 的乘積結果，因為不用完全計算靈敏度矩陣，大大節(jié)約了計算時間和存儲空間，如式(10)所示：

反演就是通過一定的策略來估計符合要求的模型參數。令模擬電阻網絡離散化后的地電模型 m= (ρ1, ρ2, …, ρi×j×k)T，對于存在n 個視電阻率觀測數據dob=(ρs1, ρs2, …, ρsn)T，反演問題即為尋找適當的模型m 以滿足

采用共軛梯度法快速反演的具體迭代算法如下：

1) 設定最大反演和共軛梯度迭代次數imax和jmax，設定允許的擬合差ε，給定初始模型m0與參考模型m0；

3) 如果 idΔ ≤ε 則跳轉到8)；

4) 計算 bi= ATΔdi；

5) 共軛梯度循環(huán)for j=1 to jmax：

式中：f(m)為模型m 的正演計算數據；ε 為允許的擬合差。由于正演模型與實際物理模型的差異以及觀測數據的噪音影響等因素，這類問題往往沒有確定解或存在多解。地球物理反演問題的不適定性在引入Tikhonov 正則化后，可以得到適當的改善，求解穩(wěn) 定[21]，如式(7)所示：

最終，可以獲得反演問題的最小二乘迭代計算方程，如式(9)所示：

式中：A 為靈敏度矩陣； dΔ 為觀測數據dob與正演計

4月26日上午，十一屆全國人大常委會第二十六次會議在人民大會堂舉行聯(lián)組會議，專題詢問國務院關于農田水利建設工作情況。受國務院委托，水利部、發(fā)展改革委、財政部、國土資源部、農業(yè)部、銀監(jiān)會等六部委負責人到會聽取意見，回答詢問（本期“特別關注”全文刊發(fā)）。這是全國人大常委會2012年舉行的第一次專題詢問。受吳邦國委員長委托，烏云其木格副委員長主持會議。

7) 如沒有達到最大反演次數imax則跳轉到2)；

8) 反演結束。

3 計算模擬與實驗測試

3.1 模擬數據測試

首先采用模擬合成數據對反演算法進行了測試。設三維模型xy 平面由8×10 個網格組成，z 方向為深度分8層，共計640 個網格，每個網格尺寸均為1.1 m×1.1 m×0.46 m。圖2(a)所示為模擬數據模型1。模型1 的第1 層x 與y 方向的第1、2 網格電阻率為10 ?·m，其余網格電阻率均為200 ?·m，低阻網格用來模擬地表存在的污染源。圖2(b)所示為模擬數據模型2。模型2 的第1 層與第2 層x 方向第1、2 網格以及y 方向的第1、2、3 網格電阻率為10 ?·m，模擬地表污染擴散至這些區(qū)域，其余網格電阻率為200 ?·m。

設視電阻率觀測采用二極裝置，電極在地表xy平面布置，共布置5×6 個電極，電極距為2.1 m。然后，利用模擬阻抗網絡數值方法進行正演模擬，計算出各電極點位置的二級裝置電壓與供電電流的比值。利用計算獲得的模擬觀測數據對三維反演算法進行測試。由于采用模擬阻抗網絡進行正演模擬以及采用共軛梯度反演算法，網格內無需插值，反演計算時間復雜度為一階o(n)，耗時較少，640 個網格的計算在普通個人計算機上(intel 2.1 GHz CPU)上的計算耗時僅為采用有限單元法耗時的一半。圖3 所示為反演計算測試結果，圖3(a)和(b)分別為模型1 和模型2 模擬數據的反演三維視電阻率結果的水平切片。圖3 中只給出了深度分別為0、0.5、0.9、1.4 和1.8 m 共5 個切片，重點研究反演三維視電阻率結果對模型電阻率發(fā)生變化區(qū)域的響應情況。

反演結果能很好的反映模型1 與模型2 之間的電阻率變化情況，并且對電阻率變化的區(qū)域與電阻率數值大小均有比較精確的響應，表明該反演方法可以有效實現(xiàn)對直流電阻率數據的快速三維反演成像計算。

3.2 實測數據反演

為進一步驗證反演方法的有效性，還開展了實測數據的計算測試。實測數據來自某場地進行鹽水注射過程中的視電阻率監(jiān)測實驗。實驗在鉆孔中注入高濃度的鹽水來模擬污染物在地下介質中的擴散情況。視電阻率監(jiān)測采用AGI 公司生產的Sting R1 型電法儀器，通過自動電極轉換器同時布置30 個電極，開展二極裝置視電阻率成像觀測。試驗場地監(jiān)測區(qū)域大小約為10 m×12 m。圖4 所示為視電阻率監(jiān)測電極布置。x 和y 方向電極距均為2.1 m。二級裝置的一個供電電極和一個測量電極均布置在遠離監(jiān)測區(qū)域的“無窮遠”處，圖4 中的各測點依次供電進行循環(huán)掃描觀測。

圖2 模擬數據模型 Fig.2 Models of synthetic data: (a) Model 1; (b) Model 2

圖3 模擬數據反演視電阻率結果 Fig. 3 Results of inversing synthetic resistivity data: (a) Model 1; (b) Model 2

圖4 視電阻率監(jiān)測電極布置 Fig. 4 Survey array of apparent resistivity monitoring

鹽水注射孔深度1 m，位于場地一側2 個電極之間(見圖4)，在注射鹽水前先進行了一次二極法視電阻率掃描觀測，然后開始注射鹽水。注射鹽水后的第2天開始，每天進行一次視電阻率掃描觀測，監(jiān)測用來 模擬污染物的鹽水在地下介質中的擴散情況和影響范圍。對每次采集到的視電阻率觀測數據進行三維反演計算，獲得一系列的三維反演成像結果，便可以直觀地了解鹽水連續(xù)擴散的具體情況。圖5 所示為3 次不同時間監(jiān)測數據進行三維反演成像的計算結果。圖 5(a)為開始注射鹽水前觀測數據的三維視電阻率反演結果水平切片，深度分別為0、0.9、1.8、2.7 和3.6 m。反演結果表明該場地淺地表具有較明顯的電性分層，表層厚度約1 m 左右電阻率稍高，隨著深度增加電阻率逐步降低。圖5(b)和(c)分別為注射鹽水后第3 天和第6 天觀測數據的三維視電阻率反演結果。從圖5(b)和(c)中可以明顯看出，由于鹽水注射導致的低電阻率區(qū)域。對比圖5(a)～(c)可以很直觀地了解鹽水在地下的滲透擴散情況，由于鹽水滲透擴散引起的低電阻率區(qū)域主要沿注射孔底部向四周擴散，并且在水平y(tǒng) 方向的擴散速度高于其他方向的速度。測試表明：反演算法對實際監(jiān)測數據的反演計算也具有良好的應用效果，可以提供一系列直觀的三維視電阻率影像，快捷地實現(xiàn)對污染物視電阻率監(jiān)測數據的動態(tài)實時解譯。

圖5 實測數據反演視電阻率結果 Fig. 5 Results of inverting field resistivity data at different time: (a) Before injection; (b) 3 days later after injection; (c) 6 days later after injection

4 結論

1) 視電阻率快速三維反演方法能有效地實現(xiàn)對污染監(jiān)測數據的快速有效反演，使得三維反演可在普通計算機上短時內完成，有利于提高監(jiān)測數據解釋的時效性。

2) 采用模擬阻抗網絡模型對地下污染地電模型進行數值模擬，在保證模擬精度的同時減少正演計算量，為快速反演提供良好的模型基礎。

3) 共軛梯度迭代過程中只需要計算靈敏度矩陣或其轉置矩陣與向量的乘積，而不用完全計算靈敏度矩陣，大大節(jié)約計算時間和存儲空間，使得反演算法可以快速實現(xiàn)對監(jiān)測數據的反演。

4) 針對污染監(jiān)測視電阻率數據的快速反演方法，能實時直觀地提供一系列準確反映地下污染擴散動態(tài)范圍的三維視電阻率分布，可為地下污染的防控與處理及時提供可靠的依據。

致謝：

感謝劉瀾波教授提供場地實驗及相關監(jiān)測數據。

[1] 劉 月, 董穎博, 林 海, 陳月芳, 于明利. 云南某典型錫礦選礦廠重金屬污染特征[J]. 中國有色金屬學報, 2014, 24(4): 1084-1090. LIU Yue, DONG Ying-bo, LIN Hai, CHEN Yue-fang, YU Ming-li. Pollution characteristics of heavy metals in typical tin mineral processing plant of Yunnan, China[J]. The Chinese Journal of Nonferrous Metals, 2014, 24(4): 1084-1090.

[2] 程業(yè)勛, 楊 進, 趙章元. 環(huán)境地球物理學的現(xiàn)狀與發(fā)展[J]. 地球物理學進展, 2007, 22(4): 1364-1369. CHENG Ye-xun, YANG Jin, ZHAO Zhang-yuan. Status and development of environmental geophysics[J]. Progress in Geophysics, 2007, 22(4): 1364-1369.

[3] 鄭軍衛(wèi), 張志強, 董連成. 環(huán)境地球物理學及其現(xiàn)狀與進展[J]. 地球科學進展, 2000, 15(1): 40-47. ZHENG Jun-wei, ZHANG Zhi-qiang, DONG Lian-cheng. Environmental geophysics and its advance[J]. Advance in Earth sciences, 2000, 15(1): 40-47.

[4] 王 迪, 尉 斌, 閆永利, 馬曉冰, 王顯祥, 楊京勤, 魯笑穎. 垃圾填埋場地下環(huán)境污染檢測方法技術研究[J]. 地球物理學報, 2012, 55(6): 2115-2119. WANG Di, WEI Bin, YAN Yong-Li, MA Xiao-Bing, WANG Xian-xiang, YANG Jing-qin, LU Xiao-ying. Study of detecting technology for landfill underground environmental pollution[J]. Chinese Journal of Geophysics, 2012, 55(6): 2115-2119.

[5] VAUDELET P, SCHMUTZ M, PESSEL M, FRANCESCHI M, GUéRIN R, ATTEIA O, BLONDEL A, NGOMSEU C, GALAUP S, REJIBA F, BéGASSAT P. Mapping of contaminant plumes with geo-electrical methods. A case study in urban context[J]. Journal of Applied Geophysics, 2011, 75(4): 738-751.

[6] GALLAS J D F, TAIOLI F, MALAGUTTI FILHO W. Induced polarization, resistivity, and self-potential: A case history of contamination evaluation due to landfill leakage[J]. Environmental Earth Sciences, 2011, 63(2):251-261.

[7] ABDULLAHI N K, OSAZUWA I B. Geophysical imaging of municipal solid waste contaminant pathways[J]. Environmental Earth Sciences, 2011, 62(6): 1173-1181.

[8] METWALY M, KHALIL M A, AL-SAYED E S, EL-KENAWY A. Tracing subsurface oil pollution leakage using 2D electrical resistivity tomography[J]. Arabian Journal of Geosciences, 2013, 6(9): 3527-3533.

[9] OLOWOFELA J A, AKINYEMI O D, OGUNGBE A S. Imaging and detecting underground contaminants in landfill sites using electrical impedance tomography (EIT): A case study of Lagos, southwestern, Nigeria[J]. Research Journal of Environmental and Earth Sciences, 2012, 4(3): 270-281.

[10] 尚 浩, 郭秀軍, 王仙麗, 孟慶生. 垃圾填埋場滲漏污染三維分布式電學監(jiān)測系統(tǒng)研制及應用[J]. 地球物理學進展, 2010, 25(4): 1485-1491. SHANG Hao, GUO Xiu-jun, WANG Xian-li, MENG Qing-sheng. Design and experiment of the 3-D electrical resistivity detecting system for monitoring underground soil contaminated by landfill leachate[J]. Progress in Geophysics, 2010, 25(4): 1485-1491.

[11] 郭秀軍, 魏 麗, 賈永剛, 黃瀟雨. 垃圾填埋場滲濾液污染地下含水層及修復過程的三維動態(tài)監(jiān)測實驗[J]. 地球物理學進展, 2006, 21(2): 637-642. GUO Xiu-jun, WEI Li, JIA Yong-gang, HUANG Xiao-yu. Experiment of monitoring underground aquifers contaminated by leachate of landfill and remation with 3-D electrical method[J]. Progress in Geophysics, 2006, 21(2): 637-642.

[12] GASPERIKOVA E, HUBBARD S S, WATSON D B, BAKER G S, PETERSON J E, KOWALSKY M B, SMITH M, BROOKS S. Long-term electrical resistivity monitoring of recharge-induced contaminant plume behavior[J]. Journal of Contaminant Hydrology, 2012, 142/143(11): 33-49.

[13] 潘克家, 湯井田, 胡宏伶, 陳傳淼. 直流電阻率法2.5 維正演的外推瀑布式多重網格法[J]. 地球物理學報, 2012, 55(8): 2769-2778. PAN Ke-jia, TANG Jing-tian, HU Hong-ling, CHEN Chuan-miao. Extrapolation cascadic multigrid method method for 2.5D direct current resistivity modeling[J]. Chinese Journal of Geophysics, 2012, 55(8): 2769-2778.

[14] 黃俊革, 阮百堯, 鮑光淑. 基于有限單元法的三維地電斷面電阻率反演[J]. 中南大學學報(自然科學版), 2004, 35(2): 295-299. HUANG Jun-ge, RUAN Bai-yao, BAO Guang-shu. Resistivity inversion on 3-D section based on FEM[J]. Journal of Central South University (Natural Science), 2004, 35(2): 295-299.

[15] ZHANG J, MACKIE R L, MADDEN T R. 3-D resistivity forward modeling and inversion using conjugate gradients[J]. Geophysics, 1995, 60(5): 1313-1325.

[16] 戴前偉, 江沸菠. 基于混沌振蕩PSO-BP算法的電阻率層析成像非線性反演[J]. 中國有色金屬學報, 2013, 23(10): 2897-2904. DAI Qian-wei, JIANG Fei-bo. Nonlinear inversion for electrical resistivity tomography based on chaotic oscillation PSO-BP algorithm[J]. The Chinese Journal of Nonferrous Metals, 2013, 23(10): 2897-2904.

[17] LI Shu-cai, NIE Li-chao, LIU Bin, SONG Jie, LIU Zheng-yu, SU Mao-xin, XU Lei. 3D electrical resistivity inversion using prior spatial shape constraints[J]. Applied Geophysics, 2013, 10(4): 361-372.

[18] 劉 斌, 李術才, 李樹忱, 聶利超, 鐘世航, 李利平, 宋 杰, 劉征宇. 基于不等式約束的最小二乘法三維電阻率反演及其算法優(yōu)化[J]. 地球物理學報, 2012, 55(1): 260-268. LIU Bin, LI Shu-cai, LI Shu-chen, NIE Li-chao, ZHONG Shi-hang, LI Li-ping, SONG Jie, LIU Zheng-yu. 3D electrical resistivity inversion with least-squares method based on inequality constraint and its computation efficiency optimization[J]. Chinese Journal of Geophysics, 2012, 55(1): 260-268.

[19] 吳小平, 徐果明. 利用共軛梯度法的電阻率三維反演研究[J]. 地球物理學報, 2000, 43(3): 420-427. WU Xiao-ping, XU Guo-ming. Study on 3-D resistivity inversion using conjugate gradient method[J]. Chinese Journal of Geophysics, 2000, 43(3): 420-427.

[20] 黃俊革, 阮百堯, 王家林. 坑道直流電阻率法超前探測的快速反演[J]. 地球物理學報, 2007, 50(2): 619-624. HUANG Jun-ge, RUAN Bai-yao, WANG Jia-lin. The fast inversion for advanced detection using DC resistivity in tunnel[J]. Chinese Journal of Geophysics, 2007, 50(2): 619-624.

[21] 朱自強, 曹書錦, 魯光銀. 基于混合正則化的重力場約束反演[J]. 中國有色金屬學報, 2014, 24(10): 2601-2608. ZHU Zi-qiang, CAO Shu-jin, LU Guang-yin. 3D gravity inversion with bound constraint based on hyper-parameter regularization[J]. The Chinese Journal of Nonferrous Metals, 2014, 24(10): 2601-2608.