2016, Vol. 31
2. 中南大学 有色资源与地质灾害探查湖南省重点实验室, 长沙 410083
2. Hunan Key Laboratory of Non-ferrous Resources and Geological Hazard Detection, Central South University, Changsha 410083, China
自然电场法是一种快速有效的地球物理方法，已经广泛地应用于矿产油气勘查(杨春梅等, 2005; 宋子齐等, 2009; 李虎等, 2013)、地震研究(范莹莹等, 2010; 谭大诚等, 2012, 2014)地热火山研究(张贵宾等, 1999; Aizawa, 2008; Mauri et al., 2010)、岩石学(张朋等, 2013)等领域.同时在环境工程和水利工程领域得到了有效的应用，在渗漏检测(郭秀军等, 2007)、地下水调查(Jardani et al., 2009; 刘加文等, 2009)、水位监测(Martínez-Pagán et al., 2010)和孔洞检测(Jardani et al., 2007)的应用效果都很优异.由于自然电场法对与污染伴生的渗流扩散、氧化还原以及微生物活动信号相当敏感(Rittgers et al., 2013)，近年来，在污染检测方面的研究应用越来越广泛.尚浩等(2010)通过自然电场法监测了垃圾的泄露情况，Edmundo通过自然电场法成功地监测了尾砂库废液的泄漏泄漏情况(Placencia-Gómez et al., 2010)，Che-Alota则利用其监测了烃类污染物的情况(Che-Alota et al., 2009).
因此展开自然电场异常的反演研究是有必要的.最初，被用来解释自然电位数据的方法有特征点法、曲线拟合法，特征点法是通过少数异常点来比较，曲线拟合法则是需要与理论曲线相比较，此类方法都比较繁琐且效果并不理想.随着研究的发展，二阶移动平均法(Abdelrahman et al., 2009)，小波分析(Sailhac and Marquis, 2001)，直接差分方法(Skianis et al., 2006), 最小二乘法(Abdelrahman et al., 2008)等方法被用来解释自然电位数据，但此类方法都具有局限性或只能定性分析.近年来，一些非线性算法如(El-Kaliouby and Al-Garni, 2009)神经网络算法、遗传算法、模拟退火算法(Göktürkler and Balkaya, 2012)、粒子群优化算法(Monteiro Santos, 2010; Pekşen et al., 2011; 朱肖雄等, 2015)被应用于自然电场的反演，但此类方法多应用于简单模型组合的解析解，不能满足复杂的模型，而在反演参数较多时收敛很慢，精度不高.总之，自然电位数据的反演解释依旧有待研究.因此，本文作者开展相关研究，把目标区域内的多种自然电位源统一处理，在采用有限单元进行高精度正演模拟的基础上，通过最小二乘正则化反演对自然电位异常进行二维反演成像.
1 自然电位源模型自然电位是由电子与非电子之间的相对流动和驱动力的耦合作用所产生.电流密度j=q1[A·m-2]的产生可以包括四种耦合作用，如公式(1)所示，其中jk(x)表示由流动力产生的流动电流密度，jd(x)表示由化学动力产生的扩散电流密度，jt(x)表示由热力产生的热电流密度.jc(x)表示由电场力产生的传导电流密度.公式(1)为
|
(1) |
通常自然电场的形成过程都非常缓慢的.因此，自然电场的问题适用于准静态条件下的麦克斯韦方程，对安培定律的两端取散度得公式(2), 再代入高斯定律公式，在准静态的条件下，电荷密度的时间导数项可以被忽略，得到了电流守恒公式，如公式(3)所示.其中表示E电场强度，H表示磁场强度，ε[F·m-1]表示介电常数，ρq[C·m-3]表示电荷密度，表示总电流密度.公式(2)和(3)为
|
(2) |
|
(3) |
通常自然电场的产生含有一个或多个耦合现象，因此可把这些耦合现象当成一个整体来处理.把公式(1)代入公式(3)中，移项等式两端取散度可得自然电位的微分方程，如公式(4)所示，其中σ表示电导率，φ表示电位，js(x)=jk(x)+jd(x)+jt(x)，P(x)[A·m-3]表示源项.公式(4)为
|
(4) |
对于二维的地电构造，设z轴平行于构造走向，即电导率σ在z方向保持不变，用一个足够大的空间代替半无限空间，并在地表Γs和Γ∞无穷远处添加边界条件，可以把自然电位的边值问题归结为公式(5).其中u表示电位，δ(x)表示狄克函数，l是点源大小，n是Γ∞处的外法向方向，c是比例系数，r′是点源至无穷远边界Γ∞的距离，xi, yi, zi表示场源在空间的位置.公式(5)为
|
(5) |
通过傅里叶变换可以将它变成二维的问题，如公式(6)所示.其中r表示的是二维平面内点源至无穷远边界Γ∞的距离，U表示k表示傅里叶变化系数，K0和K1分别表示的是第二类零阶和一阶修正贝塞尔函数.公式(6)为
|
(6) |
有限单元法解边值问题需要把问题转化成相应的变分问题，对区域进行网格剖分，再通过解相应的高阶线性方程组得到网格节点处的电位.构造泛函，带入边界条件，公式(6)的边值问题可转换成下列变分问题, 公式为
|
(7) |
有限元是通过有界函数逼近模拟，为了提高模拟精度，减少地表测量值受到点源奇异性和边界问题的影响，网格剖分如图 1所示，点源只存在于目标区域.采取矩形单元，双线性插值，通过单元分析将公式(7)化为各小单元的积分得到公式(8), 其中K 1θ, K 2θ, K 3θ为与电导率、网格长度相关的系数矩阵.
|
图 1 自然电场二维正反演网格剖分示意图 Figure 1 Model of the finite element method network |
把单元节点扩展成全体节点组成的矩阵，全部单元相加, 并求变分令其等于零，可以得到公式(9).其中表示合成后得扩展矩阵，P表示源项:P=(0, 0, …, 0, I/2, 0, …)T，通过解方程组可以得到傅里叶空间的U.本文采用了变带宽矩阵储存节约储存空间及高效的LDL T分解法求解方程组.公式(8)为
|
(8) |
|
(9) |
解得的U值是对应于特定傅里叶变换系数kj的，因此还需要通过傅里叶反变换把U值转换成三维空间的电位u值, 如公式(10)所示，其中gj表示傅里叶反变换系数.本文中采用傅里叶变换系数为k1=0.004758, k2=0.0407011, k3=0.1408855, k4=0.393225, k5=1.088038。g1=0.0099472, g2=0.0381619, g3=0.0980327, g4=0.2511531, g5=0.7260814.公式(10)为
|
(10) |
自然电位的产生通常地下含有多个点源，为了进一步提高模拟精度，减少测量值受到点源奇异性和边界的影响，我们对产生的自然电位进行校正.校正公式如公式(11)所示.其中ucorr表示校正后的电位, uH和uh分别表示在σ=σ
|
(11) |
最小二乘正则化反演是一种线性迭代的方法，总目标函数如式(12)所示，其中φ(m)为数据目标函数.α是正则化因子，s(m)是模型约束目标函数，公式为
|
(12) |
在二维自然电位的反演问题中，目标函数如式(13)所示.式中dobs为实测数据，F(m)为正演响应函数. W m为光滑度矩阵，m ref为先验模型，在没有先验模型的情况下m ref取为零.公式(13)为
|
(13) |
通过泰勒公式把F(m)展开，略去二次项以上的多项式，如式(14)所示，其中m k为模型的第k次迭代值，F(m k)=d k，Δ m=m k+1-m k，J k是雅克比灵敏度矩阵.公式(14)为
|
(14) |
本文中采用伴随矩阵的方法来计算灵敏度矩阵，J k的定义如公式所示.将公式(9)和公式(10)带入公式(15)推导可得公式(16)为
|
(15) |
|
(16) |
把式(14)代入式(13)，对Δ m进行求导并等于零，可推得Δm的表达公式(17)，通过公式(18)迭代可以获得最终解.本文采用的是最光滑模型约束.公式(18)为
|
(17) |
|
(18) |
采用模拟合成数据对反演算法进行了试算，模型均采用28×15的网格剖分，其中目标区域为14×7，目标区域单元网格的边长为5 m，反演的自然电场场源的节点为15×8.模型一，偶极源模型：ρ1=10 Ω·m, h1=30 m, ρ2=100 Ω·m均匀两层介质中，在h0=10 m处有一对正负点异常源I=±2 mA，其中两点源间距为35 m，模型示意图如图 2a所示，地表自然电位异常如图 2b所示.模型二，多源模型：ρ1=30 Ω·m, h1=30 m, ρ2=100 Ω·m均匀两层介质中，左侧h0=10 m处有异常源I=-3 mA，右侧有一组斜排列的异常源I=+1 mA，模型示意图如图 3a所示，地表自然电位异常如图 3b所示.
|
图 2 模型一 (a)模型一示意图; (b)模型一模拟异常曲线. Figure 2 Model 1 (a) Model 1;(b) Synthetic SP anomaly of model 1. |
|
图 3 模型二 (a)模型二示意图; (b)模型二模拟异常曲线. Figure 3 Model 2 (a) Model 2;(b) Synthetic SP anomaly of model 2. |
对模型一进行反演试算，通过地表自然电位异常反演所得电流拟断面等值线图如图 4a所示，由图 4a可知, 根据电流拟断面正负等值线的高低可以推断自然电场场源的位置和深度，并且与模型一的吻合程度很高.然后对分别加入了10%和20%高斯随机误差的自然电位进行反演试算，反演结果如图 4b和图 4c所示，电流等值线图依然能够较为准确地反应自然电场场源的极性，位置和深度稍有偏差.说明在偶极源模型中，反演算法具有良好的稳定性和抗噪声能力.
|
图 4 模型一反演电流等值线图 (a)理想数据反演结果; (b)10%高斯误差数据反演结果; (b)20%高斯误差数据反演结果.其中白色五角星标注了异常源的正确位置. Figure 4 Inversion results of model 1 (a) Synthetic data without error; (b) Synthetic data with 10% error; (c) Synthetic data with 20% error. The true source locations are depicted by white pentagram. |
对模型二进行反演试算，反演所得电流拟断面等值线图如图 5所示，在理想数据和10%误差数据的情况下，能准确地反应左侧负电流源的位置和深度，而在20%误差数据的情况下，左侧负电流源的位置和深度稍有偏差.而三种数据的反演都能反应出右侧正电流源排列的位置和趋势.说明在多源的情况下，依旧具有良好的反演效果和抗噪声能力.
|
图 5 模型二反演电流等值线图 (a)理想数据反演结果; (b)10%高斯误差数据反演结果; (b)20%高斯误差数据反演结果.其中白色五角星标注了异常源的正确位置. Figure 5 Inversion results of model 2 (a) Synthetic data without error; (b) Synthetic data with 10% error; (c) Synthetic data with 20% error. The true source locations are depicted by white pentagram. |
通过把目标区域的自然电位源进行统一处理，通过基于奇异边界矫正的有限单元法实现二维自然电位的高精度模拟.
5.2基于最小二乘正则化的自然电位二维反演能够有效地对自然电位异常源进行反演，能对异常源的位置进行准确地定位，反演算法收敛快速稳定，具有良好的抗噪声能力.在多源的情况下反演效果依旧很好，能反映出异常源分布的趋势.
5.3下一步的研究重点是通过实测数据检验算法的实用性.
致谢 感谢审稿专家和编辑部的大力支持.| [] | Abdelrahman E M, El-Araby T M, Essa K S .2009. Shape and depth determinations from second moving average residual self-potential anomalies[J]. Journal of Geophysics and Engineering, 6 (1) : 43–52. DOI:10.1088/1742-2132/6/1/005 |
| [] | Abdelrahman E S M, Essa K S, Abo-Ezz E R, et al .2008. New least-squares algorithm for model parameters estimation using self-potential anomalies[J]. Computers & Geosciences, 34 (11) : 1569–1576. |
| [] | Aizawa K .2008. Classification of self-potential anomalies on volcanoes and possible interpretations for their subsurface structure[J]. Journal of Volcanology and Geothermal Research, 175 (3) : 253–268. DOI:10.1016/j.jvolgeores.2008.03.011 |
| [] | Che-Alota V, Atekwana E A, Atekwana E A, et al .2009. Temporal geophysical signatures from contaminant-mass remediation[J]. Geophysics, 74 (4) : B113–B123. DOI:10.1190/1.3139769 |
| [] | El-Kaliouby H M, Al-Garni M A .2009. Inversion of self-potential anomalies caused by 2D inclined sheets using neural networks[J]. Journal of Geophysics and Engineering, 6 (1) : 29–34. DOI:10.1088/1742-2132/6/1/003 |
| [] | Fan Y Y, Du X B, Zlotnicki J, et al .2010. The electromagnetic phenomena before the MS 8.0 Wenchuan earth quake[J]. Chinese Journal of Geophysics, 53 (12) : 2887–2898. DOI:10.3969/j.issn.0001-5733.2010.12.012 |
| [] | Göktürkler G, Balkaya Ç .2012. Inversion of self-potential anomalies caused by simple-geometry bodies using global optimization algorithms[J]. Journal of Geophysics and Engineering, 9 (5) : 498–507. DOI:10.1088/1742-2132/9/5/498 |
| [] | Guo X J, Meng Q S, Wang J C, et al .2007. Leakage detection of underground pressureless industrial waste water pipeline using GPR、electrical and SP method[J]. Progress in Geophysics, 22 (1) : 279–282. DOI:10.3969/j.issn.1004-2903.2007.01.041 |
| [] | Jardani A, Revil A, Santos F, et al .2007. Detection of preferential infiltration pathways in sinkholes using joint inversion of self-potential and EM-34 conductivity data[J]. Geophysical Prospecting, 55 (5) : 749–760. DOI:10.1111/gpr.2007.55.issue-5 |
| [] | Jardani A, Revil A, Barrash W, et al .2009. Reconstruction of the water table from self-potential data:A bayesian approach[J]. Ground Water, 47 (2) : 213–227. DOI:10.1111/gwat.2009.47.issue-2 |
| [] | Li H, Fan Y R, Hu Y Y, et al .2013. Joint inversion of HDIL and SP with a five-parameter model for estimation of connate water resistivity[J]. Chinese Journal of Geophysics, 56 (2) : 688–695. DOI:10.6038/cjg20130233 |
| [] | Liu J W, Wang Z J, Du Z W .2009. Application of natural electric field in site groundwater exploration[J]. Chinese Journal of Engineering Geophysics, 6 (5) : 612–615. |
| [] | Martínez-Pagán P, Jardani A, Revil A, et al .2010. Self-potential monitoring of a salt plume[J]. Geophysics, 75 (4) : WA17–WA25. DOI:10.1190/1.3475533 |
| [] | Mauri G, Williams-Jones G, Saracco G .2010. Depth determinations of shallow hydrothermal systems by self-potential and multi-scale wavelet tomography[J]. Journal of Volcanology and Geothermal Research, 191 (3-4) : 233–244. DOI:10.1016/j.jvolgeores.2010.02.004 |
| [] | Monteiro Santos F A .2010. Inversion of self-potential of idealized bodies' anomalies using particle swarm optimization[J]. Computers & Geosciences, 36 (9) : 1185–1190. |
| [] | Pekşen E, Yas T, Kayman A Y, et al .2011. Application of particle swarm optimization on self-potential data[J]. Journal of Applied Geophysics, 75 (2) : 305–318. DOI:10.1016/j.jappgeo.2011.07.013 |
| [] | Placencia-Gómez E, Parviainen A, Hokkanen T, et al .2010. Integrated geophysical and geochemical study on AMD generation at the Haveri Au-Cu mine tailings, SW Finland[J]. Environmental Earth Sciences, 61 (7) : 1435–1447. DOI:10.1007/s12665-010-0459-9 |
| [] | Qin X Y, Xiao L Z, Zhang Y Z .2005. Methods of natural gas reservoir identification and evaluation of Erdos Basin[J]. Progress in Geophysics, 20 (4) : 1099–1107. DOI:10.3969/j.issn.1004-2903.2005.04.034 |
| [] | Rittgers J B, Revil A, Karaoulis M, et al .2013. Self-potential signals generated by the corrosion of buried metallic objects with application to contaminant plumes[J]. Geophysics, 78 (5) : EN65–EN82. DOI:10.1190/geo2013-0033.1 |
| [] | Sailhac P, Marquis G .2001. Analytic potentials for the forward and inverse modeling of SP anomalies caused by subsurface fluid flow[J]. Geophysical Research Letters, 28 (9) : 1851–1854. DOI:10.1029/2000GL012457 |
| [] | Shang H, Guo X J, Wang X L, et al .2010. Design and experiment of the 3-D electrical resistivity detecting system for monitoring underground soil contaminated by landfill leachate[J]. Progress in Geophysics, 25 (4) : 1485–1491. DOI:10.3969/j.issn.1004-2903.2010.04.042 |
| [] | Skianis G A, Papadopulos T D, Vaiopoulos D A .2006. Direct interpretation of self-potential anomalies produced by a vertical dipole[J]. Journal of Applied Geophysics, 58 (2) : 130–143. DOI:10.1016/j.jappgeo.2005.05.003 |
| [] | Song Z Q, Li W F, Tang C J, et al .2009. Dividing sedimentary facies and reservoir distributions by using natural potential and natural gamma ray logging curves[J]. Progress in Geophysics, 24 (2) : 651–656. DOI:10.3969/j.issn.1004-2903.2009.02.038 |
| [] | Tan D C, Zhao J L, Liu X F, et al .2014. Features of regional variations of the spontaneous field[J]. Chinese Journal of Geophysics, 57 (5) : 1588–1598. DOI:10.6038/cjg20140522 |
| [] | Tan D C, Zhao J L, Xi J L, et al .2012. The variation of waveform and analysis of composition for the geoelectrical field before moderate or strong earthquakes in Qinghai-Tibetan plateau regions[J]. Chinese Journal of Geophysics, 55 (3) : 875–885. DOI:10.6038/j.issn.0001-5733.2012.03.017 |
| [] | Yang C M, Li H Q, Lu D W, et al .2005. Electrokinetic phenomena in oilfield development[J]. Progress in Geophysics, 20 (4) : 1140–1144. DOI:10.3969/j.issn.1004-2903.2005.04.038 |
| [] | Zhang G B, Aubert M, Qi J L .1999. A study of self-potential method in hydrogeological exploration of volcanic area[J]. Chinese Journal of Geophysics, 42 (S) : 250–256. |
| [] | Zhang P, Liu S D, Xu J P, et al .2013. Experimental on geoelectric field of water seepage of complex rock mass[J]. Progress in Geophysics, 28 (3) : 1557–1562. DOI:10.6038/pg20130352 |
| [] | Zhu X X, Cui Y A, Li X Y, et al .2015. Inversion of self-potential anomalies based on particle swarm optimization[J]. Journal of Central South University (Science and Technology), 46 (2) : 579–585. |
| [] | 范莹莹, 杜学彬, ZlotnickiJ, 等.2010. 汶川MS 8.0大震前的电磁现象[J]. 地球物理学报, 53 (12) : 2887–2898. DOI:10.3969/j.issn.0001-5733.2010.12.012 |
| [] | 郭秀军, 孟庆生, 王基成, 等.2007. 地球物理方法在含油工业污水管道渗漏探测中的应用[J]. 地球物理学进展, 22 (1) : 279–282. DOI:10.3969/j.issn.1004-2903.2007.01.041 |
| [] | 李虎, 范宜仁, 胡云云, 等.2013. 基于阵列感应与自然电位联合反演地层水电阻率[J]. 地球物理学报, 56 (2) : 688–695. DOI:10.6038/cjg20130233 |
| [] | 刘加文, 王治军, 杜志, 伟.2009. 自然电场法在场地地下水勘查中的应用[J]. 工程地球物理学报, 6 (5) : 612–615. |
| [] | 秦绪英, 肖立志, 张元中.2005. 鄂尔多斯盆地天然气有效储层识别与评价方法[J]. 地球物理学进展, 20 (4) : 1099–1107. DOI:10.3969/j.issn.1004-2903.2005.04.034 |
| [] | 尚浩, 郭秀军, 王仙丽, 等.2010. 垃圾填埋场渗漏污染三维分布式电学监测系统研制及应用[J]. 地球物理学进展, 25 (4) : 1485–1491. DOI:10.3969/j.issn.1004-2903.2010.04.042 |
| [] | 宋子齐, 李伟峰, 唐长久, 等.2009. 利用自然电位与自然伽马测井曲线划分沉积相带及储层分布[J]. 地球物理学进展, 24 (2) : 651–656. DOI:10.3969/j.issn.1004-2903.2009.02.038 |
| [] | 谭大诚, 赵家骝, 刘小凤, 等.2014. 自然电场的区域性变化特征[J]. 地球物理学报, 57 (5) : 1588–1598. DOI:10.6038/cjg20140522 |
| [] | 谭大诚, 赵家骝, 席继楼, 等.2012. 青藏高原中强地震前的地电场变异及构成解析[J]. 地球物理学报, 55 (3) : 875–885. DOI:10.6038/j.issn.0001-5733.2012.03.017 |
| [] | 杨春梅, 李洪奇, 陆大卫, 等.2005. 油田开发过程中的动电现象研究[J]. 地球物理学进展, 20 (4) : 1140–1144. DOI:10.3969/j.issn.1004-2903.2005.04.038 |
| [] | 张贵宾, AubertM, 齐剑玲.1999. 自然电场法勘探火山岩地区地下水的研究[J]. 地球物理学报, 42 (S) : 250–256. |
| [] | 张朋, 刘盛东, 许进鹏, 等.2013. 复杂岩体水渗流地电场响应试验分析[J]. 地球物理学进展, 28 (3) : 1557–1562. DOI:10.6038/pg20130352 |
| [] | 朱肖雄, 崔益安, 李溪阳, 等.2015. 基于粒子群优化的自然电场数据反演[J]. 中南大学学报(自然科学版), 46 (2) : 579–585. |