缺失聲波條件下的頁巖儲層地應力測井解釋方法

時 賢 程遠方 常 鑫 蔣 恕 王 欣

1.中國石油大學（華東）石油工程學院 2.美國猶他大學能源與地球科學研究院 3.中國石油勘探開發研究院廊坊分院

北美頁巖氣商業化開發的巨大突破受益于水平分段壓裂等技術的快速進步，衍生出旨在創造復雜縫網的壓裂增產理論［1－3］。地應力是分析水力裂縫擴展規律，進行射孔方案及壓裂決策與設計的核心參數［4－7］。對頁巖儲層而言，因層理面發育等地質特征，造成頁巖層表現出較強的非均質性和各向異性，傳統地應力模型存在較多不適應性。國內外學者Hornby、Schatz、Cicotti、Xu、Barree等通過實驗和理論論證了各向異性引 起 的 地 層 應 力 改 變［8－14］。Jaeger、Deenadayalu、Khan、George、Higgins、Olson等則分析了頁巖各向異性力學參數在壓裂設計中的作用，并將儲層力學參數和地應力各向異性看作是復雜裂縫網絡形成的原因之一［15－21］。利用聲波測井數據是計算地應力的常用方法，較單點深度的室內巖石力學實驗而言更具應用價值，但因聲波測井費用昂貴，實際工程中進行雙極聲波測井的井數還不到總施工井數的1%，所以地應力計算數據資料十分有限［22］。

以上客觀條件為頁巖儲層的地應力解釋帶來巨大不便，所以需要在缺失聲波時差條件下，形成一套現場實用的頁巖層各向異性地應力解釋方法。基于權重分析思想，可參照臨井聲波數據，并通過中子、密度、伽馬等測井資料來反求未進行雙極測井油氣井的相關聲波數據，表現出較好應用價值［22－23］。筆者利用該方法，實現了對缺失聲波數據的回歸，同時基于該數據分析了某區塊頁巖儲層的各向異性，并分別通過H法（Higgnis）和B法（Blanton）計算了考慮各向異性的地應力大小。筆者旨在系統介紹一種缺失聲波測井數據條件下如何進行各向異性地應力計算的實用方法，更好的服務頁巖氣增產開發工作。

1 基于綜合測井資料的聲波回歸解釋方法

通過聲波測井解釋地應力避免了取心對地層資料原始性質破壞，連續性好、可靠性較高、被廣泛用于現場。通過現場應用發現，雖然一個區域內僅在幾口井進行了聲波測井，但這些聲波測井數據和臨井其他測井資料，如電阻率、伽馬、中子孔隙度和人工神經網絡回歸的孔隙度等其他數據存在一定關系，所以，可以以此作為基礎，實現對缺失聲波資料的油氣井聲波數據的反演解釋［22］。以縱波數據為例，如果存在傳導率、中子孔隙度和密度孔隙度等測井數據，則需要首先繪制傳導率、中子孔隙度和密度孔隙度等倒數與已有臨井縱波時差的交匯圖，再線性擬合出回歸方程，然后將各數據帶入回歸方程中，即可獲得初始聲波數據。為了進一步提高精度，需要對初始聲波時差與臨井聲波時差取平均值，再根據平均加權后獲得權重因子，從而建立起精度更好的鄰井縱波時差的解釋模型。同樣，橫波時差數據也可通過相同計算步驟獲取。該方法的優勢在于在有限聲波測井資料的條件下，可以擬合出高精度的聲波數據。并且，提出的權重因子可使回歸聲波速度結果避免井眼沖蝕和氣層存在的影響。下面以縱波時差為例介紹該模型的計算公式，在已知中子孔隙度，電阻率及人工神經網絡回歸孔隙度等測井資料的條件下，參照鄰井實測聲波時差（dvp）計算缺失井縱波時差（dvpcom）的表達式為：

2 頁巖各向異性成因及本構特征

