趙德明，龐 銳，王海波

(中國(guó)石油化工股份有限公司石油物探技術(shù)研究院，江蘇 南京 211103)

0 引 言

微地震監(jiān)測(cè)是指利用井下或者地表布置的一個(gè)個(gè)監(jiān)測(cè)檢波器接受非常規(guī)油氣生產(chǎn)過程中的微小的破裂或錯(cuò)斷信號(hào)，并通過一系列處理方法反推出信號(hào)發(fā)生位置、延伸方向等信息，最后根據(jù)這些信息反過來優(yōu)化非常規(guī)壓裂設(shè)計(jì)方案，降低成本提高出產(chǎn)率[1]。隨著中國(guó)非常規(guī)油氣儲(chǔ)量不斷探明，豐富的儲(chǔ)量對(duì)微地震監(jiān)測(cè)的需求越來越強(qiáng)烈，而監(jiān)測(cè)的重要目的之一就是精確定位微地震震源。震源位置反演就是利用檢波器收集的地震波信息確定事件發(fā)生位置，即確定壓裂導(dǎo)致的裂點(diǎn)坐標(biāo)，進(jìn)而確定壓裂形成的裂縫延伸到角度和長(zhǎng)度[2]。所以，如何更好更快速地反演出微地震的震源位置，是微地震監(jiān)測(cè)的關(guān)鍵問題之一。

很多學(xué)者對(duì)微地震震源位置反演方法進(jìn)行了大量有效的研究。均勻介質(zhì)微地震位置反演方法主要有縱橫波時(shí)差法和同型波時(shí)差法。非均勻介質(zhì)反演方法主要有蒙特卡洛法、模擬退火法、遺傳算法及牛頓法、共軛梯度法等[3]。為了解決更復(fù)雜的反演問題，有的還提出了基于正演模型的反演方法、多個(gè)參數(shù)聯(lián)合反演[4]等。文獻(xiàn)[5]針對(duì)微地震信號(hào)信噪比低、淺地表干擾多、不易拾取初值等問題，在傳統(tǒng)的線性及非線性反演方法都不能給出可靠結(jié)果的情況下提出了貝葉斯反演定位方法。Kao和Shan提出了震源掃描算法(source scanning algorithm，SSA)，該算法不需要預(yù)先知道斷層信息，就能利用波形數(shù)據(jù)系統(tǒng)地掃描設(shè)定的時(shí)間及位置是否有震源以發(fā)現(xiàn)整個(gè)震源序列分布，對(duì)波形速度模型的抗干擾能力強(qiáng)，同時(shí)由于掃描過程不需要人為干預(yù)，能保證結(jié)果的客觀性和準(zhǔn)確性[6]。文獻(xiàn)[7]利用面波的特點(diǎn)對(duì)SSA的分辨率進(jìn)行了改進(jìn)。文獻(xiàn)[8-9]對(duì)SSA算法進(jìn)行了詳細(xì)闡述，并用實(shí)際工區(qū)資料驗(yàn)證了該算法對(duì)微小地震事件能很好的定位。

但隨著勘探的深入，遇到的問題越來越復(fù)雜，實(shí)際生產(chǎn)應(yīng)用的數(shù)據(jù)規(guī)模也隨之增加。為了不斷適應(yīng)這種變化，算法的復(fù)雜性和計(jì)算量呈指數(shù)性增加。在前面的文獻(xiàn)中不難看出，大部分研究側(cè)重于定位精度的改進(jìn)，但實(shí)際生產(chǎn)過程中往往需要實(shí)時(shí)給出壓裂破裂點(diǎn)及延伸方向。由于SSA算法的優(yōu)越性，這里以該定位算法為例進(jìn)行優(yōu)化研究。為了能快速地給出定位結(jié)果，通過研究掃描算法特點(diǎn)和并行解決方案并結(jié)合現(xiàn)場(chǎng)處理?xiàng)l件限制，采用并行化變網(wǎng)格方法加快震源定位速度。

