低滲透致密砂巖氣層產(chǎn)能預(yù)測(cè)方法研究

范宜仁, 宋巖, 張海濤, 張鵬, 李思

(1.中國(guó)石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院, 山東 青島 266580; 2.中國(guó)石油大學(xué)(華東)CNPC測(cè)井重點(diǎn)實(shí)驗(yàn)室, 山東 青島 266580; 3.中國(guó)石油長(zhǎng)慶油田分公司勘探開發(fā)研究院, 陜西 西安 710018; 4.中國(guó)石油長(zhǎng)慶油田分公司, 陜西 西安 710018; 5.中國(guó)石油集團(tuán)測(cè)井有限公司華北事業(yè)部, 河北 任丘 062552)

0 引 言

儲(chǔ)層產(chǎn)能評(píng)價(jià)與預(yù)測(cè)結(jié)果可以檢驗(yàn)油氣勘探的成果,為油氣田開發(fā)提供最基本的依據(jù)[1-2]。產(chǎn)能預(yù)測(cè)主要有3種方法:①通過平面徑向流模型與測(cè)井參數(shù)結(jié)合求得油氣產(chǎn)量適用于中高孔隙度滲透率儲(chǔ)層,但對(duì)于滲流機(jī)理復(fù)雜、物性參數(shù)難以準(zhǔn)確計(jì)算的低滲透儲(chǔ)層精度不高[3];②根據(jù)地區(qū)經(jīng)驗(yàn)提取產(chǎn)能敏感參數(shù),在產(chǎn)能分級(jí)的基礎(chǔ)上應(yīng)用產(chǎn)能指數(shù)法、測(cè)試法、解釋圖版法等建立相應(yīng)的產(chǎn)能預(yù)測(cè)模型[4];③通過主成分分析、灰色建模理論、神經(jīng)網(wǎng)絡(luò)等建立智能產(chǎn)能估算模型,實(shí)現(xiàn)產(chǎn)能預(yù)測(cè)[5]。上述方法對(duì)儲(chǔ)層物性參數(shù)計(jì)算精度有較高要求,在低孔隙度低滲透率致密儲(chǔ)層中應(yīng)用效果較差。

針對(duì)致密砂巖含氣儲(chǔ)層孔滲關(guān)系復(fù)雜,滲透率計(jì)算精度低,導(dǎo)致傳統(tǒng)產(chǎn)能預(yù)測(cè)方法效果較差的問題,本文以鄂爾多斯盆地東部上古生界致密砂巖氣層為例,從巖石物理實(shí)驗(yàn)出發(fā),結(jié)合Biot-Gassmann理論和Xu-White模型,分析儲(chǔ)層含氣性與彈性參數(shù)的關(guān)系,提取產(chǎn)能敏感因子,探索從陣列聲波測(cè)井?dāng)?shù)據(jù)獲取表征儲(chǔ)層含氣量的地層彈性參數(shù),從而建立產(chǎn)能預(yù)測(cè)模型,以期提高致密砂巖氣層產(chǎn)能預(yù)測(cè)的準(zhǔn)確性。

1 基于常規(guī)測(cè)井資料的氣層產(chǎn)能預(yù)測(cè)

1.1 基于孔隙結(jié)構(gòu)的儲(chǔ)層分類

鄂爾多斯盆地東部多層系低滲透砂巖氣藏屬典型的巖性氣藏,受物源、沉積環(huán)境的影響,巖石礦物組分多變,純石英砂巖和巖屑砂巖并存;孔隙結(jié)構(gòu)復(fù)雜,儲(chǔ)集空間類型多樣,原生、次生孔隙并存,非均質(zhì)性強(qiáng)。對(duì)于孔隙結(jié)構(gòu)越好的儲(chǔ)層,其產(chǎn)氣能力越高,反之亦然。

