郭來功 戴廣龍 楊本才 薛俊華 陳本良

(①安徽理工大學電信學院，安徽淮南 232001； ②安徽理工大學能源學院，安徽淮南 232001；③深部煤炭開采與環(huán)境保護國家重點實驗室，安徽淮南 232001)

1 引言

直流電法是一種重要的地球物理勘探方法，在礦井突水、地質(zhì)滑坡、隧道施工等巖土工程領(lǐng)域，地下水勘探、環(huán)境地質(zhì)監(jiān)測等環(huán)境工程領(lǐng)域，以及礦產(chǎn)資源調(diào)查等資源勘探領(lǐng)域均有廣泛應用[1-3]。直流電法三維成像存在固有的多解性問題，會出現(xiàn)異常體定位不準確的現(xiàn)象。另外，反演數(shù)據(jù)量大、求解效率不高，其在時效性要求較高的應用中易受到限制，因此提高反演精度和計算效率是直流電三維勘探需要解決的主要問題[4,5]。

在三維反演中引入更多的先驗信息進行約束，是提高成像效果的主要手段之一。劉征宇等[6]提出幾種有效方法，一是引入距離加權(quán)約束，構(gòu)建距離加權(quán)因子分布模型，對于不同的電源點距離，施加指數(shù)規(guī)律的權(quán)重因子，對不同單元賦予差異權(quán)重，反演的異常體形態(tài)得到有效改善；二是引入不等式約束，將介質(zhì)電阻率的取值范圍作為先驗約束條件，提高反演精度[7]；三是引入?yún)⒖寄Ｐ图s束，將已經(jīng)獲取的異常體結(jié)構(gòu)信息作為先驗信息施加到約束條件中，有效去除反演成像中的假異常[8]。劉鑫明等[9]采用光滑約束構(gòu)造三維反演函數(shù)中的粗糙度矩陣，采用合理的網(wǎng)格剖分，在注漿檢測中取得良好的效果。吳小平等[10]采用不完全Gauss-Newton法，應用非結(jié)構(gòu)網(wǎng)格實現(xiàn)任意偶極—偶極視電阻率的三維反演，較好地展現(xiàn)了地下真實電阻率結(jié)構(gòu)。Johnson等[11]在正演建模時通過埋設導電件，建立精準模型，提高反演成像效果。Johnson 等[12]還采用并行分布式計算方法，提高了反演效率。Farquharson[13]針對最小構(gòu)造反演，提出一種分段常數(shù)的地質(zhì)模型構(gòu)造方法，采用L1范數(shù)準則下的迭代重加權(quán)最小二乘方法，成像結(jié)果與實際地質(zhì)結(jié)構(gòu)高度接近。在反演方法上，基于全局優(yōu)化的模擬退火方法、人工神經(jīng)網(wǎng)絡方法等受到計算機計算能力的限制，目前還無法全面應用到實踐，當前的工程應用以基于最小二乘準則的線性反演為主[14-18]。

實際應用中，仍有諸多問題需要解決。關(guān)于反演算法，當目標函數(shù)確定后，先驗信息無法靈活地施加到約束中，例如引入距離加權(quán)約束后，無法再施加不等式約束信息；另外，反演耗時過長，普通計算機無法實現(xiàn)快速反演。本文采用分布式并行計算方法[19]，在最小二乘準則的基礎(chǔ)上，在重加權(quán)函數(shù)中增加結(jié)構(gòu)度量、權(quán)重均值和標準偏差等參數(shù)，對模型的目標區(qū)域施加不同約束，并通過合成數(shù)據(jù)進行直流電三維反演算例分析。

2 目標函數(shù)的建立

三維有限元模型空間中，模型參數(shù)用向量m表示，模型響應為f，測量數(shù)據(jù)用列向量d表示。反演過程中，網(wǎng)格固定不變，其目的是確定模型參數(shù)。反演通過構(gòu)造一個包含模型參數(shù)和測量數(shù)據(jù)的目標函數(shù)，使其最小化。根據(jù)Farquharson[13]提出的一般準則，目標函數(shù)Φ可表示為

Φ=Φd(ud)+βΦm(um)

(1)

式中：下標“m”和“d”分別代表模型和數(shù)據(jù)；Φd(ud)是度量模型響應數(shù)據(jù)與測量數(shù)據(jù)d之間擬合差的函數(shù)；Φm(um)是度量模型參數(shù)變化的函數(shù)；β是正則化參數(shù)，控制擬合數(shù)據(jù)和模型向量變化對目標函數(shù)的影響。其中ud的表達式為

