專利名稱：基于圖形處理器計算起伏地表直接疊前逆時偏移的方法
技術(shù)領域：
本發(fā)明涉及計算計算起伏地表直接疊前逆時偏移的方法，尤其涉及基于圖形處理器計算計算起伏地表直接疊前逆時偏移的方法。
背景技術(shù)：
隨著地震勘探技術(shù)的不斷發(fā)展，油氣地震勘探的重點正轉(zhuǎn)向起伏地表條件和復雜地質(zhì)條件的區(qū)域，如山地、灘海和沼澤地區(qū)等，這對地震勘探工作及資料處理提出了新的挑戰(zhàn)。傳統(tǒng)地震成像方法在復雜地表、復雜構(gòu)造地區(qū)主要存在兩方面的問題一是地表起伏大，表層速度結(jié)構(gòu)十分復雜；二是地下構(gòu)造復雜，如褶皺強烈、斷層發(fā)育、構(gòu)造陡峻、地層產(chǎn)狀變化大等。這種雙復雜結(jié)構(gòu)使得傳統(tǒng)地震成像方法在這類地區(qū)難以準確成像。起伏地表地區(qū)地震資料成像主要采用兩種方法一種是先進行表層波場校正，再偏移成像，該法在實際生產(chǎn)中占主導地位；另一種是直接從起伏地表進行深度域的偏移成像。表層波場校正一直是實際地震資料處理中一項關(guān)鍵處理技術(shù)，通常采用高程基準面靜校正將地震數(shù)據(jù)校正到固定基準面或者浮動基準面上的辦法解決地形起伏的影響。這種方法隱含著一個明顯的基本假設就是地表一致性假設，即在地表起伏不大，低速帶橫向速度變化緩慢的地區(qū)，地下淺、中、深層的反射經(jīng)過低速帶時，幾乎遵循同一路徑近乎垂直入射至地表，這時它們的靜校正量基本相等。在地形平緩的地區(qū)，近地表速度比地下速度慢很多，且射線路徑出射角較小，采用上述校正方法比較合適，但在復雜地表和復雜地下地區(qū)， 近地表速度和地下速度相差不大，且射線路徑出射角較大，因此，簡單的校正可能會扭曲波場，降低地震成像的質(zhì)量，此時常規(guī)處理不能產(chǎn)生正確的成像。波動方程基準面延拓方法為表層波場校正提供了一種更精確的技術(shù)手段。只要能夠得到比較準確的近地表速度模型， 以“波動方程基準面延拓”替換常規(guī)“時移靜校正”，對提高山前帶復雜構(gòu)造成像精度非常重要。直接從起伏地表進行深度域偏移成像的方法，將野外靜校正隱含地包括在其中，這種校正不但包含了旅行時的縱向分量，同時還包含了旅行時的橫向分量，因此，基于起伏地表的疊前深度偏移可以同時解決上文中提到的復雜地表和復雜地下構(gòu)造這種地質(zhì)上的雙重復雜性。起伏地表直接疊前偏移對山前帶復雜構(gòu)造成像的價值在于一方面，在近地表反射信息較好的地區(qū)，可基于起伏地表疊前深度偏移進行速度分析，自地表逐層向下優(yōu)化深度一速度模型；另一方面，如果前期工作能夠提供滿足構(gòu)造成像要求的宏觀速度模型，起伏地表疊前深度一步法偏移則可提供地表以下全面、精細的構(gòu)造圖像。接下來的問題是選擇何種偏移算法進行起伏地表直接偏移。KirchhofT積分法簡便、高效，適應任意觀測方式，且易于對局部目標成像。但Kirchhoff積分法是一種高頻近似，不能對復雜波場中的焦散、干涉現(xiàn)象進行很好的處理，在強橫向非均勻介質(zhì)中成像精度不高；與Kirchhoff偏移法相比， 起伏地表單程波波動方程疊前深度偏移法對復雜地區(qū)構(gòu)造成像精度高，且能處理多波至問題，對陡傾界面、斷層和斷面的刻畫更清晰。然而，單程波方法同樣是一種近似方法，成像角度限制在90度以內(nèi)，直接求解雙程波波動方程的逆時偏移方法突破了成像傾角的限制，是最精確地成像方法。因此，在復雜地表和復雜地下構(gòu)造這種雙重復雜性的地區(qū)應該開展起伏地表直接逆時偏移工作。對于起伏地表雙程波波動方程的求解，有限元法處理不規(guī)則地表比較方便，但計算效率低，計算精度難以提高，使用也不方便。有限元和其他方法相結(jié)合的離散波數(shù)和區(qū)域分裂等方法，盡管在一定程度上提高了計算效率，但計算精度低、使用不方便等問題并沒有根本解決。邊界積分或邊界元法的半解析性質(zhì)決定了該方法不能適用于地表速度變化較大情況，限制了其在勘探實際中的應用。