根據(jù)研究區(qū)20口井50塊巖心的壓汞和核磁共振實(shí)驗(yàn)數(shù)據(jù),通過曲線形態(tài)及孔喉特征參數(shù)對(duì)儲(chǔ)層孔隙結(jié)構(gòu)特征進(jìn)行了研究與分析,排驅(qū)壓力0.18～2.92 MPa,中值壓力8 MPa左右,該區(qū)大孔隙發(fā)育較少,小孔隙占優(yōu);平均孔喉半徑為0.19 μm,中值孔喉半徑分布0.01～0.34 μm,孔喉均值絕大多數(shù)小于0.4 μm,孔喉半徑整體較小,屬微細(xì)喉道,以小孔-微細(xì)喉為主;微觀均質(zhì)系數(shù)0.04～0.23,均質(zhì)程度較差,表明研究區(qū)儲(chǔ)層非均質(zhì)性較強(qiáng)。通過對(duì)壓汞參數(shù)、物性資料及測(cè)井?dāng)?shù)據(jù)的綜合分析,將該區(qū)儲(chǔ)層分為3類,各類型孔隙結(jié)構(gòu)的特征、分類標(biāo)準(zhǔn)以及常規(guī)測(cè)井識(shí)別標(biāo)準(zhǔn)等見表1。

1.2 基于常規(guī)測(cè)井資料產(chǎn)能預(yù)測(cè)模型

從達(dá)西滲流產(chǎn)量公式出發(fā),以相對(duì)滲透率與含水飽和度的函數(shù)關(guān)系為紐帶[6],導(dǎo)出油氣儲(chǔ)層產(chǎn)能與儲(chǔ)層有效孔隙度、滲透率以及電阻率之間的理論模型,進(jìn)而根據(jù)常規(guī)測(cè)井資料計(jì)算各層產(chǎn)能。在鄂爾多斯盆地多層合試的試氣工藝條件下,總產(chǎn)氣量則通過各段氣層產(chǎn)能的累加獲得。

在油氣田開采過程中,油井投產(chǎn)后穩(wěn)定生產(chǎn)和壓差符合平面徑向流產(chǎn)量公式

(1)

式中,Qg為氣層日產(chǎn)氣量,m3/d;C為單位換算系數(shù);K0為滲透率,×10-3μm2;H為氣層厚度,m;Tsc為地面標(biāo)準(zhǔn)溫度,K;T為氣層溫度,K;pe為氣層壓力,MPa;pwf為流動(dòng)壓力, MPa;

psc為地層標(biāo)準(zhǔn)壓

*非法定計(jì)量單位,1 mD=9.87×10-4μm2,下同

力,MPa;μg為氣體黏度,mPa·s;Z為氣體偏差系數(shù);re為供氣半徑,m;rwf為井徑,m;S為視表皮系數(shù)。

效滲透率K0與絕對(duì)滲透率K以及相對(duì)滲透率Kro之間關(guān)系為

K0=KKro

(2)

相對(duì)滲透率與飽和度之間的相關(guān)公式為

(3)

式中,e和l分別為對(duì)應(yīng)常數(shù)項(xiàng)和指數(shù)項(xiàng)系數(shù),將式(2)和式(3)帶入式(1),可得儲(chǔ)層每米產(chǎn)氣指數(shù)的表達(dá)式

(4)

(5)

影響儲(chǔ)層產(chǎn)能的主要因素為有效孔隙度、滲透率以及儲(chǔ)層電阻率,它們之間為復(fù)雜的非線性關(guān)系。對(duì)于合試氣層則采用累加法進(jìn)行計(jì)算,見式(5)。

系數(shù)A0、A1、A2可采用最小二乘法、最優(yōu)化法、神經(jīng)網(wǎng)絡(luò)等確定,這里采用麥夸特優(yōu)化法,該法能克服多項(xiàng)式插值帶來的數(shù)值振蕩,收斂速度快且能有效減小擬合結(jié)果的誤差。從研究區(qū)34個(gè)試油層中選取20個(gè)層代表性樣本點(diǎn),確定待定系數(shù),14個(gè)層作為檢驗(yàn)樣本,計(jì)算不同類型儲(chǔ)層對(duì)應(yīng)的A0、A1、A2,結(jié)果見表2。實(shí)際預(yù)測(cè)Ⅲ類儲(chǔ)層產(chǎn)能時(shí)應(yīng)用Ⅱ類儲(chǔ)層模型。

表2 基于常規(guī)測(cè)井資料產(chǎn)能預(yù)測(cè)模型系數(shù)表

應(yīng)用上述模型對(duì)單層產(chǎn)能計(jì)算,再通過累加法求得合試層產(chǎn)能,計(jì)算效果見圖1,其相關(guān)系數(shù)為0.81?；诔Ｒ?guī)測(cè)井資料產(chǎn)能預(yù)測(cè)方法具有一定效果,但由于致密砂巖氣層的低滲透性,儲(chǔ)層滲透率計(jì)算精度低,導(dǎo)致預(yù)測(cè)效果偏差較大。