ud=Wd(d-f)

(2)

式中Wd是數(shù)據(jù)加權(quán)對角矩陣，其元素是測量數(shù)據(jù)中誤差(噪聲)的標準偏差σi的倒數(shù)，即1/σi。如果測量時沒有考慮誤差，則需對誤差做合理的估計[20]。假設誤差由測量電壓誤差ΔU和一個固定百分比誤差(q%)構(gòu)成，則第i個數(shù)據(jù)的標準偏差為

(3)

式中：Ui為第i個測量電壓； ΔUi為第i個數(shù)據(jù)的測量電壓誤差。目標函數(shù)Φ中模型參數(shù)向量um可寫作

um=Wm(m-mref)

(4)

式中：Wm是模型參數(shù)權(quán)重矩陣；mref表示參考模型。在一種簡單的情況下，考慮三維空間的各向異性，有

(5)

式中：Cx、Cy和Cz分別是x、y、z方向的平滑矩陣；αx、αy和αz分別是x、y、z方向的相對權(quán)重。平滑矩陣C的一般形式是一階差分矩陣[21]。

3 基于迭代重加權(quán)最小二乘法的三維反演

采用最小二乘法時，反演方程的求解采用將非線性問題線性化的方法，經(jīng)過多次迭代修正模型參數(shù)，使估計解不斷逼近真實解。迭代重加權(quán)的核心思想是通過計算一系列權(quán)系數(shù)并不斷更新，得到目標函數(shù)的最優(yōu)解。假設在第n次迭代時，得到最佳參數(shù)模型m(n)，使目標函數(shù)最小化。由式(1)，第n次迭代的目標函數(shù)可以表示為

(6)

設模型參數(shù)向量的變化量為Δm=m(n)-m(n-1)，則第n次迭代時，得到模型m(n-1)的正演結(jié)果

f(n)≈f(n-1)+JΔm

(7)

式中J為偏導數(shù)的靈敏度雅克比矩陣，其元素為

(8)

表示第j個模型參數(shù)的變化導致的第i個模型響應的變化，因此，由式(2)和式(4)可得

ud=Wd[d-f(n-1)-JΔm]

(9)

um=Wm[m(n-1)+Δm-mref]

(10)

式(6)是關(guān)于Δm的函數(shù)，最小化式(6)，令其偏導為零，結(jié)合式(9)、式(10)，應用Gauss-Newton法求解，可得

(11)

式中：R是對角矩陣。Rd的元素為

Rd,i,i=

(12)

式中：p表示所使用范數(shù)的階次；ε是一個極小的數(shù)，用以保證di-fi-1-Ji,jΔmi→0； 矩陣Rd是非奇異的；Rm與Rd具有相似的形式，其元素為

Rm,i,i=

(13)

式(11)寫成通用形式

(14)

在求解地球物理反演問題時，以Lp范數(shù)(p≥1)作為測度進行解的估計，根據(jù)觀測數(shù)據(jù)和模型參數(shù)的統(tǒng)計特征，使用不同的范數(shù)。如果實際地下電阻率平滑地漸變，測量誤差服從高斯分布，L2范數(shù)是最合適的方法，它對每一個元素進行平方，將向量的L2范數(shù)最小化，在數(shù)據(jù)擬合時，剔除相鄰單元的異常值，空間區(qū)域之間無明顯邊界，排除了空間結(jié)構(gòu)的突變。但在實際地質(zhì)結(jié)構(gòu)中，如煤層、金屬礦山等的邊界會出現(xiàn)明顯的突變，此時，需要采用L1范數(shù)求解，它對模型電阻率值的絕對值變化最小化，對誤差大的數(shù)據(jù)不敏感，保留了相鄰元素之間的異常值。

4 約束條件

對于迭代重加權(quán)最小二乘法，每次迭代的權(quán)函數(shù)均不相同，為了在重加權(quán)函數(shù)中施加先驗約束信息，定義結(jié)構(gòu)度量方程[19]

X=|mi-mt|

(15)

式中，與L1范數(shù)不同，mt在不同約束信息下代表不同的模型信息，例如相鄰單元的電導率、常數(shù)電導率等。引入關(guān)于X的重加權(quán)函數(shù)Wf，根據(jù)不同的約束要求，Wf可以用不同的函數(shù)替代。本文選擇互補誤差函數(shù)Wf1或誤差函數(shù)Wf2

