基于FCM聚類算法的二維RMT電阻率與介電常數(shù)聯(lián)合反演

易柯， 張志勇,2*， 李曼， 原源， 郭一豪， 周峰,3

1 東華理工大學地球物理與測控技術學院， 南昌 330013 2 東華理工大學核資源與環(huán)境國家重點實驗室， 南昌 330013 3 中南大學有色金屬成礦預測與地質環(huán)境監(jiān)測教育部重點實驗室， 長沙 410083

0 引言

射頻大地電磁(Radio-magnetotelluric, RMT)通過接收正交的高頻(10kHz～1MHz)電場與磁場分量進行勘探，相較于(音頻)大地電磁((A)MT)測深，RMT對百米內的地下介質具有更好的分辨率.該方法在歐洲已應用于地下水(Turberg et al., 1994; Pedersen et al., 2005; Ismail et al., 2011)、環(huán)境工程(Tezkan et al., 2000, 2005)、地質構造(Linde and Pedersen, 2004a,b)、地質災害(Shan et al., 2014, 2016; Wang et al., 2016)等勘探，但在我國RMT方面的研究相對較少(原源等, 2015, 2019).

RMT初期采用與甚低頻電阻率(VLF-R)觀測方法相似的標量形式(Benson et al., 1997; Turberg et al., 1994)；隨后發(fā)展了10 kHz～250 kHz頻段的張量Enviro-MT儀器(Bastani, 2001)；當前儀器頻段得到進一步拓展，研發(fā)了10 kHz～1 MHz的張量RMT-F系統(tǒng)(Tezkan and Saraev, 2008)；得益于高頻電磁勘探觀測時間短的優(yōu)點，發(fā)展了連續(xù)觀測的RMT系統(tǒng)，方便了水域勘探(Bastani et al., 2015).此外，研發(fā)的人工場源射頻大地電磁(CSRMT)技術有效提高了觀測數(shù)據(jù)的信噪比(Pedersen et al., 2005; Bastani et al., 2009; Ismail et al., 2011; Saraev et al., 2017).

隨RMT觀測頻段不斷提高，位移電流的影響不容忽視.理論計算表明(圖1)，對于相對介電常數(shù)εr為5的均勻介質，電阻率ρ大于100 Ωm且頻率超過100 kHz時，位移電流占比超過1%，電磁波的趨膚深度在高頻、高阻時受頻率影響變?。浑S著相對介電常數(shù)的增加，位移電流在相對更低頻和低阻處即產(chǎn)生較大的影響.高阻介質(ρ≥10000 Ωm)考慮位移電流的二維RMT正演研究表明，受介電常數(shù)(εr=5)的影響，視電阻率和相位值從10 kHz開始降低(Kalscheuer et al., 2008)；位移電流受地形影響嚴重，采用非結構化網(wǎng)格進行帶地形正演可有效提高計算精度(原源等, 2015).試算表明對RMT數(shù)據(jù)同時反演電阻率和介電常數(shù)是可行的，但缺少聯(lián)合約束時，待解參數(shù)的增加將直接影響反演結果的精確性(原源等, 2019).此外，RMT與直流聯(lián)合反演電阻率，提高了資料解釋精度(Seher and Tezkan, 2007; Bastani et al., 2012; Yogeshwar et al., 2012; Wang et al., 2018).

圖1 趨膚深度和位移電流與傳導電流之比左側為趨膚深度，右側為位移電流與傳導電流之比，(a)、(c)，(b)、(d)對應的相對介電常數(shù)εr分別為5、81.Fig.1 Skin depth and the ratio of displacement current to conduction currentLeft two figures represent the skin depth, right two figures represent the ratio of displacement current to conduction current, the relative permittivity of (a) and (c), (b) and (d) are 5 and 81 respectively.