發(fā)明內(nèi)容
針對上述傳統(tǒng)方法的缺點，本發(fā)明提供了基于圖形處理器計算起伏地表直接疊前逆時偏移的方法。在本發(fā)明中包括下述兩個技術(shù)概念，基于圖形處理器的計算，所采用的計算方法是有限差分方法。申請人在此作出解釋
采用有限差分(FD)方法計算起伏地表直接疊前逆時偏移，此方法計算效率高，使用方便。然而，受算法穩(wěn)定性和頻散關(guān)系的影響，有限差分法波場模擬計算量巨大。本發(fā)明中波場傳播采用的是時間域二階差分、空間域高階差分格式。高階差分格式表明，波場傳播和應用成像條件時空間中每個網(wǎng)格點都是解耦的，即獨立的。此種方法下所有網(wǎng)格點可以并行計算，并行粒度很小。針對該問題，由于有限差分每個網(wǎng)格點并行粒度很小，而圖形處理器多由大量地計算核心組成，因此其處理運算時每個線程只計算一個網(wǎng)格點或者幾個網(wǎng)格點的值，對這種小粒度并行的問題相比傳統(tǒng)的中央處理器更有優(yōu)勢。如附圖1所示，顯示了本發(fā)明的方法算實現(xiàn)逆時偏移的原理流程圖，計算密集的波場傳播與成像部分借助圖形處理器并行計算，中間結(jié)果需要借助通訊保存在中央處理器內(nèi)存或硬盤上。由于運用空間高階差分法需要大量的內(nèi)存讀寫，以三維8階差分格式為例， 每計算一個網(wǎng)格點的值需要讀取網(wǎng)格點周圍25個網(wǎng)格點的數(shù)據(jù)，內(nèi)存讀取冗余度很高。針對這一問題，與傳統(tǒng)的計算方法直接使用圖形處理器上的內(nèi)存相比，本發(fā)明的方法進一步的通過圖形處理器共享存儲器(share memory)機制重用內(nèi)存數(shù)據(jù)，可以大幅度降低內(nèi)存讀取冗余度，提高算法效率。在上述基礎上，有限差分法波場模擬的一個缺陷是處理起伏地表比較困難自由邊界處的網(wǎng)格點需要特殊處理，即需要對起伏地表網(wǎng)格點進行細致的分類，在程序?qū)崿F(xiàn)上需要大量的邏輯判斷，這將大大影響圖形處理器的運行效率。本發(fā)明的方法通過起伏地表邊界條件核心函數(shù)克服了此一問題。在本發(fā)明中，GPU的全稱是圖形處理器(Graphic Processing Unit)，俗稱“顯卡”， 在市場上目前有多家公司有研發(fā)生產(chǎn)此類硬件，包括NVIDIA、AMD (原ATI)、VIA等x86體系顯卡生產(chǎn)商、POWERVR、ARM等基于RISC體系的顯卡生產(chǎn)商。本發(fā)明公開了基于圖形處理器計算起伏地表直接疊前逆時偏移的方法，包括下述步驟
(1)、讀取偏移參數(shù)及起伏地表邊界參數(shù)；
(2)、讀入單炮數(shù)據(jù)，并根據(jù)觀測系統(tǒng)確定偏移孔徑；
(3)、確定本炮數(shù)據(jù)的時間延拓步長；
(4 )、將偏移孔徑內(nèi)的速度模型以及震源子波上傳到圖形處理器；(5)、利用有限差分的方法正演震源波場到最大接收時間，記錄不同時刻的正演波場；
(6)、將單炮數(shù)據(jù)讀入圖形處理器并進行數(shù)據(jù)規(guī)則化；
(7)、利用有限差分的方法沿時間方向反延拓一步，并運用起伏地表邊界條件；
(8)、讀取步驟(5)中存儲的相應時刻的波場并運用互相關(guān)成像條件進行成像運算；
(9)、重復步驟(7)、(8)，至反延拓到t=0時刻；
(10)、重復步驟(2)-(9)，至最后一炮數(shù)據(jù);
(11)、對偏移結(jié)果去除低頻成像噪音，輸出偏移結(jié)果，將偏移結(jié)果傳輸回CPU內(nèi)存； 其中，上述步驟(5)、(7)、(8)、(9)、(10)是在圖形處理器中進行的。 在本發(fā)明中，上述步驟(9 )、( 10 )如果在步驟(8 )完成后一次反延拓到t=0或達到最后一炮，則步驟(9)、(10)可省略，對于本領域技術(shù)人員而言，這種情況是顯而易見可以理解的。在本發(fā)明中，步驟(4)、(6)所述的偏移孔徑內(nèi)的速度模型、單炮數(shù)據(jù)、震源子波傳遞給圖形處理器，是指將偏移孔徑內(nèi)的速度模型、單炮數(shù)據(jù)、震源子波的數(shù)據(jù)傳遞到GPU的內(nèi)存中(通常稱之為顯存)。所述的GPU的內(nèi)存是直接位于GPU芯片上的內(nèi)存。在本發(fā)明中，所述圖形處理器尤其可以是近些年廣泛興起的支持通用計算功能的圖形處理器，詞類處理器通常包括基于統(tǒng)一計算設備架構(gòu)平臺CUDA的圖形處理器，此圖形處理器由顯卡廠商NVIDIA推出；基于開放通用計算模型OpenCL的圖形處理器，此類圖形處理器包括NVIDIA、AMD、Intel, VIA等均有推出；基于微軟DirectCompute標準的圖形處理器，此類圖形處理器包括NVIDIA、AMD、htel、VIA等均有推出；優(yōu)選的，基于成熟度和開發(fā)快捷的需要，本發(fā)明中所用的圖形處理器優(yōu)選為基于統(tǒng)一計算設備架構(gòu)平臺CUDA的圖形處理器。本領域技術(shù)人員根據(jù)本發(fā)明的實質(zhì)精神采用其他類型的支持通用計算功能的圖形處理器依舊屬于本發(fā)明保護范圍?；谏鲜龌A，在本發(fā)明中，在步驟(4)、(6)的偏移孔徑內(nèi)的速度模型、震源子波數(shù)據(jù)、單炮數(shù)據(jù)由中央處理器通過圖形處理器的CUDA編程接口，傳送到圖形處理器中的渲染管道，到達可編程片段處理器；所述步驟(5)、(7)、(8)、(9)、(10)在圖形處理器的多個可編程片段處理器中的多個線程并行進行，通過圖形處理器的共享存儲器為線程提供、存儲、 傳輸數(shù)據(jù)，其具體的計算過程如下
步驟(5)調(diào)用統(tǒng)一計算設備架構(gòu)平臺CUDA的有限差分核心函數(shù)； 步驟(7)調(diào)用統(tǒng)一計算設備架構(gòu)平臺CUDA的起伏地表邊界條件核心函數(shù)； 步驟(8)調(diào)用統(tǒng)一計算設備架構(gòu)平臺CUDA的互相關(guān)成像條件函數(shù)； 在上述基礎上，步驟(9)用于判斷是否反延拓到t=0時刻，若否，則重復(7)，(8)兩步； 若是，則到步驟(10);步驟(10)用于判斷是否是最后一炮數(shù)據(jù)，若否，則重復(2)- (9)步驟；若是，則將偏移結(jié)果傳輸回CPU內(nèi)存。盡管上述給出的是采用基于統(tǒng)一計算設備架構(gòu)平臺CUDA的圖形處理器實現(xiàn)本發(fā)明目的的實現(xiàn)方法，然而本領域技術(shù)人員根據(jù)上述實現(xiàn)方法的實質(zhì)精神采用其它支持通用計算功能的圖形處理器并調(diào)用其中對應的函數(shù)庫或相應函數(shù)實現(xiàn)相應功能也是可以為本領域技術(shù)人員所能知悉的，依舊屬于本發(fā)明的保護范圍。在上述所述中，所用的詞匯“渲染管道、可編程片段處理器”對于本領域技術(shù)人員是廣為理解的，其他類似的稱呼如管線、R0P、處理單元、流處理器等亦可為本領域技術(shù)人員所理解為與本發(fā)明所用詞匯含義一致，改用不同的名稱，只要其針對的是圖新處理器上相同的硬件單元結(jié)構(gòu)，依舊屬于本發(fā)明的保護范圍。其它的技術(shù)用詞如共享存儲器用于指位于圖形處理器上的share memory，在硬件架構(gòu)中與圖形處理器片上緩存Ll Cache類似，所有的共享存儲器與Cache相通，用于為Cache提供有效和高速的數(shù)據(jù)支持；由于有限差分的方法需要大量的內(nèi)存讀寫，內(nèi)存冗余度高，尤其需要高速的數(shù)據(jù)傳輸支持，本發(fā)明的方法通過使用圖形處理器共享存儲器重用內(nèi)存數(shù)據(jù)，從而大幅度降低內(nèi)存讀取冗余度。與傳統(tǒng)方法計算起伏地表逆時偏移相比，本發(fā)明的方法具有如下技術(shù)優(yōu)勢
1、運用逆時偏移方法直接從起伏地表偏移，與傳統(tǒng)的起伏地表處理方法相比較，成像效果優(yōu)勢明顯，并且近地表的層位都能夠準確的成像，因此在復雜地表和復雜地下構(gòu)造這種地質(zhì)上的雙重復雜性的區(qū)域具有廣泛的應用前景。2、針對傳統(tǒng)的雙程波波動方程求解時計算效率低，計算精度難以提高，使用不方便等缺點，本發(fā)明的方法在具體求解上采用有限差分法求解。針對有限差分法波場模擬的處理起伏地表困難，自由邊界處的網(wǎng)格點需要特殊處理這一缺陷，本發(fā)明在處理自由邊界網(wǎng)格點時，采用統(tǒng)一計算設備架構(gòu)平臺CUDA的起伏地表邊界條件核心函數(shù)將網(wǎng)格點統(tǒng)一處理，既解決了自由邊界反射問題，又避免了大量的邏輯判斷，從而使得圖形處理器得以快速應用。3、本發(fā)明的方法利用圖形處理器實現(xiàn)計算方法的加速，解決了逆時偏移計算量甚巨的問題。