(16)

(17)

式中：M是重加權(quán)的權(quán)重均值；S是權(quán)重標準偏差。重加權(quán)函數(shù)是一個基于正態(tài)的累積分布函數(shù)，其目的是在X變化時，施加不同的權(quán)重。Wf1、Wf2隨X的變化如圖1所示。

圖1 權(quán)重函數(shù)曲線

(1)L2范數(shù)平滑約束：mt代表相鄰單元的電導率，選擇Wf1，取一個較大的M值和較小的S值，當每個單元與相鄰單元的電導率差異減小且小于M+2S時，則加權(quán)函數(shù)施加權(quán)重；當電導率差異小于M-2S時，則施加滿權(quán)，使X最小化。此時的約束鼓勵兩個相鄰單元的電導率接近，等效于L2范數(shù)施加在模型空間。

(2)L1范數(shù)塊狀約束：mt代表相鄰單元的電導率，選擇Wf1，取一個較小的M值和較大的S值，此時，當每個單元與相鄰單元的電導率之差增加時，施加在相鄰單元的電導率平滑約束值減少； 當二者的差異大于M+2S時，平滑約束被移除。這樣，相鄰單元的電導率差異被保留，兩個單元的電導率值之間形成明顯邊界，等效于L1范數(shù)施加在模型空間。

(3)電導率極值約束：通過地質(zhì)勘探等手段獲取地下電導率最大值和最小值的合理范圍時，反演可以施加電導率極值約束。mt代表電導率的極限值，選擇Wf1，設定M=0和一個較小的S值。此時將每個單元的電導率與mt的差值作為約束，當X小于M-2S時，權(quán)重函數(shù)提供滿權(quán)以最小化X；當X大于M+2S時，權(quán)重函數(shù)為0，移除約束，這使得反演結(jié)果中的元素電導率的最小值為mt。同樣，如果選擇Wf2，則反演結(jié)果中的元素電導率的最大值為mt。

在已知勘探區(qū)域的電阻率極大值和極小值、區(qū)域斷面位置和地下空間結(jié)構(gòu)等信息的條件下，可將反演區(qū)域分為若干子區(qū)域，通過修改重加權(quán)函數(shù)的參數(shù)，增加、減少或改變每個子區(qū)域的約束信息，提高反演精度。

5 反演方法

共軛梯度法是一種常用的地球物理反演優(yōu)化算法，在求解式(14)時，采用并行共軛梯度最小二乘法(Parallel Conjugate-Gradient Least-Squares,PCGLS)[12]。該方法將處理器分為主處理器和多個從處理器，反演時由主處理器分配任務給不同的從處理器，自身也同時執(zhí)行計算數(shù)據(jù)量較小的任務，待從處理器完成計算任務后，主處理器再分配下一次計算任務，主、從處理器協(xié)同工作，達到提高計算效率的目的，IRLS反演流程見圖2。令式(14)的右側(cè)為b，設PCGLS的內(nèi)迭代次數(shù)為j=0，1，2，…，n，λ(j)為Δm的修正因子，γ(j)是b的修正因子，數(shù)據(jù)擬合的條件是測量數(shù)據(jù)與模擬數(shù)據(jù)之間的均方誤差小于設定的容許值，算法收斂的條件是滿足數(shù)據(jù)擬合，或者內(nèi)部迭代次數(shù)j達到預設值。

圖2 IRLS法迭代流程

6 合成數(shù)據(jù)反演實例

6.1 三維模型

為了評價最小二乘法三維電阻率反演效果，利用合成數(shù)據(jù)進行三維反演。應用有限元法，由TetGen軟件建立四面體非結(jié)構(gòu)網(wǎng)格三維模型，并采用局部加密的方式，減少網(wǎng)格數(shù)量，提高網(wǎng)格精度。三維模型如圖3a所示，四面體模型由三角形網(wǎng)格構(gòu)成，有限元空間定義為800m×800m×500m，監(jiān)測區(qū)尺寸為16m×12m×10m，內(nèi)部設有兩個2m×2m×2m的低阻異常體，其中低阻塊1的電導率為0.020S/m，低阻塊2的電導率為0.200S/m，圍巖電導率為0.002S/m，低阻塊埋深均為1m。測區(qū)布置3條測線，共48個電極，電極距為1.0m，測線間距為5.5m，電極分布和剖面位置見圖3b和圖3c。生成的合成數(shù)據(jù)引入5%高斯噪聲。