物性相關和結構耦合是構建聯(lián)合約束的兩種主要方法.當不同物性參數(shù)存在巖石物理學的耦合關系時可采用物性相關約束(Lees and VanDecar, 1991; Nielsen and Jacobson, 2000)，但該耦合關系通常只適用于局部地區(qū)(陳曉等, 2016).采取交叉梯度約束(Gallardo and Meju, 2003, 2004)無需得到參數(shù)間的顯式解析關系；然而，當物性參數(shù)變化無明顯方向性時，交叉梯度的約束效果有限(王俊等, 2013).近年來，模糊C均值(Fuzzy C-means, FCM)聚類算法開始應用于聯(lián)合反演，并取得了良好效果(Paasche and Tronicke, 2007; Carter-Mcauslan et al., 2015; Sun and Li, 2015, 2016).基于正則化理論，F(xiàn)CM聚類算法通過聚類分析得到的隸屬度可直接用于反演迭代(Paasche and Tronicke, 2007)；將巖石地球物理先驗信息以參考聚類中心的方式融入反演過程(Sun and Li, 2015)，可有效提高反演效果.此外，F(xiàn)CM聚類算法在磁法(Li and Sun, 2016)、重力(Maag and Li, 2018)等數(shù)據(jù)反演中的應用均取得了良好效果.

目前，尚無基于聯(lián)合約束的RMT電阻率和介電常數(shù)反演研究，為提高聯(lián)合反演效果，本文在經(jīng)典最小結構模型正則化約束的基礎上，采用FCM聚類算法開展了二維RMT電阻率與介電常數(shù)的聯(lián)合反演研究.研究工作首先驗證了二維有限單元正演計算的精度；隨后開展了介電常數(shù)、電阻率單獨反演，進一步驗證反演算法各環(huán)節(jié)的正確性；最后模擬實際工作情況設計了帶地形模型，進一步討論了FCM聚類聯(lián)合反演的效果以及參數(shù)選擇對反演結果的影響.

1 全電流二維RMT有限元正演

取時諧因子eiω t時的頻率域麥克斯韋方程為

(1)

(2)

二維地電模型見圖2，y為地質體走向，z方向垂直向下，磁導率μ取真空磁導率μ0，則電磁場分量與地電參數(shù)σ、ε只沿x、z方向變化.

圖2 二維地球物理模型Fig.2 Sketch of two dimensional geophysical model

二維條件下，全電流Helmholtz方程為

TE模式：

(3)

TM模式：

(4)

式(3)、(4)可統(tǒng)一表示為偏微分方程：

(5)

綜上所述，二維RMT邊界條件為

(6)

式(5)、(6)與求解下列變分問題等價：

(7)

對區(qū)域進行非結構化三角單元離散后，采用有限單元法對(7)式變分問題進行求解，得以下線性系統(tǒng)：

Bu=0,

(8)

其中B為剛度矩陣，u為節(jié)點場值，求解(8)式線性方程，得到各節(jié)點主場，并通過主場梯度計算輔助場(徐世浙, 1994)，最后可計算視電阻率和相位響應：

(9)

(10)

Zre、Zim分別為阻抗Z的實部與虛部.

為驗證正演計算精度，采用非結構化網(wǎng)格有限單元法對100 Hz～1 MHz頻段進行了數(shù)值計算，并與一維解析解進行對比分析.均勻半空間和三層水平層狀模型對比結果見圖3，其中圖3a、3b分別為均勻半空間(ρ=10000 Ωm,εr=5)視電阻率和相位響應，一維解析解與二維有限元解曲線吻合，視電阻率平均相對誤差為0.4%，相位平均相對誤差為0.1%；當頻率超過30 kHz后，由于位移電流的影響視電阻率大幅度下降，相位則從10 kHz開始明顯下降.圖3c、3d分別為三層層狀模型視電阻率和相位響應，視電阻率平均相對誤差為2.2%，相位平均相對誤差為1.1%.綜上，正演對比試驗驗證了開發(fā)的二維程序在100 Hz～1 MHz頻段計算的正確性，計算精度滿足要求.

