胡祖志 石艷玲 劉云祥 劉雪軍 孫衛(wèi)斌 何展翔

(①中國(guó)石油集團(tuán)東方地球物理公司綜合物化探處，河北涿州 072751； ②南方科技大學(xué)前沿與交叉科學(xué)研究院，廣東深圳 518055； ③南方科技大學(xué)地球與空間科學(xué)系，廣東深圳 518055)

0 引言

重、磁、電單方法解決問(wèn)題的局限性、觀測(cè)數(shù)據(jù)的有限性以及數(shù)據(jù)存在的觀測(cè)誤差等諸多因素，造成了反演結(jié)果的多解性，影響了地質(zhì)解釋的可靠性。因此，綜合利用重、磁、電等方法，從不同角度研究同一地質(zhì)體，開展各種方法之間的聯(lián)合反演與綜合解釋，降低多解性，成為眾多學(xué)者的研究重點(diǎn)。

地震資料對(duì)中淺地層有很高的分辨率，但在復(fù)雜區(qū)對(duì)深層結(jié)構(gòu)揭示不清，而重、磁、電資料可以為地震勘探提供重要的補(bǔ)充。利用已知的地震、地質(zhì)、電測(cè)井等資料進(jìn)行重、磁、電資料的約束反演是提高重、磁、電勘探分辨率的重要途徑之一[1-10]。SEG(勘探地球物理學(xué)家學(xué)會(huì))為此主辦了多類關(guān)于數(shù)據(jù)聯(lián)合反演的研討會(huì)，并且在其他很多專題報(bào)告中，約束聯(lián)合反演已成為一項(xiàng)提高分辨率的重要技術(shù)之一。Colombo等[11]通過(guò)建立密度、電阻率與速度的經(jīng)驗(yàn)關(guān)系式，通過(guò)重、電、地震聯(lián)合反演圈定了隱伏鹽丘； De Stefano等[12]則提出了一種多種地球物理數(shù)據(jù)同時(shí)聯(lián)合反演的新方法，提高了墨西哥灣地區(qū)地震偏移數(shù)據(jù)的質(zhì)量；Zhdanov等[13]提出了基于Gramian約束的廣義多種地球物理數(shù)據(jù)聯(lián)合方法； 陳華根等[14]研究了大地電磁(MT)與重力數(shù)據(jù)的模擬退火聯(lián)合反演方法，取得了較好的效果； Liu等[15]提出了基于井?dāng)?shù)據(jù)和地震解釋層位約束的三維重力反演方法，有效地提高了深部鹽巖頂部構(gòu)造的預(yù)測(cè)精度；李華東等[7]基于隨機(jī)分布共網(wǎng)格單元模型的建模技術(shù)，利用密度、電阻率、磁性參數(shù)和速度等多種物性參數(shù)建立模型，通過(guò)引入正則化思想，實(shí)現(xiàn)了共網(wǎng)格條件下的重力、磁力、電法和地震資料的同步聯(lián)合反演，較好地確定了物性界面及地質(zhì)結(jié)構(gòu)；陳曉等[10]提出了“多次建模，綜合約束，分步反演”的多方法聯(lián)合反演策略，并以大地電磁和地震數(shù)據(jù)聯(lián)合反演為例，探討了該框架在建立玄武巖模型中的適用性；潘豪杰等[16]提出在貝葉斯理論框架下基于非線性的簡(jiǎn)化三相方程和改進(jìn)Archie公式的聲波—電性數(shù)據(jù)聯(lián)合反演方法，進(jìn)行儲(chǔ)層參數(shù)預(yù)測(cè)，有效降低因敏感性和噪聲等問(wèn)題產(chǎn)生的不確定性。

本文基于井—震約束的大地電磁與重力數(shù)據(jù)的人工魚群聯(lián)合反演，利用已知的測(cè)井?dāng)?shù)據(jù)、地震剖面、地質(zhì)解釋剖面等先驗(yàn)信息約束建模，引入人工魚群的群智能反演算法，通過(guò)物性關(guān)系建立電阻率與密度之間聯(lián)系，利用MPI編程環(huán)境對(duì)反演方法進(jìn)行并行設(shè)計(jì)和編程實(shí)現(xiàn)，最后用模型數(shù)據(jù)和實(shí)測(cè)數(shù)據(jù)測(cè)試反演方法的正確性和可行性。