巖石具有各向異性的原因被認為是地層沉積、壓實、膠結等復雜物理和化學過程綜合作用的結果，其中層理發育是主因。除了巖石層理結構的因素外，天然裂縫、構造運動、應力加載歷史和應力集中效應、孔隙壓力改變等也都有可能導致巖石各向異性的產生。

對于頁巖層而言，其各向異性通常表現在垂直和水平兩個平面上的差異。所以，可以通過橫觀各向同性模型（Transverse Isotropic）對此進行描述（圖1）。

圖1 橫觀各向同性模型中各向異性的單元體圖

該模型認為頁巖在水平方向是各向同性的，并以垂向和橫向上不同的彈性模量Eh和Ev對垂向的各向異性特征進行定量描述［10，24］。根據線性材料應變與應力本構關系，需要四階張量，共81個獨立分量組成的矩陣進行表達。橫觀各向同性模型從胡克定律出發，針對各向同性特點，應力應變關系僅用5個獨立的彈性常數進行描述，模型的本構形式如（3）所示。即

對于橫向各向同性材料而言，Cij可以表達為：

在該模型當中，僅用5個獨立變量C11，C12，C13，C55，C66就可以表達彈性模量（E）、泊松比（v）及剪切模量（G）。Ev與Eh為水平和垂直楊氏模量，其表達公式為式（5）和（6）［25］：

vh與vv為水平和垂直方向的泊松比，可由下式進行計算：

其中，C11、C12和C66存在線性關系，但是C33、C44、C66需要通過縱波，橫波及管流聲波速度分別進行計算，通常需要高級的聲波解釋工具來進行測量［21］。

通過上述測井聲波解釋方法，就實現了對各向異性巖石力學參數的測量。另外，實驗室也提供了計算各向異性彈性模量等靜力學參數的測量方法，但通常式樣大小有限，費用昂貴，所以常作為核心目的層測井解釋數據的校正。實驗測試方法如圖2所示。它需要在垂直，水平和角度為45°這3個方向進行分別測試，才可取得各向異性的力學參數結果［19］。

3 頁巖各向異性地應力解釋模型

橫觀各向同性地應力模型和常規地應力模型的差異主要在水平地應力的計算上，目前有兩種計算模式：①通過提出各向異性參數（K0）的Higgins法來對經典地應力模型的進行修正，簡稱H法；②考慮巖石巖性差異引起的應變不同，并耦合溫度等其他參數的Blanton法，簡稱 B法［25］。

圖2 考慮各向異性的巖石力學實驗室取心方法圖

3.1 Higgins地應力解釋方法

各向異性參數（K0）代表橫向和縱向的彈性模量及泊松比對巖石在平行和垂直層理面的兩個方向上的對巖石力學性質的不同影響程度，即

所以，根據各向異性參數（K0）建立的最小地應力解釋方程為：

根據式（10）可以看出，對于各向異性巖石，水平彈性模量增加會導致Eh／Ev項的增加，從而使最終計算水平最小主應力結果也增加。對于各向同性巖石，Eh／Ev值為1，則與常規地應力模型相同。

Higgins等基于應變理論提出一套新的地應力解釋方法，它主要將各向異性對構造應力的影響更加準確地考慮到方程之中，其計算公式為［9］：

目前，沒有理想的方法確定εH與εh的標定關系，一般在計算中多假設應變εH為水平最小主應力方向上的應變εh的兩倍。

3.2 Blanton地應力解釋方法

在進行最小水平主應力計算時，很多時候并沒有考慮構造應力，巖石線彈性和熱變量的影響。因此，Blanton等提出一種綜合系數校正方法來計算的最小水平主應力［8］：

其中構造應變的大小需要通過對目的層進行相關測試來獲取，熱擴展系數（εT）在巖石拉伸時取負值，而壓縮時取正值。通常，αT的取值和巖性密切相關，對于砂巖來說，αT一般為5.56×10－6／F，頁巖為5.00×10－6／F，而碳酸鹽大約為4.44×10－6／F。對于ΔT而言，則可以通過地區地熱梯度和井底溫度進行反推。

