本發(fā)明屬于地震數(shù)據(jù)數(shù)值模擬和偏移成像領(lǐng)域，具體涉及一種基于GPU的波動方程反偏移方法，可用于地震數(shù)據(jù)正演模擬和最小二乘偏移數(shù)據(jù)重構(gòu)。
背景技術(shù)：
：：地震數(shù)據(jù)數(shù)值模擬是研究地震波在地下介質(zhì)中傳播特征的重要手段，也是地震成像(包括逆時(shí)偏移，最小二乘偏移、全波形反演等)實(shí)現(xiàn)過程中的一個(gè)重要步驟。基于偏移結(jié)果的反偏移是一種地震數(shù)據(jù)數(shù)值模擬方法，它是偏移成像的反過程，由偏移結(jié)果通過數(shù)學(xué)運(yùn)算得到疊前道集數(shù)據(jù)。目前常用的反偏移是基于kirchhoff積分算子的射線類反偏移方法(可參考：孫建國，Kirchhoff型真振幅偏移與反偏移[J]，勘探地球物理進(jìn)展，2002，25(6)：1-5；查樹貴，陳洪堤，龔小金，反偏移技術(shù)對疊前資料的噪聲消除[J]，中外能源，2006，11(3)：16-20；孫建國，Kirchhoff型高保真反偏移理論[J]，地球科學(xué)與環(huán)境學(xué)報(bào)，2011，33(2)：207-212，216；楊慶道，蔣兵，岳玉波，郭朝斌，基于Kirchhoff偏移/反偏移的隨機(jī)噪聲壓制方法[J]，石油物探，2013，52(5)：530-536)，可以較好地模擬一次波，但對復(fù)雜波場無法精確模擬?；陔p程波方程的逆時(shí)反偏移方法(可參考：YuZhang,LianDuanandYiXie.Astableandpracticalimplementationofleast-squaresreversetimemigration[J].ExpandedAbstractsof83thAnnualInternationalMeeting,SEG,2013,3715-3720；YuZhang,LianDuan.Predictingmultiplesusingareversetimedemigration[J].ExpandedAbstractsof82ndAnnualInternationalMeeting,SEG,2012,1-5)可以預(yù)測一次波、多次波，是一種精確的反偏移技術(shù)，然而這種方法計(jì)算量大，目前基于CPU(CentralProcessingUnit，中央處理器)的實(shí)現(xiàn)算法效率較低。近幾年，出現(xiàn)一種基于擴(kuò)展成像結(jié)果的波 動方程反偏移方法(可參考：HervéChauris1andMondherBenjemaa.Seismicwave-equationdemigration/migration[J].Geophysics,2010，75(3):S111-S119；WiktorWaldemarWeibullandArntsen.Reverse-timedemigrationusingtheextended-imagingcondition[J].Geophysics,2014,79(3):WA97-WA105)，依據(jù)爆炸反射面的思想，利用波動方程算子對每一個(gè)偏移距的成像剖面進(jìn)行反偏移。這種方法效率較快，但無法得到疊前炮集數(shù)據(jù)體，因而算不上嚴(yán)格意義的反偏移。上述現(xiàn)有反偏移方法技術(shù)在效果或效率方面無法滿足實(shí)際需求，較難得到應(yīng)用和推廣。技術(shù)實(shí)現(xiàn)要素：本發(fā)明針對現(xiàn)有反偏移方法存在的效果不佳或效率較低的問題，提供一種基于GPU(GraphicProcessingUnit，圖形處理器)的波動方程反偏移方法，輸入光滑速度模型和已有的偏移成像，逐炮按照炮點(diǎn)坐標(biāo)和接收點(diǎn)坐標(biāo)，利用時(shí)間空間域波動方程高階差分算子在GPU計(jì)算平臺對背景波場和散射波場進(jìn)行波場傳播外推計(jì)算，每個(gè)時(shí)刻在地表記錄散射波場，得到反偏移的疊前炮集數(shù)據(jù)體。本發(fā)明提供一種基于GPU的波動方程反偏移方法，其特征在于，包括以下步驟：輸入速度模型、偏移結(jié)果和控制參數(shù)；CPU獲得炮點(diǎn)位置和震源子波；CPU獲得接收點(diǎn)位置，確定炮數(shù)據(jù)范圍及波動方程計(jì)算網(wǎng)格，并從輸入的速度模型和偏移結(jié)果中獲取計(jì)算網(wǎng)格內(nèi)的速度c0和偏移值m；將CPU上的炮點(diǎn)位置、震源子波、接收點(diǎn)位置、計(jì)算網(wǎng)格內(nèi)的速度c0和偏移值m傳送到GPU；在GPU上進(jìn)行散射波場延拓并記錄反偏移波場；將GPU上記錄的反偏移波場傳送到CPU；以及CPU將記錄的反偏移波場寫入存儲器。進(jìn)一步地，在GPU上進(jìn)行散射波場延拓并記錄反偏移波場的步驟包括：在炮點(diǎn)處加載震源子波；利用計(jì)算網(wǎng)格內(nèi)的速度c0進(jìn)行背景波場u0的傳播；利用計(jì)算網(wǎng)格內(nèi)的速度c0和偏移值m進(jìn)行散射波場u1的傳播；在接收點(diǎn)處記錄該時(shí)刻的 散射波場u1的值，即為反偏移波場。進(jìn)一步地，利用計(jì)算網(wǎng)格內(nèi)的速度c0進(jìn)行背景波場u0的傳播，是基于輸入的速度模型和已知的震源位置，采用時(shí)間空間域雙程波方程高階有限差分算子對地震波場進(jìn)行正向延拓。進(jìn)一步地，利用計(jì)算網(wǎng)格內(nèi)的速度c0和偏移值m進(jìn)行散射波場u1的傳播，是基于輸入的速度模型和偏移結(jié)果，采用時(shí)間空間域雙程波方程高階有限差分算子對地震波場進(jìn)行正向延拓。進(jìn)一步地，時(shí)間空間域雙程波方程為：其中，t表示時(shí)間，c0表示計(jì)算網(wǎng)格內(nèi)的速度，m表示偏移值，u0表示背景波場，u1表示散射波場。優(yōu)選地，所述利用計(jì)算網(wǎng)格內(nèi)的速度c0進(jìn)行背景波場u0的傳播以及所述利用計(jì)算網(wǎng)格內(nèi)的速度c0和偏移值m進(jìn)行散射波場u1的傳播，是通過利用GPU并行技術(shù)進(jìn)行計(jì)算的。在時(shí)間空間域雙程波方程高階有限差分計(jì)算中，x、y，z三個(gè)方向上的網(wǎng)格數(shù)為Nx、Ny、Nz，則共需要對Nx×Ny×Nz個(gè)網(wǎng)格點(diǎn)進(jìn)行差分計(jì)算，每個(gè)網(wǎng)格點(diǎn)的計(jì)算是相互獨(dú)立的，由GPU獨(dú)立地同時(shí)進(jìn)行并行計(jì)算。GPU并行計(jì)算中采用每個(gè)線程進(jìn)行一個(gè)網(wǎng)格點(diǎn)的差分計(jì)算，采用多線程的并行計(jì)算技術(shù)。進(jìn)一步地，本發(fā)明的基于GPU的波動方程反偏移方法可以逐個(gè)炮點(diǎn)計(jì)算反偏移。進(jìn)一步地，本發(fā)明的基于GPU的波動方程反偏移方法可以逐個(gè)時(shí)刻進(jìn)行散射波場延拓并記錄反偏移波場。本發(fā)明針對現(xiàn)有反偏移技術(shù)效果不佳或效率較低的問題，采用雙程波方程進(jìn)行散射波場傳播獲得反偏移數(shù)據(jù)，可由偏移結(jié)果模擬一次波和多次波，是一種精確的反偏移方法。此外，為應(yīng)對雙程波方程高階差分計(jì)算效率低的缺點(diǎn)，采取基于GPU多線程的并行計(jì)算技術(shù)，大幅提高波動方程反偏移的計(jì)算速度。 本發(fā)明形成了高精度、高效率的反偏移方法，可用于地震數(shù)據(jù)正演模擬和最小二乘偏移數(shù)據(jù)重構(gòu)。本發(fā)明的其它特征和優(yōu)點(diǎn)將在隨后的具體實(shí)施方式部分予以詳細(xì)說明。附圖說明附圖是用來提供對本發(fā)明的進(jìn)一步理解，并且構(gòu)成說明書的一部分，與下面的具體實(shí)施方式一起用于解釋本發(fā)明，但并不構(gòu)成對本發(fā)明的限制。在附圖中：圖1是根據(jù)本發(fā)明一種實(shí)施方式的基于GPU的波動方程反偏移方法的步驟流程圖；圖2是根據(jù)本發(fā)明一種實(shí)施方式的基于GPU進(jìn)行散射波場延拓并記錄反偏移波場的步驟流程圖；圖3是在一個(gè)2D理論模型數(shù)據(jù)的反偏移計(jì)算中輸入的偏移結(jié)果；圖4是根據(jù)一個(gè)實(shí)施例的反偏移得到的一炮地震數(shù)據(jù)；圖5是根據(jù)一個(gè)實(shí)施例的原始炮集數(shù)據(jù)；圖6是根據(jù)一個(gè)實(shí)施例的伊拉克某工區(qū)一條偏移剖面；以及圖7是根據(jù)一個(gè)實(shí)施例的伊拉克某工區(qū)反偏移得到的一炮數(shù)據(jù)。具體實(shí)施方式以下結(jié)合附圖對本發(fā)明的具體實(shí)施方式進(jìn)行詳細(xì)說明。應(yīng)當(dāng)理解的是，此處所描述的具體實(shí)施方式僅用于說明和解釋本發(fā)明，并不用于限制本發(fā)明。根據(jù)本發(fā)明的一種實(shí)施方式，提供一種基于GPU的波動方程反偏移方法，其特征在于，包括以下步驟：輸入速度模型、偏移結(jié)果和控制參數(shù)；CPU獲得炮點(diǎn)位置和震源子波；CPU獲得接收點(diǎn)位置，確定炮數(shù)據(jù)范圍及波動方程計(jì)算網(wǎng)格，并從輸入的速度模型和偏移結(jié)果中獲取計(jì)算網(wǎng)格內(nèi)的速度c0和偏移值m；將CPU上的炮點(diǎn)位置、震源子波、接收點(diǎn)位置、計(jì)算網(wǎng)格內(nèi)的速度c0和偏移值m 傳送到GPU；在GPU上進(jìn)行散射波場延拓并記錄反偏移波場；將GPU上記錄的反偏移波場傳送到CPU；以及CPU將記錄的反偏移波場寫入存儲器。本實(shí)施方式的基于GPU的波動方程反偏移方法可以逐個(gè)炮點(diǎn)計(jì)算反偏移，具體地，可以在輸入速度模型、偏移結(jié)果和控制參數(shù)之前或者之后，判斷是否針對所有炮點(diǎn)計(jì)算完畢，如果計(jì)算完畢則方法流程結(jié)束，否則開始進(jìn)行下一個(gè)炮點(diǎn)的計(jì)算，計(jì)算完畢后CPU將記錄的反偏移波場寫入存儲器，再次進(jìn)行判斷，直至針對所有的炮點(diǎn)計(jì)算完畢。進(jìn)一步地，在GPU上進(jìn)行散射波場延拓并記錄反偏移波場的步驟包括：在炮點(diǎn)處加載震源子波；利用計(jì)算網(wǎng)格內(nèi)的速度c0進(jìn)行背景波場u0的傳播；利用計(jì)算網(wǎng)格內(nèi)的速度c0和偏移值m進(jìn)行散射波場u1的傳播；在接收點(diǎn)處記錄該時(shí)刻的散射波場u1的值，即為反偏移波場。利用計(jì)算網(wǎng)格內(nèi)的速度c0進(jìn)行背景波場u0的傳播，是基于輸入的速度模型和已知的震源位置，采用時(shí)間空間域雙程波方程高階有限差分算子對地震波場進(jìn)行正向延拓。利用計(jì)算網(wǎng)格內(nèi)的速度c0和偏移值m進(jìn)行散射波場u1的傳播，是基于輸入的速度模型和偏移結(jié)果，采用時(shí)間空間域雙程波方程高階有限差分算子對地震波場進(jìn)行正向延拓。本實(shí)施方式的基于GPU的波動方程反偏移方法可以逐個(gè)時(shí)刻進(jìn)行散射波場延拓并記錄反偏移波場。具體地，在開始進(jìn)行散射波場延拓時(shí)，首先判斷時(shí)間循環(huán)是否結(jié)束，如果時(shí)間循環(huán)結(jié)束，則散射波場延拓流程結(jié)束；如果時(shí)間循環(huán)沒有結(jié)束，則開始進(jìn)行下一時(shí)刻的散射波場延拓。針對一個(gè)時(shí)刻的散射波場延拓結(jié)束后，再次進(jìn)行時(shí)間循環(huán)是否結(jié)束的判斷，直至針對所有時(shí)刻的散射波場延拓結(jié)束。優(yōu)選地，時(shí)間空間域雙程波方程為：其中，t表示時(shí)間，c0表示計(jì)算網(wǎng)格內(nèi)的速度，m表示偏移值，u0表示背景波場，u1表示散射波場。優(yōu)選地，所述利用計(jì)算網(wǎng)格內(nèi)的速度c0進(jìn)行背景波場u0的傳播以及所述利用計(jì)算網(wǎng)格內(nèi)的速度c0和偏移值m進(jìn)行散射波場u1的傳播，是通過利用GPU并行 技術(shù)進(jìn)行計(jì)算的。在時(shí)間空間域雙程波方程高階有限差分計(jì)算中，x、y，z三個(gè)方向上的網(wǎng)格數(shù)為Nx、Ny、Nz，則共需要對Nx×Ny×Nz個(gè)網(wǎng)格點(diǎn)進(jìn)行差分計(jì)算，每個(gè)網(wǎng)格點(diǎn)的計(jì)算是相互獨(dú)立的，由GPU獨(dú)立地同時(shí)進(jìn)行并行計(jì)算。GPU并行計(jì)算中采用每個(gè)線程進(jìn)行一個(gè)網(wǎng)格點(diǎn)的差分計(jì)算，采用多線程的并行計(jì)算技術(shù)。下面參照附圖描述本發(fā)明的具體實(shí)施方式。圖1是根據(jù)本發(fā)明一種實(shí)施方式的基于GPU的波動方程反偏移方法的步驟流程圖。首先輸入速度模型、偏移結(jié)果及關(guān)鍵控制參數(shù)。反偏移是偏移的逆過程，偏移時(shí)輸入實(shí)際地震采集數(shù)據(jù)和由該數(shù)據(jù)分析得到的速度模型，通過運(yùn)算輸出偏移結(jié)果。反偏移是把速度模型、偏移結(jié)果，以及實(shí)際地震采集數(shù)據(jù)(利用炮點(diǎn)位置和接收點(diǎn)位置的信息)作為輸入，輸出模擬的地震采集數(shù)據(jù)。控制參數(shù)包括反偏移的炮數(shù)、反偏移子波頻率、反偏移計(jì)算步數(shù)、反偏移計(jì)算步長、速度模型x、y、z三個(gè)方向的網(wǎng)格數(shù)、速度模型x、y、z三個(gè)方向的網(wǎng)格間隔等。然后判斷是否針對所有炮點(diǎn)計(jì)算完畢，如果計(jì)算完畢則方法流程結(jié)束，否則開始進(jìn)行下一個(gè)炮點(diǎn)的計(jì)算。本實(shí)施例中逐炮進(jìn)行反偏移計(jì)算并保存反偏移結(jié)果，所有炮計(jì)算完畢則反偏移結(jié)束。接下來的流程是針對一個(gè)炮點(diǎn)的計(jì)算。1)CPU上獲得炮點(diǎn)位置和震源子波。炮點(diǎn)位置從實(shí)際地震采集數(shù)據(jù)中獲得，震源子波可以是人為給定的。2)CPU上獲得接收點(diǎn)位置，確定炮數(shù)據(jù)范圍及波動方程計(jì)算網(wǎng)格，并從輸入的速度模型和偏移結(jié)果中獲取計(jì)算網(wǎng)格內(nèi)的速度c0和偏移值m。接收點(diǎn)位置可從實(shí)際地震采集數(shù)據(jù)中獲得。炮數(shù)據(jù)范圍及波動方程計(jì)算網(wǎng)格是根據(jù)接收點(diǎn)位置確定的。通常輸入的偏移結(jié)果是整體的，一炮范圍對應(yīng)的計(jì)算網(wǎng)格是其中的一部分，由炮范圍的坐標(biāo)值可以從整體偏移結(jié)果中抽取該炮范圍內(nèi)的偏移結(jié)果。而輸入的速度模型也是整體的，一個(gè)炮范圍對應(yīng)的計(jì)算網(wǎng)格是其中的一部分，因此由炮范圍的坐標(biāo)值可以從整體速度模型中抽取該炮范圍內(nèi)的速度模型。3)將CPU上的炮點(diǎn)位置、震源子波、接收點(diǎn)位置和計(jì)算網(wǎng)格內(nèi)的速度c0和 偏移值m等數(shù)據(jù)傳送到GPU。4)在GPU上進(jìn)行散射波場延拓并記錄反偏移波場。5)將GPU上記錄的反偏移波場發(fā)送到CPU。6)CPU將反偏移波場寫入存儲器。優(yōu)選地，在GPU上進(jìn)行散射波場延拓并記錄反偏移波場的步驟4)，可參見附圖2。圖2是根據(jù)本發(fā)明一種實(shí)施方式的基于GPU進(jìn)行散射波場延拓并記錄反偏移波場的步驟流程圖。本實(shí)施方式中，逐個(gè)時(shí)刻進(jìn)行散射波場延拓和記錄偏移波場，每個(gè)時(shí)刻的實(shí)現(xiàn)方法是：41)在炮點(diǎn)處加載震源子波；42)利用速度c0進(jìn)行背景波場u0的傳播；43)利用速度c0和偏移值m進(jìn)行散射波場u1的傳播；44)在接收點(diǎn)處記錄該時(shí)刻的散射波場u1的值，即為反偏移波場。在利用速度c0進(jìn)行背景波場u0的傳播的步驟42)中，背景波場u0的正向傳播，是基于輸入的光滑速度模型和已知的震源位置，采用時(shí)間空間域雙程波方程高階有限差分算子對地震波場進(jìn)行正向延拓?？蛇x地，震源傳播雙程波方程和高階差分算子可參考“劉守偉、王華忠、陳生昌等，三維逆時(shí)偏移GPU/CPU機(jī)群實(shí)現(xiàn)方案研究[J]，地球物理學(xué)報(bào)，2013，56(10)：3487-3496”。其中震源位置可從輸入的地震采集數(shù)據(jù)中獲得。在利用速度c0和偏移值m進(jìn)行散射波場u1的傳播的步驟43)中，散射波場u1的正向傳播，是基于輸入的光滑速度模型和偏移結(jié)果，采用時(shí)間空間域雙程波方程高階有限差分算子對地震波場進(jìn)行正向延拓。上面所述的時(shí)間空間域雙程波方程為：1c02∂2u1∂2t-▿2u1=-m∂2u0∂2t]]>其中，t表示時(shí)間，c0表示速度，m表示偏移值，相當(dāng)于速度擾動，u0表示傳播的背景波場，u1表示散射波場。在接收點(diǎn)位置記錄每一個(gè)時(shí)刻的偏移波場， 得到的地震記錄即為一炮的波動方程反偏移結(jié)果。優(yōu)選地，上述步驟42)、43)是利用GPU并行技術(shù)進(jìn)行計(jì)算。在時(shí)間空間域雙程波方程高階有限差分計(jì)算中，假設(shè)x、y，z三個(gè)方向上的網(wǎng)格數(shù)為Nx、Ny、Nz，則共需要對Nx×Ny×Nz個(gè)點(diǎn)進(jìn)行差分計(jì)算，而每個(gè)點(diǎn)的計(jì)算是相互獨(dú)立的，由GPU獨(dú)立地同時(shí)進(jìn)行并行計(jì)算。GPU計(jì)算中，內(nèi)核函數(shù)是以線程網(wǎng)格(Grid)的形式組織的，每個(gè)線程網(wǎng)格由若干個(gè)線程塊(Block)組成，而每個(gè)線程塊又由多個(gè)線程(thread)組成。因此，采用每個(gè)線程進(jìn)行一個(gè)網(wǎng)格點(diǎn)的差分計(jì)算，可以實(shí)現(xiàn)并行快速運(yùn)算，提高了計(jì)算效率。下面參照圖3-圖7描述本發(fā)明的兩個(gè)具體實(shí)施例。實(shí)施例一：一個(gè)2D理論模型數(shù)據(jù)的反偏移計(jì)算圖3是輸入的偏移結(jié)果，圖4是通過反偏移得到的某炮集記錄，圖5是切除了直達(dá)波后的原始炮集記錄。對比圖4和圖5，兩者比較接近，由此說明根據(jù)本發(fā)明方法的反偏移可以很好地模擬一次波和多次波。實(shí)施例二：伊拉克某3D工區(qū)地震數(shù)據(jù)反偏移計(jì)算圖6是輸入的偏移結(jié)果中的一條剖面，由此利用本發(fā)明方法進(jìn)行反偏移計(jì)算，圖7是反偏移得到的一炮地震數(shù)據(jù)。反偏移計(jì)算網(wǎng)格為356×191×801，波場傳播時(shí)間步數(shù)為5000步，時(shí)間間隔為1ms，采用NVIDIA公司產(chǎn)品TeslaM2090型號的GPU卡進(jìn)行計(jì)算，一炮反偏移用時(shí)23分鐘。對于相同的計(jì)算量，CPU計(jì)算用時(shí)105分鐘?？梢?，本發(fā)明通過GPU多線程并行運(yùn)算，相比于現(xiàn)有方法中通過CPU運(yùn)算，具有很高的計(jì)算效率。本發(fā)明針對現(xiàn)有反偏移技術(shù)效果不佳或效率較低的問題，提出一種基于GPU的波動方程反偏移方法。通過逐炮進(jìn)行反偏移計(jì)算，對每一炮而言，CPU先將炮點(diǎn)、接收點(diǎn)信息和震源子波、速度、偏移值發(fā)送到GPU，然后在GPU上逐個(gè)時(shí)刻進(jìn)行散射波場延拓并進(jìn)行采樣記錄，記錄的散射波場即為一炮的反偏移結(jié)果，再將其發(fā)送到CPU上并寫入存儲器。通過本發(fā)明的方法反偏移精度好，計(jì)算效率高，可用于地震數(shù)據(jù)正演模擬和最小二乘偏移數(shù)據(jù)重構(gòu)。以上結(jié)合附圖詳細(xì)描述了本發(fā)明的優(yōu)選實(shí)施方式，但是，本發(fā)明并不限于上述實(shí)施方式中的具體細(xì)節(jié)，在本發(fā)明的技術(shù)構(gòu)思范圍內(nèi)，可以對本發(fā)明的技術(shù)方案進(jìn)行多種簡單變型，這些簡單變型均屬于本發(fā)明的保護(hù)范圍。另外需要說明的是，在上述具體實(shí)施方式中所描述的各個(gè)具體技術(shù)特征，在不矛盾的情況下，可以通過任何合適的方式進(jìn)行組合。為了避免不必要的重復(fù)，本發(fā)明對各種可能的組合方式不再另行說明。此外，本發(fā)明的各種不同的實(shí)施方式之間也可以進(jìn)行任意組合，只要其不違背本發(fā)明的思想，其同樣應(yīng)當(dāng)視為本發(fā)明所公開的內(nèi)容。當(dāng)前第1頁1 2 3 當(dāng)前第1頁1 2 3