張生強,韓立國,李 才,閆 濤,王玉秀,麻旭剛

(1.中海石油(中國)有限公司天津分公司,天津300452;2.吉林大學地球探測科學與技術(shù)學院,吉林長春130026)

基于高分辨率反演譜分解的儲層流體流度計算方法研究

張生強1,韓立國2,李 才1,閆 濤1,王玉秀1,麻旭剛1

(1.中海石油(中國)有限公司天津分公司,天津300452;2.吉林大學地球探測科學與技術(shù)學院,吉林長春130026)

反射地震數(shù)據(jù)中的低頻信息包含了與儲層及流體有關(guān)的豐富信息，從地震數(shù)據(jù)中提取儲層流體流度屬性可以為利用地震低頻信息進行儲層預測和流體識別提供一種新的途徑。為此，研究并提出了基于高分辨率稀疏反演譜分解的儲層流體流度計算方法。首先基于Biot孔隙介質(zhì)依賴頻率的反射系數(shù)低頻漸近分析理論，推導出了儲層流體流度屬性的計算表達式；然后利用地震數(shù)據(jù)低頻段優(yōu)勢頻率的瞬時譜振幅代替相應(yīng)頻率處的反射系數(shù)，給出了儲層流體流度屬性的直接近似計算方法，其中關(guān)于瞬時譜的計算采用了基于交替方向算法的高分辨率稀疏反演譜分解方法，該方法相對于常規(guī)譜分解方法具有更高的時間分辨率和頻率分辨率。陸上和海上二維疊后地震資料的試處理結(jié)果表明，基于高分辨率稀疏反演譜分解的儲層流體流度計算方法得到的儲層流體流度屬性剖面分辨率非常高，對于含油氣儲層顯示了良好的成像能力，降低了儲層流體識別的多解性和不確定性。

儲層流體識別；流體流度；反演譜分解；交替方向法；低頻信息

反射地震數(shù)據(jù)中的低頻信息包含了與儲層及流體有關(guān)的豐富信息，對于地震儲層預測和油氣檢測非常重要[1]。Silin等[1-5]對Biot模型含流體孔隙介質(zhì)分界面的地震反射信息進行了低頻漸近分析，通過引入一個小的無量綱參數(shù)進行推導，得到了依賴頻率的含流體孔隙介質(zhì)的反射系數(shù)漸近表達式，該反射系數(shù)正比于流體流度、地震信號的頻率和流體密度三者乘積的平方根。同時，Silin等[1-5]預測縱波反射系數(shù)在低頻端達到峰值?；贐iot孔隙介質(zhì)模型的低頻漸近分析理論可應(yīng)用于地震正演、反演以及屬性分析等方面。由于儲層流體流度的大小變化可以引起烴飽和儲層介質(zhì)的速度頻散變化[6]，因此可以利用儲層流體流度屬性進行儲層預測和流體檢測。

Korneev等[7]利用低頻漸近分析理論對儲層產(chǎn)油率進行了地震成像。Goloshubin等[8]在低頻漸近分析理論的基礎(chǔ)上提出了成像屬性的概念(其正比于流體流度的平方根)，將井的產(chǎn)油率與流體流度聯(lián)系起來，并研究了成像屬性在含油氣儲層成像、預測產(chǎn)油率和區(qū)分油藏油水界面中的作用。Goloshubin等[9]進一步根據(jù)低頻漸近分析理論進行地震屬性分析，估計儲層滲透率。Goloshubin等[10]使用低頻漸近分析理論嘗試對氣藏下方產(chǎn)生的低頻陰影進行解釋。蔡涵鵬[11]在低頻漸近分析理論的基礎(chǔ)上對成像屬性的計算、流體流度屬性計算、滲透率與成像屬性的關(guān)系、成像屬性與油氣生產(chǎn)率的關(guān)系以及優(yōu)勢頻率的確定方法進行了系統(tǒng)研究，并用實際數(shù)據(jù)進行了試算。Chen等[12-13]在低頻漸近分析理論基礎(chǔ)上，利用地震低頻信息和廣義S變換方法推導了計算儲層流體流度的方法，并使用該方法處理了實際地震資料，取得了一定的效果。

