朱肖雄，崔益安，李溪陽，佟鐵鋼，紀銅鑫

?

基于粒子群優(yōu)化的自然電場數(shù)據(jù)反演

朱肖雄1, 2，崔益安1, 2，李溪陽1, 2，佟鐵鋼1, 2，紀銅鑫1, 2

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

在分析測試粒子數(shù)、速度因子、目標函數(shù)等算法參數(shù)對粒子群優(yōu)化算法效果的影響規(guī)律的基礎(chǔ)上，設(shè)計自然電場粒子群優(yōu)化反演算法，并對加入不同程度白噪聲模擬數(shù)據(jù)進行反演試算。研究結(jié)果表明：設(shè)計的粒子群優(yōu)化算法能有效實現(xiàn)對自然電場數(shù)據(jù)的反演，算法具有收斂速度快、穩(wěn)定、反演精度較高和抗噪音能力強等優(yōu)點，可以較為準確地得到異常體的位置、形態(tài)、極化角等參數(shù)，能較好地滿足生成實際要求。

自然電場；反演；粒子群優(yōu)化；參數(shù)分析

自然電場法是一種基于巖、礦石天然產(chǎn)生電位差的觀測方法。自然電場法除了應(yīng)用在傳統(tǒng)的金屬硫化物和石墨[1]等礦產(chǎn)勘查領(lǐng)域外，近年來在環(huán)境和工程地球物理領(lǐng)域的應(yīng)用越來越廣泛，主要應(yīng)用于堤壩的滲漏檢測[2?3]、地下水調(diào)查[4]和孔洞檢測[5]等。在這些應(yīng)用中，經(jīng)常需要圈定如礦脈、滲漏點、孔洞這類異常體。這些異常體通?？梢越瞥蓡蝹€球體、圓柱體和板狀體等簡單幾何模型或多個簡單幾何模型的組合。一些圖形和數(shù)值的方法已經(jīng)被用來解釋自然電位異常，如特征點法[6]、曲線擬合法[7]、最小二乘法[8]、梯度法[9]、傅里葉分析法[10]、校正平均殘余異常法 等[11]。特征點法和曲線擬合法比較繁瑣并且反演效果不佳。后幾種方法當測量數(shù)據(jù)有干擾時，反演結(jié)果會受到較大的影響，因此，很有必要采用抗干擾能力強的方法來反演自然電場數(shù)據(jù)。粒子群優(yōu)化算法(PSO)是Kennedy等于1995年提出的一種新穎的進化算 法[12]，它是一種基于群智能的算法，通過模擬鳥集群飛行覓食的行為來達到最優(yōu)化的目的。該算法已經(jīng)成功地應(yīng)用于很多領(lǐng)域，例如信號處理[13]、圖像處理[14]以及巖土力學(xué)等[15]。近年來，粒子群優(yōu)化算法在地球物理領(lǐng)域的應(yīng)用研究[16?18]也越來越多。例如，崔益安等[19?20]將粒子群優(yōu)化算法用于激電數(shù)據(jù)的聯(lián)合反演和中梯電阻率多異常體反演的研究。Fernando等[21?22]開展了關(guān)于自然電位異常數(shù)據(jù)反演的研究，并證明了粒子群優(yōu)化算法是一種有效的反演自然電位異常數(shù)據(jù)的方法，但其只對10%噪音模擬數(shù)據(jù)進行了測試?？紤]到自然電場數(shù)據(jù)易受干擾而導(dǎo)致實際噪音干擾遠大于10%，為了實現(xiàn)抗噪音能力更強的反演算法，本文作者研究了粒子群優(yōu)化算法參數(shù)對算法效果的影響規(guī)律，以期通過更佳的參數(shù)設(shè)置來提高算法的性能和抗噪音能力。

1 正演計算

地下極化球體模型如圖1所示。假設(shè)在均勻充滿電阻率為1的介質(zhì)里，有1個電阻率為2、半徑為0的球體。當球體被均勻極化時，球體表面形成不均勻的(異常)雙電層，其電位差Δ(近似看作是偶電層的電位躍變值)隨極化方向的坐標呈線性變化。均勻極化球體內(nèi)、外電位的分布以極化軸為對稱軸，與方位角無關(guān)，故滿足球坐標系的拉普拉斯方程：