2 RMT反演理論

2.1 目標函數(shù)

根據(jù)正則化理論(Tong et al., 2018)構建RMT單獨反演介電常數(shù)和電阻率的目標函數(shù)分別為

P(mε)=φd(mε)+λεφm(mε)

(11)

P(mρ)=φd(mρ)+λρφm(mρ)

(12)

其中，φd為數(shù)據(jù)擬合泛函，φm為模型穩(wěn)定泛函，λε、λρ為模型穩(wěn)定泛函的正則化因子，Wd為數(shù)據(jù)協(xié)方差矩陣，可表示為

其中，s為極小正數(shù)，nd為觀測數(shù)據(jù)總個數(shù)，i為第i個觀測數(shù)據(jù)的噪音水平.Wε、Wρ為最小結構模型加權矩陣，dobs為觀測數(shù)據(jù)，A(mε,mρ)為正演響應，mε、mρ為待解參數(shù)，為參考模型.本文研究中，采用最小二乘法計算非結構化網(wǎng)格模型粗糙度(Lelièvre and Farquharson, 2013)，取參考模型與模型粗糙度之比為1‰(αε 0=αρ 0=0.001)，模型粗糙度可表示為

(13)

其中，αx、αz為比例系數(shù)，取αx=αz=1.

圖3 100 Hz～1 MHz正演二維有限元解與一維解析解對比(a)、(b) 分別為均勻半空間(ρ=10000 Ωm, εr=5)視電阻率、相位響應; (c)、(d) 分別為三層水平層狀模型視電阻率、相位響應.Fig.3 Comparison of 2D FEM solution with 1D analytic solution from 100 Hz to 1 MHz(a) and (b) show the calculated apparent resistivity and phase of homogenous model with ρ=10000 Ωm and εr=5 respectively, (c) and (d) show the calculated apparent resistivity and phase of a three-layered model, respectively.

雙參數(shù)聯(lián)合反演目標函數(shù)為

(14)

(15)

其中C為聚類中心個數(shù)，M為模型單元總個數(shù)；mj=(mεj,mρj)T為第j個模型單元物性值，其中mεj、mρj分別為第j個單元的介電常數(shù)與電阻率值；vk=(vεk,vρk)T為反演的第k個聚類中心，其中vεk、vρk分別為第k個介電常數(shù)與電阻率的聚類中心；ujk為第j個模型單元物性值對第k個聚類中心的隸屬度，q為模糊化參數(shù)，本文取q=2.0；tk=(tεk,tρk)T為第k個參考聚類中心，tεk、tρk分別為第k個介電常數(shù)與電阻率的參考聚類中心；ηk為第k個參考聚類中心tk的權重因子，表示第k個參考聚類中心的置信度.

vk、ujk可表示為(Sun and Li, 2016)

(16)

(17)

聚類項φFCM中第一項可表示為

(18)

結合式(14)、(15)、(18)可將雙參數(shù)聯(lián)合反演目標函數(shù)表示為

(19)

2.2 目標函數(shù)優(yōu)化求解

對于雙參數(shù)mε、mρ的反演，采用迭代的高斯-牛頓算法(Zaslavsky et al., 2013)優(yōu)化求解目標函數(shù)(19)的最小值.聚類項中vk和ujk均采用迭代法求解(Sun and Li, 2016)，其第n+1次迭代形式為

(20)

(21)

目標函數(shù)(19)的第n+1次迭代形式為

(22)

式(22)中J=(Jε,Jρ)T為靈敏度矩陣，Jε、Jρ分別為介電常數(shù)與電阻率的靈敏度矩陣，為減少計算工作量，采用互換定理(Rodi, 1976)進行計算；Δm=(Δmε,Δmρ)T為模型改變量，Δmε、Δmρ分別為介電常數(shù)與電阻率的模型改變量.

