沈驥千

（江蘇煤炭地質物測隊，南京 210046）

三維彈性波逆時偏移中多卡GPU的應用

沈驥千

（江蘇煤炭地質物測隊，南京 210046）

地震勘探對象的日益復雜化，對成像效果提出了更高的要求。逆時偏移在理論上可以對復雜構造準確成像，且二維逆時偏移的研究也取得了一系列成果。而逆時偏移巨大的存儲、計算量以及野外單炮數(shù)據(jù)量的龐大，已成為制約其在三維多分量資料中應用的瓶頸。針對有效邊界存儲在三維逆時偏移中存儲較大的問題，采用隨機邊界條件以降低逆時偏移存儲量并提高偏移效率。該邊界條件無需保存各分量在邊界處的波場值，即可實現(xiàn)正向延拓波場的逆時重構；考慮到三維逆時偏移需開辟較大數(shù)組空間用于計算，提出應用區(qū)域分解技術將計算數(shù)據(jù)分配到不同節(jié)點，實現(xiàn)基于MPI+CUDA的協(xié)同并行。SEG/EAGE鹽丘模型算例證明了上述方法的有效性。

三維彈性波逆時偏移；隨機邊界；區(qū)域分解；MPI+CUDA協(xié)同并行

0 引言

基于標量波理論的反射縱波勘探技術在以往的能源開發(fā)過程中發(fā)揮了重要作用，但隨著油氣開發(fā)對地震勘探精度要求的不斷提高和勘探目標構造和巖性復雜程度的不斷增加，縱波勘探精度的提高越來越受制于其理論上的局限與波場信息單一等缺陷，在碳酸鹽巖裂縫性儲層、煤層氣和頁巖氣儲層等的勘探方面難以獲得滿意效果。近年來，以矢量波動理論為基礎的多波地震勘探技術引起了業(yè)界的廣泛重視，同常規(guī)縱波勘探技術相比，多波地震技術的理論假設與實際情況更為接近，且橫波在復雜介質中具有獨特的傳播機理，故理論上多波地震技術解決實際問題的能力更強，更有利于提高勘探精度或實現(xiàn)某些特殊目標的精確勘探。

疊前深度偏移是多波地震技術的重要研究內容，對多波地震資料進行疊前深度偏移處理的意義在于：一方面，地下構造的精確成像需要借助高精度的疊前偏移技術來實現(xiàn)；另一方面，多波偏移剖面中所包含的地震波動力學特征還可用于巖性和流體識別；同時偏移得到的共成像點道集還可為振幅隨偏移距變化（AVO）、振幅隨入射角變化（AVA）或其它疊前反演工作提供輸入數(shù)據(jù)。因此，多波地震資料偏移應當實現(xiàn)兩個目標：（1）地下構造的縱、橫波準確成像；（2）獲取保真的縱橫波疊前反演道集和疊加剖面。

基于彈性波理論的多分量地震資料聯(lián)合逆時偏移技術是實現(xiàn)多波地震縱橫波偏移成像的有效工具[1]。一方面，彈性波逆時偏移技術以雙程波動方程為理論出發(fā)點，它無傾角限制，理論上可實現(xiàn)任意復雜構造的成像[2-5]；另一方面，相較于標量波動方程逆時偏移，多分量聯(lián)合逆時偏移保留了地震波的彈性、矢量特征，可以對氣云帶等特殊目標獲得良好成像效果[6]。

逆時偏移思想提出之初，在疊后地震資料偏移成像中獲得了良好的成像剖面[3]；但受限于計算機設備的發(fā)展，其巨大的計算量和存儲I/O（輸入/輸出）成為其工業(yè)化應用的巨大瓶頸。隨著計算機硬件設備的革新和GPU（Graphic Processing Unit）以及其編程架構CUDA（Computing Unified Device Archi?tecture）的快速發(fā)展，標量波動方程逆時偏移并行化算法已漸趨成熟[7-10]。但多分量地震資料聯(lián)合逆時偏移的并行加速尚存在諸多問題，主要表現(xiàn)在：①多分量地震聯(lián)合逆時偏移存儲量和計算量遠大于標量波動方程逆時偏移；②基于CPU/GPU協(xié)同并行的標量波動方程逆時偏移相關配套技術應用于彈性波逆時偏移尚需改進，如有效邊界存儲等。