由此可見，從地震數(shù)據(jù)的低頻信息中提取儲層流體流度屬性可以為利用地震低頻信息進行儲層預測和流體識別提供一種新的途徑，具有重要的研究意義。然而，現(xiàn)有的關(guān)于儲層流體流度的計算方法都存在相同的缺點，即分辨率比較低。針對該問題，本文應(yīng)用Biot模型孔隙介質(zhì)依賴頻率的反射系數(shù)低頻漸近分析理論，推導出基于高分辨率稀疏反演譜分解的儲層流體流度計算方法，并進一步測試了利用儲層流體流度屬性進行儲層預測和流體識別的效果。實際資料處理結(jié)果表明，該方法計算得到的儲層流體流度信息分辨率非常高，對于含油氣儲層顯示了良好的成像能力。

1 方法原理

1.1 儲層流體流度的計算

Silin等[5]對Biot孔隙彈性介質(zhì)模型進行了低頻漸近分析，取得了一系列成果：①在低頻范圍內(nèi)獲得了有效的漸近簡諧波解；②在平面縱波入射到滲透地層界面的情況下，獲得了反射系數(shù)和透射系數(shù)相對簡單的顯式表達式；③獲得了來自滲透地層的快縱波反射共振頻率的顯式表達式。本文主要根據(jù)滲透地層的快縱波反射共振頻率表達式和流體飽和孔隙介質(zhì)在低頻域中的地震反射系數(shù)漸近表達式推導儲層流體流度的計算表達式。

Silin等[5]推導的滲透地層的快縱波反射共振頻率表達式為：

(1)

Silin等[5]推導的流體飽和孔隙介質(zhì)在低頻域中的地震反射系數(shù)漸近表達式可表示為[7]：

(2)

式中：R0和R1是實系數(shù)，它們是巖石和孔隙流體的力學特性(包括密度、孔隙度和彈性系數(shù))的無量綱函數(shù)；i是虛數(shù)單位；ω是地震信號的角頻率；ρf是儲層的流體密度。從(2)式可以看出：含流體孔隙介質(zhì)儲層依賴頻率的地震反射系數(shù)R正比于儲層流體流度、孔隙流體密度和地震信號頻率三者乘積的平方根。

在滲透力學中，儲層的流體流度F被定義為儲層的滲透率與流體的粘滯系數(shù)之比，即F=κ/η。從(2)式可知，地震反射系數(shù)R是儲層流體流度F,地震信號的角頻率ω和流體密度ρf的復函數(shù)。所以，為了從地震數(shù)據(jù)中提取流體流度屬性信息，利用(2)式對頻率求導，即：

(3)

(4)

式中，系數(shù)C是巖石骨架和孔隙流體的彈性性質(zhì)以及地震信號角頻率的復函數(shù)，可通過測井數(shù)據(jù)估算[8]。

由(4)式可得儲層流體流度屬性的計算表達式

(5)

從(5)式可知，儲層流體流度與反射系數(shù)對頻率的一階導數(shù)的平方成正比。對地震信號進行時頻分解后，其單頻瞬時譜振幅或能量能夠準確刻畫該頻率的地震反射振幅或能量。所以，在實際應(yīng)用計算中，我們用頻率為ω的瞬時譜振幅a(ω)代替相應(yīng)頻率處的反射系數(shù)R，則有

(6)

由于相對于傳統(tǒng)的時頻分析方法，即Gabor變換、連續(xù)小波變換(CWT)、S變換(或廣義S變換)等，基于稀疏反演的譜分解方法具有更高的分辨率和更好的聚焦性[14]，所以本文采用基于稀疏反演的譜分解方法計算(6)式中的a(ω)。在理論計算中，還需要利用(1)式確定計算儲層流體流度的優(yōu)勢頻率ω。Silin等[1,5]指出，當巖石和孔隙流體的特性在合理范圍時，依賴頻率的地震反射系數(shù)在地震低頻端達到峰值，這個結(jié)論與在實際數(shù)據(jù)中觀測到的低頻陰影一致[8,15-18]。因此，我們根據(jù)(6)式使用地震信號的低頻成分計算儲層的流體流度屬性是合理的。

(6)式將儲層流體流度表示為拋物線函數(shù)的形式。所以，在有已知井的產(chǎn)能或滲透率參數(shù)的情況下，利用各井處計算得到的[?a(ω)/?ω]2，即可通過拋物線擬合求得參數(shù)C。而在勘探區(qū)內(nèi)少井或無井的情況下，由于C相當于一個比例系數(shù)，仍可在不依賴井的情況下直接計算[?a(ω)/?ω]2來獲取儲層流體流度的相對值[12-13]。