附圖1為本發(fā)明基于圖形處理器計算起伏地表直接疊前逆時偏移的方法原理圖； 附圖2為本發(fā)明基于圖形處理器計算起伏地表直接疊前逆時偏移的方法流程圖； 附圖3為SEG起伏地表速度模型；
附圖4為采用本發(fā)明的方法計算基于起伏地表直接疊前逆時偏移結(jié)果。
具體實施例方式下面結(jié)合附圖結(jié)合實施例對發(fā)明進行進一步闡述，本領域技術(shù)人員應知，本發(fā)明也可以有多種不同形式或使用不同的圖形處理器或者不同的通用計算函數(shù)實施，因此不應認為它局限于說明書列出的實施例。本發(fā)明的實際應用效果與具體的圖形處理器硬件相關(guān)，在一定范圍內(nèi)，圖形處理器所擁有的流處理器越多、速度越快，速度提升的效果將越好。結(jié)合附圖1、2，可以理解本發(fā)明的基于圖形處理器計算起伏地表直接疊前逆時偏移的方法的原理和流程，即利用圖形處理器可以同時處理眾多并發(fā)線程的原理實現(xiàn)高速計算，具體的包含下述步驟
(1)、讀取偏移參數(shù)及起伏地表邊界參數(shù)；
(2)、讀入單炮數(shù)據(jù)，并根據(jù)觀測系統(tǒng)確定偏移孔徑；
(3)、確定本炮數(shù)據(jù)的時間延拓步長；
(4 )、將偏移孔徑內(nèi)的速度模型以及震源子波上傳到圖形處理器；
(5)、利用有限差分的方法正演震源波場到最大接收時間，記錄不同時刻的正演波場；
(6)、將單炮數(shù)據(jù)讀入圖形處理器并進行數(shù)據(jù)規(guī)則化；(7)、利用有限差分的方法沿時間方向反延拓一步，并運用起伏地表邊界條件；
(8)、讀取步驟(5)中存儲的相應時刻的波場并運用互相關(guān)成像條件進行成像運算；
(9)、重復步驟(7)、(8)，至反延拓到t=0時刻；
(10)、重復步驟(2)-(9)，至最后一炮數(shù)據(jù);
(11)、對偏移結(jié)果去除低頻成像噪音，輸出偏移結(jié)果，將偏移結(jié)果傳輸回CPU內(nèi)存； 其中，上述步驟(5)、(7)、(8)、(9)、(10)是在圖形處理器中進行的。附圖3為國際上標準的SEG起伏地表模型(Amoco和BP公司設計的加拿大起伏地表逆掩斷層模型)數(shù)據(jù)，其地質(zhì)原型是加拿大英屬哥倫比亞(British Columbia)東北部的逆掩斷層構(gòu)造。該數(shù)據(jù)共277炮，每炮480道，道距15 m，中間激發(fā)，記錄長度8s，采樣率 4ms,偏移距-3600nT3600m，地形最大高差1527m，速度360(T6000m/s。速度模型剖面長度 25 km，深度 10 km。附圖4為采用本發(fā)明的方法對原始數(shù)據(jù)進行起伏地表直接逆時偏移的偏移結(jié)果， 如圖4所示，逆掩斷層成像效果優(yōu)勢明顯，并且近地表的層位都能夠準確的成像。并且與單純采用中央處理器的方法相比，速度提升了 100倍以上，具有顯著的優(yōu)勢。
權(quán)利要求
1.基于圖形處理器計算起伏地表直接疊前逆時偏移的方法，包括下述步驟(1)、讀取偏移參數(shù)及起伏地表邊界參數(shù)；(2)、讀入單炮數(shù)據(jù)，并根據(jù)觀測系統(tǒng)確定偏移孔徑；(3)、確定本炮數(shù)據(jù)的時間延拓步長；(4 )、將偏移孔徑內(nèi)的速度模型以及震源子波上傳到圖形處理器；(5)、利用有限差分的方法正演震源波場到最大接收時間，記錄不同時刻的正演波場；(6)、將單炮數(shù)據(jù)讀入圖形處理器并進行數(shù)據(jù)規(guī)則化；(7)、利用有限差分的方法沿時間方向反延拓一步，并運用起伏地表邊界條件；(8)、讀取步驟(5)中存儲的相應時刻的波場并運用互相關(guān)成像條件進行成像運算；(9)、重復步驟(7)、(8)，至反延拓到t=0時刻；(10)、重復步驟(2)-(9)，至最后一炮數(shù)據(jù)；(11)、對偏移結(jié)果去除低頻成像噪音，輸出偏移結(jié)果，將偏移結(jié)果傳輸回中央處理器內(nèi)存；上述步驟(5)、(7)、(8)、(9)、(10)在圖形處理器中進行。
2.根據(jù)權(quán)利要求1所述的基于圖形處理器計算起伏地表直接疊前逆時偏移的方法，其特征在于步驟(4)的偏移孔徑內(nèi)的速度模型以及震源子波上傳于圖形處理器的內(nèi)存中、步驟(6)的單炮數(shù)據(jù)讀入圖形處理器的內(nèi)存中。
3.根據(jù)權(quán)利要求1所述的基于圖形處理器計算起伏地表直接疊前逆時偏移的方法，其特征在于所述圖形處理器是基于統(tǒng)一計算設備架構(gòu)平臺CUDA的圖形處理器。
4.根據(jù)權(quán)利要求3所述的基于圖形處理器計算起伏地表直接疊前逆時偏移的方法，其特征在于步驟(4)、(6)的偏移孔徑內(nèi)的速度模型、震源子波數(shù)據(jù)、單炮數(shù)據(jù)由中央處理器通過圖形處理器的CUDA編程接口，傳送到圖形處理器中的渲染管道，到達可編程片段處理器；所述步驟(5)、(7)、(8)、(9)、(10)在圖形處理器的多個可編程片段處理器的多個線程并行運算，通過圖形處理器共享內(nèi)存為線程提供數(shù)據(jù)。
5.根據(jù)權(quán)利要求3所述的基于圖形處理器計算起伏地表直接疊前逆時偏移的方法，其特征在于所述步驟(5)調(diào)用統(tǒng)一計算設備架構(gòu)平臺CUDA的有限差分核心函數(shù)。
6.根據(jù)權(quán)利要求3所述的基于圖形處理器計算起伏地表直接疊前逆時偏移的方法，其特征在于所述步驟(7)調(diào)用統(tǒng)一計算設備架構(gòu)平臺CUDA的起伏地表邊界條件核心函數(shù)。
7.根據(jù)權(quán)利要求3所述的基于圖形處理器計算起伏地表直接疊前逆時偏移的方法，其特征在于所述步驟(8)調(diào)用統(tǒng)一計算設備架構(gòu)平臺CUDA的互相關(guān)成像條件函數(shù)。
8.根據(jù)權(quán)利要求1所述的基于圖形處理器計算起伏地表直接疊前逆時偏移的方法，其特征在于采用有限差分法計算起伏地表直接疊前逆時偏移，其中波場傳播采用的是時間域二階差分、空間域高階差分格式，波場傳播和應用成像條件時空間中每個網(wǎng)格點都是解耦的，獨立的。
全文摘要
本發(fā)明公開了基于圖形處理器計算起伏地表直接疊前逆時偏移的方法，通過使用運用逆時偏移方法直接從起伏地表偏移，并采用有限差分法求解，利用圖形處理器的并行計算特性顯著提高了計算速度和效率。
文檔編號G01V1/28GK102353988SQ201110190520
公開日2012年2月15日 申請日期2011年7月8日 優(yōu)先權(quán)日2011年7月8日
發(fā)明者佟小龍, 劉洪 , 劉紅偉, 劉欽, 李博 申請人:中國科學院地質(zhì)與地球物理研究所, 北京吉星吉達科技有限公司