對目標函數(shù)(19)的第n+1次迭代形式求偏導并令其等于0：

(23)

可得高斯-牛頓迭代方程：

(24)

采用雙共軛梯度穩(wěn)定算法(BICGSTAB)求解式(24)可得新的模型參數(shù)：

m(n+1)=m(n)+γTΔm(n),

(25)

其中γ=(γε,γρ)T為沿改進量Δm的搜索步長(Zhdanov, 2002)，γε、γρ分別為介電常數(shù)與電阻率的搜索步長.

3 模型試算

3.1 電阻率與介電常數(shù)反演

為驗證算法的正確性，首先進行了相對介電常數(shù)異常和電阻率異常單獨反演模型試驗.設計模型如下：在相對介電常數(shù)為5、電阻率為3000 Ωm的均勻介質中，設置長寬分別為400 m、50 m的矩形異常體，異常體上頂埋深20 m.當進行介電常數(shù)反演時，取矩形異常體電阻率為3000 Ωm，相對介電常數(shù)為80，即只有介電常數(shù)異常(圖4a)；當進行電阻率反演時，取矩形異常體電阻率為100 Ωm，相對介電常數(shù)為5，即只有電阻率異常(圖4c).

在10 kHz～400 kHz頻段對數(shù)等間距取13個頻點，點距10 m，共100個測點，采用有限單元法進行正演計算，將視電阻率數(shù)據(jù)加入5%的隨機噪聲作為反演數(shù)據(jù)，進行反演試算.反演取初始正則化因子λε和λρ均為1000，后期根據(jù)數(shù)據(jù)誤差、模型誤差等信息進行經(jīng)驗選?。环囱莸某跏寄Ｐ蜑楸尘熬鶆蚰Ｐ?，在迭代72次后得到反演的相對介電常數(shù)模型圖4b與電阻率模型圖4d.反演結果表明矩形相對介電常數(shù)和電阻率異常體的物性值和輪廓均得到較為準確的恢復，表明本文反演程序正確.

圖4 相對介電常數(shù)和電阻率單獨反演結果(a)、(c) 分別為設計的相對介電常數(shù)和電阻率模型; (b)、(d) 分別為反演的相對介電常數(shù)和電阻率模型.Fig.4 Separate inverted model of relative permittivity and resistivity(a) and (c) are designed relative permittivity and resistivity model respectively, (b) and (d) are inverted relative permittivity and resistivity model respectively.

3.2 FCM聚類雙參數(shù)聯(lián)合反演

為研究FCM聚類對介電常數(shù)和電阻率聯(lián)合約束的效果，設計了帶地形的多異常模型見圖5，背景電阻率為1000 Ωm，相對介電常數(shù)為10，其中加入三個矩形異常體編號(1)、(2)、(3)，其屬性見表1.

圖5 設計的帶地形模型(a) 相對介電常數(shù)模型; (b) 電阻率模型.Fig.5 Sketch of designed model with terrain(a) Relative permittivity model; (b) Resistivity model.

表1 地下異常體屬性Table 1 Properties of the underground blocks

通過正演計算圖5所示模型的響應，測點距、測點數(shù)、頻點與模型試算3.1節(jié)相同.分別討論參考聚類中心權重因子ηk、聚類項權重因子β、最小結構正則化因子λε與λρ對反演的影響.

首先討論參考聚類中心權重因子ηk對反演結果的影響，根據(jù)位移電流與傳導電流的比值至少為1‰，介電常數(shù)能產(chǎn)生明顯影響，設置λε=1、λρ=1000；聚類項的權重因子β為1固定不變；參考聚類中心為t1[10000 Ωm，3]、t2[10 Ωm，50]、t3[20 Ωm，30]、t4[1000 Ωm，10]，方括號內第1個值為電阻率，第2個值為相對介電常數(shù).在本文算例中認為四個參考聚類中心的置信度相同，后文中將ηk統(tǒng)一表示為η，分別對參考聚類中心權重因子η為1、100、10000的情況進行試算，同時對比不施加FCM聚類聯(lián)合約束(β=0)的雙參數(shù)反演結果，試算結果見圖6，左側、右側分別為相對介電常數(shù)與電阻率反演結果，(a)、(e)為無聚類約束的結果，(b)、(f)，(c)、(g)，(d)、(h)分別對應η=1、100、10000.

