一種利用gpu提高波場延拓計算效率的方法
【專利摘要】本發(fā)明提供了一種利用GPU提高波場延拓計算效率的方法，屬于地震資料處理領(lǐng)域。所述方法首先擴(kuò)大震源波場和記錄波場的數(shù)據(jù)網(wǎng)格，然后利用GPU進(jìn)行多線程并行計算，并應(yīng)用CUDA提供的二維傅里葉變換函數(shù)實(shí)現(xiàn)波場延拓。本發(fā)明的方法簡單實(shí)用且效率高，利用本發(fā)明得到的成像效果與CPU單程波動方程疊前深度偏移得到的成像效果一致，但是計算效率提高了30倍以上。
【專利說明】一種利用GPU提高波場延拓計算效率的方法
【技術(shù)領(lǐng)域】
[0001]本發(fā)明屬于地震資料處理領(lǐng)域，具體涉及一種利用GPU提高波場延拓計算效率的方法。
【背景技術(shù)】
[0002]波動方程疊前深度偏移最核心的工作是地震波場的深度外推。深度外推過程必須依賴于波場外推算子，而波場外推算子實(shí)質(zhì)上是求解單向波動方程。單炮道集波動方程疊前深度偏移成像要分別對炮點(diǎn)波場和檢波點(diǎn)波場進(jìn)行下行波和上行波外推，利用激發(fā)時間成像條件提取每一外推層的成像值。
[0003]波動方程疊前深度偏移成像計算是以單炮為一個計算單元，單炮的波場延拓，每外推一步，需要兩次二維傅里葉正變換，一次相移，兩次二維傅里葉反變換，一次時移，其中二維傅里葉變換計算量最大，因此需要提高其計算效率。

【發(fā)明內(nèi)容】

[0004]本發(fā)明的目的在于解決上述現(xiàn)有技術(shù)中存在的難題，提供一種利用GPU提高波場延拓計算效率的方法，根據(jù)GPU (Graphic Processing Unit圖形處理器)的計算特點(diǎn),應(yīng)用CUDA (Compute Unified Device Architecture,統(tǒng)一計算設(shè)備架構(gòu))提供的傅里葉變換函數(shù)，有效的提高單程波動方程疊前深度偏移成像的計算效率，縮短單程波動方程疊前深度偏移成像的處理周期，快速地進(jìn)行波場延拓計算，提高波場延拓的計算效率，提高炮域疊前深度偏移的效率。
[0005]本發(fā)明是通過以下技術(shù)方案實(shí)現(xiàn)的:
[0006]一種利用GPU提高波場延拓計算效率的方法，所述方法首先擴(kuò)大震源波場和記錄波場的數(shù)據(jù)網(wǎng)格，然后利用GPU進(jìn)行多線程并行計算，并應(yīng)用CUDA提供的二維傅里葉變換函數(shù)實(shí)現(xiàn)波場延拓。
[0007]所述方法包括以下步驟:
[0008](I)擴(kuò)大震源波場和記錄波場的數(shù)據(jù)網(wǎng)格得到處理后的震源波場和記錄波場；
[0009](2)應(yīng)用裂步傅里葉法(SSF)對所述處理后的震源波場和記錄波場進(jìn)行波場延拓；
[0010](3)設(shè)延拓深度iz = I ;
[0011](4)利用CUDA提供的二維傅里葉變換函數(shù)對震源波場進(jìn)行二維正傅里葉變換；
[0012](5)利用CUDA提供的二維傅里葉變換函數(shù)對記錄波場進(jìn)行二維正傅里葉變換；
[0013](6)參考速度場的相移計算:利用步驟用(4)和步驟(5)的結(jié)果進(jìn)行相移計算，參考速度場是平均速度，即輸入速度場的平均值(請參考Stoffa P L.Split-step Fouriermigration.Geophysics,1990,55(2) ;410 ?421)；
[0014](7)利用CUDA提供的二維傅里葉變換函數(shù)對震源波場進(jìn)行二維反傅里葉變換；
[0015](8)利用CUDA提供的二維傅里葉變換函數(shù)對記錄波場進(jìn)行二維反傅里葉變換；[0016](9)擾動速度場的波場計算與相關(guān)成像:利用步驟(7)和步驟(8)的結(jié)果進(jìn)行相移計算，擾動速度場是輸入速度場與參考速度場的差值(請參考StofTa P L.Split-stepFourier migration.Geophysics,1990,55(2) ;410 ?421)；
[0017](10)判斷iz+1是否大于NZ，所述NZ =延拓深度/延拓步長，如果是，則轉(zhuǎn)入步驟
(11)，如果否，則返回步驟⑷；
[0018](11)結(jié)束。
[0019]所述步驟(I)是這樣實(shí)現(xiàn)的:對震源波場和記錄波場的長度進(jìn)行鑲邊處理，即當(dāng)震源波場和記錄波場的長度小于1000時，將其長度鑲邊處理為128的倍數(shù)，對于不足的數(shù)據(jù)處填充零；當(dāng)震源波場和記錄波場的長度大于1000而小于2048時，將其長度鑲邊處理成256的倍數(shù)，對于不足的數(shù)據(jù)處填充零。
[0020]所述步驟⑵至(11)都是在GPU上進(jìn)行的。
[0021]與現(xiàn)有技術(shù)相比，本發(fā)明的有益效果是:本發(fā)明的方法簡單實(shí)用且效率高，利用本發(fā)明得到的成像效果與CPU單程波動方程疊前深度偏移得到的成像效果一致，但是計算效率提高了 30倍以上。
【專利附圖】