圖1 常規(guī)方法產(chǎn)能預(yù)測(cè)效果圖

2 基于陣列聲波測(cè)井的氣層產(chǎn)能預(yù)測(cè)方法

陣列聲波測(cè)井能夠提供反映地層真實(shí)骨架和流體特征的聲學(xué)信息,對(duì)儲(chǔ)層含氣性具有一定的敏感性[8],可利用聲學(xué)特性實(shí)驗(yàn)分析,提取產(chǎn)能敏感因子,建立低滲透致密砂巖氣層的產(chǎn)能預(yù)測(cè)模型。

2.1 致密砂巖聲學(xué)特性實(shí)驗(yàn)分析

2.1.1 含氣飽和度和縱橫波速的關(guān)系

根據(jù)威利平均時(shí)間公式(Wyllie Time-average Equation)

Δt=(1-φ)Δtma+φΔtf

(6)

式中,Δt、Δtma、Δtf分別為純巖石、骨架和孔隙流體的聲波時(shí)差,μs/m;φ為純巖石的孔隙度,小數(shù)。巖性和孔隙度一定時(shí),巖樣的縱波時(shí)差主要取決于孔隙流體,與含水飽和度存在正相關(guān)關(guān)系。

圖2 36塊巖樣平均縱橫波速度與含氣飽和度關(guān)系

圖2為常溫常壓不同飽和度下該地區(qū)巖樣縱橫波速度測(cè)試結(jié)果。由圖2可以看出氣水兩相條件下隨含氣飽和度的增加,縱波速度由高變低,干燥砂巖縱波速度最低,飽含水時(shí)最高,且含氣飽和度在10%左右時(shí)縱波速度變化幅度最大。橫波速度基本不受含氣性的影響。在獲得實(shí)驗(yàn)數(shù)據(jù)的同時(shí),基于Krief理論建立石英骨架砂巖體積模型,進(jìn)行了不同含氣飽和度下聲波速度的數(shù)值模擬計(jì)算,模型中骨架密度為2.65 g/cm3,體積模量38.2 GPa,切變模量42.7 GPa,孔隙度為10%,孔隙流體為氣和水,其體積模量分別為0.07 GPa、2.73 GPa。數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果相吻合(見圖3)。

圖3 縱橫波速度與含氣飽和度數(shù)值計(jì)算結(jié)果

2.1.2 含氣飽和度和彈性模量的關(guān)系

圖4 巖心彈性模量與含氣飽和度的關(guān)系

巖石的彈性性質(zhì)主要受巖石的孔隙結(jié)構(gòu)、固體相及孔隙流體的影響,反之巖石的彈性性質(zhì)也能反映巖石的孔隙結(jié)構(gòu)及孔隙流體的性質(zhì)[9]。巖石的彈性參數(shù)主要有泊松比、彈性模量等。實(shí)驗(yàn)結(jié)果(見圖4)表明,巖心含氣飽和度對(duì)其體積模量的影響最顯著,體積模量隨含氣飽和度增大而減小;彈性模量隨含氣飽和度增加略有減小,不如體積模量的變化明顯;巖心含氣飽和度對(duì)切變模量的影響很小,巖心切變模量隨含氣飽和度的增大基本不變。體積模量對(duì)儲(chǔ)層含氣性較為敏感,能更好地反映氣層的產(chǎn)能高低。

2.2 構(gòu)建產(chǎn)能敏感參數(shù)2.2.1 構(gòu)建產(chǎn)能敏感參數(shù)——體積模量的差值ΔK

基于對(duì)聲波實(shí)驗(yàn)的研究與分析,當(dāng)?shù)貙?00%飽含水時(shí),通過Biot-Gassmann方程預(yù)測(cè)的100%飽含水儲(chǔ)層的體積模量與通過陣列聲波測(cè)井資料計(jì)算的儲(chǔ)層體積模量基本相等,其差值ΔK=0;當(dāng)儲(chǔ)層含氣時(shí),預(yù)測(cè)其100%飽含水的體積模量與實(shí)際陣列聲波測(cè)井資料計(jì)算的體積模量的差值ΔK>0,ΔK越大,說明含氣越多,產(chǎn)能Qi越高。