1 人工魚群約束反演

人工魚群算法是一種群智能最優(yōu)化算法。其基本思想是在一片水域中，魚生存數(shù)目最多的地方一般就是該水域中富含營(yíng)養(yǎng)物質(zhì)最多的地方，數(shù)學(xué)上根據(jù)這個(gè)特點(diǎn)模仿魚群的覓食行為、聚群行為和追尾行為，從而實(shí)現(xiàn)全局尋優(yōu)。因此，人工魚群算法主要包括以下四個(gè)步驟：模型選取、覓食行為、聚群行為、追尾行為，人工魚群的約束反演主要通過(guò)指定模型參數(shù)變化空間來(lái)實(shí)現(xiàn)[17]。

1.1 約束模型選取

假設(shè)人工魚群反演模型矢量M的約束范圍為[Mmin,Mmax]，其模型參數(shù)可以是電阻率、密度或者厚度。在模型空間中隨機(jī)產(chǎn)生初始模型，模擬人工魚群的隨機(jī)行為，并通過(guò)均勻概率分布產(chǎn)生第i條人工魚的初始值

(1)

1.2 覓食行為

(2)

(3)

1.3 聚群行為

(4)

1.4 追尾行為

(5)

式中u4是[0,1]的隨機(jī)數(shù)。

2 二維電磁—重力數(shù)據(jù)并行約束聯(lián)合反演

本文基于電磁與重磁數(shù)據(jù)之間的關(guān)系，并利用已知的測(cè)井?dāng)?shù)據(jù)、地震解釋剖面、地質(zhì)解釋剖面等先驗(yàn)信息約束建模，采用人工魚群智能反演算法和MPI并行設(shè)計(jì)，實(shí)現(xiàn)MT—重力數(shù)據(jù)的并行約束聯(lián)合反演。

二維MT—重力數(shù)據(jù)聯(lián)合反演的目標(biāo)函數(shù)為

(6)

式中：α和β分別為大地電磁和重力反演目標(biāo)函數(shù)的加權(quán)因子； 下標(biāo)“M”和“G”分別代表大地電磁和重力數(shù)據(jù)； 上標(biāo)“TE”和“TM”分別代表TE模式和TM模式； 上標(biāo)“cal”和“obs”分別代表計(jì)算值與觀測(cè)值；Ms和Ns分別表示二維MT測(cè)點(diǎn)數(shù)和重力測(cè)點(diǎn)數(shù)；Mf為二維MT測(cè)點(diǎn)的頻點(diǎn)數(shù)；ρ和φ分別代表視電阻率和相位；g代表重力值。

通過(guò)統(tǒng)計(jì)各個(gè)地層的密度、電阻率變化范圍，建立層密度—電阻率之間的物性關(guān)系，在迭代反演過(guò)程可相互轉(zhuǎn)換；分別計(jì)算大地電磁和重力的目標(biāo)函數(shù)以及各自的加權(quán)因子α和β，獲得聯(lián)合反演的擬合差；再不斷進(jìn)行人工魚群迭代反演，使目標(biāo)函數(shù)達(dá)到期望值。聯(lián)合反演過(guò)程中，要進(jìn)行多次二維大地電磁和二維重力數(shù)據(jù)的正演計(jì)算，筆者利用MPI并行設(shè)計(jì)方法，分別對(duì)二維大地電磁和二維重力數(shù)據(jù)的正演模塊進(jìn)行并行計(jì)算，以期提高反演效率。

3 模型數(shù)據(jù)和實(shí)測(cè)數(shù)據(jù)試算

3.1 理論模型

