基于泊松碟采樣的小生境遺傳算法剩余靜校正

周 強(qiáng) 曹琳昱 陳愛萍

(①中國石油集團(tuán)東方地球物理公司物探技術(shù)研究中心成都分中心，四川成都 610213； ②中國石油集團(tuán)東方地球物理公司研究院西南分院，四川成都 610213； ③中國石油集團(tuán)山地地震技術(shù)試驗(yàn)基地，四川成都 610213)

1 引言

山地復(fù)雜構(gòu)造是油氣資源的富集區(qū)之一。在其形成過程中，經(jīng)歷了多期構(gòu)造運(yùn)動，使地表形態(tài)和地下構(gòu)造復(fù)雜，靜校正問題嚴(yán)重。無論采用折射或?qū)游龇椒?gòu)建的近地表結(jié)構(gòu)，都只是對速度、厚度的一個近似解，不可避免存在剩余靜校正問題[1-4]。這將直接影響地震資料成像效果，從而影響對地下構(gòu)造認(rèn)識的精度和可信度。

目前主要的反射剩余靜校正方法有三種。第一種是高斯賽德爾分解方法，以能量最大確定延遲時，建立反射時間方程組，采用高斯賽德爾迭代求取剩余靜校正量，計(jì)算速度快?，F(xiàn)今主流商業(yè)軟件大多配備此類模塊，該類方法適用于絕大多數(shù)地震資料。但由于采用線性方法求解剩余靜校正這種典型的非線性問題[5，6]，在處理山地低信噪比資料時，多輪之后會出現(xiàn)剩余靜校正量不收斂的現(xiàn)象。第二種是能量迭代法，它對低信噪比資料有一定適應(yīng)性，但因同一物理點(diǎn)對應(yīng)的靜校正量須滿足隨機(jī)性假設(shè)，否則存在累計(jì)誤差，進(jìn)而影響計(jì)算結(jié)果。第三種是非線性方法，包括模擬退火、遺傳算法等。常采用不同搜索策略開展全局尋優(yōu)[7-11]，但存在一個通病就是運(yùn)算速度極慢，難以滿足實(shí)際數(shù)據(jù)處理中的高效要求。因此，亟待研發(fā)一種能全局尋優(yōu)、快速解決山地低信噪比資料剩余靜校正問題的、可工業(yè)化實(shí)施的非線性剩余靜校正方法。

本文基于常規(guī)遺傳算法，利用泊松碟采樣生成初始解； 在進(jìn)化過程中通過多線程構(gòu)建多種群小生境，采用計(jì)算多核并行搜索最優(yōu)剩余靜校正量； 最后應(yīng)用理論模型和實(shí)際資料從效果和效率兩方面驗(yàn)證了該方法的有效性。

2 方法技術(shù)

Wilson等[12，13]率先將遺傳算法應(yīng)用于求解地震信號剩余靜校正量。以疊加能量最大為準(zhǔn)則，對多個剩余靜校正量組合進(jìn)行交換、選擇、變異操作，最終優(yōu)選能量最大的靜校正量組合，有效地避免了線性方法易陷入局部極值的問題。

常規(guī)遺傳算法剩余靜校正在編代碼實(shí)現(xiàn)時，因C語言自帶的隨機(jī)函數(shù)生成初始種群的隨機(jī)性不夠，采用單一種群容易產(chǎn)生“早熟”現(xiàn)象； 另外，計(jì)算效率低也是其明顯缺陷。針對常規(guī)遺傳算法進(jìn)行剩余靜校正存在的弊端，本文設(shè)計(jì)出一種基于泊松碟采樣的小生境遺傳算法(圖1)，用于更準(zhǔn)確地求取地震資料剩余靜校正量。

圖1 本文方法概要

2.1 泊松碟采樣均勻初始解

隨機(jī)數(shù)通常以系統(tǒng)時間為種子生成，隨機(jī)剩余靜校正量基本在同一時間完成，使隨機(jī)性不夠，導(dǎo)致初始解空間不均勻，包含最優(yōu)解的可能性低，對后期的遺傳進(jìn)化依賴性強(qiáng)，進(jìn)而影響遺傳算法的全局尋優(yōu)，且搜索耗時長[14]。針對該問題，擬采用泊松碟采樣思路。