ΔK=Ksat-K

(7)

(8)

(9)

式中,K為陣列測(cè)井聲波資料計(jì)算體積模量;ρ為密度測(cè)井?dāng)?shù)據(jù);vp為縱波波速;vs為橫波波速;Ksat為飽含流體巖石的體積模量;Kdry為干巖石的體積模量;Km為巖石基質(zhì)的體積模量;Kfl為流體的體積模量。Ksat通過低頻條件下的Biot-Gassmann理論進(jìn)行預(yù)測(cè)[10]見式(9);Kdry則應(yīng)用徐懷特模型結(jié)合Gassmann方程、Kuster-Toksoz理論以及有效差分介質(zhì)(DEM)理論進(jìn)行求解[11-12]。

2.2.2 構(gòu)建產(chǎn)能敏感參數(shù):含氣飽和度指示因子Ig

根據(jù)陣列聲波測(cè)井資料計(jì)算得到的體積模量,結(jié)合Xu-White模型計(jì)算的干巖石骨架體積模量Km,可以求得孔隙流體體積模量Kfl如式(10),進(jìn)而應(yīng)用含氣飽和度聲波巖石物理模型A.Brie公式計(jì)算近似含水飽和度Sw

(10)

(11)

式中,e為地區(qū)經(jīng)驗(yàn)參數(shù),通常為2～5。通過體積模量構(gòu)建含氣飽和度指示因子Ig,Ig越大,則氣層產(chǎn)量Q0越高。

(12)

通過構(gòu)建的參數(shù)ΔK和Ig即可建立對(duì)應(yīng)的產(chǎn)能預(yù)測(cè)模型,準(zhǔn)確建模的前提是求準(zhǔn)干巖石骨架體積模量Kdry。

2.2.3 應(yīng)用Xu-White模型計(jì)算干巖石骨架的體積模量

Xu和White[13-14]基于Gassmann方程、Kuster-Toksoz理論以及有效差分介質(zhì)(DEM)理論提出了一種利用孔隙度和泥質(zhì)含量估算砂泥巖縱橫波速度的混合介質(zhì)模型,即Xu-White模型。該模型假定巖石的總孔隙空間由與砂巖顆粒和泥巖顆粒(包含束縛水)相關(guān)的孔隙組成(見圖5)。

圖5 Xu-White模型示意圖

Xu-White模型綜合考慮了巖石基質(zhì)性質(zhì)、孔隙大小、孔隙形狀以及孔隙流體性質(zhì)的影響,其本質(zhì)是認(rèn)為這2種孔隙空間的幾何形狀(孔隙縱橫比)差異較大,砂巖具有較大的孔隙縱橫比,而泥巖具有較小的孔隙縱橫比,因而對(duì)彈性模量的影響存在差異。孔隙空間按照礦物、黏土體積各自所占的百分比為

(13)

式中,φ為總孔隙度;φs為砂巖孔隙度;φc為泥質(zhì)(黏土)孔隙度;νs和νc分別為砂巖和泥巖體積百分比,其中νs=1-φ-νc。

為簡(jiǎn)化計(jì)算,Key和Xu[15-16]假定干巖石骨架的泊松比近似不隨孔隙度變化,將上述求解問題轉(zhuǎn)換成一個(gè)線性常微分方程組的求解,干巖石骨架彈性模量近似表達(dá)為

Kdry=Km(1-φ)p

(14)

(15)

2.3 基于陣列聲波測(cè)井的產(chǎn)能預(yù)測(cè)模型

應(yīng)用上述模型計(jì)算構(gòu)建的產(chǎn)能敏感參數(shù)ΔK與Ig所選取的一些參數(shù)及初始值見表3,其中巖石骨架的體積模量、骨架密度及骨架聲波為通過巖石物理實(shí)驗(yàn)獲得,氣體、泥巖和純水的體積模量、流體密度、流體聲波、e值及泥巖截止值為地區(qū)經(jīng)驗(yàn)值,基于陣列聲波測(cè)井資料,對(duì)產(chǎn)能敏感參數(shù)進(jìn)行求取。

表3 計(jì)算產(chǎn)能敏感參數(shù)時(shí)采用的參數(shù)及初值表