圖1 地下極化球體模型

有邊界條件為

式中：1和2分別為球體的外部和內(nèi)部電位；1n和2n為球體內(nèi)外法向的電流密度。

地下異常源模型如圖2所示。通過分離變量法解得空間內(nèi)任意一點的自然電位異常為：

(a) 球體模型；(b) 水平圓柱體模型；(c) 垂直圓柱體模型

當極化軸偏轉(zhuǎn)與負半軸呈時，如圖2(a)所示，可推出地表的自然電位異常公式為：

類似地，可推出地下水平圓柱體和垂直圓柱體的自然電位異常公式分別為：

綜合式(5)~(8)可得類球體的地表自然電位異常公式：

其中：為埋藏源的深度；為極化角度；為形狀因子，0為異常源的坐標，為電偶極矩。形狀因子取值為1.5，1.0，0.5時分別表示球體、水平圓柱體和垂直圓柱體。在異常源地表處(=0)的異常值為

地下異板狀模型如圖3所示。設(shè)在均勻充滿電阻率為的介質(zhì)里，有一傾斜板狀體?？臻g中電位分布為：

其中：1和2分別為傾斜板狀體兩端到空間某一點的距離；為單位長度上的電流；′為電偶極矩。

(a) 模型示意圖；(b) 模型坐標圖

可推出地表上的自然電位異常式為

其中：為傾斜板狀體的傾斜角度；為傾斜板狀體長度的1/2；0為異常源的坐標；為傾斜板狀體中心的深度。

2 粒子群優(yōu)化算法的設(shè)計

粒子群優(yōu)化算法的基本思想是隨機初始化一群粒子，將每個粒子視為優(yōu)化問題的1個可行解，粒子的好壞由一個事先設(shè)定的目標函數(shù)來確定。每個粒子將在可行解空間中運動，并由1個速度變量確定其方向和距離。通常粒子將追尋當前的最優(yōu)粒子，經(jīng)過多次迭代得到最優(yōu)解。

自然電位異常的粒子群優(yōu)化反演算法具體的設(shè)計流程如下。

1) 設(shè)定搜索空間，即反演參數(shù)的上、下界。設(shè)置迭代過程中的最大迭代次數(shù)max，粒子更新速度范圍為[min，max]，

式中：為比例系數(shù)；||為搜索范圍的上界。

2) 初始化粒子群位置和更新速度。標記初始粒子群位置為粒子群最優(yōu)位置Best，通過目標函數(shù)算評價得到粒子群中的最優(yōu)粒子。

3) 粒子的狀態(tài)更新，通過式(15)和式(16)更新每一個粒子的速度和位置。

式中：W為慣性權(quán)重；start取0.9，end取0.4；，為速度因子，根據(jù)速度因子策略更新；Best為當前粒子群的最優(yōu)位置；為迭代次數(shù)；1和2為在(0,1)之間的隨機數(shù)。

4) 對粒子進行評價，通過目標函數(shù)評價粒子群X+1，更新粒子群最優(yōu)位置和最優(yōu)粒子。

5) 迭代至預(yù)設(shè)的最大迭代次數(shù)max。輸出反演 結(jié)果。

根據(jù)反演結(jié)果，通過單一變量法對粒子數(shù)、速度比例系數(shù)、速度因子策略、和目標函數(shù)進行分析比對，如圖4所示。

(a) 粒子數(shù)分析；(b) 速度比例因子分析；(c) 學(xué)習(xí)因子策略分析；(d) 目標函數(shù)分析

通過測試，在迭代次數(shù)max取200，粒子數(shù)取50，速度比例系數(shù)取0.1，線性遞減速度因子策略，目標函數(shù)選擇式(18)時，基于粒子群優(yōu)化的自然電場反演的效果較好。

3 反演示例