泊松碟均勻初始解空間的方法是每個剩余靜校正量隨機(jī)生成，任意兩個值的差值大于閾值，以該差值為半徑的圓內(nèi)不再生成其他值。具體做法如下。

統(tǒng)計(jì)疊前地震數(shù)據(jù)體中炮點(diǎn)總數(shù)(Ns)及檢波點(diǎn)總數(shù)(Nr)，結(jié)合用戶設(shè)定的靜校正量范圍[a，b]，并采用實(shí)數(shù)編碼，獲取初始解空間

ri,j=U[a,b]

i∈[1,N]j∈[1,Ns+Nr]

(1)

式中：ri，j為初始解，且i為個體序號，j為基因位置；N為一個種群的個體數(shù)；U為解空間；a、b分別為最小、最大剩余靜校正量。

為了增強(qiáng)初始解遍布解空間的均勻度，采用泊松碟采樣方法在全部解空間中抽取相應(yīng)的個體。具體方法為：首先計(jì)算解空間的規(guī)模，將解空間的最小范數(shù)解與最大范數(shù)解之間的距離作為解空間的規(guī)模參數(shù)

(2)

式中γ、κ分別為解空間中的最小范數(shù)和最大范數(shù)解。采用下列閾值

(3)

若該閾值不滿足，則重復(fù)計(jì)算；通過篩選迭代，最終獲得下列初始種群矩陣

(4)

2.2 種群內(nèi)遺傳進(jìn)化

種群內(nèi)采用“精英+賭輪”選擇、余下個體與精英兩兩交叉、小擾動變異策略進(jìn)行遺傳進(jìn)化，實(shí)現(xiàn)全局尋優(yōu)。

計(jì)算種群中每個個體對應(yīng)的目標(biāo)函數(shù)值，即為種群中個體的適應(yīng)度值。以剖面能量為例，目標(biāo)函數(shù)計(jì)算公式為

]2

si,ns,ri,nr∈R

(5)

式中：Ei為疊加剖面能量；s、r分別為炮點(diǎn)和檢波點(diǎn)的剩余靜校正量；d為疊前地震信號幅值；nr、ns分別為第i個個體對應(yīng)的炮點(diǎn)、檢波點(diǎn)號；i為個體序號；y為CDP號；h為炮檢距；t為時間。根據(jù)個體適應(yīng)度值，保留k個適應(yīng)度最高的個體作為精英，對余下個體計(jì)算選擇概率

(6)

式中：E為個體對應(yīng)的目標(biāo)函數(shù)；P為選擇概率；R為個體。之后選擇直接進(jìn)入子代的個體，生成隨機(jī)數(shù)i，若滿足

(7)

則選擇第i個父代個體進(jìn)入子代，其余個體依據(jù)

r∈U[0,1]

(8)

實(shí)現(xiàn)交換策略。若

r

(9)

則有

(10)

否則不交換；若交換，即由種群中最優(yōu)個體與其交換產(chǎn)生新個體。采取兩點(diǎn)交換策略，即隨機(jī)產(chǎn)生兩個交換位置，將兩個待交換個體中兩個交換位置之間的基因段進(jìn)行交換，產(chǎn)生子代新個體； 否則，父代個體直接進(jìn)入子代成為子代個體。

對交換后種群，保留k個適應(yīng)度最高的個體作為精英，其余個體依據(jù)變異概率Pm判斷是否進(jìn)行變異，即產(chǎn)生[0，1]間隨機(jī)數(shù)r。若r

(11)

式中：k為常數(shù)，可據(jù)實(shí)際設(shè)定；E0為初始目標(biāo)函數(shù)；Eiter為當(dāng)前代最優(yōu)解的目標(biāo)函數(shù)； 靜校正量范圍[a，b]滿足

[ri,j-c,ri,j+c]∩[a,b]

(12)

2.3 種群間聯(lián)姻

常規(guī)遺傳算法具有極大概率出現(xiàn)早熟現(xiàn)象，即種群中某單個優(yōu)秀個體的基因會漸漸“統(tǒng)治”其他基因，造成整個種群無法跳出該個體基因的“統(tǒng)治”而陷入局部極值，進(jìn)而導(dǎo)致算法無法得到最優(yōu)解。