由圖6可見電阻率反演均取得了較好的效果，并且隨參考聚類中心權重因子η的增加邊界略有改善.不施加聯(lián)合約束得到的相對介電常數(shù)反演模型(圖6a)中，(3)號高阻異常體的介電常數(shù)得到較好恢復，但兩個低阻異常體的介電常數(shù)并未恢復，這與低阻異常體內位移電流不足有關；圖6b中低阻體的介電常數(shù)仍未得到有效恢復；圖6c和圖6d中低阻體的介電常數(shù)得到一定恢復，但邊界、物性值與真實模型均有較大差距，總體上圖6c和圖6d的反演模型類似，通過加大參考聚類中心權重的影響已經(jīng)到達極限.

根據(jù)圖6的反演結果可知，基于FCM聚類的電阻率與介電常數(shù)聯(lián)合反演具有可行性.為改善反演效果，進一步討論聚類項權重因子β、最小結構正則化因子λε與λρ的影響.當η=100時聚類約束對異常的恢復有一定改善，但更大的η對反演結果沒有改善，所以取η=100、λε=1、λρ=1000討論β=10、100時的反演結果，見圖7.圖7左側、右側分別為相對介電常數(shù)與電阻率反演結果，(a)、(c)，(b)、(d)分別對應β=10、100.

圖6 不同參考聚類中心權重因子η對反演結果的影響(λε=1、λρ=1000，β=1)左側為相對介電常數(shù)反演模型，右側為電阻率反演模型，其中(a)、(e)為無聚類約束的結果，(b)、(f)，(c)、(g)，(d)、(h)分別對應η=1，η=100，η=104.Fig.6 Test for target cluster centers weight factor η with λε=1, λρ=1000 and β=1Left four figures are inverted relative permittivity models and right four figures are inverted resistivity models, (a) and (e) show the recovered models without FCM clustering constraint, and the test conditions for (b) and (f), (c) and (g), (d) and (h) are η=1, η=100, η=104 respectively.

對比分析圖6和圖7可以發(fā)現(xiàn)，隨著β值增大，反演的相對介電常數(shù)異常數(shù)值有明顯改善，當β增大到10，圖7a中(1)號低阻異常體的相對介電常數(shù)值接近真實值30，但邊界未得到明顯改善；當β增大到100時(圖7b)，相對介電常數(shù)異常數(shù)值有所改善，但異常體邊界變差，而(2)號異常體右下方還出現(xiàn)假異常；電阻率反演模型，β=10時(圖7c)高阻異常體電阻率值更接近真值，當β增大到100時(圖7d)，(2)號異常體電阻率模型出現(xiàn)了類似于介電常數(shù)模型的假異常.

圖7 不同F(xiàn)CM聚類約束權重因子β對反演的影響(λε=1、λρ=1000，η=100)左側為相對介電常數(shù)反演模型，右側為電阻率反演模型，其中(a)、(c)，(b)、(d)分別對應β=10，β=100.Fig.7 Test for FCM clustering constraint weight factor β with λε=1、λρ=1000，η=100Left two figures are inverted relative permittivity models and right two figures are inverted resistivity models, the test conditions for (a) and (c), (b) and (d) are β=10, β=100 respectively.

為分析FCM聚類聯(lián)合約束效果，繪制了無聯(lián)合約束反演結果與施加聯(lián)合約束(β=10，η=100)相對介電常數(shù)與電阻率交會圖，分別見圖8(a、b)，橙色點為真實模型，藍色點為反演模型.無FCM聚類聯(lián)合約束(圖8a)中反演物性值未形成有效的聚類，而施加聯(lián)合約束(圖8b)反演物性在真實聚類中心附近聚集.