1 SSA算法描述

震源掃描疊加方法(SSA)是地震學(xué)中進(jìn)行微小能量地震震源定位的方法。水力壓裂微地震監(jiān)測(cè)信號(hào)具有與其相類似的特點(diǎn)，因此，可將震源掃描疊加方法應(yīng)用于微地震震源定位的實(shí)際應(yīng)用中[10]。

震源掃描算法是一種可以處理分布的震源信號(hào)進(jìn)而成像定位的方法，這種新方法不用預(yù)先獲得斷層的產(chǎn)狀和規(guī)模，可以根據(jù)相對(duì)振幅和初至等波形數(shù)據(jù)，判斷出某個(gè)知道的時(shí)間和坐標(biāo)有沒有發(fā)生地震震源。它會(huì)利用序列化掃過某個(gè)試探性的位置范圍和起震時(shí)刻，發(fā)現(xiàn)整個(gè)震源排序分布而不需要拾取初至，也不用計(jì)算理論地震圖[11]。

SSA算法是由Kao和Shan在2004年提出的，他們用該方法對(duì)2003年3月北卡斯卡迪亞俯沖帶的慢地震資料處理后得到了意外的結(jié)果。在該資料一段時(shí)間里發(fā)生了幾個(gè)微地震信號(hào)，信號(hào)延續(xù)范圍比一般地震事件長(zhǎng)很多，其振幅水平卻又與噪聲一致，傳統(tǒng)的定位方法無法得出同一震源的位置，由于不能找到準(zhǔn)確的震源，信號(hào)波形反演也無效。最后利用SSA算法的“亮度”函數(shù)處理，準(zhǔn)確得出了微地震事件的震源坐標(biāo)和發(fā)震時(shí)刻[7]，其基本原理如圖1所示。

圖1 震源掃描原理

假設(shè)一個(gè)微地震信號(hào)由N個(gè)檢波器接收到(圖中的1、2、3點(diǎn))，第一步分別對(duì)每個(gè)檢波器記錄rn歸一化，接著計(jì)算某個(gè)點(diǎn)(c)某個(gè)時(shí)刻(t0)的“亮度”函數(shù)。

(1)

式中，N表示檢波器個(gè)數(shù)，t0表示開始檢測(cè)的時(shí)間，rn表示檢波器n的微地震信號(hào)，tcn表示某個(gè)信號(hào)由空間c傳到檢波器n經(jīng)過的時(shí)間。地層中的某一點(diǎn)c在t0時(shí)刻的“亮度”能通過對(duì)應(yīng)檢波器理論到達(dá)時(shí)刻(即t0加上走時(shí)t1c、t2c、t3c)推算得出，如果該處“亮度”較大(空心☆)，表示此處有事件，如果“亮度”較小(實(shí)心★)，說明此處沒有事件。假設(shè)監(jiān)測(cè)的全部振幅在空間c處和t0時(shí)間產(chǎn)生，li(c,t0)=1。如果li(c,t0)=0.5，表示在這點(diǎn)有事件的概率為50%。一般來講，li(c,t0)等于1的概率很低，所以就以li函數(shù)值大小為衡量標(biāo)準(zhǔn)確定微地震事件位置。

但在生產(chǎn)應(yīng)用過程中，計(jì)算走時(shí)的實(shí)際速度模型與理論速度模型有可能存在一些偏差，進(jìn)而影響結(jié)果的準(zhǔn)確性。公式中tcn僅考慮到直達(dá)波的理想情況，沒有考慮信號(hào)傳播過程中帶來的反射、折射等一系列噪音信號(hào)的干擾。因此Kao等人于2007年對(duì)上述算法進(jìn)行了改進(jìn)，提高了算法的抗噪能力和準(zhǔn)確性[2]。首先，在采樣點(diǎn)上下開一個(gè)時(shí)窗，用時(shí)窗內(nèi)的所有點(diǎn)的振幅替換該時(shí)刻的振幅。其次，添加權(quán)重值Q可以隨速度模型準(zhǔn)備與否靈活調(diào)整，改進(jìn)后公式如下：