目前的野外采集普遍向高密度、三維方向發(fā)展，為確保復雜構造觀測波場的連續(xù)性，多采用小炮距、小接收線距、長排列接收，部分探區(qū)單炮記錄過萬道[11]。這導致三維彈性波逆時偏移過程中開辟的數(shù)組空間所占內存達10 G，依靠單節(jié)點GPU的加速不能滿足存儲需求。另一方面，在二維波動方程逆時偏移中為解決逆時延拓與正向延拓時間序列不一致性，常采用有效邊界[12]存儲正向延拓邊界波場，實現(xiàn)正演波場的重構來保證與逆時延拓的時間一致性。但是在三維多分量地震資料聯(lián)合逆時偏移中，需要保存的邊界波場比之于二維情況增加數(shù)個數(shù)量級，加重了存儲負擔和存取I/O。

本文針對有效邊界存儲在三維情況下存儲過大的情況，提出應用隨機邊界條件，在不保存邊界波場的情況下，實現(xiàn)正向延拓波場的重構；考慮三維逆時偏移情況下，開辟數(shù)組空間內存較大，提出應用區(qū)域分解技術將計算數(shù)據(jù)分配到不同節(jié)點，實現(xiàn)基于MPI+CUDA的協(xié)同并行。

1 彈性波逆時偏移的基本原理

彈性波逆時偏移將多分量地震記錄當作彈性波場的邊值問題[13]來處理，由計算機實現(xiàn)炮點波場的正向延拓和接收點波場的逆時延拓，并在延拓過程中采用相應的算法實現(xiàn)縱橫波的成像。其輸入數(shù)據(jù)一般為未經坐標旋轉的三分量炮集，輸出數(shù)據(jù)一般包括反射縱波成像結果和轉換橫波成像結果。

彈性波逆時偏移的基本原理為：假定地下介質由一系列繞射點組成，地面、海底或井中接收到的三分量記錄為各繞射點產生的縱橫繞射波場在接收點處的疊加響應，其中各繞射點縱橫繞射波的產生時間為炮點波場的主能量到達該點的時間[14]，繞射點產生的縱橫繞射波振幅取決于該點的反射系數(shù)[15]，由此可通過求取各繞射點產生縱橫繞射波的時間和該時刻的縱橫波振幅來確定地下各繞射點的空間位置和反射系數(shù)，實現(xiàn)地下介質的縱橫波成像。實際成像過程一般由以下幾部分組成：①炮點波場的正向延拓和重構；②接收點波場的逆時延拓；③炮、檢波場延拓過程中的縱橫波解耦與橫波的標量化處理；④縱橫波成像。

1.1 一階速度——應力彈性波方程及其高階延拓差分格式

根據(jù)各向同性介質中，位移-應變關系、應力-應變關系、應力-位移關系[16]可推導得到各向同性介質中的一階應力-速度彈性波方程：

其中，υx、υy、υz為質點振動速度；σxx、σyy、σzz、σxy、σxz、σyz為應力分量；ρ為密度；υp為縱波速度，υs為橫波速度；x、y、z分別為空間坐標；t為時間。

在交錯網格空間中對式（1）進行差分離散可得到時間2階，空間任意偶數(shù)階的差分格式，以式（1）中第一個方程為例：

其中，Δx、Δy、Δz分別表示坐標系中x、y、z方向上的網格間距；Δt表示時間采樣間隔；i，j，k分別代表x、y、z方向上離散點的序號；n表示時間上的離散點的序號。且ρ、μ半網格點上的值由同一方向上相鄰整網格點求取算術平均得到。2N為差分階數(shù)，D(mN)為2N階交錯網格差分系數(shù)，其計算公式為[16]：

1.2 穩(wěn)定性條件與吸收邊界條件

彈性波方程交錯網格有限差分格式屬于顯示差分型遞推格式，計算當前時刻波場值時，需利用上一時刻波場值。由于數(shù)值計算截斷誤差的影響，隨著時間遞推誤差會累積，破壞數(shù)值穩(wěn)定性甚至引起計算溢出[17]。

通過對三維各向同性介質一階雙曲型彈性波方程進行平面諧波分析[16]可得穩(wěn)定性條件為：