【附圖說明】
[0022]圖1是本發(fā)明方法的波場延拓的示意圖。
[0023]圖2是本發(fā)明方法的GPU運(yùn)算的流程圖。
[0024]圖3-1是本發(fā)明方法實(shí)施例中利用CPU計算得到的鹽丘模型偏移剖面。
[0025]圖3-2是本發(fā)明方法實(shí)施例中利用GPU計算得到的鹽丘模型偏移剖面。
[0026]圖4-1是本發(fā)明方法實(shí)施例中利用CPU計算得到的負(fù)向構(gòu)造偏移剖面。
[0027]圖4-2是本發(fā)明方法實(shí)施例中利用GPU計算得到的負(fù)向構(gòu)造偏移剖面。
[0028]圖5是利用CPU和GPU進(jìn)行波場延拓的計算效率對比圖。
[0029]圖6是本發(fā)明方法的波場延拓的程序框圖。
【具體實(shí)施方式】
[0030]下面結(jié)合附圖對本發(fā)明作進(jìn)一步詳細(xì)描述:
[0031]如圖6所示(圖6中的NW表示頻率片)，本發(fā)明利用GPU進(jìn)行多線程并行計算，應(yīng)用CUDA (CUDA是一種將GPU作為數(shù)據(jù)并行計算設(shè)備的軟硬件體系，采用類C語言進(jìn)行軟件開發(fā)。)提供的二維(2D)傅里葉變換函數(shù)，結(jié)合波場延拓的計算特點(diǎn)(即二維傅立葉計算密集度大)，采用擴(kuò)大波場數(shù)據(jù)網(wǎng)格和多信號同時變換的策略，從而達(dá)到提高計算效率的目的，該方法簡單實(shí)用，效率高。
[0032]具體來說，在GPU計算中，內(nèi)核函數(shù)是以線程網(wǎng)格(Grid)的形式組織的，每個線程網(wǎng)格由若干個線程塊(Block)組成，而每個線程塊又由多個線程(thread)組成。在每次調(diào)用GPU計算時，都要指定線程的網(wǎng)格結(jié)構(gòu)〈〈〈dimGrid，dimBloCk>>>。在相移和時移的計算中，根據(jù)炮點(diǎn)波場和檢波點(diǎn)波場數(shù)據(jù)體分配線程結(jié)構(gòu):每個線程計算三維體中的一個點(diǎn)，所有的點(diǎn)可以在由一個線程網(wǎng)格計算完成。因?yàn)樘岣吡怂惴ǖ牟⑿卸?，所以計算效率大幅度提升。同時，將擾動速度場(時移)計算和相關(guān)成像計算進(jìn)行合并計算，可以減少數(shù)據(jù)一次讀寫和一次線程調(diào)用，進(jìn)一步提高了計算效率。[0033]本發(fā)明對炮域疊前深度偏移算法中廣泛使用的二維傅里葉變換CUFFT (由CUDA數(shù)學(xué)函數(shù)庫提供)進(jìn)行了深入的測試分析，發(fā)現(xiàn)CUFFT對信號長度比較敏感，當(dāng)信號的長度為奇數(shù)時，運(yùn)算速度較慢；當(dāng)信號的長度為偶數(shù)時，速度較快；當(dāng)信號長度是長度為2的冪次時，速度最快。
[0034]CUFFT的計算速度非常快，計算時間與信號長度基本上呈正比關(guān)系。從這個結(jié)果出發(fā)，本發(fā)明對CUFFT的信號長度做鑲邊處理，以提高CUFFT的計算效率。但是，如果將信號長度都鑲邊成2n，那么在提高CUFFT計算效率的同時，也會增加數(shù)據(jù)的存儲量，因此要找一個折中點(diǎn)，也就是在稍微增加數(shù)據(jù)量的基礎(chǔ)上大幅提高CUFFT的計算效率；當(dāng)信號長度小于1000時，信號長度為128的倍數(shù)的有較優(yōu)的計算效率，當(dāng)信號長度在2048以內(nèi)，信號長度為256的倍數(shù)的有較優(yōu)的計算效率。因此本發(fā)明提出擴(kuò)大波場數(shù)據(jù)網(wǎng)格以提高CUFFT的計算效率。
[0035]在CUFFT中，一次計算NW個2DFFT要比NW次計算2DFFT快一倍，因此可以將原來的波場延拓的NW次計算震源波場2D正FFT，記錄波場2D正FFT，相移計算，震源波場2D反FFT,記錄波場2D反FFT，擾動速度場(時移)計算，相關(guān)成像七步計算，更改為一次震源波場2D正FFT (NW)，記錄波場2D正FFT (NW)，相移計算，震源波場2D反FFT (NW)，記錄波場2D反FFT (NW)，擾動速度場(時移)計算+相關(guān)成像六步計算，從而提高波場延拓的計算速度。
[0036]用CUDA提供的cufftPlanMany (CUDA傅立葉變換批量處理)進(jìn)行多信號同時變換(是在(4) (5) (7) (8)這四個步驟中實(shí)現(xiàn)的)，具體如下:
[0037]cufftPlanMany (&plan，2，dims, NULL, 1,0, NULL, 1,0, CUFFT C2C, nw)；
[0038]cufftExecC2C(plan, d_signal_s, d_signal_s,CUFFT_F0RWARD)；
[0039]cufftExecC2C(plan, d_signal_r, d_signal_r, CUFFT_F0RWARD)；
[0040]cufftExecC2C(plan, d_signal_s, d_signal_s, CUFFT_INVERSE)；
[0041]cufftExecC2C (plan, d_signal_r, d_signal_r, CUFFT_INVERSE);所述
[0042]圖1給出的是波場延拓示意圖，單炮道集波動方程疊前深度偏移成像要分別對炮點(diǎn)波場和檢波點(diǎn)波場進(jìn)行下行波和上行波外推，利用激發(fā)時間成像條件提取每一外推層的成像值。
[0043]圖2給出的是GPU運(yùn)算最終流程，具體如下:在每個波場延拓的過程中，依次進(jìn)行以下計算正變換，相移計算，fft反變換，時移+相關(guān)成像。
[0044]圖3-1是本發(fā)明方法實(shí)施例中利用CPU計算得到的鹽丘模型偏移剖面，圖3-2是本發(fā)明方法實(shí)施例中利用GPU計算得到的鹽丘模型偏移剖面，對比圖3-1和圖3-2可以看出利用CPU和GPU得到的成像結(jié)果是一致的。圖4-1是本發(fā)明方法實(shí)施例中利用CPU計算得到的負(fù)向構(gòu)造偏移剖面，圖4-2是本發(fā)明方法實(shí)施例中利用GPU計算得到的負(fù)向構(gòu)造偏移剖面，對比圖4-1和圖4-2可以看出利用CPU和GPU得到的成像結(jié)果是一致的。
[0045]圖5表示的是圖2的每一步分別在CPU和GPU上的計算時間對比，從圖5中可以看出，利用GPU后在計算效率提高了 30倍以上。
[0046]上述技術(shù)方案只是本發(fā)明的一種實(shí)施方式，對于本領(lǐng)域內(nèi)的技術(shù)人員而言，在本發(fā)明公開了應(yīng)用方法和原理的基礎(chǔ)上，很容易做出各種類型的改進(jìn)或變形，而不僅限于本發(fā)明上述【具體實(shí)施方式】所描述的方法，因此前面描述的方式只是優(yōu)選的，而并不具有限制性的意義。
【權(quán)利要求】
1.一種利用GPU提高波場延拓計算效率的方法，其特征在于:所述方法首先擴(kuò)大震源波場和記錄波場的數(shù)據(jù)網(wǎng)格，然后利用GPU進(jìn)行多線程并行計算，并應(yīng)用CUDA提供的二維傅里葉變換函數(shù)實(shí)現(xiàn)波場延拓。
2.根據(jù)權(quán)利要求1所述的利用GPU提高波場延拓計算效率的方法，其特征在于:所述方法包括以下步驟: (1)擴(kuò)大震源波場和記錄波場的數(shù)據(jù)網(wǎng)格得到處理后的震源波場和記錄波場； (2)應(yīng)用裂步傅里葉法對所述處理后的震源波場和記錄波場進(jìn)行波場延拓； (3)設(shè)延拓深度iz= I ; (4)利用CUDA提供的二維傅里葉變換函數(shù)對震源波場進(jìn)行二維正傅里葉變換； (5)利用CUDA提供的二維傅里葉變換函數(shù)對記錄波場進(jìn)行二維正傅里葉變換； (6)參考速度場的相移計算:利用步驟用⑷和步驟(5)的結(jié)果進(jìn)行相移計算，參考速度場是平均速度，即輸入速度場的平均值； (7)利用CUDA提供的二維傅里葉變換函數(shù)對震源波場進(jìn)行二維反傅里葉變換； (8)利用CUDA提供的二維傅里葉變換函數(shù)對記錄波場進(jìn)行二維反傅里葉變換； (9)擾動速度場的波場計算與相關(guān)成像:利用步驟(7)和步驟(8)的結(jié)果進(jìn)行相移計算，擾動速度場是輸入速度場與參考速度場的差值； (10)判斷iz+1是否大于NZ,所述NZ=延拓深度/延拓步長，如果是，則轉(zhuǎn)入步驟(11)，如果否，則返回步驟(4)； (11)結(jié)束。
3.根據(jù)權(quán)利要求2所述的利用GPU提高波場延拓計算效率的方法，其特征在于:所述步驟(I)是這樣實(shí)現(xiàn)的:對震源波場和記錄波場的長度進(jìn)行鑲邊處理，即當(dāng)震源波場和記錄波場的長度小于1000時，將其長度鑲邊處理為128的倍數(shù)，對于不足的數(shù)據(jù)處填充零；當(dāng)震源波場和記錄波場的長度大于1000而小于2048時，將其長度鑲邊處理成256的倍數(shù)，對于不足的數(shù)據(jù)處填充零。
4.根據(jù)權(quán)利要求2所述的利用GPU提高波場延拓計算效率的方法，其特征在于:所述步驟⑵至(11)都是在GPU上進(jìn)行的。
【文檔編號】G01V1/28GK103675895SQ201210315284
【公開日】2014年3月26日 申請日期:2012年8月30日 優(yōu)先權(quán)日:2012年8月30日
【發(fā)明者】孔祥寧, 張慧宇, 段心標(biāo), 徐兆濤, 張兵, 孫武亮 申請人:中國石油化工股份有限公司, 中國石油化工股份有限公司石油物探技術(shù)研究院