圖3 MB井橫、縱聲波時差與中子孔隙度和密度孔隙度的交匯圖

圖4 MC井橫、縱聲波時差與中子孔隙度和密度孔隙度的交匯圖

4 實例分析及模型應用效果分析

以美國某頁巖區塊為例介紹在缺失聲波資料的情況下如何回歸巖石力學參數和計算地應力的基本過程，并分析該地區的各向異性程度。該頁巖層位于3 200～3 900m的白堊系 Mancos地層，屬于海相沉積頁巖，油氣藏上部由粉砂巖，泥巖及頁巖夾層組成，中間有大段較純的黑色頁巖層，其中頁巖有機質含量介于3%～10%，孔隙度介于4%～8%，滲透率介于（5～20）×10－6mD，屬于極低滲透的頁巖氣藏。所選區域目前共有3口頁巖氣井 MA、MB、MC，其中僅有1口氣井MA具有裸眼的雙極聲波測井數據，井MB和MC同時具有中子孔隙度和密度孔隙度測井資料，MC另外具有自然伽馬資料，MB具有部分電阻率測井資料。由于只有中子孔隙度和密度孔隙度測井資料。所以，首先將MB和MC的中子孔隙度，密度孔隙度測井資料與井MA的縱波和橫波數據分別進行回歸解釋。圖3（a～d）分別表示MB井橫波時差和縱波時差分別和密度孔隙度、中子孔隙度的交匯圖，圖4（a～d）分別表示MC井橫波時差和縱波時差分別和密度孔隙度、中子孔隙度的交匯圖。

在獲取上述交匯圖后，可以回歸出相應的聲波時差。根據加權平均的思想，通過計算求得MB和MC井縱橫波中子孔隙度和密度孔隙度權重因子（表1）。

表1 MB井及MC井權重因子統計表

圖5 MB井與MC井垂直與水平彈性模量與泊松比對比圖

表2 定井深條件下各向異性系數與巖性解釋結果對比表

根據已知A井聲波時差資料，可以回歸出MB井及MC井的聲波數據，然后結合測井解釋的其他各向異性參數，可以對動態楊氏模量和泊松比進行分別計算（圖5）。表2是通過QEMSCAN的巖性解釋結果與巖石力學參數的對比。

運用美國猶他大學的QEMSCAN技術對取心處巖樣的巖石礦物組分進行了定量分析，顯示了礦物的具體分布結構，實驗結果發現所取巖樣的膠結物成分主要為伊利石、高嶺石、石英以及伊／蒙混層等，圖6－a中顯示的黃色區域為黃鐵礦層，并清晰顯示出頁巖的層理，圖6－b中的層理較不發育，與定井深點的各向異性參數比較如表2所示。通過巖石力學參數的解釋結果比較發現，石英含量較高，脆性指數較高的頁巖層具有較低的泊松比，而黏土含量高的頁巖泊松比較高，主要由其軸向和橫向的變形差異較小導致。通過回歸的垂直與水平楊氏模量和泊松比的曲線發現，整體上看，水平彈性模量要略高于垂直彈性模量，垂直彈性模量更符合垂向的巖性變化，并且頁巖層的各向異性相對較為突出。對于泊松比而言，垂向泊松比值比橫向泊松比略大，但整體差異較小。

根據解釋的力學參數結果，分別應用傳統方法和兩種考慮頁巖各向異性的方法對地應力進行計算分析，H法主要根據公式（11），而B法則根據公式（14）對最小水平主應力進行計算，計算結果如圖7所示。構造應力的計算需要通過小型壓裂獲取的閉合應力結果進行校正，同時構造應力的大小也和巖性相關，需要根據層位進行分別標定。所以，對于B法中的應變值，一般隨巖性不同而變化，對于砂巖來說，一般取0.000 27，而頁巖一般取0.000 24。