結(jié)合試氣資料,首先分析ΔK、Ig與Q0的相關(guān)關(guān)系如圖6、圖7,其相關(guān)關(guān)系較好,產(chǎn)能與ΔK為線性關(guān)系,與Ig為指數(shù)關(guān)系,其中Q0是通過單層試氣產(chǎn)量與單層有效厚度作對(duì)比獲得。

圖6 每米日產(chǎn)氣和體積模量差值交會(huì)圖

綜合2個(gè)敏感因子,從研究區(qū)17個(gè)有陣列聲波測(cè)井資料的試氣層中選取12個(gè)單層試氣層進(jìn)行建模,預(yù)測(cè)模型見表4。Ⅰ類儲(chǔ)層與Ⅱ類儲(chǔ)層預(yù)測(cè)模型都具有較高相關(guān)系數(shù),分別為0.94和0.95。

應(yīng)用上述預(yù)測(cè)模型,對(duì)研究區(qū)全部22個(gè)有陣列聲波測(cè)井資料的試氣層進(jìn)行預(yù)測(cè),合試氣層采用累加法求取見式(5),計(jì)算效果圖如圖8。實(shí)測(cè)產(chǎn)量與預(yù)測(cè)產(chǎn)量具有較好的相關(guān)性,相對(duì)于常規(guī)預(yù)測(cè)方法,誤差更小,精度更好。

圖8 實(shí)測(cè)產(chǎn)能與陣列聲波預(yù)測(cè)產(chǎn)能交會(huì)圖

圖9 A井產(chǎn)能預(yù)測(cè)效果圖

3 效果分析

分別應(yīng)用常規(guī)產(chǎn)能預(yù)測(cè)模型與基于陣列聲波資料的產(chǎn)能預(yù)測(cè)模型對(duì)研究區(qū)塊2口井資料進(jìn)行處理。圖9為A井2 529～2 560 m井氣層產(chǎn)能預(yù)測(cè)效果圖。該井段滲透率計(jì)算誤差較大,其原因?yàn)榭紫督Y(jié)構(gòu)復(fù)雜,呈雙峰或三峰結(jié)構(gòu);該井段巖屑含量較高,非均質(zhì)性強(qiáng),孔滲關(guān)系復(fù)雜,導(dǎo)致應(yīng)用孔滲回歸計(jì)算滲透率精度低。因此,基于常規(guī)測(cè)井參數(shù)的產(chǎn)能預(yù)測(cè)模型效果差。該井16號(hào)氣層與17號(hào)氣層為Ⅰ類儲(chǔ)層,常規(guī)預(yù)測(cè)模型誤差分別為40.3%;基于陣列聲波測(cè)井預(yù)測(cè)模型誤差分別為14.9%,準(zhǔn)確度大大提高,說明該方法能夠降低參數(shù)計(jì)算精度低導(dǎo)致產(chǎn)能預(yù)測(cè)偏差大的問題,在致密砂巖氣層有較好的應(yīng)用效果(見表5)。

圖10為B井2 530～2 600 m井段氣層產(chǎn)能預(yù)測(cè)效果圖。13號(hào)氣層與15號(hào)氣層為Ⅰ類儲(chǔ)層,16號(hào)干層與17號(hào)差氣層合試,分別為Ⅲ類儲(chǔ)層和Ⅱ類儲(chǔ)層,常規(guī)預(yù)測(cè)模型誤差分別為25.2%、26.8%和17.3%;基于陣列聲波測(cè)井預(yù)測(cè)模型誤差分別為13.8%、6.0%和4.7%,預(yù)測(cè)效果優(yōu)于常規(guī)方法(見表5)。

表5 不同方法產(chǎn)能預(yù)測(cè)結(jié)果對(duì)比表

圖10 B井產(chǎn)能預(yù)測(cè)效果圖

4 結(jié) 論

(1) 在儲(chǔ)層分類的基礎(chǔ)上,利用常規(guī)測(cè)井資料,從滲流理論出發(fā),建立產(chǎn)能預(yù)測(cè)模型,可實(shí)現(xiàn)產(chǎn)能快速預(yù)測(cè),但對(duì)于致密砂巖氣層,常規(guī)測(cè)井資料分辨率低,導(dǎo)致預(yù)測(cè)效果偏差較大。