(2)

式中，N表示檢波器個(gè)數(shù)，W表示時(shí)間窗口范圍，Δt表示信號(hào)采集間隔，t0和tcn同上，Un表示歸一化記錄結(jié)果，Q表示理論到達(dá)時(shí)刻與檢測(cè)到達(dá)時(shí)刻差值確定的權(quán)重值。

2 并行解決方案

目前比較常用的編程方法有MPI(message passing interface)消息傳遞方式、基于集群的OpenMP多線程方式、基于CUDA的異構(gòu)編程方式及混合異構(gòu)方式。由于實(shí)際應(yīng)用的復(fù)雜性，又產(chǎn)生了多級(jí)聯(lián)合編譯方式：MPI和OpenMP混合編譯方式，MPI和OpenMP及CUDA混合編譯方式[12]。另外，基于硬件GPU的多種編程方式也有廣泛應(yīng)用。按數(shù)據(jù)消息通信方式，以上所有方式大體可以分為數(shù)據(jù)并行和消息傳遞兩種編程方式[10]。數(shù)據(jù)并行方式編程相對(duì)簡(jiǎn)單，但對(duì)數(shù)據(jù)可并行化程度有很高的要求，適用范圍比較小。消息傳遞方式編程對(duì)數(shù)據(jù)要求不高，但需要考慮每一步運(yùn)算的數(shù)據(jù)間通信，容易產(chǎn)生數(shù)據(jù)不同步、死鎖等問題，編程太過復(fù)雜。為了能在編程時(shí)兼具以上兩種方式的優(yōu)點(diǎn)，谷歌公司于2004年提出了MapReduce[11]編程模型。這種方式能讓編程人員不需要考慮復(fù)雜的底層庫實(shí)現(xiàn)，只需要利用其提供的Map、Reduce和Partition等函數(shù)輕松并發(fā)處理海量數(shù)據(jù)[13]。利用MapReduce可以實(shí)現(xiàn)自動(dòng)并發(fā)處理數(shù)據(jù)，可以很方便地利用其提供的接口實(shí)現(xiàn)任務(wù)調(diào)度、負(fù)載均衡、容錯(cuò)和一致性管理等[12]，執(zhí)行過程如圖2所示。

圖2 MapReduce執(zhí)行過程

(1)啟動(dòng)MapReduce運(yùn)行時(shí)庫后，程序員可以手動(dòng)或根據(jù)可用工作單元負(fù)載情況自動(dòng)分割數(shù)據(jù)，并給每份數(shù)據(jù)分配一個(gè)key/value組，同時(shí)向即將運(yùn)算的工作單元復(fù)制一份并行程序；

(2)只要有Map運(yùn)算任務(wù)，工作單元會(huì)獲得對(duì)應(yīng)的key/value組，并把其作為數(shù)據(jù)的關(guān)鍵字存儲(chǔ)在內(nèi)存中；

(3)程序會(huì)定期將內(nèi)存中的key/value組寫入本地存儲(chǔ)空間，程序員定義的劃分函數(shù)將其細(xì)分為不同分組，任務(wù)調(diào)度器會(huì)根據(jù)Reduce任務(wù)執(zhí)行情況將本地存儲(chǔ)空間key/value組的物理地址傳給Reduce任務(wù)工作單元；

(4)通過RPC(remote procedure call)，Reduce任務(wù)執(zhí)行單元讀取上面獲取的Map執(zhí)行單元的key/value組；

(5)當(dāng)Reduce工作單元獲得所有需要的key/value組后，將具有相同key的key/value組排序形成key/values組，作為Reduce函數(shù)的輸入；