圖6 高精度（微米級精度）QEMSCAN掃描結果圖

圖7 考慮各向異性的MB井和MC井地應力剖面與常規方法地應力剖面關系圖

表3 小型壓裂數據與C井回歸地應力數據對比表

將所求地層最小水平主應力和小型壓裂閉合應力數據進行比較，并做誤差分析，結果如表3所示。通過應力結果比較發現，H法和B法預測的結果都比常規方法的地應力值要高，其中B法的值最大，主要在于模型中泊松比各向系數最小。H法和B法計算地應力值較高的原因認為是模型中考慮了不同巖性導致的構造應力差異。通常在構造板塊活躍的區域，構造應力的作用忽略會導致計算的最小水平主應力的增加。對于常規方法計算的應力剖面而言，可以看出，如果考慮構造應力項，只能在最終應力剖面上增加一條平滑的曲線，忽略了地層的差異性。考慮各向異性的方法由于將巖石力學性質改變導致的應力改變施加到不同的井深當中，所以結果將更加精確。從壓裂施工的角度來看，考慮各向異性計算出最小水平主應力值偏高。因此將更加有利射孔層位的選擇及對水力裂縫縫高延伸的約束。

5 結論

1）通過綜合測井資料和臨井聲波資料，可以實現缺失聲波數據的反演，大大擴展了聲波數據的利用廣度。從技術操作上看，該種方法應用簡單且精確度較高，可為力學參數與地應力的計算提供技術參數。

2）進行考慮各向異性頁巖的力學參數解釋時，需要分別計算垂向和橫向楊氏模量值和泊松比值。一般來說，水平楊氏模量值比垂直彈性模量值高，并且垂直楊氏模量值更符合巖性規律，而垂向與橫向泊松比值相差不多，但垂向泊松比值要比橫向泊松比值略大。

3）由于傳統地應力模型在構造應力的處理上僅僅做了簡單的處理，所以附加在整個深度總應力剖面上的構造應力為平滑直線，無法體現出巖石力學性質對應力的影響作用，同時也不利于區分各層的應力差異。橫觀各向同性地應力模型考慮了巖石力學性質改變對每個深度上地應力的影響，其中泊松比值越小則最小主應力越大，其中B法預測的最小地應力最大正是因為泊松比系數較低所導致，H法主要受橫向泊松比影響，所以計算的最小水平主地應力較大。使用橫觀各向同性模型將更有助于準確計算設計的裂縫幾何形態，但對于各向同性地層，橫觀各向同性模型和常規模型的結果則基本相同。另外，通過成像測井等技術可對所求地應力結果做進一步驗證。

致謝：本研究得到了美國猶他大學能源與地學研究院（EGI）提供的軟件及資料支持，另外對中國留學基金委（CSC）的海外留學資助也表示衷心感謝。

符 號 說 明

Ed為動態楊氏模量，MPa；vd為泊松比；vp為縱波速度，m／μs；vs為橫波速度，m／μs；ρ 為密度，g／cm3；Z 為地層深度，m；ρ（z）為上覆巖層密度函數，kg／m3；g 為重力加速度，m／s2；σv為上覆巖層壓力，MPa；pp為孔隙壓力，MPa；α為Biot系數，無量綱；αT為構造應力，MPa；vs為巖石靜態泊松比，無量綱；dvphin為中子孔隙度回歸聲波時差，μs／m；dvp為實測聲波時差，μs／m；W1～W4為縱波聲波時差權重因子，無量綱；dvpcom為中子孔隙度回歸聲波時差，μs／m；dvs為實測聲波時差，μs／m；M1～M4為橫波聲波時差權重因子，無量綱；σij、εkl為二階應力與應變張量；Cijkl為四階剛度張量；ξ為線彈性常數；εH與εh為最大水平主應力和最小水平主應力方向的應變；C1為平面應變模量；C2為巖石楊氏模量。