式（4）中，Δx、Δy、Δz分別表示坐標系中x、y、z方向上的網格間距；Δt表示時間采樣間隔；υmax為背景縱、橫波速度的極大值。對于三維各向同性彈性介質，地震波傳播速度與傳播方向無關，所以，在三個方向上的穩(wěn)定條件相同。

采用無分裂式PML邊界條件解決階段邊界的偽反射問題[18]，依據(jù)無分裂PML邊界條件的實現(xiàn)思路，可得到式（1）基于無分裂PML邊界條件的控制方程；以式（1）中第一個方程為例：

上式中Ωxx、Ωxy、Ωxz為引入的中間變量，其所對應的控制方程為：

其中，d（x）、d（y）、d（z）為三個方向的衰減因子；根據(jù)公式（5）、（6）可依此類推基于無分裂PML邊界條件的彈性波控制方程，其計算步驟為：（1）對包含PML邊界的整個區(qū)域進行常規(guī)的彈性波波場求解；（2）在PML邊界區(qū)域將得到的常規(guī)波場減去衰減量，然后按照時間軸方向更新輔助變量。由控制方程和計算步驟可知，無分裂的PML邊界條件避免了對彈性波波場進行分裂，降低了存儲。由于采用同一套方程也不會出現(xiàn)在邊界處兩套方程耦合的編程難題。

1.3 互相關成像條件

成像條件決定著逆時偏移成像剖面的質量。目前業(yè)界廣泛應用互相關成像條件，其基本原理為：利用炮點波場和接收點逆時延拓波場做零延時互相關，在反射點位置的波場將得到加強，非反射點處的波場將削弱，最后所有時刻波場相加得到偏移剖面。互相關成像條件的實現(xiàn)算法為：

其中，Ipp（x，z）為PP波偏移剖面，Ips（x，z）為PS波偏移剖面，Isp（x，z）為SP波偏移剖面，Iss（x，z）為SS波偏移剖面，S_P（x，z）、S_S（x，z）為炮點縱、橫波波場，R_P（x，z）、R_S（x，z）為接收點縱波和橫波波場，純縱橫波波場的獲得可參考文獻[19，20]，tmax為地震記錄長度。

2 基于多卡GPU的三維彈性波逆時偏移及其相關技術

上述彈性波逆時偏移基本流程中，在利用互相關成像條件時需同時應用相同時刻的檢波點波場和炮點波場，上述兩個波場的時間序列相反，因此需要重構炮點波場以實現(xiàn)縱橫波成像；本文采用隨機邊界條件散射人工邊界的偽反射，實現(xiàn)波場的重構。為了解決三維情況下單卡GPU顯存不足的難題，利用區(qū)域分解技術將數(shù)據(jù)劃分，實現(xiàn)大規(guī)模數(shù)據(jù)的三維彈性波逆時偏移。

2.1 隨機邊界條件的應用

在解決炮點、檢波點波場時間序列不一致的問題中，國內外學者先后提出Checkpointing方法、有效邊界存儲策略[12]，來降低存儲量。在二維彈性波逆時偏移中，有效邊界存儲策略已獲得業(yè)界廣泛認可；有效邊界存儲只需保存PML區(qū)域差分階數(shù)層的邊界值，即可實現(xiàn)波場的準確重構。但在三維情況下，其存存儲量仍相當可觀，假設三維模型網格點個數(shù)為200 m×300 m×400 m，采樣點數(shù)為5 000，有限差分階數(shù)為12階，若采用有效邊界存儲其存儲量達60G左右，這顯然難以滿足生成要求。因此，提出使用隨機邊界對波場進行散射[21]，因隨機邊界不涉及波場的衰減，其過程可逆，學者們提出“以計算換存儲”，利用最后時刻波場即可重構得到每一時刻震源波場值。

隨機邊界條件實質為在計算區(qū)域外部隨機賦速度值，利用隨機速度的無規(guī)律性，將波場值散射，這樣在進行互相關成像和疊加時，無規(guī)律的波場值能量得到減弱，反射點成像值相對得到增強。所以，隨機邊界條件的關鍵之處在于隨機邊界層內隨機速度的選擇，常用隨機速度賦值公式為：

式中，Vrandom（x，y，z）為邊界層內（x，y，z）處應取的隨機速度，v（x，y，z）為對應的背景速度，rand是隨機函數(shù)求得的隨機數(shù)，d為隨機邊界層內的點（x，y，z）與區(qū)域邊界的距離，v0（x，y，z）為最大穩(wěn)定速度。