圖3 合成模型(a)三維模型示意圖； (b)模型水平截面圖； (c)模型y=0切面圖

為了對比混合約束條件下先驗信息的效果，分別采用平滑約束、塊狀約束和混合約束對該模型進行反演，比較不同約束信息下的反演結(jié)果。計算機CPU為 Intel i5-3470，主頻為3.2G，內(nèi)存為16GB，模型網(wǎng)格信息見表1。

有限元正演過程采用預條件共軛梯度法計算剛度矩陣，求解空間的電場分布[7,23]。

表1 模型參數(shù)信息表

6.2 不同先驗信息約束的反演

6.2.1 平滑約束反演

對測區(qū)使用平滑約束反演。結(jié)構(gòu)度量中，mt代表相鄰單元的電導率，選擇Wf1，設M=10，標準偏差S=0.1，測區(qū)及低阻塊采用同樣的約束規(guī)則。反演結(jié)果如圖4所示，可見在平滑約束下，反演結(jié)果中低阻異常體的位置基本正確，但反演精度較低，電導率從低阻的中心位置向圍巖逐漸降低，地質(zhì)結(jié)構(gòu)具有平滑過渡的特征；在低阻塊2中，反演低阻區(qū)域形態(tài)與模型不完全一致，低阻區(qū)向測線位置延伸。

圖4 平滑約束反演結(jié)果(a)三維切片； (b)y=0切片； (c)y=1m切片； (d)y=2m切片

6.2.2 塊狀約束反演

對測區(qū)及低阻塊施加塊狀約束反演(非連續(xù)地質(zhì)結(jié)構(gòu)反演)，仍采用函數(shù)Wf1，結(jié)構(gòu)度量中mt代表相鄰單元的電導率，約束信息的施加通過減小權(quán)重函數(shù)的均值實現(xiàn)。取值M=0.1、S=0.01時的反演結(jié)果如圖5所示，可以看出，M值減小后，反演時會移除更多的相似約束，導致相鄰單元的電導率呈現(xiàn)較大差異，反演結(jié)果中低阻異常體位置基本準確，但精度不夠。與地質(zhì)原型相比，形態(tài)出現(xiàn)較大偏差，低阻塊2表現(xiàn)得尤為明顯，同時低阻塊1和低阻塊2向測線位置偏移。

從圖4和圖5的反演結(jié)果可以看出，單純采用平滑約束(L2范數(shù))或者塊狀約束(L1范數(shù))的反演精度均不夠理想。在工程實踐中，由于測量誤差和人工操作的影響，反演結(jié)果的精度會更低。

6.2.3 混合約束1反演

對測區(qū)和低阻塊同時施加平滑約束、最大值約束、最小值約束進行反演(稱為混合約束1)，對測區(qū)和低阻塊分別施加約束值?？紤]實際地質(zhì)結(jié)構(gòu)中，同一塊地質(zhì)單元電導率可能會有變化，甚至可能出現(xiàn)數(shù)倍的差異，因此設定某區(qū)域電導率最大值和最小值時，保留較大的余量。對每個區(qū)域同時施加三個約束條件，以測區(qū)為例，施加約束信息為：①對于平滑約束，mt代表相鄰單元的電導率，選擇Wf1，設定M=10，S=0.1； ②對于最大電導率值約束，設定mt=0.01S/m，選擇Wf1，M=0，S=0.1；③對于最小電導率值約束，設定mt=0.0004S/m，選擇Wf2，M=0，S=0.1。對低阻塊1和低阻塊2，同樣施加平滑約束和電導率極值約束信息，其中低阻塊1最大電導率值mt=0.1S/m，最小電導率值mt=0.004S/m；對于低阻塊2，設定最大電導率值mt=1.0S/m，最小電導率值mt=0.04S/m。三種約束信息同時施加的反演結(jié)果如圖6所示，可以看出，施加多個約束信息后，反演結(jié)果的精度有較大的提高，反演的測區(qū)電導率值較接近于理論模型，其中最大電導率σmax=0.03705S/m，最小電導率σmin=0.001209S/m，低阻塊的位置形態(tài)與理論模型基本一致。需要說明的是，對低阻塊2施加電導率極值約束(0.04S/m≤σ≤1.00S/m)后，實際上該區(qū)域電導率接近0.04S/m，在圖例中接近紅色區(qū)域，因此視覺上像是一個紅色方塊。

