基于解析雅克比矩陣的E-Ex廣域電磁法一維并行約束反演

索光運， 李帝銓， 胡艷芳

(中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長沙 410083)

0 引言

廣域電磁法采用電磁場單分量進行測量并定義具有全區(qū)特性的廣域視電阻率，其主要優(yōu)點是可以用較小的收發(fā)距獲得很大的有效勘探測深度，裝備輕便，野外工作效率高，在油氣、礦產(chǎn)勘探中有著廣泛的應(yīng)用前景[1-2]。目前，投入實際應(yīng)用的主要是E-Ex形式的廣域電磁法，其采用電性源供電，并測量與電性源平行的電場分量。

關(guān)于電磁法正反演的研究現(xiàn)在正由1D向2D、3D轉(zhuǎn)變[3-6]，但在野外實測數(shù)據(jù)的反演解釋工作中，由于計算時間精度的限制，使用最多的是1D反演技術(shù)。Routh[7]進行了CSAMT一維全資料反演，反演結(jié)果要比做過近場校正后的MT方法要好；王若等[8]聯(lián)合使用網(wǎng)格參數(shù)法和剝皮法進行了不做校正的CSAMT一維全資料反演，該方法能較好地再現(xiàn)模型的層狀結(jié)構(gòu)，并對噪聲有一定的壓制效果；李帝銓等[9]利用遺傳算法完成了CSAMT的最小構(gòu)造反演，該方法具有不用考慮近場校正、初始模型依賴小、反演結(jié)果更真實等優(yōu)點；湯井田等[10]實現(xiàn)了基于頻點的最小構(gòu)造反演，使得初始模型的建立更加簡單、靈活；周俊杰等[11]提出用Bostick半定量結(jié)果作為一維反演的初始模型，結(jié)合多種先驗信息優(yōu)化層參數(shù)并控制反演過程的迭代方向，反演結(jié)果更加貼近真實模型；余云春[12]采用最小二乘法實現(xiàn)了廣域電磁法的一維反演，不過其方法對初值依賴性較大。

除了少數(shù)幾種完全非線性反演方法以外，在多數(shù)反演方法的計算過程中都需要計算雅可比矩陣。在一維反演中，傳統(tǒng)上習(xí)慣使用擾動法來計算雅克比矩陣，即給單個自變量一個微小的擾動，再次正演獲得擾動后的電磁響應(yīng)，從而計算擾動前后的一階差分來近似代替雅克比矩陣。如果對于N層模型同時進行電阻率和厚度兩個參數(shù)的反演，則一次反演迭代，就需要進行2N-1次正演。當(dāng)層數(shù)較多時，計算時間將會大大增加。而使用解析法時，雅克比矩陣將直接由公式給出[13]，其計算時間大致與1次～2次正演相同。

筆者使用探測深度和等對數(shù)域剖分的方式來剖分層厚，先使用MT反演來建立初始模型，再采用MATLAB自帶的最優(yōu)化工具箱來進行反演最優(yōu)化計算，并且使用了解析法來計算雅克比矩陣，替代原函數(shù)中的擾動法，從而大大加快了反演的計算速度。另外，還使用了parfor循環(huán)來對反演過程進行并行計算。對于實際工區(qū)資料，分別進行了無約束和井震約束時的反演計算，并將兩種方法的結(jié)果進行了對比。

1 解析法求靈敏度矩陣

1.1 E-Ex廣域視電阻率對層參數(shù)的偏導(dǎo)數(shù)

E-Ex廣域視電阻率是將所測得的|Ex|看作是均勻半空間的電場響應(yīng)，通過迭代法或逆樣條插值法求解出的該均勻半空間的電阻率。其對層參數(shù)的偏導(dǎo)數(shù)可由式(1)與式(2)給出。

(1)

(2)