圖1所示為根據(jù)某實(shí)際地震解釋剖面設(shè)計(jì)的一個(gè)復(fù)雜二維電阻率—密度模型。反演目的層的電阻率為400Ω·m、剩余密度為-0.024g/cm3、埋深為3000～8000m。目的是檢驗(yàn)本文反演方法對(duì)埋藏深度大、深度變化范圍也很大的層位是否能夠有效識(shí)別。

二維MT和重力數(shù)據(jù)模擬的點(diǎn)距為200m，測(cè)點(diǎn)為99個(gè)，正演剖分網(wǎng)格均為115×111。MT數(shù)據(jù)采用二維有限差分進(jìn)行正演，重力數(shù)據(jù)采用經(jīng)典的二度體多邊形正演算法。

對(duì)正演MT和重力數(shù)據(jù)分別進(jìn)行二維單獨(dú)反演和二維并行聯(lián)合反演，以約束異常體的形態(tài)、反演異常體的電阻率和剩余密度。反演的模型空間變化范圍為真實(shí)模型的±30%，最大迭代次數(shù)為30，人工魚的數(shù)量為10，目標(biāo)函數(shù)擬合終止條件是目標(biāo)函數(shù)ΔE<1.0×10-6，MT和重力二維反演的模型剖分網(wǎng)格與正演相同，都為115×111。經(jīng)過(guò)30次迭代反演之后，程序都到達(dá)最大迭代次數(shù)而正常結(jié)束。

圖1 復(fù)雜二維電阻率(左)與剩余密度(右)模型

圖2為二維MT數(shù)據(jù)單獨(dú)反演及MT—重力數(shù)據(jù)聯(lián)合反演的電阻率剖面。由圖可見，兩種方法都較好地恢復(fù)了模型的真實(shí)值。圖3為二維MT數(shù)據(jù)單獨(dú)反演及MT—重力數(shù)據(jù)聯(lián)合反演的異常體電阻率分布柱狀圖。如果反演的電阻率異常值柱狀分布越集中于實(shí)際值，說(shuō)明反演的效果越好，因此通過(guò)圖3右可進(jìn)一步分析、評(píng)價(jià)反演效果。從圖3可以看出，聯(lián)合反演的電阻率柱狀圖的中心更為集中，并且中心值對(duì)應(yīng)的最大數(shù)量大于300，而單獨(dú)反演的最大數(shù)量約為200，說(shuō)明聯(lián)合反演結(jié)果更接近于真實(shí)值，反演效果亦更好。圖4為二維重力數(shù)據(jù)單獨(dú)反演及MT—重力數(shù)據(jù)聯(lián)合反演的剩余密度剖面，圖5是與之對(duì)應(yīng)的剩余密度分布柱狀圖。從圖5可見，相對(duì)于較單獨(dú)用重力數(shù)據(jù)進(jìn)行反演，MT—重力數(shù)據(jù)聯(lián)合反演的異常體剩余密度值分布中心更集中。

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

實(shí)測(cè)數(shù)據(jù)來(lái)自川中高石梯—龍女寺地區(qū)的大地電磁—重力勘探項(xiàng)目。以往的鉆井及地震資料表明，川中地區(qū)震旦系油氣的分布與深部裂谷的發(fā)育有一定關(guān)系[8]。為了進(jìn)一步了解川中高石梯—龍女寺地區(qū)震旦系裂谷發(fā)育情況，在該區(qū)部署了重力和MT測(cè)線，重力測(cè)點(diǎn)與MT測(cè)點(diǎn)重合，圖6所示為測(cè)線分布圖，本文以L1測(cè)線為例分析本文方法的效果。L1測(cè)線共有299個(gè)測(cè)點(diǎn)，點(diǎn)距約為500m，長(zhǎng)度為161km。