(6)Reduce工作單元將函數(shù)運(yùn)算結(jié)果發(fā)送到指定文件中。當(dāng)所有的Map、Reduce任務(wù)都接受后，任務(wù)調(diào)度器將結(jié)果和代碼返回到主程序[14]。

QT編程模型中QtConcurrent機(jī)制提供了MapReduce的編程接口，簡(jiǎn)單編程就可以將任務(wù)分發(fā)到處理器所有的核上。QtConcurrent命名空間里的run()函數(shù)可以實(shí)現(xiàn)多線程的并發(fā)執(zhí)行，QtConcurrent提供了map()函數(shù)可以自動(dòng)分配硬件線程去執(zhí)行各個(gè)任務(wù),直到所有任務(wù)都執(zhí)行完。MapReduce和多線程執(zhí)行代碼完全被隱藏在QtConcurrent框架下，所以可以直接使用而不用考慮細(xì)節(jié)[14]。

3 微地震震源掃描算法優(yōu)化

為了達(dá)到實(shí)時(shí)處理，該項(xiàng)目針對(duì)微地震震源掃描算法進(jìn)行了優(yōu)化?，F(xiàn)有的微地震震源掃描算法的掃描網(wǎng)格是固定的，該文在已有的定網(wǎng)格震源掃描定位處理的基礎(chǔ)上，進(jìn)行網(wǎng)格插值，將其大范圍粗網(wǎng)格掃描范圍通過網(wǎng)格插值的方法插值實(shí)現(xiàn)小區(qū)域的細(xì)網(wǎng)格進(jìn)行掃描，并對(duì)細(xì)網(wǎng)格進(jìn)行并行化計(jì)算，最終實(shí)現(xiàn)粗、細(xì)雙重網(wǎng)格并行化掃描。其效果上，數(shù)據(jù)計(jì)算量大幅度減少，有效提高了計(jì)算效率[15]。

二次插值法亦是用于一元函數(shù)f(α)在確定的初始區(qū)間內(nèi)搜索極小點(diǎn)的一種方法。二次插值的基本原理：在求解一元函數(shù)f(α)的極小點(diǎn)時(shí)，常常利用一個(gè)比函數(shù)低次的插值多項(xiàng)式p(α)來逼近原目標(biāo)函數(shù)，然后求這個(gè)多項(xiàng)式的極小點(diǎn)(低次多項(xiàng)式的極小點(diǎn)比較容易計(jì)算)，并以此作為目標(biāo)函數(shù)f(α)的近似極小點(diǎn)。如果其近似的程度尚未達(dá)到所要求的精度時(shí)，可以反復(fù)使用此法，逐次擬合，直到滿足給定的精度時(shí)為止。

插值多項(xiàng)式p(α)如果為二次多項(xiàng)式則叫二次插值法，如果是三次多項(xiàng)式則叫三次插值法。文中采用的是二次插值法，這里以此為例進(jìn)行介紹。

假定目標(biāo)函數(shù)在初始搜索區(qū)間[a,b]中有三點(diǎn)α(1)、α(2)、和α(3)(a≤α(1)<α(2)<α(3)≤b)，其函數(shù)值分別為f1、f2和f3(見圖3)，且滿足f1>f2，f2>f3，即滿足函數(shù)值為兩頭大中間小的性質(zhì)。利用這三點(diǎn)及相應(yīng)的函數(shù)值作一條二次曲線，其函數(shù)p(α)為一個(gè)二次多項(xiàng)式。

圖3 二次插值示意圖

p(α)=a0+a1α+a2α2

(3)

式中，a0、a1、a2為待定系數(shù)。

根據(jù)插值條件，插值函數(shù)p(α)與原函數(shù)f(α)在插值節(jié)點(diǎn)P1、P2、P3處函數(shù)值相等，得：

(4)

p'(α)=a1+2a2α=0

(5)