1.2 基于稀疏反演的譜分解方法

反演譜分解的主要思想是[19]：首先將譜分解描述為一個線性反演問題，然后使用稀疏反演算法求解該線性反演問題，從而得到一個稀疏的時頻譜。反演譜分解可以作為一個基追蹤問題或基追蹤去噪問題[20]來求解。

1.2.1 反演譜分解的數(shù)學模型

非平穩(wěn)地震褶積模型[14]表達式為

(7)

式中：s(t)表示地震記錄；wi(t)表示以下角標i對應(yīng)的頻率為主頻的子波；ri(t)表示與wi(t)對應(yīng)的依賴頻率的反射系數(shù)；N表示參與計算的反射系數(shù)或子波向量的數(shù)量；n(t)表示隨機噪聲，“*”代表褶積運算。

(7)式的物理意義是指地震記錄s(t)是由一系列不同頻率的子波wi(t)和與之對應(yīng)的依賴頻率的地層偽反射系數(shù)(這里ri(t)稱為偽反射系數(shù)是為了區(qū)別于真實的地下介質(zhì)反射系數(shù))ri(t)褶積后相加，再加上隨機噪聲n(t)組成的。根據(jù)線性代數(shù)理論，可以將(7)式寫成矩陣與向量相乘的形式：

(8)

或

b=Ax+n

(9)

式中：s=b表示地震記錄；Wi表示子波wi(t)對應(yīng)的褶積矩陣；A=(W1W2…WN)代表子波褶積矩陣庫；x=(r1r2…rN)T代表依賴頻率的地層偽反射系數(shù)矩陣。

由(9)式可知，當已知子波褶積矩陣庫A和地震記錄b時，便可以通過求解線性反演問題得到依賴頻率的偽反射系數(shù)矩陣x。數(shù)據(jù)矩陣(r1r2…rN)的行和列分別代表時間和頻率，所以將x轉(zhuǎn)化為(r1r2…rN)的形式，即可以認為是反演得到的時頻譜。通常在地球物理反演中，由(9)式定義的反演譜分解問題是一個欠定問題，為了降低問題(9)的不確定性和多解性，得到稀疏的時頻譜，需要對時頻譜x進行稀疏約束，進而將問題(9)轉(zhuǎn)化為基追蹤去噪問題，即：

(10)

1.2.2 交替方向優(yōu)化算法

最近幾年，學者們發(fā)展了多種求解基追蹤去噪問題的快速算法，如譜投影梯度法(Spectral Projected Gradient，簡稱SPGL1)[21]、快速迭代軟閾值法(Fast Iterative Shrinkage-Thresholding Algorithm，簡稱FISTA)[22]、分塊坐標梯度下降法(Block Coordinate Gradient Descent，簡稱為CGD)[23]和交替方向法(Alternating Direction Method，簡稱ADM)[24]等。相對于SPGL1，F(xiàn)ISTA和CGD等算法，Yang等[24]證明了ADM算法是一種高效、穩(wěn)健的重構(gòu)算法，具有更好的數(shù)值計算性能。因此，本文采用ADM算法求解無約束基追蹤去噪問題(10)，實現(xiàn)反演譜分解。

應(yīng)用ADM算法求解問題(10)時，首先要引入輔助變量r∈Cm，將問題(10)轉(zhuǎn)換為如下等價形式：

(11)

(11)式的增廣Lagrangian子問題為：

(12)

式中：Re(·)表示對復數(shù)取實部運算；y∈Cm是Lagrangian乘子；H表示共軛轉(zhuǎn)置運算；罰參數(shù)β>0。如果給定(xk,yk)，則可以通過求解交替極小化問題(12)的方法得到(rk+1,xk+1,yk+1)。求解原始問題(10)的ADM算法迭代框架如下[24]：

(13)

式中：Shrink(·,·)表示一維收縮算子；gkAH·(Axk+rk+1-b-yk/β)；γ>0是一常數(shù)；τ>0是近端參數(shù)(proximal parameter)。

以上所述算法即為ADM算法，該算法具有很好的數(shù)值計算性能，并且適合處理復數(shù)，所以適用于稀疏反演譜分解。