［1］賈承造，鄭民，張永峰.中國非常規油氣資源與勘探開發前景［J］.石油勘探與開發，2012，39（2）：129－136.JIA Chengzao，ZHENG Min，ZHANG Yongfeng.Unconventional hydrocarbon resources in China and the prospect of exploration and development［J］.Petroleum Exploration and Development，2012，39（2）：129－136.

［2］林臘梅，張金川，唐玄，等.中國陸相頁巖氣的形成條件［J］.天然氣工業，2013，33（1）：35－39.LIN Lamei，ZHANG Jinchuan，TANG Xuan，et al.Conditions of continental shale gas accumulation in China［J］.Natural Gas Industry，2013，33（1）：35－39.

［3］程遠方，李友志，時賢，等.頁巖氣體積壓裂縫網模型分析及應用［J］.天然氣工業，2014，34（9）：53－58.CHENG Yuanfang，LI Youzhi，SHI Xian，et al.Analysis and application of fracture network models of volume fracturing in shale gas reservoirs［J］.Natural Gas Industry，2014，34（9）：53－58.

［4］胡永全，賈鎖剛，趙金洲，等.縫網壓裂控制條件研究［J］.西南石油大學學報：自然科學版，2013，35（4）：126－132.HU Yongquan，JIA Suogang，ZHAO Jinzhou，et al.Study on controlling conditions in network hydraulic fracturing［J］.Journal of Southwest Petroleum University：Science ＆Technology Edition，2013，35（4）：126－132.

［5］李勇明，彭瑤，王中澤.頁巖氣壓裂增產機理與施工技術分析［J］.西南石油大學學報：自然科學版，2013，35（2）：90－96.LI Yongming，PENG Yao，WANG Zhongze.Analysis of shale gas fracture stimulation mechanism and operating techniques［J］.Journal of Southwest Petroleum University：Science ＆ Technology Edition，2013，35（2）：90－96.

［6］趙金洲，任嵐，胡永全.頁巖儲層壓裂縫成網延伸的受控因素分析［J］.西南石油大學學報：自然科學版，2013，35（1）：1－9.ZHAO Jinzhou，REN Lan，HU Yongquan.Controlling factors of hydraulic fractures extending into network in shale formations［J］.Journal of Southwest Petroleum University：Science ＆ Technology Edition，2013，35（1）：1－9.

［7］王香增，吳金橋，張軍濤.陸相頁巖氣層的CO2壓裂技術應用探討［J］.天然氣工業，2014，34（1）：64－67.WANG Xiangzeng，WU Jinqiao，ZHANG Juntao.Application of CO2fracturing technology for terrestrial shale gas reservoirs［J］.Natural Gas Industry，2014，34（1）：64－67.

［8］BLANTON T L，OLSON J E.Stress magnitudes from logs：Effects of tectonic strains and temperature［J］.SPE Reservoir Evaluation ＆ Engineering，1999，2（1）：62－68.

［9］HIGGINS S M，GOODWIN S A，BRATTON T R，et al.Anisotropic stress models improve completion design in the Baxter Shale［C］∥paper 115736－MS presented at the SPE Annual Technical Conference and Exhibition，21－24September 2008，Denver，Colorado，USA.New York：SPE，2008.

［10］SCHOENBERG M，MUIR F，SAYERS C.Introducing ANNIE：A simple and three－parameter anisotropic velocity model for shales［J］.Journal of Seismic Exploration，1996，5（1）：35－50.

［11］VERNIK L，LIU X.Velocity anisotropy in shales：A petrophysical study［J］.Exploration Geophysicists，1997，62（2）：521－532.

［12］BULLER D，HUGHES S，MARKET J.Petrophysical evaluation for enhancing hydraulic stimulation in horizontal shale gas wells［C］∥paper 132990－MS presented at the SPE Annual Technical Conference and Exhibition，19－22 September 2010，Florence，Italy.New York：SPE，2010.