為了驗證隨機邊界條件的適用性，本文采用二維均勻速度模型，驗證隨機邊界的重構和散射效果。二維均勻速度模型大小1500 m×1500 m，縱波速度為2500 m/s，橫波速度1800 m/s，炮點位置位于（750 m，750 m），縱波源激發(fā)，記錄不同時刻vz分量的正演和重構波場快照。圖1即為均勻介質隨機邊界速度場，圖2、圖3、圖4即為不同時刻正向延拓和逆時重構的波場值，由此可驗證隨機邊界的有效性及波場重構的正確性。

2.2 區(qū)域分解技術及基于MPI+CUDA的彈性波逆時偏移流程

圖1 二維均勻介質隨機邊界速度場Figure1 2D homogeneous medium stochastic boundary velocity field

圖2 二維均勻介質隨機邊界正向延拓（250 ms）及重構波場快照Figure 2 2D homogeneous medium stochastic boundary forward continuation（250 ms）and wave field reconstruction snapshot

GPU中有大量晶體管，是專為執(zhí)行復雜的數(shù)學和幾何計算而設計的。當今的GPU已不再將這些計算局限于圖形渲染，其通用計算技術的發(fā)展和天生的并行處理機制在計算方面提供了高于CPU數(shù)十倍乃至上百倍于CPU的性能。業(yè)界已經利用GPU的并行特性在二維逆時偏移取得了一系列突破[22-24]，然而對于大規(guī)模數(shù)據(jù)處理通常需要較大的計算機內存（顯存），GPU目前的顯存一般為6～8 G，對于模型較大和實際數(shù)據(jù)的三維逆時偏移這顯然是不能滿足計算要求，需要進行特殊情況考慮。本文的解決策略為將地震數(shù)據(jù)進行區(qū)域分解[10]，形成多GPU協(xié)同處理單炮地震數(shù)據(jù)的快速計算方法；此方法在理論上可以對任意規(guī)模的三維地震資料進行逆時偏移成像。

圖3 二維均勻介質隨機邊界正向延拓（300 ms）及重構波場快照Figure 3 2D homogeneous medium stochastic boundary forward continuation（300 ms）and wave field reconstruction snapshot

圖4 二維均勻介質隨機邊界正向延拓（350 ms）及重構波場快照Figure 4 2D homogeneous medium stochastic boundary forward continuation（350 ms）and wave field reconstruction snapshot

在對彈性波方程的差分格式，進行有限差分計算編程時，將一個維度上（以y方向為例）的區(qū)域進行劃分，按照實際資料的大小和顯卡顯存的限制來決定劃分的數(shù)量。利用區(qū)域分解進行差分計算，需要解決的核心問題在于如何隱藏卡與卡間的數(shù)據(jù)交換。如圖5所示，為了保證差分計算的正確性，在每一個時間步需要在虛線區(qū)域進行數(shù)據(jù)交換。

圖5 三維數(shù)據(jù)體區(qū)域分解示意圖Figure 5 A schematic diagram of 3D data volume domain decomposition

以其中的兩塊卡為例，說明分解區(qū)域后每塊卡的計算與數(shù)據(jù)傳遞。如圖6所示，每塊卡可以分為三部分（圖中分別用A、B和C表示），C區(qū)域為兩塊卡之間的重疊部分，厚度由差分格式來決定，B區(qū)域為傳遞部分與C區(qū)域保持一樣的厚度，A區(qū)域為顯卡本身的計算區(qū)域。在正演和延拓計算的每一個時間步內，C區(qū)域的網格處在邊界而無法進行差分運算，所以需要將B區(qū)域的數(shù)據(jù)傳遞到C。這樣計算過程就分為了如下步驟：首先是每塊卡計算B區(qū)域；然后，每塊卡計算A區(qū)域的同時，傳遞B區(qū)域數(shù)據(jù)到相鄰顯卡。

圖6 區(qū)域分解數(shù)據(jù)交換示意圖Figure 6 A schematic diagram of domain decomposition data exchange