1.2.3 合成地震數(shù)據(jù)測試分析

圖1a是由不同頻率Ricker子波合成的地震信號，用于測試譜分解方法。第1個子波的中心時間為50ms，0相位，主頻為60Hz；第2個子波的中心時間為150ms，-90°相位，主頻為40Hz；第3個子波是由相同子波中心時間(270ms)、相同相位(-180°)、不同頻率(20Hz和60Hz)的兩個子波組成的；第4個子波是由兩個中心時間間隔很短(間隔40ms)的子波構(gòu)成的，其中一個子波中心時間為380ms，0度相位，主頻為30Hz，另一個子波中心時間為420ms，0度相位，主頻為30Hz。圖1b 為圖1a所示合成地震信號的理論時頻譜，圖1c 展示的是CWT時頻譜，圖1d展示的是S變換時頻譜，圖1e展示的是基于ADM算法的稀疏反演時頻譜。

對比圖1給出的3種譜分解方法所得結(jié)果可以看出：相對于常規(guī)譜分解方法CWT(圖1c)和S變換方法(圖1d)，基于ADM算法的反演譜分解方法(圖1e)具有更高的時間分辨率和頻率分辨率，可以分辨出多頻率成份重疊的兩個子波(第3個子波)及時間間隔很近的兩個子波(第4個子波)；而常規(guī)譜分解方法受時頻分辨率相互制約的限制很難將其分開。

圖1 合成地震信號及其基于ADM算法的稀疏反演譜分解與常規(guī)譜分解方法所得結(jié)果對比

2 實際地震數(shù)據(jù)試處理

將本文提出的基于高分辨率反演譜分解的儲層流體流度計算方法分別應(yīng)用于陸上和海上二維實際地震數(shù)據(jù)試處理，以驗證該方法在儲層預測及流體識別中的有效性。

2.1 陸上地震數(shù)據(jù)試處理

首先應(yīng)用儲層流體流度屬性檢測陸上實際地震數(shù)據(jù)的含氣儲層。圖2為某地區(qū)陸上二維實際地震數(shù)據(jù)的疊后剖面，圖中紅色箭頭指示的是目標儲層位置，該目標儲層已經(jīng)被勘探井證實是含氣儲層。經(jīng)頻譜分析可知該疊后地震數(shù)據(jù)的主頻在45Hz左右。

分別采用基于CWT譜分解和基于高分辨率稀疏反演譜分解的儲層流體流度計算方法處理圖2 對應(yīng)的實際地震數(shù)據(jù)，得到的儲層流體流度屬性剖面(紅色表示強能量異常)分別如圖3和圖4所示。

圖2 某地區(qū)陸上二維實際地震數(shù)據(jù)疊后剖面

圖3 陸上二維地震數(shù)據(jù)基于CWT譜分解的儲層流體流度計算方法得到的儲層流體流度剖面

圖4 陸上二維地震數(shù)據(jù)基于高分辨率稀疏反演譜分解的儲層流體流度計算方法得到的儲層流體流度剖面

從圖3和圖4中可以看到，含氣目標儲層在各自的儲層流體流度屬性剖面中均顯示為明顯的強能量異常(紅色箭頭標注)。通過對比可以看出，基于高分辨率稀疏反演譜分解方法計算得到的儲層流體流度屬性剖面(圖4)具有更高的分辨率和可識別性，背景干擾較少。同時，相對于圖2所示的疊后地震剖面，基于高分辨率稀疏反演譜分解方法計算得到的儲層流體流度屬性剖面能夠更好地反映含氣儲層的展布。

2.2 海上地震數(shù)據(jù)試處理

圖5為某海區(qū)過井A和井B的實際地震數(shù)據(jù)疊后剖面，黃色虛線標記的是目標儲層位置，綠色橢圓標注的是有效目標儲層位置，該目標儲層已經(jīng)被勘探井證實是含油儲層。井A在有效目標儲層位置鉆遇油層，而井B在目標儲層位置處鉆遇到水層。經(jīng)頻譜分析可知，該疊后地震數(shù)據(jù)的主頻約為55Hz。

分別采用基于CWT譜分解和基于高分辨率稀疏反演譜分解的儲層流體流度計算方法處理圖5 對應(yīng)的實際地震數(shù)據(jù)，得到的儲層流體流度屬性剖面(紅色表示強能量異常)分別如圖6和圖7所示。