［13］FRANQUET J A，RODRIGUEZ E F.Orthotropic horizontal stress characterization from logging and core derived acoustic anisotropies［C］∥paper ARMA－2012－644 presented at the 46thU.S.Rock Mechanics／Geomechanics Symposium，24－27June 2012，Chicago，Illinois，USA.New York：SPE，2012.

［14］RAJAPPA B，CONLET J T.In－situ sress profiling in the Agency Draw Field，Utah［C］∥paper 107818－MS presented at the Rocky Mountain Oil ＆ Gas Technology Symposium，16－18April 2007，Colorado，USA.New York：SPE 2007.

［15］SONG L，HARELAND G.Minimum horizontal stress profile from logging data for Montney Formation of North East British Columbia［C］∥paper 162233－MS presented at the SPE Canadian Unconventional Resources Conference，30October －1November 2012，Calgary，Alberta，Canada.New York：SPE，2012.

［16］CIPPLLA C L，LIU Diabin，KYTE D G.Practical application of in－situ stress profiles［C］∥paper 28607－MS presented at the SPE Annual Technical Conference and Exhibition，25－28September 1994，New Orleans，Louisiana，USA.New York：SPE，1994.

［17］KHAN S，ANSARI S A，HAN Hongxue，et al.Importance of shale anisotropy in estimating in－situ stresses and wellbore stability analysis in Horn River Basin［C］∥paper 149433－MS presented at the Canadian Unconventional Resources Conference，15－17November 2001，Alberta，Canada.New York：SPE，2011.

［18］WATERS G A，LEWIS R E，BENTLEY D.The effect of mechanical properties anisotropy in the generation of hydraulic Fractures in organic shales［C］∥paper 146776－MS presented at the SPE Annual Technical Conference and Exhibition，30October －2November 2011，Denver，Colorado，USA.New York：SPE，2011.

［19］BARREE R D，GILBERT J V，CONWAYA M W.Stress and rock property profiling for unconventional reservoir stimulation［C］∥paper 118703－MS presented at the SPE Hydraulic Fracturing Technology Conference，19－21January 2009，the Woodlands，Texas，USA.New York：SPE，2009.

［20］LARKIN S D，BROWN E K，BAZAN L W.Stimulation design and post fracture production analysis：A tight gas sand case history［C］∥paper 74361－MS presented at the SPE International Petroleum Conference and Exhibition in Mexico，10－12February 2002，Villahermosa，Mexico.New York：SPE，2002.

［21］WALSH J，SINHA B，DONALD A.Formation anisotropy parameters using borehole sonic data［C］∥paper 2006－TT presented at the SPWLA 47thAnnual Logging Symposium，4－7June 2006，Veracruz，Mexico.

［22］MULLEN M J，ROUNDTREE R，TURK G A.A composite determination of mechanical rock properties for stimulation design（what to do when you don＇t have a sonic log）［C］∥paper 108139－MS presented at the Rocky Mountain Oil ＆ Gas Technology Symposium，16－18April 2007，Denver，Colorado，USA.New York：SPE，2007.

［23］RICKMA R，MULLEN M J.A practical use of shale petrophysics for stimulation design optimization：All shale plays are not clones of the Barnett Shale［C］∥paper 115258－MS presented at the SPE Annual Technical Conference and Exhibition，21－24September 2008，Denver，Colorado，USA.New York：SPE，2008.

［24］AMADEI B.Importance of anisotropy when estimating and measuring in situ stresses in rock.［J］.International Journal of Rock Mechanics and Mining Sciences ＆ Geomechanics，1996，33（3）：293－325.

［25］AMADEI B，SAVAGE W Z，SWOLFS H S.Gravitational stresses in anisotropic rock masses［J］.International Journal of Rock Mechanics and Mining Sciences ＆ Geomechanics，1987，24（1）：5－14.