首先對單個異常體模型的正演計算數(shù)據(jù)進行反演測試。為了便于比較分析反演效果，采用了與文獻[22]相同的地電模型參數(shù)，具體參數(shù)為：垂直圓柱體(=4 m，=0.5，=40°，=?500 mV，0=?321.39 mV，0=0 m)和傾斜板狀體模型(=10 m，=5 m，=35°，0=8 m，′=150 mV)。為了測試抗噪音能力，對模擬數(shù)據(jù)加入10%~30%的隨機噪音再進行反演。其中反演結(jié)果誤差的評價表達式為

分別在理想情況下和加入噪音的狀態(tài)下對單個異常體的模擬數(shù)據(jù)進行反演，結(jié)果如圖5和圖6所示。由圖5(a)和圖6(a)可見：粒子群優(yōu)化反演曲線與正演模擬的曲線完全相同。說明了粒子群優(yōu)化反演和有很高的反演精度。在加入了10%~30%的噪音情況下，具體數(shù)值如表1和表2所示。粒子群優(yōu)化反演仍然獲得了很好的反演效果，說明該算法有較強的抗噪音能力。

(a) 理想模擬數(shù)據(jù)反演結(jié)果；(b) 帶噪音的模擬數(shù)據(jù)反演結(jié)果

(a) 理想模擬數(shù)據(jù)反演結(jié)果；(b) 帶噪音的模擬數(shù)據(jù)反演結(jié)果

表1 垂直圓柱體理想數(shù)據(jù)及噪音數(shù)據(jù)的PSO反演結(jié)果

表2 傾斜板狀體理想數(shù)據(jù)及噪音數(shù)據(jù)的PSO反演結(jié)果

同時，用傳統(tǒng)的高斯算法對模擬數(shù)據(jù)進行反演，反演結(jié)果如表3和表4所示。由反演結(jié)果可知：無噪音和噪音較小時，高斯反演的參數(shù)精度很高；隨著噪音加大，高斯法反演的參數(shù)精度要比粒子群反演的精度低，且高斯反演嚴重依賴初始值設(shè)置，在目標函數(shù)具有多極值的情況下，很難設(shè)置合適的初值，而若設(shè)置稍有不當，則會發(fā)散，反演不出結(jié)果。

表3 垂直圓柱體理想數(shù)據(jù)及噪音數(shù)據(jù)的高斯反演結(jié)果

表4 傾斜板狀體理想數(shù)據(jù)及噪音數(shù)據(jù)的高斯反演結(jié)果

對多個異常模型的正演模擬數(shù)據(jù)進行反演測試。測試的模型由1個球體模型(=2.6 m，=1.5，=40°，=?1 615 mV，0=?75 mV，0=?55 m)和1個傾斜板狀體模型(=5 m，=3 m，=30°，0=6 m，′=102 mV)組成。并對加入了10%~30%的隨機噪音的模擬數(shù)據(jù)進行測試。在理想情況下和帶噪應(yīng)條件下對多個異常體的模擬數(shù)據(jù)反演結(jié)果，對比如圖8所示。由圖8可見：反演曲線與正演曲線完全相同。說明反演算法具有同時反演多個目標體的能力，且反演精度很高。在加入10%~30%噪音的情況下對多個異常體的模擬數(shù)據(jù)進行反演，粒子群優(yōu)化反演具有很強的抗噪音能力。反演效果完全能滿足對異常體的位置、形態(tài)、極化角、埋深參數(shù)準確估計的要求。

(a) 理想模擬數(shù)據(jù)反演結(jié)果；(b) 帶噪音的模擬數(shù)據(jù)反演結(jié)果

4 結(jié)論

1) 粒子群優(yōu)化反演對單個和多個異常體的反演效果都很好，反演精度高，收斂速度快，抗噪音干擾能力強。且由于正演計算很快，在迭代200次的情況下，耗時可忽略不計，反應(yīng)效率很高。相比于傳統(tǒng)的線性算法，粒子群優(yōu)化算法的優(yōu)點在于不用設(shè)置初值，前期搜索全局，后期搜索局部，不易發(fā)散。這種快速的高精度反演具有重要實際意義。

2) 下一步的研究重點并通過對典型實測數(shù)據(jù)進行反演，檢驗算法實用性。