(2) 基于巖石聲學(xué)物理實(shí)驗(yàn),研究了對(duì)儲(chǔ)層含氣性較為敏感的彈性參數(shù),提出基于陣列聲波測(cè)井資料的產(chǎn)能預(yù)測(cè)模型,大大提高了產(chǎn)能預(yù)測(cè)精度,在鄂爾多斯盆地東部上古生界致密砂巖氣層產(chǎn)能預(yù)測(cè)取得較好的應(yīng)用效果,并為含氣性儲(chǔ)層產(chǎn)能評(píng)價(jià)和預(yù)測(cè)的深入研究提供了一定的依據(jù)。

(3) 產(chǎn)能敏感參數(shù)的求取具有區(qū)域性和經(jīng)驗(yàn)性,結(jié)合研究區(qū)特征,探索了基于Xu-White模型和Biot-Gassmann方程的敏感參數(shù)計(jì)算方法,效果較好。

參考文獻(xiàn):

[1] 鄒才能, 陶士振, 袁選俊, 等. 連續(xù)型油氣藏形成條件與分布特征 [J]. 石油學(xué)報(bào), 2009, 30(3): 324-331.

[2] 楊華, 付金華, 劉新社, 等. 鄂爾多斯盆地上古生界致密氣成藏條件與勘探開發(fā) [J]. 石油勘探與開發(fā), 2012, 39(3): 295-303.

[3] 韓雪, 潘保芝. 利用測(cè)井資料預(yù)測(cè)儲(chǔ)層產(chǎn)能 [J]. 吉林大學(xué)學(xué)報(bào): 地球科學(xué)版, 2010, 40: 102-105.

[4] 孫建孟, 運(yùn)華云, 馮春珍. 測(cè)井產(chǎn)能預(yù)測(cè)方法與實(shí)例 [J]. 測(cè)井技術(shù), 2012, 36(6): 628-634.

[5] 譚成仟, 馬娜蕊, 蘇超. 儲(chǔ)層油氣產(chǎn)能的預(yù)測(cè)模型和方法 [J]. 地球科學(xué)與環(huán)境學(xué)報(bào), 2004, 26(2): 42-46.

[6] 成志剛, 張蕾, 趙建武, 等. 利用巖石聲學(xué)特性評(píng)價(jià)致密砂巖儲(chǔ)層含氣性 [J]. 測(cè)井技術(shù), 2013, 37(3): 253-257.

[7] 毛志強(qiáng), 李進(jìn)福. 油氣層產(chǎn)能預(yù)測(cè)方法及模型 [J]. 石油學(xué)報(bào), 2000, 21(5): 58-61.

[8] 張海濤, 石玉江, 張鵬, 等. 基于偶極橫波測(cè)井的低

滲透砂巖氣層識(shí)別方法 [J]. 測(cè)井技術(shù), 2015, 39(5): 591-595.

[9] 范宜仁, 邢東輝, 鄧少貴, 等. 低滲透巖石聲學(xué)特征及在含氣性預(yù)測(cè)中的應(yīng)用 [J]. 西南石油大學(xué)學(xué)報(bào): 自然科學(xué)版, 2015, 37: 64-70.

[10] Mavko G, Mukerji T, Dvorkin J. The Rock Physics Handbook: Tools for Seismic Analysis of Porous Media [M]. Cambridge: Cambridge University Press, 2009.

[11] Gassmann, F. Uber die Elastizitat Poroser Medien, Vierteljahrsschrift der Naturforschenden Gesellschaft in Zurich, 96: 1-23.

[12] Russell B H, Hedlin K, Hilterman F J, et al. Fluid-property Discrimination with AVO: A Biot-Gassmann Perspective [J]. Geophysics, 2003, 68(1): 29-39.

[13] Xu S, White R E. A New Velocity Model for Clay-sand Mixtures [J]. Geophysical Prospecting, 1995, 43(1): 91-118.

[14] Xu S, White R. Poro-elasticity of Clastic Rocks: A Unified Model [C]∥The SPWLA 36th Annual Logging Symposium, 1995.

[15] Key R G, Xu S. An Approximation for the Xu-White Velocity Model [J]. Geophysical, 1996, 67(5): 1406-1414.

[16] Biot M A. General Theory of Three-dimensional Consolidation [J]. Journal of Applied Physics, 1941, 12(2): 155-164.