計算能力高于2.0的NVIDIA顯卡支持在同節(jié)點上進行顯卡與顯卡之間的數(shù)據(jù)交換。而本文提出的多卡并行是基于多節(jié)點的并行，在數(shù)據(jù)傳遞時需要將傳遞數(shù)據(jù)通過內存中轉傳遞，由MPI的每一個進程控制每塊顯卡，在MPI進程組內進行非阻塞的數(shù)據(jù)交換。整個的實現(xiàn)流程圖如圖7所示：①由主進程讀取到炮記錄和速度模型并分發(fā)于進程；②各進程控制GPU顯卡實現(xiàn)每一時間步的圖4B區(qū)域波場計算；③GPU回傳數(shù)據(jù)并在內存中交換，同時實現(xiàn)圖4A區(qū)域的波場計算；④循環(huán)迭代波場，最后由主進程規(guī)約得到最終的偏移結果。

3 數(shù)值算例

SEG/EAEG鹽丘模型是一個國際上通用測試三維復雜構造成像效果的地質模型。本文采用SEG/ EAEG鹽丘速度模型作為縱波模型，通過固定的縱橫波速度比構造橫波模型，炮集記錄由一階彈性波方程有限差分正演獲得。模型參數(shù)如下：模型大小3380 m×3380 m×1050 m，網格大小Δx=Δy=Δz=5 m，時間采樣間隔0.35 ms，記錄長度2.8 s，震源為35 Hz雷克子波。觀測系統(tǒng)參數(shù)為：總炮數(shù)4690炮，67條測線，炮線間距50 m，炮間距60 m，接收線間距10 m，接收道間距5 m，每條測線70炮，每炮30400道接收；偏移算法采用本文所提出的基于MPI+CUDA的三維彈性波逆時偏移并行算法。圖8為三維鹽丘模型的縱橫波速度模型以及多分量逆時偏移剖面。

另外，PP、PS的成像效果較SP、SS的成像效果更好。原因為：PP，PS波是在縱波傾斜入射到彈性界面時產生能量較強，而SP、SS是由縱波傾斜入射到彈性界面產生的轉換橫波再次入射得到的地震記錄，該情況涉及波型的多次轉換，情況復雜。因此PP成像結果主要代表了P波的反射系數(shù)，PS主要成像結果代表了PS波的反射系數(shù)；而SP、SS成像結果并不對應一種波的反射系數(shù)，多次轉換后并不清楚它的能量分布。PP波和PS波均可以對鹽丘模型構造得到準確的偏移成像結果。PP波主頻比PS波主頻低，其原因在于：模型資料的縱橫波記錄是由同一縱波源激發(fā)得到，且轉換橫波在空間域中的波長小于縱波波長，所以在模型資料的深度偏移成像結果中PS成像結果分辨率高于PP成像結果。

圖7 三維彈性波逆時偏移協(xié)同并行流程圖Figure 7 3D elastic wave RTM cooperative parallel flow chart

4 結論

隨機邊界條件在不保存邊界波場的情況下，可以實現(xiàn)正向延拓波場的逆時重構，避免了三維彈性波逆時偏移大規(guī)模的邊界存儲問題；區(qū)域分解技術的應用使得大規(guī)模數(shù)據(jù)資料的逆時偏移處理成為可能；在上述基礎上，本文提出的基于MPI+CUDA協(xié)同并行的三維彈性波逆時偏移方法可以對模型資料準確成像。

[1]Gaiser J，Strudley A.Acquisition and application of multicomponent vector wavefields：are they practical[J].first break，2005，23(6)：61-67.

[2]Whitmore，N.D.Iterative depth migration by backward time propaga?tion.1983 SEG Annual Meeting.Society of Exploration Geophysicists，1983.

[3]McMechan G A.Migration by extrapolation of time-dependent bound?ary values.Geophysical Prospecting，1983，31(3):413-420.

[4]Baysal E，Kosloff D D，Sherwood J W C.Reverse time migration.Geo?physics，1983，48(11):1514-1524.

[5]Loewenthal D，Mufti I R.Reverse time migration in the spatial fre?quency domain.Geophysics，1983，48(5):627-635.

[6]楊佳佳.多分量地震波逆時偏移的關鍵技術研究[D].中國海洋大學，2015.

[7]何兵壽，張會星，韓月.雙程聲波方程疊前逆時深度偏移及其并行算法.煤炭學報，2010，35(3)：458-462.

[8]李博，劉國峰，劉洪.地震疊前逆時偏移的一種圖形處理器提速實現(xiàn)方法[J].地球物理學報，2009，52(12)：245-252.