從圖6和圖7中可以看到，井A處含油目標儲層在兩種方法計算的儲層流體流度屬性剖面中都同樣顯示為明顯的強能量異常(綠色橢圓標注)，而在井B處的目標儲層段無能量異常特征。與陸上地震數(shù)據(jù)的處理結(jié)果類似，基于高分辨率稀疏反演譜分解方法計算得到的儲層流體流度屬性剖面同樣具有更高的分辨率和可識別性，背景干擾少，且檢測結(jié)果與鉆井結(jié)果完全吻合。同時，相對于圖5 所示的疊后地震剖面，基于高分辨率稀疏反演譜分解方法計算得到的流體流度屬性剖面能夠更好地反映含油儲層的分布。

圖5 某海區(qū)過井A和井B的二維實際地震數(shù)據(jù)疊后剖面

圖6 海上二維地震數(shù)據(jù)基于CWT譜分解的儲層流體流度計算方法得到的儲層流體流度剖面

圖7 海上二維地震數(shù)據(jù)基于高分辨率稀疏反演譜分解的儲層流體流度計算方法得到的儲層流體流度剖面

3 結(jié)束語

本文應(yīng)用含流體孔隙介質(zhì)依賴頻率的反射系數(shù)漸近分析理論，推導出了基于高分辨率稀疏反演譜分解的儲層流體流度計算方法。使用該方法計算得到的儲層流體流度信息對于含油氣儲層具有良好的成像能力，且分辨率很高，為利用地震資料的低頻信息進行儲層預測和流體識別提供了一種新的途徑。同時，該方法同樣適用于無井的情況(特別是海上油氣勘探)，通過將流體流度計算公式中的比例系數(shù)C設(shè)定為某一常數(shù)，計算的儲層相對流體流度可以用來進行儲層流體的定性分析。

需要指出的是，由于儲層流體流度屬性的計算基于地震數(shù)據(jù)的低頻信息，所以在地震數(shù)據(jù)采集到處理的整個過程中都要注意保留地震信號的低頻成分，以提高所提取儲層流體流度屬性信息的可靠性。

[1] Silin D B,Korneev V A,Goloshubin G M,et al.A hydrologic view on Biot’s theory of poroelasticity[R].Berkeley:Lawrence Berkeley National Laboratory,2004:1-23

[2] Silin D B,Korneev V A,Goloshubin G M,et al.Low-frequency asymptotic analysis of seismic reflection from a fluid-saturated medium[J].Transport in Porous Media,2006,62(3):283-305

[3] Silin D,Goloshubin G.Seismic wave reflection from a permeable layer:low-frequency asymptotic analysis[C]∥ASME 2008 International Mechanical Engineering Congress and Exposition.Boston：American Society of Mechanical Engineers，2008:47-56

[4] Silin D.A low-frequency asymptotic model of seismic reflection from a high-permeability layer[R].Berkeley:Lawrence Berkeley National Laboratory,2009:1-28

[5] Silin D,Goloshubin G.An asymptotic model of seismic reflection from a permeable layer[J].Transport in Porous Media,2010,83(1):233-256

[6] Batzle M L,Han D H,Hofmann R.Fluid mobility and frequency-dependent seismic velocity-direct measurements[J].Geophysics,2006,71(1):N1-N9

[7] Korneev V A,Silin D,Goloshubin G M,et al.Seismic imaging of oil production rate[J].Expanded Abstracts of 74thAnnual Internat SEG Mtg,2004,1476-1479

[8] Goloshubin G,van Schuyver C,Korneev V,et al.Reservoir imaging using low frequencies of seismic reflections[J].The Leading Edge,2006,25(5):527-531

[9] Goloshubin G,Silin D,Vingalov V,et al.Reservoir permeability from seismic attribute analysis[J].The Leading Edge,2008,27(3):376-381

[10] Goloshubin G,Chabyshova E.A possible explanation of low frequency shadows beneath gas reservoirs[J].Expanded Abstracts of 82ndAnnual Internat SEG Mtg,2012,1-5

[11] 蔡涵鵬.基于地震資料低頻信息的儲層流體識別[D].成都:成都理工大學,2012 Cai H P.Reservoir fluid identification from low-frequency seismic wave[D].Chengdu:Chengdu University of Technology,2012