解式(5)即求得插值函數(shù)的極小點(diǎn)。

(6)

式(6)中要確定的系數(shù)a1,a2可在方程組(2)中利用相鄰兩個(gè)方程消去a0而得:

a1=

(7)

a2=

(8)

(9)

為便于計(jì)算，可將式(9)改寫為：

(10)

式中：

(11)

(12)

圖4 變網(wǎng)格震源掃描程序流程

變網(wǎng)格震源掃描算法實(shí)現(xiàn)步驟為：(1)輸入微地震有效事件記錄；(2)對(duì)微地震記錄進(jìn)行預(yù)處理(壞道剔除、強(qiáng)噪音壓制等)；(3)進(jìn)行地表高程校正；(4)讀入速度模型，計(jì)算粗網(wǎng)格炮點(diǎn)在各個(gè)網(wǎng)格點(diǎn)上的旅行時(shí)文件；(5)進(jìn)行粗網(wǎng)格震源掃描定位計(jì)算每個(gè)炮點(diǎn)的亮點(diǎn)值；(6)求得粗網(wǎng)格的亮點(diǎn)值；(7)讀取亮點(diǎn)值所在粗網(wǎng)格點(diǎn)的周圍網(wǎng)格上的旅行時(shí)值；(8)通過插值算法，求得小范圍的細(xì)網(wǎng)格的旅行時(shí)文件；(9)通過震源定位算法進(jìn)行細(xì)網(wǎng)格掃描；(10)對(duì)整個(gè)微地震事件進(jìn)行分時(shí)窗局部數(shù)據(jù)疊加；(11)求得細(xì)網(wǎng)格范圍的掃描亮點(diǎn)值；(12)掃描程序結(jié)束。

該文采用Qt的run()多線程機(jī)制和map()機(jī)制相結(jié)合的并行方法對(duì)變網(wǎng)格震源掃描方法并行化實(shí)現(xiàn)，具體步驟如下：

圖5 變網(wǎng)格微震震源掃描并行化流程

(1)創(chuàng)建震源掃描類，在震源掃描類中創(chuàng)建網(wǎng)格參數(shù)獲取函數(shù)獲取震源定位處理流程中sx、sy、sz的最大最小值及震源掃描的視窗大小、旅行時(shí)文件、速度文件等參數(shù)和文件路徑。

(2)在震源掃描類中創(chuàng)建多線程run()函數(shù)，獲取校正量，讀文件頭參數(shù)并檢查，如果旅行時(shí)文件震源參數(shù)和掃描參數(shù)不符則返回。讀取旅行時(shí)文件數(shù)據(jù)到內(nèi)存中， 利用事先創(chuàng)建好的震源掃描任務(wù)類為每個(gè)點(diǎn)構(gòu)建一個(gè)邏輯任務(wù),把任務(wù)的數(shù)據(jù)粒度細(xì)化到一個(gè)點(diǎn),這樣可以保證震源掃描任務(wù)類的數(shù)據(jù)結(jié)構(gòu)盡量簡(jiǎn)單明了。如果微地震監(jiān)測(cè)數(shù)據(jù)的輔助道上沒有同相軸，則不進(jìn)行震源掃描。對(duì)本次要掃描的監(jiān)測(cè)數(shù)據(jù)疊加道進(jìn)行相關(guān)處理獲得震源掃描初始數(shù)據(jù)，將數(shù)據(jù)與上面構(gòu)建的邏輯任務(wù)進(jìn)行關(guān)聯(lián)。