[1] 李林, 解立發(fā). 自然電場法在晶質(zhì)石墨資源勘查中的應(yīng)用[J]. 中國非金屬礦工業(yè)導(dǎo)刊, 2006, 52: 57?60.
LI Lin, XIE Lifa. Self-potential method in the application of graphite exploration[J]. China Non-metallic Mining Industry Herald, 2006, 52: 57?60.

[2] 王紹彪, 湯浩. 綜合物探方法在探測基坑圍堰滲漏中的應(yīng)用[J]. 人民珠江, 2011, 32(Z1): 52?53.
WANG Shaobiao, TANG Hao. Comprehensive Geophysical method in the application of detecting leaks on foundation pits and cofferdam[J]. Pearl River, 2011, 32(Z1): 52?53.

[3] 孫士安. 伏漫水庫主壩滲漏通道勘察方法[J]. 廣西水利水電, 2012(1): 12?15.
SUN Shian, Reconnaissance method of detecting leaks in passages on main dam of Fuman reservoir[J]. Guangxi Water Resources & Hydropower Engineering, 2012(1): 12?15.

[4] 劉加文, 王治軍, 杜志偉. 自然電場法在場地地下水勘查中的應(yīng)用[J]. 工程地球物理學(xué)報, 2009, 6(5): 612?615.
LIU Jiawen, WANG Zhijun, DU Zhiwei. Application of natural electric field in site groundwater exploration[J]. Chinese Journal of Engineering Geophysics, 2009, 6(5): 612?615.

[5] Jardani A, Revil A, Santos F A M, et al. Detection of preferential infiltration pathways in sinkholes using joint inversion of self-potential and EM-34 conductivity data[J]. Geophysical Prospecting, 2007, 55(5): 749?760.

[6] Sharma S P, Biswas A. Interpretation of self-potential anomaly over a 2D inclined structure using very fast simulated-annealing global optimization - An insight about ambiguity[J]. Geophysics, 2013, 78(3): WB3?WB15.

[7] El-Araby H M. A new method for complete quantitative interpretation of self-potential anomalies[J]. Journal of Applied Geophysics, 2004, 55(3/4): 211?224

[8] Abdelrahman E M, Essa K S, Abo-Ezz E R, et al. New least-square algorithm for model parameters estimation using self- potential anomalies[J]. Computers & Geosciences, 2008, 34: 1569?1576.

[9] Abdelrahman E M, El-Araby H M, El-Araby T M, et al. New methods for shape and depth determinations from SP data[J]. Geophysics, 2003, 68(4): 1202?1210.

[10] Fedi M, Abbas M A. A fast interpretation of self-potential data using the depth from extreme points method[J]. Geophysics, 2013, 78(2): E107?E116.

[11] Abdelrahman E M, El-Araby T M, Essa K S. Shape and depth determinations from second moving average residual self-potential anomalies[J]. Journal of Geophysics and Engineering, 2009(6): 43?52.

[12] JIE Jing, ZENG Jianchao, HAN Chongzhao, et al. Knowledge- based cooperative particle swarm optimization[J]. Applied Mathematics and Computation, 2008, 205(2): 861–873

[13] Jun S, Wei F, Wenbo X. A quantum-behaved particle swarm optimization with diversity- guided mutation for the design of two-dimensional IIR digital filters[J]. IEEE Transactions on Circuits and Systems II: Express Briefs, 2010, 57(2): 141?145.

[14] Hao G. Multilevel thresholding for image segmentation through an improved quantum- behaved particle swarm algorithm[J]. IEEE Transactions on Instrumentation and Measurement, 2010, 59(4): 934?946.

[15] 楊文東, 張強勇, 李術(shù)才, 等. 粒子群算法在時效變形參數(shù)反演中的應(yīng)用[J]. 中南大學(xué)學(xué)報(自然科學(xué)版), 2013, 44(1): 282?288.
YANG Wendong, ZHANG Qiangyong, LI Shucai, et al. Application of particle swarm optimization in time-dependent parameters[J]. Journal of Central South University (Science and Technology), 2013, 44(1): 282?288.

[16] SONG Xianhai, TANG Li, Lü Xiaochun, et al. Application of particle swarm optimization to interpret Rayleigh wave dispersion curves[J]. Journal of Applied Geophysics, 2012, 84: 1?13.