其中，式(2)是兩個復(fù)數(shù)表達式的乘積。在實際計算時，由于計算誤差的存在，其乘積往往并不是實數(shù)，此時無論是取乘積的實部或者模，計算結(jié)果都與真實值存在一定差別。因此，這里采用式(1)來計算E-Ex廣域視電阻率對層參數(shù)的偏導(dǎo)數(shù)。計算時要注意一點，復(fù)變函數(shù)的模的導(dǎo)數(shù)并不等于該函數(shù)的導(dǎo)數(shù)的模。

(3)

正確的表達式應(yīng)該是

(4)

1.2 E-Ex廣域視電阻率對|Ex|的偏導(dǎo)數(shù)

關(guān)于均勻半空間表面電偶極子的電場響應(yīng)，當(dāng)采用不同的電流源(正諧或負諧)時，其表達式有所不同。這里采用文獻[14]中的表達式：

(5)

(6)

Ex=Ercosφ-Eφsinφ

(7)

其中:k2=iωμ/ρ。對Ex廣域視電阻率求關(guān)于|Ex|的偏導(dǎo)數(shù)，有

(8)

其中，

(9)

(10)

(11)

將式(9)～式(11)代入到式(8)中，便可以計算出廣域視電阻率對|Ex|的偏導(dǎo)數(shù)。

1.3 |Ex|對層參數(shù)的偏導(dǎo)數(shù)

對于N層水平層狀介質(zhì)中的電場響應(yīng)，仍然采用文獻[14]中的公式：

(12)

(13)

其中，

(14)

空間頻率特性函數(shù)R*和R的迭代格式可寫成

(15)

(16)

式中：I為諧變電流的幅值；ω為諧變電流的圓頻率；μ為磁導(dǎo)率；φ為電偶極源中點到接收點的連線和電偶極源的夾角；k為波數(shù)；m為空間頻率；ρl和hl分別是第l層的電阻率和厚度，l=1、2、…、N。

因此，|Ex|對層參數(shù)的偏導(dǎo)數(shù)可表示為

(17)

(18)

這樣問題變成了如何計算Er和Eφ關(guān)于層參數(shù)的偏導(dǎo)數(shù)值。

假設(shè)有

再對式(12)和式(13)稍作變量代換并求偏導(dǎo)，可以得到

(19)

(20)

分別對A、B、C求關(guān)于層參數(shù)的偏導(dǎo)，由于n1和k1里面含有ρ1，在計算對ρ1的偏導(dǎo)時，與計算對其他層參數(shù)的偏導(dǎo)不同，需要單獨進行。

(21)

(22)

(23)

(24)

(25)

(26)

計算R*和R分別關(guān)于地層電阻率和厚度的偏導(dǎo)數(shù)的公式和推導(dǎo)過程可以參考文獻[13]。需要注意，文獻[13]中的(26)式可能是印刷錯誤，正確表達式應(yīng)該為

(27)

這樣|Ex|對層參數(shù)的偏導(dǎo)數(shù)已計算完成。

1.4 算例驗證

為了驗證解析法計算雅克比矩陣的正確性，將根據(jù)上面公式所編制的程序的計算結(jié)果與使用擾動法的計算結(jié)果進行了對比。

采用三層模型，由淺入深電阻率分別為100 Ω·m、1 000 Ω·m、100 Ω·m，前兩層的厚度都是300 m，第三層為均勻半空間。電偶極子長度為1 000 m，供電電流20 A，收發(fā)距6.4 km，接收點相對于發(fā)射偶極的角度為90°。頻率范圍為0.25 Hz ～8 192 Hz，按2的整指數(shù)冪變化，共16個頻點。擾動法的擾動量為1e-3乘以當(dāng)前值。圖1是兩種方法的計算結(jié)果。由圖1可以看到，由兩種方法計算出的偏導(dǎo)數(shù)值完全重合，從而驗證了所推導(dǎo)的解析法公式的準(zhǔn)確性，也為下一步廣域電磁法的快速反演打下了基礎(chǔ)。

2 E-Ex廣域電磁法一維約束反演

2.1 建立初始無約束模型