(3)調(diào)用震源掃描類中的sourcescanxy()函數(shù)，利用MapReduce框架進(jìn)行并行化處理，把任務(wù)列表中的任務(wù)map出去，qtconcurrent::map()會(huì)自動(dòng)分配硬件線程去執(zhí)行各個(gè)任務(wù)，自動(dòng)調(diào)用震源掃描任務(wù)類任務(wù)類中的sourcescantaskmapfunction()函數(shù)，直到所有任務(wù)都執(zhí)行完畢。本來應(yīng)該等所有任務(wù)都結(jié)束了,在這里對(duì)所有任務(wù)的結(jié)果進(jìn)行歸約操作，但是如果每個(gè)任務(wù)都保存結(jié)果，則會(huì)占用很大內(nèi)存，因此把歸約提前，提前到每個(gè)任務(wù)中，在每個(gè)任務(wù)中得到結(jié)果后就立即進(jìn)行歸約操作，然后釋放結(jié)果內(nèi)存，可以減少內(nèi)存占用。

(4)每個(gè)工作單元的CPU核在接收到任務(wù)后執(zhí)行震源掃描任務(wù)類中的runcomputation()函數(shù)，調(diào)用震源定位算法moveoutstack()和snr_()，并立即執(zhí)行Reduce操作。為防止數(shù)據(jù)不一致，在獲取Reduce操作需要的數(shù)據(jù)時(shí)進(jìn)行加鎖。等所有的Map和Reduce操作執(zhí)行完畢，將定位到的微地震壓裂事件數(shù)據(jù)寫回到數(shù)據(jù)文件中。

4 實(shí)際資料測(cè)試與分析

在一個(gè)八核CPU的移動(dòng)工作站上對(duì)某頁巖氣實(shí)際微地震監(jiān)測(cè)資料進(jìn)行處理，震源掃描定位算法的輸入是有事件拾取的微地震剖面，這里每個(gè)剖面的時(shí)間跨度為三秒，如圖6所示。

圖6 拾取事件的微地震剖面

為了驗(yàn)證文中優(yōu)化方法的性能，分別用優(yōu)化前和優(yōu)化后的震源掃描算法對(duì)不同數(shù)據(jù)量的拾取剖面進(jìn)行測(cè)試，并計(jì)算加速比，結(jié)果如表1所示。

表1 實(shí)驗(yàn)結(jié)果

從表1可以看出，該文采用的優(yōu)化方法可以達(dá)到較好的加速比。在八核處理器的工作站上，加速比約等于7。同時(shí)，拾取剖面數(shù)量增加，加速比也略微增加，說明該優(yōu)化方法具有較好的穩(wěn)定性和實(shí)用性，能很好地處理微地震剖面數(shù)據(jù)。

通常一段頁巖氣微地震壓裂監(jiān)測(cè)持續(xù)5～6個(gè)小時(shí)，而可以監(jiān)測(cè)到的微地震事件大約可以達(dá)到五百個(gè)，也就是五百個(gè)拾取事件的剖面。從表1可以估算出在對(duì)震源掃描算法優(yōu)化之前，一段壓裂監(jiān)測(cè)數(shù)據(jù)的定位時(shí)間需要接近兩個(gè)小時(shí)，優(yōu)化后執(zhí)行時(shí)間卻可以控制在半個(gè)小時(shí)以內(nèi)，完全達(dá)到現(xiàn)場(chǎng)處理的要求，可以對(duì)微地震壓裂方案提供實(shí)時(shí)數(shù)據(jù)指導(dǎo)。

5 結(jié)束語

大多學(xué)者對(duì)微地震震源定位方法的研究集中在定位準(zhǔn)確性的改進(jìn)上，該文則考慮到現(xiàn)場(chǎng)處理實(shí)時(shí)性的需求，采用變網(wǎng)格方法和MapReduce并行編程模型，對(duì)定位方法中的震源掃描算法(SSA)進(jìn)行了優(yōu)化。實(shí)際資料測(cè)試證明，優(yōu)化后的震源掃描算法獲得了較好的加速比，完全符合現(xiàn)場(chǎng)實(shí)時(shí)處理的需求。今后的工作主要對(duì)影響加速比的因素進(jìn)行研究，如工作站的內(nèi)存大小、硬盤容量、CPU處理核數(shù)等。