為了增加基因段的多樣性，同時減少早熟現(xiàn)象的發(fā)生，采取多種群同時進(jìn)行遺傳操作。如圖2所示，在上一輪遺傳進(jìn)化之后，進(jìn)行聯(lián)姻操作。具體做法是將每個種群中的適應(yīng)度函數(shù)值最大的個體的大部分基因交換至另一種群適應(yīng)度函數(shù)值中等的個體，再開始下一輪自由進(jìn)化。該方式下若干不同特征的種群獨(dú)立進(jìn)行遺傳操作，獲得所有種群的一個最優(yōu)解后，在各種群間進(jìn)行通婚操作，使得全部種群都可分享優(yōu)良基因。

圖2 多種群小生境遺傳算法操作示意圖

為進(jìn)一步提高計(jì)算效率，調(diào)用omp_set_num_ threads函數(shù)設(shè)置進(jìn)程個數(shù)。該進(jìn)程個數(shù)取決于設(shè)置的小生境種群個數(shù)，通常小于CPU的核數(shù)，但一般要求不小于4，以確保種群良性發(fā)展。每個線程處理一個種群的遺傳操作，包括每個個體的基因選擇、交換機(jī)變異。采用多線程并行方式處理小生境遺傳，能充分利用計(jì)算機(jī)多核更快地尋優(yōu)。

如此迭代計(jì)算，直至達(dá)到預(yù)設(shè)的迭代次數(shù)或誤差門檻值。幾個種群中疊加能量大的個體，即為估算的炮點(diǎn)、檢波點(diǎn)剩余靜校正量。

3 理論模型測試

構(gòu)建一個“雙平層+單斜層”理論模型。其二維理論數(shù)據(jù)采用30Hz雷克子波合成，1ms采樣，采樣長度為1000ms； 觀測系統(tǒng)為1785-15-30-15-1785，120道接收，檢波點(diǎn)、炮點(diǎn)樁號范圍分別是1～318、61～259，滿覆蓋次數(shù)為30。

圖3a和圖3c分別是炮點(diǎn)、檢波點(diǎn)靜校正量：黃線為期望值，范圍是-20～20ms，有較強(qiáng)隨機(jī)性； 藍(lán)線是常規(guī)剩余靜校正量； 粉線是本文方法計(jì)算的小生境遺傳算法剩余靜校正量。三種顏色曲線整體上吻合較好，說明大部分剩余靜校正問題已得到解決。但從炮點(diǎn)、檢波點(diǎn)靜校正量誤差值(圖3b、圖3d)可明顯看出，藍(lán)線有個別異常值，說明常規(guī)剩余靜校正方法存在局部解，而本文提出的小生境遺傳算法剩余靜校正能全局尋優(yōu)，未陷入局部解。

從疊加剖面來看，常規(guī)剩余靜校正的局部解導(dǎo)致剖面上同相軸出現(xiàn)錯斷及虛假同相軸(圖4a)現(xiàn)象，而本文剩余靜校正方法除了小號附加段因覆蓋參數(shù)不足出現(xiàn)稍許彎曲外，其余剖面都恢復(fù)了構(gòu)造形態(tài)，且能量得到很好聚焦(圖4b)。

進(jìn)一步分析共中心點(diǎn)道集(CMP1281)。將加入了隨機(jī)擾動量的道集(圖5a)、常規(guī)剩余靜校正后道集(圖5b)與用本文方法的剩余靜校正后道集(圖5c)進(jìn)行對比，可見常規(guī)剩余靜校正的兩個局部解造成兩道地震道發(fā)生了較大時移，在疊加剖面上表現(xiàn)為同相軸的錯斷。

圖6a展示的是小生境遺傳算法和常規(guī)遺傳算法計(jì)算的疊加能量(縱軸)隨大循環(huán)(即小生境算法中的聯(lián)姻)次數(shù)收斂的情況，可見小生境遺傳算法最終的疊加能量更大，而常規(guī)算法收斂于某一能量值就不再上升，即所謂的“早熟”。