圖8 有無FCM聚類約束時反演的相對介電常數(shù)和電阻率交會圖對比(a) 無聚類約束的反演物性交會圖; (b) FCM聚類聯(lián)合反演物性交會圖，橙色點為真實模型物性，藍色點為反演模型物性.Fig.8 Cross-plot comparison of relative permittivity and resistivity inversion with or without FCM clustering constraint(a) is the cross-plot of inversion without clustering joint constraint, and (b) is the cross-plot of inversion with FCM clustering joint constraint, the orange dots represent true model values and the blue dots represent inverted model values.

為了進一步改進反演效果，在β=10、η=100條件下，繼續(xù)討論λε、λρ對反演效果的影響.分別采用λε=10、λρ=104和λε=100、λρ=105兩組參數(shù)進行反演，反演結果見圖9.圖9(a、b)為相對介電常數(shù)反演模型，圖9(c、d)為對應的電阻率反演模型，其中(a、c)，(b、d)分別對應λε=10、λρ=104和λε=100、λρ=105.由圖9可見當λε=10、λρ=104時無論電阻率還是介電常數(shù)都較λε=100、λρ=105時更好.

圖9 不同最小結構正則化因子λε、λρ對反演的影響(β=10，η=100)左側為相對介電常數(shù)反演模型，右側為電阻率反演模型，其中(a)、(c)，(b)、(d)分別對應λε=10、λρ=104，λε=100、λρ=105.Fig.9 Test for minimum structure regularized factors λε and λρ with β=10, η=100Left two figures are inverted relative permittivity models and right two figures are inverted resistivity models, the test conditions for (a) and (c), (b) and (d) are λε=10、λρ=104, λε=100、λρ=105 respectively.

圖6、圖7、圖9綜合對比分析表明，10～400 kHz頻段RMT數(shù)據(jù)基于FCM聚類聯(lián)合約束的反演算法可以實現(xiàn)介電常數(shù)與電阻率雙參數(shù)反演，尤其是在參考聚類中心tk為真值的前提下，合理選擇參考聚類中心權重因子η、聚類項權重因子β、最小結構正則化因子λε與λρ等參數(shù)可取得理想反演效果.

4 結論與展望

本文在同時考慮傳導電流與位移電流影響的條件下，開展了二維RMT數(shù)據(jù)電阻率和介電常數(shù)的聯(lián)合反演研究.在經(jīng)典最小結構模型正則化約束的基礎上，采用FCM聚類算法進行聯(lián)合約束，提高了電阻率與介電常數(shù)的反演效果.研究工作的成果與不足如下.

(1)二維RMT數(shù)據(jù)，同時反演電阻率和介電常數(shù)具有可行性，但對位移電流占比較小的低阻異常體，介電常數(shù)反演結果不理想；

(2)采用FCM聚類約束進行電阻率與介電常數(shù)聯(lián)合反演，可有效改善雙參數(shù)同時反演的效果；當參考聚類中心為已知值時，可有效地引導反演方向，從而在一定程度上改善了低阻體介電常數(shù)反演效果不理想的問題，但聚類項參數(shù)的選擇很關鍵；

(3)引入FCM聚類約束后，增加了正則化目標函數(shù)的不確定參數(shù)，研究工作尚未研發(fā)出自動選擇相關參數(shù)的算法；

(4)當前研究頻段(10～400 kHz)，低阻介質位移電流占比不足，后續(xù)將進一步開展更高頻段的電阻率與介電常數(shù)反演研究工作.

致謝感謝課題組成員在本文研究工作及寫作過程中的辛苦勞動，特別感謝王順國博士對論文修改提出的寶貴意見，并向匿名審稿人提出的寶貴修改意見表示衷心的感謝.