實際地層的層數(shù)很多，同時進行電阻率和層厚的反演，既容易陷入局部極小值，也容易導(dǎo)致反演時間過長。可以將地層剖分為很多層，每層的層厚都較小，只進行電阻率單參數(shù)反演。另外，考慮到分辨率隨深度增大成指數(shù)減小的特點，可以使用對數(shù)域剖分方法來剖分地層。首先，根據(jù)視電阻率曲線來確定最小、最大探測深度Dmin和Dmax。dmin、dmax可以是人為給出的最小、最大深度的參考值。然后再根據(jù)下式來剖分層厚h，其中，diff代表前向差分，logspace代表對數(shù)剖分。層數(shù)NL可以剖分的多一些，但是不能超過頻點數(shù)。

基里爾表示，BPC特別重視中國市場，BPC每年很大一部分氯化鉀出口量是針對中國市場。近幾年，國際市場和中國市場變得比較快，所以BPC必須每年做出相應(yīng)的計劃調(diào)整?！敖衲晡覀兊哪繕?biāo)是200萬噸。但這并不是最高的目標(biāo)，去年按照中國鉀肥進口數(shù)據(jù)，從白俄羅斯的鉀肥進口量約有160萬噸。但是，2015年BPC在中國的銷售量高達240萬噸?！被餇栂壬a充說。BPC每年向100個多個國家出口鉀肥，今年的出口量超過了1050萬噸，其中中國和印度市場各占18%左右。

圖1 三層模型計算結(jié)果Fig.1 Results of 3-layered model(a)dρ/dρ1曲線;(b)dρ/dρ2曲線;(c)dρ/dρ3曲線;(d)dρ/dh1曲線;(e)dρ/dh2曲線

圖2 處理約束層時的幾種情況Fig.2 Conditions of constrained layer(a)單層-深部；( b)單層-淺部；(c)多層-深部；(d)多層-淺部

(28)

h=diff(logspace(Dmin,Dmax,NL))

(29)

2.2 建立初始約束模型

建立約束模型的過程與上面類似，只是在約束層如何處理上有所不同。針對不同問題，分為單層和多層兩種方法。單層是所剖分的約束層會作為完整的一層存在于最終的模型數(shù)據(jù)中，對于深度大、厚度小的層狀體(油氣)，建議使用此種方式。多層即按照給定的層數(shù)和最大深度，在對數(shù)域上計算好分界點埋深后確定分界點處的地層。此種方式，所剖分的地層在最終的模型數(shù)據(jù)中通常會占據(jù)多層，但也可能一層都沒有。對于埋深較小的塊狀體(金屬礦)，建議使用此種方式。幾種可能情況如圖2所示。

2.3 反演算法

由于廣域視電阻率和MT卡尼亞電阻率在曲線形式上的相似性[11]，可以在進行帶源反演之前先使用MT反演來建立一個更為合理的初始模型。

MATLAB工具箱中提供了lsqnonlin函數(shù)來求解非線性最小二乘優(yōu)化問題[15]，其調(diào)用格式為X=lsqnonlin(fun,x0,lb,ub,options)。其中x0為初值，options為優(yōu)化選項，這里選用Levenberg-Marquardt法來進行最優(yōu)化計算。另外，ub和lb分別為自變量的上下限，反演時可以對地層電阻率的范圍進行約束。電磁法反演很容易受到多解性的影響，而根據(jù)先驗信息來對反演參數(shù)進行約束，可以有效地降低多解性，提高反演效果。

此外，現(xiàn)有的計算機都是多核多線程的，利用parfor循環(huán)能充分利用計算機的計算資源[16]，通過并行進一步縮短反演時間。

3 反演算例驗證

3.1 理論模型反演

選用四層KH模型來驗證反演算法的有效性，其電阻率由淺到深分別為100 Ω·m、1 000 Ω·m、100 Ω·m、1 000 Ω·m，前三層厚度為100 m、400 m、500 m。裝置參數(shù)與圖1基本相同，收發(fā)距變?yōu)?2 km，頻率變?yōu)?-6:0.5:13，共計39個頻點。模型剖分了30層，最大迭代次數(shù)為30次。圖3為視電阻率擬合曲線和反演結(jié)果。