地層電阻率統(tǒng)計(jì)以區(qū)內(nèi)及工區(qū)周邊電阻率數(shù)據(jù)均較全的7口鉆井資料為主。根據(jù)電阻率統(tǒng)計(jì)數(shù)據(jù)，研究區(qū)上三疊統(tǒng)以上地層呈低阻特征；下三疊統(tǒng)嘉陵江組—中三疊統(tǒng)雷口坡組為一套次高阻層；下三疊統(tǒng)飛仙關(guān)組—中二疊統(tǒng)呈低阻特征；下二疊統(tǒng)地層呈高阻特征；寒武系—奧陶系地層呈中—低阻特征；上震旦統(tǒng)呈高阻特征。統(tǒng)計(jì)川中地區(qū)各地層的電阻率和密度資料，發(fā)現(xiàn)研究區(qū)地層的密度與電阻率比較穩(wěn)定，有一定的相關(guān)性。通過(guò)分層建立電阻率與密度的關(guān)系式，為MT與重力的聯(lián)合反演提供物性關(guān)系的紐帶。

圖2 二維MT數(shù)據(jù)單獨(dú)反演(左)及MT—重力數(shù)據(jù)聯(lián)合反演(右)的電阻率剖面

圖3 二維MT數(shù)據(jù)單獨(dú)反演(左)及MT—重力數(shù)據(jù)聯(lián)合反演(右)的目的層電阻率柱狀圖

圖4 二維重力數(shù)據(jù)單獨(dú)反演(左)及MT—重力數(shù)據(jù)聯(lián)合反演(右)的剩余密度剖面

圖5 二維重力數(shù)據(jù)單獨(dú)反演(左)及MT—重力數(shù)據(jù)聯(lián)合反演(右)的目的層剩余密度柱狀圖

圖6 研究區(qū)測(cè)線位置圖

對(duì)測(cè)線L1分別采用二維MT反演、二維重力反演和二維MT—重力聯(lián)合反演。MT和重力反演的剖分網(wǎng)格大小都為318×121，MT數(shù)據(jù)參與反演的頻點(diǎn)為19個(gè)，頻率范圍為0.001～320Hz。因?yàn)槌Ｒ?guī)無(wú)旋轉(zhuǎn)的TM數(shù)據(jù)反演結(jié)果與已知測(cè)井信息更加吻合，因此選擇TM模式數(shù)據(jù)進(jìn)行聯(lián)合反演。利用測(cè)井和地震解釋剖面進(jìn)行建模約束，根據(jù)統(tǒng)計(jì)的地層電阻率和密度變化范圍，確定各層的物性反演初值及變化范圍； 聯(lián)合反演時(shí)，采用統(tǒng)計(jì)的地層電阻率—密度關(guān)系式。反演最大迭代次數(shù)為30， 人工魚的數(shù)量為10， 目標(biāo)函數(shù)擬合終止條件是ΔE<1.0×10-4，都采用10個(gè)CPU并行計(jì)算。經(jīng)過(guò)30次迭代反演，結(jié)果都到達(dá)最大迭代次數(shù)而正常結(jié)束。

圖7為L(zhǎng)1線常規(guī)二維MT反演、井震約束的二維MT—重力資料聯(lián)合反演電阻率和密度剖面。距離該測(cè)線3km的東北方向有一鉆井M1，其電阻率測(cè)井曲線如圖7所示。由圖7a可見，反演電阻率與M1井電阻率測(cè)井曲線一致性較好。基于該電阻率剖面，結(jié)合電測(cè)井、地震解釋及物性資料，對(duì)地質(zhì)層位以及深部裂谷進(jìn)行了解釋。常規(guī)的二維反演電阻率在幅值上與電測(cè)井有較大差異，且對(duì)深部薄層的分辨能力較差，因此需要聯(lián)合地震和鉆井資料，進(jìn)行約束下的單獨(dú)反演和聯(lián)合反演，以提高重磁電勘探的分辨率。