圖5 非連續(xù)地質(zhì)結(jié)構(gòu)約束反演結(jié)果(a)三維切片； (b)y=0切片； (c)y=1m切片； (d)y=2m切片

圖6 混合約束1反演結(jié)果(a)三維切片； (b)y=0切片； (c)y=1m切片； (d)y=2m切片

6.2.4 混合約束2反演

對低阻塊1施加平滑約束和電導率極值約束，對低阻塊2施加塊狀約束和電導率極值約束(稱為混合約束2)，方法同混合約束1，電導率極值不變。反演結(jié)果如圖7所示，可以看出，對低阻塊2的約束信息改為塊狀約束后，成像質(zhì)量較好，低阻區(qū)域的位置和形態(tài)與模型高度一致。整個測區(qū)最大電導率σmax=0.04149S/m，最小電導率σmin=0.001688S/m。由于在低阻塊2中施加了塊狀約束和電導率極值約束，反演結(jié)果中低阻塊2邊界處的電導率允許出現(xiàn)不連續(xù)。而在低阻塊2內(nèi)部，電導率極值的約束條件為0.04S/m≤σ≤1.0S/m，反演得到該區(qū)域電導率為0.04～0.04149S/m，因此低阻塊2呈現(xiàn)單獨的塊狀形態(tài)。與低阻塊1所施加的平滑約束相比，在不連續(xù)地質(zhì)結(jié)構(gòu)中，塊狀約束可以更好地反映地質(zhì)體的結(jié)構(gòu)變化，這在礦產(chǎn)資源的勘探和開采監(jiān)測中尤其重要。

圖7 混合約束2反演結(jié)果(a)三維切片； (b)y=0m切片； (c)y=0.99m切片； (d)y=1.1m切片

6.3 不同約束下PCGLS法計算效率

施加不同約束信息時，PCGLS法反演性能數(shù)據(jù)統(tǒng)計見表2。合成算例中PCGLS算法計算時間在數(shù)分鐘至數(shù)十分鐘以內(nèi)，與文獻[16](約160min)、文獻[5](約86min)的方法相比，計算速度得到較大提高。不同約束信息對迭代次數(shù)和運行時間有較大的影響，從表2可見，混合約束2反演效果優(yōu)于混合約束1，且迭代次數(shù)和運行時間更少，因此在施加反演約束信息時，應根據(jù)地質(zhì)條件合理施加。圖8給出了迭代過程中數(shù)據(jù)失配數(shù)量變化，迭代中第一個值為起始模型的數(shù)據(jù)失配數(shù)量。迭代過程中，最大數(shù)據(jù)失配數(shù)量為54個，不超過總數(shù)據(jù)量的2.6%。擬合均方根誤差(RMSE)為最后一次迭代后的結(jié)果，平滑約束和混合約束1由于采用了平滑約束，地質(zhì)結(jié)構(gòu)更加光滑，反演結(jié)果中(自然對數(shù))電導率數(shù)據(jù)與測量數(shù)據(jù)偏差較大，均方根誤差大于1.5；而采用塊狀約束時，反演允許相鄰單元的電導率出現(xiàn)較大差異，保留較多的原始數(shù)據(jù)，因此擬合均方根誤差相對較小。

表2 不同約束信息的迭代次數(shù)和運行時間統(tǒng)計表

圖8 迭代過程中的數(shù)據(jù)失配數(shù)量

7 結(jié)論

基于迭代重加權(quán)最小二乘法實現(xiàn)三維電阻率反演時，將數(shù)據(jù)加權(quán)矩陣中的重加權(quán)函數(shù)修改為誤差函數(shù)和互補誤差函數(shù)的形式，通過對重加權(quán)函數(shù)參數(shù)做簡單改動，達到不同先驗信息約束的效果；反演過程中，主、從處理器多線程協(xié)同工作，共同完成反演迭代流程。合成數(shù)據(jù)的成像結(jié)果證明該方法的有效性，與單一信息約束的反演對比說明，合理施加多信息約束可以提高成像結(jié)果中異常體位置和形態(tài)的準確性；與同類型算法相比，反演耗時較短，有助于提高三維電阻率層析成像的時效。

結(jié)合其他勘探手段，獲得更多的地質(zhì)結(jié)構(gòu)信息，是提高直流電法三維成像效果的重要途徑，對所述方法進行工程驗證，并進一步挖掘測量數(shù)據(jù)的價值，是后續(xù)研究的重點。