[12] Chen X H,He Z H,Zhu S X,et al.Seismic low-frequency-based calculation of reservoir fluid mobility and its applications[J].Applied Geophysics,2012,9(3):326-332

[13] Chen X H,He Z H,Pei X G,et al.Numerical simulation of frequency-dependent seismic response and gas reservoir delineation in turbidites:a case study from China[J].Journal of Applied Geophysics,2013,94(4):22-30

[14] 韓利.高分辨率全譜分解方法研究[D].長春:吉林大學,2013 Han L.Research on the methods of high-resolution full spectrum decomposition[D].Changchun:Jilin University,2013

[15] Castagna J P,Sun S,Siegfried R W.Instantaneous spectral analysis:detection of low-frequency shadows associated with hydrocarbons[J].The Leading Edge,2003,22(2):120-127

[16] Liu J,Marfurt K J.Instantaneous spectral attributes to detect channels[J].Geophysics,2007,72(2):P23-P31

[17] Sinha S,Routh P S,Anno P D,et al.Spectral decomposition of seismic data with continuous-wavelet transform[J].Geophysics,2005,70(6):P19-P25

[18] Wang Y.Seismic time-frequency spectral decomposition by matching pursuit[J].Geophysics,2007,72(1):V13-V20

[19] Portniaguine O,Castagna J.Inverse spectral decomposition[J].Expanded Abstracts of 74thAnnual Internat SEG Mtg,2004,1786-1789

[20] Chen S S,Donoho D L,Saunders M A.Atomic decomposition by basis pursuit[J].SIAM Journal on Scientific Computing,1998,20(1):33-61

[21] Van Den Berg E,Friedlander M P.Probing the Pareto frontier for basis pursuit solutions[J].SIAM Journal on Scientific Computing,2008,30(2):890-912

[22] Beck A,Teboulle M.A fast iterative shrinkage-thresholding algorithm for linear inverse problems[J].SIAM Journal on Imaging Sciences,2009,2(1):183-202

[23] Yun S,Toh K C.A coordinate gradient descent method for L1-regularized convex minimization[J].Computational Optimization and Applications,2011,48(2):273-307

[24] Yang J F,Zhang Y.Alternating direction algorithms for L1-problems in compressive sensing[J].SIAM Journal on Scientific Computing,2011,33(1):250-278

(編輯：戴春秋)

Computation method for reservoir fluid mobility based on high-resolution inversion spectral decomposition

Zhang Shengqiang1,Han Liguo2,Li Cai1,Yan Tao1,Wang Yuxiu1,Ma Xugang1

(1.TianjinBranchofCNOOCLtd.,Tianjin300452,China;2.CollegeofGeo-ExplorationScienceandTechnology,JilinUniversity,Changchun130026,China)

The low-frequency information of seismic reflection data contains abundant information related to reservoir and fluid.Extracting reservoir fluid mobility property from seismic data provides a new means for reservoir prediction and fluid recognition.Therefore,we study and propose a calculation method of reservoir fluid mobility based on high resolution inversion spectral decomposition.Firstly,the computation formula for reservoir fluid mobility is deduced based on asymptotic analysis theory of frequency-dependent reflection coefficient in saturate porous medium of Biot model.Then,the instantaneous spectral amplitude of the dominant frequency at low frequency replaces its reflection coefficient,and the direct approximate calculation method for reservoir fluid mobility is derived.It is emphasized that the related instantaneous spectrum is calculated by using the high resolution sparse inversion spectral decomposition based on alternating direction algorithm,which has higher time-resolution and frequency-resolution than the conventional spectral decomposition method.Finally,the reservoir fluid mobility calculation method is applied to process the land and marine post-stack seismic data.The actual data processing results show that the reservoir fluid mobility section based on high-resolution sparse inversion spectral decomposition has high-resolution and well imaging capability to oil and gas reservoir,which reduces the multi-solutions and uncertainty of reservoir fluid identification.

reservoir fluid identification,fluid mobility,inverse spectral decomposition,alternating direction method (ADM),low-frequency information

2014-04-27;改回日期：2014-10-12。

張生強(1987—)，男，博士研究生，助理工程師，主要從事儲層預測、地震資料綜合解釋方面的研究。

國家科技重大專項(2011ZX05025-001-07)資助。

P631

A

1000-1441(2015)02-0142-08

10.3969/j.issn.1000-1441.2015.02.004