[17] 師學(xué)明, 肖敏, 范建柯, 等. 大地電磁阻尼粒子群優(yōu)化反演法研究[J]. 地球物理學(xué)報, 2009, 52(4): 1114?1120.
SHI Xueming, XIAO Min, FAN Jianke, et al. The damped PSO algorithm and its application for magnetotelluric sounding data inversion[J]. Chinese Journal of Geophysics, 2009, 52(4): 1114?1120.

[18] 蔡涵鵬, 賀振華, 高剛, 等. 基于混合優(yōu)化算法的地震數(shù)據(jù)匹配追蹤分解[J]. 中南大學(xué)學(xué)報(自然科學(xué)版), 2013, 44(2): 687?694.
CAI Hanpeng, HE Zhenghua, GAO Gang. Seismic data matching pursuit using hybrid optimization algorithm and its application[J]. Journal of Central South University (Science and Technology), 2013, 44(2): 687?694.

[19] 崔益安, 李溪陽, 向恩明, 等. 基于粒子群優(yōu)化的雙頻激電數(shù)據(jù)聯(lián)合反演[J]. 中國有色金屬學(xué)報, 2013, 23(9): 2498?2505.
CUI Yian, LI Xiyang, XIANG Enming, et al. Joint inversion of dual frequency IP data using PSO[J], The Chinese Journal of Nonferrous Metals, 2013, 23(9): 2498?2505.

[20] 崔益安, 紀銅鑫, 李溪陽, 等. 基于粒子群優(yōu)化的多目標體中梯電阻率異常反演[J]. 地球物理學(xué)進展, 2013, 28(4): 2164?2170.
CUI Yian, Ji Tongxin, Li Xiyang, et al. Inversion of multi-anomalies in resistivity profiling based on particle swarm optimization[J]. Progress in Geophysics, 2013, 28(4): 2164?2170.

[21] Fernando A, Monteiro Santos. Inversion of self-potential of idealized bodies’ anomalies using particle swarm optimization[J]. Computers & Geosciences, 2010, 36: 1185?1190.

[22] Peksen E, Türker Y,A Yekta K, et al. Application of particle swarm optimization on self-potential data[J]. Journal of Applied Geophysics, 2011, 75(2): 305?318.

Inversion of self-potential anomalies based on particle swarm optimization

ZHU Xiaoxiong1, 2, CUI Yian1, 2, LI Xiyang1, 2, TONG Tiegang1, 2, JI Tongxin1, 2

(1. School of Geosciences and Info-Physics, Central South University, Changsha 410083, China;2. Hunan Key Laboratory of Non-ferrous Resources and Geological Hazard Detection,Central South University, Changsha 410083, China)

Based on testing and analyzing relevant parameters including particle quantity, rate scale factor, objective function etc., the particle swarm optimization (PSO) was used to design inversion algorithm for self-potential data. Through adding different degrees of Gauss noise, the synthetic data was used to test the designed inversion algorithm. The results show that the PSO algorithm can effectively realize the inversion of self-potential data with fast and stable convergence, high inversion accuracy and high anti-noise capability. Through the designed algorithm, the parameters contain origin of the anomaly, shape, polarization angle, etc, can relativey accurately be obtained and can meet the demands of engineering investigation and mineral exploration.

self-potential; inversion; particle swarm optimization; parameter analysis

P631

A

1672?7207(2015)02?0579?07

2014?03?06；

2014?06?20

國家自然科學(xué)基金資助項目(41274122，41374119)；教育部博士點基金資助項目(20110162130008)；國家科技支撐計劃項目(2011BAB04B08)；國家科技基礎(chǔ)性工作專項(2013FY110800)(Projects (41274122, 41374119) supported by the National Natural Science Foundation of China; Project (20110162130008) supported by the Doctoral Fund of Ministry of Education of China; Project (2011BAB04B08) supported by the Special Basic Scientific Project of China; Project (2013FY110800) supported by the National Special Basic Scientific Program of China)

崔益安，博士，副教授，從事物探方法與數(shù)據(jù)處理研究；E-mail：cuiyian@csu.edu.cn

10.11817/j.issn.1672-7207.2015.02.028

(編輯 趙俊)