[9]劉紅偉，李博，劉洪等.地震疊前逆時偏移高階有限差分算法及GPU實現(xiàn)[J].地球物理學報，2010，53(7)：1725-1733.

[10]劉守偉，王華忠，陳生昌，等.三維逆時偏移CPU/GPU機群實現(xiàn)方案研究[J].地球物理學報，2013，56(10).

[11]唐祥功，匡斌，杜繼修.多GPU協(xié)同三維疊前逆時偏移方法研究與應用：石油地球物理勘探，2013，48（6）：910-914

[12]王保利，高靜懷，陳文超等.地震疊前逆時偏移的有效邊界存儲策略[J].地球物理學報，2012，55(7):2412-2421.

[13]Chang，W.F，McMechan G A.Elastic reverse-time migration[J].Geo?physics，1987，52(10)：1365-1375.

圖8 三維彈性波速度模型及逆時偏移成像結果Figure 8 3D elastic wave velocity model and RTM imaging result

[14]Claerbout，J.F.Toward a unified theory of reflector mapping[J].Geo?physics，1971，36(3):467-481.

[15]Zhi L，Chen S，Li X.Joint AVO inversion of PP and PS waves using exact Zoeppritz equation[M]//SEG Technical Program Expanded Ab?stracts 2013.Society of Exploration Geophysicists，2013：457-461.

[16]牟永光，裴正林.三維復雜介質地震數(shù)值模擬[M].北京:石油工業(yè)出版社，2005，173-176.

[17]董良國，馬在田，曹景忠.一階彈性波方程交錯網格高階差分解法穩(wěn)定性研究[J].地球物理學報，2000，06:856-864.

[18]王鵬飛，何兵壽.聲波方程逆時偏移中的無分裂PML吸收邊界條件[J].工程地球物理學報，2015，12(5):583-590.

[19]Dellinger J，Etgen J.Wave-field separation in two-dimensional anisotropic media.Geophysics，1990，55(7)：914-919.

[20]Du Q Z，Gong X F，Zhang M，et al.3D PS-wave imaging with elastic reverse-time migration.Geophysics，2014，79(5)：S173-S184.

[21]Robert G.Clapp.Reverse time migration with random boundar?ies.79th Annual International Meeting，SEG Expanded Abstracts，2009，28：2809-2813.

[22]Micikevicius P，3D finite difference computation on GPUs using CU?DA.//Processing of 2th Workshop on General Purpose Processing on Graphics Processing Units，Expanded Abstracts，2008:79-84.

[23]Michea D，Komatitsch D，Mahovsky J，et al.Industrial-scale reverse time migration on GPU hardware.SEG Expand Abstracts，2009:2789-2793.

[24]Sun X Y，Suh S.Maxing throughput for high performance TTI-RTM: From CPU-GPU.SEG Expanded Abstracts，2011，182(1):389-402.

Application of Multi-card GPU in 3D Elastic Wave Reverse-time Migration(RTM)

Shen Jiqian
(Jiangsu Provincial Coal Geological Exploration,Geophysical Prospecting and Surveying Team,Nanjing,Jiangsu 210046)

With the increasing complexity of seismic exploration objects,higher requirements are put forward on imaging results.The RTM can be used to describe complex structures in theory.A series of research findings have been acquired in 2D RTM;however,the huge amount of storage,computational workload and field single-shot data in RTM have become the bottleneck to restrict its applica?tion in 3D multi-component data.In allusion to the issue of larger efficient boundary storage in 3D RTM,the stochastic boundary ap?proach is used to lower down RTM storage and improve migration efficiency.This boundary condition needs not to keep wave field val?ue of components at boundary and realize reverse-time reconstruction.Whereas 3D RTM needs larger array space in computation,the author applies the domain decomposition technique to distribute computational data to different nodes and realize cooperative parallel based on MPI+CUDA.The SEG/EAGE salt dome example has proved the effectiveness of the method.

3D elastic wave RTM;stochastic boundary;domain decomposition;MPI+CUDA cooperative parallel

P631.4

A

10.3969/j.issn.1674-1803.2017.03.14

1674-1803(2017)03-0065-07

國家自然科學基金（41674118），國家重大科技專項（2016ZX05027-002）聯(lián)合資助

沈驥千（1970—），男，高級工程師，主要從事地震勘探資料處理工作。

2017-02-10

責任編輯：孫常長