圖7b為井震約束的二維MT與重力數(shù)據(jù)聯(lián)合反演的電阻率剖面。由圖可見，通過(guò)建模約束反演之后，電性層分布規(guī)律與電測(cè)井一致性很強(qiáng)，電阻率剖面對(duì)層位的刻畫更加精細(xì)。在海拔-1km和-4km左右的兩套高阻層連續(xù)性更好，與下伏地層的界面更清晰，裂谷形態(tài)在局部有微小的變化，電阻率也低于上覆和下伏地層，表明聯(lián)合反演結(jié)果對(duì)深部地層電性刻畫更精細(xì)。與MT數(shù)據(jù)的反演結(jié)果(圖7a)類似，聯(lián)合反演的密度剖面(圖7c)上有四套高密度層，海拔-1km和-4km左右的兩套高密度層橫向連續(xù)性更好； 海拔-6km以下裂谷的密度分布不均，密度為2.71～2.74g/cm3，符合裂谷的相對(duì)低密特征。

圖7 L1線反演剖面

圖8為實(shí)測(cè)TM模式視電阻率剖面與聯(lián)合反演擬合計(jì)算的TM視電阻率斷面圖?？梢钥闯?，總體上來(lái)說(shuō)，聯(lián)合反演剖面較好地?cái)M合了實(shí)測(cè)數(shù)據(jù)，僅局部區(qū)域擬合效果還有待提高，比如剖面上10～20km段的中低頻段(圖中橢圓所示)，一致性較差，還需進(jìn)行精細(xì)的地質(zhì)建模，調(diào)整約束參數(shù)范圍。圖9為聯(lián)合反演擬合計(jì)算與實(shí)測(cè)重力異常對(duì)比曲線，可見除了在5km～40km之間擬合略微欠缺，其他地方擬合程度較好。圖9為聯(lián)合擬合差隨反演迭代次數(shù)變化的對(duì)比圖，聯(lián)合反演的數(shù)據(jù)擬合差在30次迭代后達(dá)到0.0204，聯(lián)合反演的擬合程度較好，說(shuō)明反演結(jié)果比較可靠。

2016年在研究區(qū)內(nèi)新鉆Y1井(圖6)，該井日產(chǎn)115萬(wàn)方工業(yè)氣流。Y1井臨近L1測(cè)線的裂谷翼端，臨近南華紀(jì)裂谷內(nèi)斷陷區(qū)，是油氣的有利指向區(qū)。根據(jù)前述反演結(jié)果并結(jié)合測(cè)井和地震資料，綜合推測(cè)此處為一級(jí)有利區(qū)，Y1井的試氣結(jié)果證實(shí)了反演結(jié)果的可靠性。

圖8 L1線實(shí)測(cè)TM模式視電阻率剖面(a)和聯(lián)合反演擬合TM模式視電阻率剖面(b)

圖9 實(shí)測(cè)及聯(lián)合反演擬合的重力異常曲線對(duì)比(左)及擬合差曲線(右)

4 結(jié)論

本文提出了基于人工魚群算法的非線性MT與重力數(shù)據(jù)的并行約束反演方法。理論模型計(jì)算結(jié)果證明了該方法的有效性；通過(guò)對(duì)實(shí)測(cè)資料應(yīng)用本文方法進(jìn)行反演，發(fā)現(xiàn)了川中深層裂谷及南華紀(jì)有利目標(biāo)，解釋結(jié)論得到后續(xù)鉆井的證實(shí)。

(1)根據(jù)測(cè)井、巖心和野外實(shí)測(cè)資料分析了研究區(qū)的地層電阻率和密度特征，并提出分層建立該工區(qū)電阻率與密度的關(guān)系式，為大地電磁與重力數(shù)據(jù)的聯(lián)合反演提供物性關(guān)系的紐帶。

(2)對(duì)研究區(qū)的兩條測(cè)線分別進(jìn)行二維大地電磁反演和二維大地電磁—重力數(shù)據(jù)聯(lián)合反演。結(jié)果表明，通過(guò)井震建模約束反演，電性層分布規(guī)律與電測(cè)井信息一致性很好，二維大地電磁與重力數(shù)據(jù)的聯(lián)合反演對(duì)層位的刻畫比大地電磁單獨(dú)反演的結(jié)果更加精細(xì)，層位更加連續(xù)。

(3)實(shí)際數(shù)據(jù)處理結(jié)果證實(shí)了本文提出的非線性約束聯(lián)合反演有很強(qiáng)的實(shí)用性。