微震震源定位算法優化研究及應用

趙德明，龐 銳，王海波

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

0 引 言

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

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

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

1 SSA算法描述

震源掃描疊加方法(SSA)是地震學中進行微小能量地震震源定位的方法。水力壓裂微地震監測信號具有與其相類似的特點，因此，可將震源掃描疊加方法應用于微地震震源定位的實際應用中[10]。

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

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

圖1 震源掃描原理

假設一個微地震信號由N個檢波器接收到(圖中的1、2、3點)，第一步分別對每個檢波器記錄rn歸一化，接著計算某個點(c)某個時刻(t0)的“亮度”函數。

(1)

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

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

(2)

式中，N表示檢波器個數，W表示時間窗口范圍，Δt表示信號采集間隔，t0和tcn同上，Un表示歸一化記錄結果，Q表示理論到達時刻與檢測到達時刻差值確定的權重值。

2 并行解決方案

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

圖2 MapReduce執行過程

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

(2)只要有Map運算任務，工作單元會獲得對應的key/value組，并把其作為數據的關鍵字存儲在內存中；

(3)程序會定期將內存中的key/value組寫入本地存儲空間，程序員定義的劃分函數將其細分為不同分組，任務調度器會根據Reduce任務執行情況將本地存儲空間key/value組的物理地址傳給Reduce任務工作單元；

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

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

(6)Reduce工作單元將函數運算結果發送到指定文件中。當所有的Map、Reduce任務都接受后，任務調度器將結果和代碼返回到主程序[14]。

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

3 微地震震源掃描算法優化

為了達到實時處理，該項目針對微地震震源掃描算法進行了優化。現有的微地震震源掃描算法的掃描網格是固定的，該文在已有的定網格震源掃描定位處理的基礎上，進行網格插值，將其大范圍粗網格掃描范圍通過網格插值的方法插值實現小區域的細網格進行掃描，并對細網格進行并行化計算，最終實現粗、細雙重網格并行化掃描。其效果上，數據計算量大幅度減少，有效提高了計算效率[15]。

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

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

假定目標函數在初始搜索區間[a,b]中有三點α(1)、α(2)、和α(3)(a≤α(1)<α(2)<α(3)≤b)，其函數值分別為f1、f2和f3(見圖3)，且滿足f1>f2，f2>f3，即滿足函數值為兩頭大中間小的性質。利用這三點及相應的函數值作一條二次曲線，其函數p(α)為一個二次多項式。

圖3 二次插值示意圖

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

(3)

式中，a0、a1、a2為待定系數。

根據插值條件，插值函數p(α)與原函數f(α)在插值節點P1、P2、P3處函數值相等，得：

(4)

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

(5)

解式(5)即求得插值函數的極小點。

(6)

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

a1=

(7)

a2=

(8)

(9)

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

(10)

式中：

(11)

(12)

圖4 變網格震源掃描程序流程

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

該文采用Qt的run()多線程機制和map()機制相結合的并行方法對變網格震源掃描方法并行化實現，具體步驟如下：

圖5 變網格微震震源掃描并行化流程

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

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

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

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

4 實際資料測試與分析

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

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

為了驗證文中優化方法的性能，分別用優化前和優化后的震源掃描算法對不同數據量的拾取剖面進行測試，并計算加速比，結果如表1所示。

表1 實驗結果

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

通常一段頁巖氣微地震壓裂監測持續5～6個小時，而可以監測到的微地震事件大約可以達到五百個，也就是五百個拾取事件的剖面。從表1可以估算出在對震源掃描算法優化之前，一段壓裂監測數據的定位時間需要接近兩個小時，優化后執行時間卻可以控制在半個小時以內，完全達到現場處理的要求，可以對微地震壓裂方案提供實時數據指導。

5 結束語

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