圖3 四層KH理論模型反演結(jié)果Fig.3 Inversion result of 4-layered KH model(a)視電阻率與反演擬合曲線;(b)理論與反演模型電阻率深度圖

圖4 JSB測線兩種反演方法效果對比Fig.4 Results of inversion(a)地震剖面；(b)自由模型反演剖面；(c)井震約束反演剖面

3.2 實測數(shù)據(jù)反演

為驗證本文算法對實際工區(qū)數(shù)據(jù)的反演效果，選取某地區(qū)實測資料，分別進行無約束自由反演和井震約束反演。根據(jù)已有地質(zhì)資料，在埋深1 000 m～3 000 m范圍內(nèi)，地層大致分為3層，分別是二疊系下統(tǒng)、奧陶上統(tǒng)+志留系中統(tǒng)、二疊系下統(tǒng)，巖性分別為多厚層灰?guī)r和瀝青質(zhì)灰?guī)r、粉砂質(zhì)頁巖和炭質(zhì)頁巖、厚層灰?guī)r，整體電阻率由淺入深呈高-低-高變化趨勢。并且測線范圍內(nèi)，各地層基本成水平層狀分布，構(gòu)造主體平緩。圖5上是地震剖面圖，可以看到垂向上地層分布較為明顯，與地質(zhì)資料較為吻合。圖5中為無約束模型自由反演剖面圖。在1 000 m～3 000 m埋深范圍內(nèi)，電阻率變化只是呈低-高變化趨勢，對于頁巖層和其上覆的高阻層反映均不明顯，橫向上也不能明顯地分辨地層的分界面。圖5下是根據(jù)已有的地震剖面和測井資料，進行了關(guān)于地層大致位置和電阻率范圍的約束建模反演的剖面結(jié)果?？梢钥吹降貙臃謱虞^為明顯，頁巖層和其上覆的高阻層在橫向上也較為連續(xù)。并且15 000～17 000處的電阻率抬升趨勢與地震資料中劃分的斷裂構(gòu)造也較為吻合。此外，無約束自由反演的擬合誤差為4.95%，井震約束反演的擬合誤差為5.17%。在擬合誤差基本一致的情況下，結(jié)合了已有信息的約束反演，能更好反映出地下真實的地層分布。

在約束反演過程中，還分別統(tǒng)計了使用擾動法和解析法來計算雅克比矩陣時，整個反演過程所耗費的時間，如表1所示。本條測線共有測點342個，觀測資料中頻點數(shù)34，模型剖分層數(shù)也為34，最大迭代次數(shù)為30。并且在反演過程中使用了parfor并行循環(huán)，設(shè)置核數(shù)為4，計算機配置為4核Intel-i7處理器，主頻3.70 GHz，8 G內(nèi)存。

表1 擾動法和解析法的計算時間對比

兩種方法的計算結(jié)果一致，只是運行時間有所不同。解析法的計算時間只有擾動法的1/35，計算效率有了大幅度提高。

4 結(jié)論與建議

1)對于層厚的劃分，采用探測深度和對數(shù)域剖分的方法，考慮了實測視電阻率和分辨率隨深度成指數(shù)降低的影響，使反演模型更加貼近真實的地層情況。

2)對于N層模型的單參數(shù)反演，使用擾動法，一次反演迭代就需要進行N次正演；而使用解析法時，其計算時間大致只與1次～2次正演相同。當(dāng)層數(shù)和迭代次數(shù)較多時，解析法的計算速度將會遠遠高于擾動法。

3)通過parfor循環(huán)可以有效縮短反演計算時間，不過這只是基于處理器的并行運算，如果能實現(xiàn)基于GPU的并行計算，反演速度將會更快。

4)對于實際地電模型，由于充分考慮到已有的先驗信息，根據(jù)已有的測井或地震資料進行約束反演可以有效降低多解性，使反演結(jié)果更加接近真實的地層情況。不過，現(xiàn)有的約束反演方法中先驗約束對于反演結(jié)果的影響較大，如何平衡實測資料和先驗信息的權(quán)重，將會是下一步研究的重點。