索光運(yùn)， 李帝銓， 胡艷芳

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

0 引言

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

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

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

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

1 解析法求靈敏度矩陣

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

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

(1)

(2)

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

(3)

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

(4)

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

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

(5)

(6)

Ex=Ercosφ-Eφsinφ

(7)

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

(8)

其中，

(9)

(10)

(11)

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

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

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

(12)

(13)

其中，

(14)

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

(15)

(16)

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

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

(17)

(18)

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

假設(shè)有

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

(19)

(20)

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

(21)

(22)

(23)

(24)

(25)

(26)

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

(27)

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

1.4 算例驗(yàn)證

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

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

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

2.1 建立初始無(wú)約束模型

實(shí)際地層的層數(shù)很多，同時(shí)進(jìn)行電阻率和層厚的反演，既容易陷入局部極小值，也容易導(dǎo)致反演時(shí)間過(guò)長(zhǎng)?？梢詫⒌貙悠史譃楹芏鄬樱繉拥膶雍穸驾^小，只進(jìn)行電阻率單參數(shù)反演。另外，考慮到分辨率隨深度增大成指數(shù)減小的特點(diǎn)，可以使用對(duì)數(shù)域剖分方法來(lái)剖分地層。首先，根據(jù)視電阻率曲線來(lái)確定最小、最大探測(cè)深度Dmin和Dmax。dmin、dmax可以是人為給出的最小、最大深度的參考值。然后再根據(jù)下式來(lái)剖分層厚h，其中，diff代表前向差分，logspace代表對(duì)數(shù)剖分。層數(shù)NL可以剖分的多一些，但是不能超過(guò)頻點(diǎn)數(shù)。

基里爾表示，BPC特別重視中國(guó)市場(chǎng)，BPC每年很大一部分氯化鉀出口量是針對(duì)中國(guó)市場(chǎng)。近幾年，國(guó)際市場(chǎng)和中國(guó)市場(chǎng)變得比較快，所以BPC必須每年做出相應(yīng)的計(jì)劃調(diào)整。“今年我們的目標(biāo)是200萬(wàn)噸。但這并不是最高的目標(biāo)，去年按照中國(guó)鉀肥進(jìn)口數(shù)據(jù)，從白俄羅斯的鉀肥進(jìn)口量約有160萬(wàn)噸。但是，2015年BPC在中國(guó)的銷(xiāo)售量高達(dá)240萬(wàn)噸。”基里爾先生補(bǔ)充說(shuō)。BPC每年向100個(gè)多個(gè)國(guó)家出口鉀肥，今年的出口量超過(guò)了1050萬(wàn)噸，其中中國(guó)和印度市場(chǎng)各占18%左右。

圖1 三層模型計(jì)算結(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 處理約束層時(shí)的幾種情況Fig.2 Conditions of constrained layer(a)單層-深部；( b)單層-淺部；(c)多層-深部；(d)多層-淺部

(28)

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

(29)

2.2 建立初始約束模型

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

2.3 反演算法

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

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

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

3 反演算例驗(yàn)證

3.1 理論模型反演

選用四層KH模型來(lái)驗(yàn)證反演算法的有效性，其電阻率由淺到深分別為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，共計(jì)39個(gè)頻點(diǎn)。模型剖分了30層，最大迭代次數(shù)為30次。圖3為視電阻率擬合曲線和反演結(jié)果。

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

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

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

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

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

表1 擾動(dòng)法和解析法的計(jì)算時(shí)間對(duì)比

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

4 結(jié)論與建議

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

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

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

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