圖6b展示的是在第一次循環(huán)時，即還沒發(fā)生聯(lián)姻前，四個種群各自選擇、交叉、變異。其進(jìn)化過程各不相同，即疊加能量的具體變化不同，而最終的疊加能量是A種群最大，則A 種群中能量最大的精英基因?qū)⑻娲渌齻€種群能量處于中間值的個體的部分基因； 之后，四個種群進(jìn)入下一次大的進(jìn)化循環(huán)。

圖3 剩余靜校正量及其誤差

圖4 理論模型成像效果對比

圖5 共中心點(diǎn)道集

圖6 能量收斂圖

4 實(shí)際應(yīng)用及效果

新疆某山地地區(qū)斷裂發(fā)育，地層傾角較大，地下構(gòu)造復(fù)雜；近地表結(jié)構(gòu)復(fù)雜，縱橫向速度、厚度變化劇烈。無論采用初至折射還是層析折射靜校正方法，都只能取得對近地表的近似成像效果，不可避免地存在較嚴(yán)重的剩余靜校正問題。

基于此，在進(jìn)行多輪常規(guī)剩余靜校正和速度迭代更新后，成像品質(zhì)仍不能很好地滿足地質(zhì)任務(wù)要求。應(yīng)用本文提出的基于泊松碟采樣的小生境遺傳算法剩余靜校正方法處理后，成像效果(同相軸連續(xù)性、波組特征)有了明顯改善。

圖7a是三維地震M測線經(jīng)多輪常規(guī)剩余靜校正和速度迭代更新后的疊加剖面，存在的問題是構(gòu)造主體部位同相軸連續(xù)性不好，地震反射特征雜亂，波組特征不夠清晰。應(yīng)用本文方法后，剖面整體成像效果(圖7b)得到改善：同相軸變得連續(xù)；尤其是通過細(xì)節(jié)對比，構(gòu)造頂部波組特征變得清晰，兩翼與頂部的銜接自然連續(xù)(圖8)。

水平疊加時，剩余靜校正和速度更新都需反復(fù)迭代，耗費(fèi)大量時間。而本文方法實(shí)現(xiàn)了多線程編程，應(yīng)用計(jì)算機(jī)多核優(yōu)勢縮短了估算剩余靜校正量的時間。通過效率對比(表1)發(fā)現(xiàn)，達(dá)到能量增加率終止條件，常規(guī)遺傳算法需608min，而本文方法只需49min，其處理效率得到大幅度提高。

表1 效率對比表

圖7 應(yīng)用本方法前(a)、后(b)三維地震M測線成像效果對比

圖8 應(yīng)用本文方法前(a)、后(b)疊加剖面構(gòu)造頂部細(xì)節(jié)對比

5 結(jié)論與討論

本文提出的估算剩余靜校正量的方法，目標(biāo)函數(shù)為反射疊加能量，求解對象為炮點(diǎn)、檢波點(diǎn)剩余靜校正量； 在此基礎(chǔ)上，應(yīng)用泊松碟采樣對求解的靜校正量實(shí)數(shù)編碼，增加初始解遍布解空間的均勻度； 構(gòu)建多種群小生境，種群內(nèi)采用“精英+賭輪”選擇、余下個體與精英兩兩交叉、小擾動變異策略進(jìn)行遺傳進(jìn)化，種群間采用多線程同時進(jìn)化，選出幾個種群中適應(yīng)度最大個體與其他種群聯(lián)姻，增加基因片段的多樣性，避免“早熟”，提高尋優(yōu)的效率與精度。理論模型和實(shí)際數(shù)據(jù)測試表明，該方法可顯著且快速改善剖面成像效果，為復(fù)雜地區(qū)(尤其是低信噪比)資料提供了一種解決高頻靜校正問題的實(shí)用方法。

本文方法采用了多線程方式，效率提高了約10倍。但隨著對地震資料品質(zhì)要求的提高，高覆蓋次數(shù)及新老資料的連片處理使得地震數(shù)據(jù)量急劇增加，要使本文方法具有更強(qiáng)適應(yīng)性，還需基于地表一致性原則構(gòu)建MPI框架，以增強(qiáng)超大數(shù)據(jù)處理能力； 同時應(yīng)用多節(jié)點(diǎn)、多線程對剩余靜校正量進(jìn)行估算，進(jìn)一步提高計(jì)算速度及處理效率。