本發明屬于頁巖氣開采技術領域，尤其涉及一種頁巖儲層地應力預測建模方法。
背景技術：
目前，常規能源的開發和利用已逐漸不能滿足現代工業的發展需求，勘探開發的重點投向了難以開采的頁巖氣。頁巖是一種沉積巖，成分復雜，但都具有薄頁狀或薄片狀的節理，主要是由黏土沉積經壓力和溫度形成的巖石，但其中混雜有石英、長石的碎屑以及其他礦物。頁巖儲層是一種低孔隙度和特低滲透率的致密儲層，其生產完全依賴于壓裂效果。對于頁巖進行壓裂施工必須對頁巖起裂的規律有足夠的認識，其中地應力模型是指導頁巖壓裂施工的極為重要的理論基礎。目前地應力預測模型大部分建立在地層均質各項同性的基礎上，然而頁巖儲層具有較強的非均質性和各項性。因此，其計算結果并不十分準確，對地應力的描述并不全面，不能完全的描述橫觀各向同性地層的地應力的特征，且建模的過程中，智能化程度較低，對探測孔數據的讀取與測量多為工作人員操作，采集數據種類較少，數據不全面，人為誤差大，效率較低，不利于工程施工進度。
因此，如何能夠對頁巖儲層地應力進行有效地預測分析，以及如何提高預測準確率和預測效率是現今應解決的問題。
技術實現要素：
本發明的目的在于提供一種頁巖儲層地應力預測建模方法，旨在解決目前地應力預測模型大部分建立在地層均質各項同性的基礎上，然而頁巖儲層具有較強的非均質性和各項性。因此，其計算結果并不十分準確，對地應力的描述并不全面，不能完全的描述橫觀各向同性地層的地應力的特征，且建模的過程中，智能化程度較低，對探測孔數據的讀取與測量多為工作人員操作，采集數據種類較少，數據不全面，人為誤差大，效率較低，不利于工程施工進度的問題。
本發明是這樣實現的，一種頁巖儲層地應力預測建模系統，包括單片機控制器和綜合數據采集器，所述單片機控制器分別與綜合數據采集器、數據處理器、RAM存儲器、ROM存儲器、數據庫、計時器和無線射頻收發器電性連接；
所述單片機控制器的輸入端與供電電源的輸出端電性連接，所述單片機控制器的輸出端與應力模擬器的輸入端電性連接；
所述綜合數據采集器的輸入端分別與第一數據采集器、第二數據采集器、第三數據采集器和第四數據采集器的輸出端電性連接；
所述無線射頻收發器通過GPRS網絡與外部設備連接，所述第一數據采集器的輸入端與距離傳感器的輸出端電性連接；
所述第二數據采集器的輸入端與孔隙壓力傳感器的輸出端電性連接，所述第三數據采集器的輸入端與應力傳感器的輸出端電性連接；
所述第四數據采集器的輸入端與應變式傳感器的輸出端電性連接。
所述距離傳感器、孔隙壓力傳感器、應力傳感器和應變式傳感器均至少為一個，所述單片機控制器的輸入端與GPS定位器的輸出端電性連接；
所述頁巖儲層地應力預測建模系統的頁巖儲層地應力預測建模方法包括：
步驟一、測量油氣藏的巖石力學參數，統計所述油氣藏的巖體力學性質影響參數；
步驟二、將頁巖試件放在抽真空飽和裝置中，利用飽和液體進行飽和；
步驟三、在所述頁巖試件上安裝軸向變形傳感器和徑向變形傳感器，并將所述巖石試件放入巖石三軸試驗機的三軸室內，對所述巖體試樣施加與所述巖體試樣所處地層深度相對應的圍壓；
步驟四、按照多個溫度值由高到低的順序，依次將高壓腔內的溫度降低到多個溫度值下，按照預設加載速率對巖體試樣加載軸壓，得到巖體試樣在各溫度值下的力學參數；
步驟五、通過測井復雜巖性分析程序進行測井解釋，獲得礦物體積、巖石孔隙度、流體飽和度、干酪根體積，進而求得礦物在巖石中的體積分數，所述礦物包括粘土、石英、方解石；使用交會分析獲得每個礦物的彈性模量和密度，使用已知的經驗值作為孔隙中流體的彈性模量和密度，所述彈性模量指體積模量、剪切模量；
步驟六、取2倍洞徑以外典型的原巖巖餅至少3塊，要求各巖餅的孔深相差在1m以內，巖餅形態類似，記錄每個巖餅發生處的埋深和孔深；
步驟七、確定每個巖餅三個主應力的方向，做巖石的巴西圓盤劈裂試驗，確定每個巖餅的巖石抗拉強度St，確定每個巖餅三個主應力的大小；
步驟八、從陣列聲波測井的波形資料中提取縱波時差Δtc和橫波時差Δts；
步驟九、根據縱波時差Δtc和橫波時差Δts，結合陣列聲波測井資料中的密度ρ，計算動態楊氏模量Ed和動態泊松比PRd；
步驟十、統計所述油氣藏巖體力學性質影響參數，包括內部影響參數和外部影響參數，內部影響參數包括巖石成分特征參數，巖石結構特征參數，巖石構造特征參數以及巖石孔隙裂縫特征參數，外部影響參數包括溫度特征參數，圍壓特征參數，孔隙流體特征參數以及孔隙壓力特征參數；
步驟十一、根據巖石力學參數采集的數據以及統計的巖體力學性質影響參數建立油氣藏的頁巖儲層地應力物理模型；
巖體力學性質影響參數包括內部影響參數和外部影響參數；所述內部影響參數包括巖石成分特征參數，巖石結構特征參數，巖石構造特征參數以及巖石孔隙裂縫特征參數；所述外部影響參數包括溫度特征參數，圍壓特征參數，孔隙流體特征參數以及孔隙壓力特征參數；
按照如下公式計算巖石摩擦系數：
其中，σ1為殘余強度；
σ2為圍壓；
p0為孔隙壓力；
σ為裂縫面上的正應力；
τ為裂縫面上的剪應力；
β為裂縫面和端面的夾角；
μ為巖石摩擦系數。
巖石力學參數包括巖石彈性參數，所述巖石彈性參數包括楊氏模量、泊松比和Biot系數，所述Biot系數用于表征孔隙彈性參數；所述巖石力學參數的各向異性強度包括楊氏模量各向異性強度、泊松比各向異性強度和Biot系數各向異性強度；
所述巖石力學參數的各向異性強度通過如下公式計算獲得：
其中，Ek為楊氏模量各向異性強度；Eh為水平方向楊氏模量；Ev為垂直方向楊氏模量；μk為泊松比各向異性強度；μh為水平方向泊松比；μv為垂直方向泊松比；
算動態楊氏模量Ed，按照如下公式進行計算：
式中：
Ed—動態楊氏模量，
Δts—橫波時差，
Δtc—縱波時差，
ρ—密度。
進一步，所述無線射頻收發器設置有信號相位計算模塊，所述相位計算模塊的相位計算方法包括：
第k通道第n路信號Sk.n(t)的表達式為：
Sk.n(t)＝expj{ω0(t+τk+nTs)+1/2μ(t+τk+nTs)2},k＝0,1,...；n＝0,1,2,...
其中，ω0是輸出波形的初始角速度，t是時間，n表示每通道中路數的序號，μ表示調頻斜率，Ts是采樣周期，τk表示第k通道信號起始相位對應的時間差；
第k通道第n路信號的相位作如下變換：
令Pk.n＝μ(τk+nTs)、Qk.n＝ω0(τk+nTs)+1/2μ(τk+nTs)2，則上式可化為：
其中為基準相位，當信號屬性參數固定時，Pk.n、Qk.n為定值；其他路的波形輸出均看做在基本相位的基礎上增加一個偏移相位得到。
進一步，所述綜合數據采集器設置有聯合估計單元，所述聯合估計單元的聯合估計方法包括：目標信號DOA估計、波達信號時延TD估計和多譜勒頻率fd估計三者聯合估計方法；
首先，對將陣列天線的接收到的能量衰減得非常厲害的非目標反射信號，經過小波分析的分解、自適應閾值選擇以及重構從而檢測和提取弱信號，所述小波分析將接收到的弱信號分成高頻部分和低頻部分，通過多層分解，重構時將每層的高頻部分按照不同的閾值進行取舍，超出閾值的部分取閾值，低于閾值的部分取實際的值；
第二步，對多普勒頻率進行估計：通過前期的弱信號檢測，對接收到的提高了信號的信噪比的目標反射信號中，含有的在頻率上有差別的多譜勒頻率信號，通過本地的正弦和余弦信號進行同步輸出并經過同步檢測電路接收信號進行正交運算，提取多普勒頻率；
具體為：經過正交推導以及電路濾波，多普勒頻率部分公式如(1)，
式中，fd＝(ω-ω0)/2π，包含多普勒頻率成分；
式中m表示信號在整個電路系統傳輸過程中的衰減因子，βj表示在空氣中的傳輸因子，ψ(θj)表示含有目標信號方向角的函數，ω表示接收到的含有多普勒頻率分量的目標反射信號的載波信號頻率，ω0表示載波頻率；
第三步，對DOA和時延進行聯合提取：利用已經提取的多普勒頻率以及ESPRIT算法中多普勒頻率與DOA、信號時延的聯合性進行矩陣變換和計算，從而提取出DOA和波達信號時延；
具體為：聯合估計DOA和TD中，利用的ESPRIT算法，最后推導公式如(2)，
Rxx-γRxy＝APAH-γAPφHAH＝AP(I-γφH)AH (2)
公式中Rxx表示陣列自協方差，Rxy表陣列互協方差，γ表示陣列信號對應的特征值，I表示單位矩陣，A表示方向矩陣，AH表示方向矩陣的希爾伯特矩陣，P表示信源部分的協方差矩陣，φH表示含有多普勒頻率成分的方向矩陣；
通過得到的多普勒頻率，然后利用ESPRIT算法進行參數配對即可得到DOA和TD。
本發明提供的頁巖儲層地應力預測建模方法，智能化程度較高，外部設備通過GPRS網絡可對單片機控制器遠程控制，實現了操作的便利性，多種傳感器可對探測孔的各項數據進行詳細的采集，單片機控制器可對采集的數據綜合處理分析，通過應力模擬器可對頁巖儲層的應力情況直觀的模擬，分析結果可視化，與傳統預測建模方法比較，避免了人為誤差，數據采集具有多樣性，保證了應力預測的準確率，極大地提高了預測效率，有助于加快工程進度。
附圖說明
圖1是本發明實施例提供的頁巖儲層地應力預測建模系統結構示意圖；
圖中：1、單片機控制器；2、綜合數據采集器；3、數據處理器；4、RAM存儲器；5、ROM存儲器；6、數據庫；7、計時器；8、無線射頻收發器；9、供電電源；10、應力模擬器；11、第一數據采集器；12、第二數據采集器；13、第三數據采集器；14、第四數據采集器；15、GPRS網絡；16、外部設備；17、距離傳感器；18、孔隙壓力傳感器；19、應力傳感器；20、應變式傳感器；21、GPS定位器。
具體實施方式
為了使本發明的目的、技術方案及優點更加清楚明白，以下結合實施例，對本發明進行進一步詳細說明。應當理解，此處所描述的具體實施例僅僅用以解釋本發明，并不用于限定本發明。
下面結合附圖對本發明的應用原理作詳細的描述。
如圖1所示，本發明實施例提供的頁巖儲層地應力預測建模系統包括單片機控制器1和綜合數據采集器2，所述單片機控制器1分別與綜合數據采集器2、數據處理器3、RAM存儲器4、ROM存儲器5、數據庫6、計時器7和無線射頻收發器8電性連接，所述單片機控制器1的輸入端與供電電源9的輸出端電性連接，所述單片機控制器1的輸出端與應力模擬器10的輸入端電性連接，所述綜合數據采集器2的輸入端分別與第一數據采集器11、第二數據采集器12、第三數據采集器13和第四數據采集器14的輸出端電性連接，所述無線射頻收發器8通過GPRS網絡15與外部設備16連接，所述第一數據采集器11的輸入端與距離傳感器17的輸出端電性連接，所述第二數據采集器12的輸入端與孔隙壓力傳感器18的輸出端電性連接，所述第三數據采集器13的輸入端與應力傳感器19的輸出端電性連接，所述第四數據采集器14的輸入端與應變式傳感器20的輸出端電性連接。
進一步，所述外部設備16為電腦、手機等具有網絡連接功能的電子產品。
進一步，所述距離傳感器17、孔隙壓力傳感器18、應力傳感器19和應變式傳感器20均至少為一個。
進一步，所述單片機控制器1的輸入端與GPS定位器21的輸出端電性連接。
進一步，所述無線射頻收發器設置有信號相位計算模塊，所述相位計算模塊的相位計算方法包括：
第k通道第n路信號Sk.n(t)的表達式為：
Sk.n(t)＝expj{ω0(t+τk+nTs)+1/2μ(t+τk+nTs)2},k＝0,1,...；n＝0,1,2,...
其中，ω0是輸出波形的初始角速度，t是時間，n表示每通道中路數的序號，μ表示調頻斜率，Ts是采樣周期，τk表示第k通道信號起始相位對應的時間差；
第k通道第n路信號的相位作如下變換：
令Pk.n＝μ(τk+nTs)、Qk.n＝ω0(τk+nTs)+1/2μ(τk+nTs)2，則上式可化為：
其中為基準相位，當信號屬性參數固定時，Pk.n、Qk.n為定值；其他路的波形輸出均看做在基本相位的基礎上增加一個偏移相位得到。
進一步，所述綜合數據采集器設置有聯合估計單元，所述聯合估計單元的聯合估計方法包括：目標信號DOA估計、波達信號時延TD估計和多譜勒頻率fd估計三者聯合估計方法；
首先，對將陣列天線的接收到的能量衰減得非常厲害的非目標反射信號，經過小波分析的分解、自適應閾值選擇以及重構從而檢測和提取弱信號，所述小波分析將接收到的弱信號分成高頻部分和低頻部分，通過多層分解，重構時將每層的高頻部分按照不同的閾值進行取舍，超出閾值的部分取閾值，低于閾值的部分取實際的值；
第二步，對多普勒頻率進行估計：通過前期的弱信號檢測，對接收到的提高了信號的信噪比的目標反射信號中，含有的在頻率上有差別的多譜勒頻率信號，通過本地的正弦和余弦信號進行同步輸出并經過同步檢測電路接收信號進行正交運算，提取多普勒頻率；
具體為：經過正交推導以及電路濾波，多普勒頻率部分公式如(1)，
式中，fd＝(ω-ω0)/2π，包含多普勒頻率成分；
式中m表示信號在整個電路系統傳輸過程中的衰減因子，βj表示在空氣中的傳輸因子，ψ(θj)表示含有目標信號方向角的函數，ω表示接收到的含有多普勒頻率分量的目標反射信號的載波信號頻率，ω0表示載波頻率；
第三步，對DOA和時延進行聯合提取：利用已經提取的多普勒頻率以及ESPRIT算法中多普勒頻率與DOA、信號時延的聯合性進行矩陣變換和計算，從而提取出DOA和波達信號時延；
具體為：聯合估計DOA和TD中，利用的ESPRIT算法，最后推導公式如(2)，
Rxx-γRxy＝APAH-γAPφHAH＝AP(I-γφH)AH (2)
公式中Rxx表示陣列自協方差，Rxy表陣列互協方差，γ表示陣列信號對應的特征值，I表示單位矩陣，A表示方向矩陣，AH表示方向矩陣的希爾伯特矩陣，P表示信源部分的協方差矩陣，φH表示含有多普勒頻率成分的方向矩陣；
通過得到的多普勒頻率，然后利用ESPRIT算法進行參數配對即可得到DOA和TD。
本發明實施例提供的頁巖儲層地應力預測建模系統，該系統通過以下方法實現頁巖儲層地應力預測：
步驟一、測量油氣藏的巖石力學參數，統計所述油氣藏的巖體力學性質影響參數；
步驟二、將頁巖試件放在抽真空飽和裝置中，利用飽和液體進行飽和；
步驟三、在所述頁巖試件上安裝軸向變形傳感器和徑向變形傳感器，并將所述巖石試件放入巖石三軸試驗機的三軸室內，對所述巖體試樣施加與所述巖體試樣所處地層深度相對應的圍壓；
步驟四、按照多個溫度值由高到低的順序，依次將高壓腔內的溫度降低到多個溫度值下，按照預設加載速率對巖體試樣加載軸壓，得到巖體試樣在各溫度值下的力學參數；
步驟五、通過測井復雜巖性分析程序進行測井解釋，獲得礦物體積、巖石孔隙度、流體飽和度、干酪根體積，進而求得礦物在巖石中的體積分數，所述礦物包括粘土、石英、方解石；使用交會分析獲得每個礦物的彈性模量和密度，使用已知文獻中的經驗值作為孔隙中流體的彈性模量和密度，所述彈性模量指體積模量、剪切模量；
步驟六、取2倍洞徑以外典型的原巖巖餅至少3塊，要求各巖餅的孔深相差在1m以內，巖餅形態類似，記錄每個巖餅發生處的埋深和孔深；
步驟七、確定每個巖餅三個主應力的方向，做巖石的巴西圓盤劈裂試驗，確定每個巖餅的巖石抗拉強度St，確定每個巖餅三個主應力的大小；
步驟八、從陣列聲波測井的波形資料中提取縱波時差Δtc和橫波時差Δts；
步驟九、根據縱波時差Δtc和橫波時差Δts，結合陣列聲波測井資料中的密度ρ，計算動態楊氏模量Ed和動態泊松比PRd；
步驟十、統計所述油氣藏巖體力學性質影響參數，包括內部影響參數和外部影響參數，內部影響參數包括巖石成分特征參數，巖石結構特征參數，巖石構造特征參數以及巖石孔隙裂縫特征參數，外部影響參數包括溫度特征參數，圍壓特征參數，孔隙流體特征參數以及孔隙壓力特征參數；
步驟十一、根據巖石力學參數采集的數據以及統計的巖體力學性質影響參數建立油氣藏的頁巖儲層地應力物理模型。
進一步，巖體力學性質影響參數包括內部影響參數和外部影響參數；所述內部影響參數包括巖石成分特征參數，巖石結構特征參數，巖石構造特征參數以及巖石孔隙裂縫特征參數；所述外部影響參數包括溫度特征參數，圍壓特征參數，孔隙流體特征參數以及孔隙壓力特征參數。
進一步，按照如下公式計算巖石摩擦系數：
其中，σ1為殘余強度；
σ2為圍壓；
p0為孔隙壓力；
σ為裂縫面上的正應力；
τ為裂縫面上的剪應力；
β為裂縫面和端面的夾角；
μ為巖石摩擦系數。
進一步，巖石力學參數包括巖石彈性參數，所述巖石彈性參數包括楊氏模量、泊松比和Biot系數，所述Biot系數用于表征孔隙彈性參數；所述巖石力學參數的各向異性強度包括楊氏模量各向異性強度、泊松比各向異性強度和Biot系數各向異性強度；
所述巖石力學參數的各向異性強度通過如下公式計算獲得：
其中，Ek為楊氏模量各向異性強度；Eh為水平方向楊氏模量；Ev為垂直方向楊氏模量；
μk為泊松比各向異性強度；μh為水平方向泊松比；μv為垂直方向泊松比。
進一步，算動態楊氏模量Ed，按照如下公式進行計算：
式中：
Ed—動態楊氏模量；
Δts—橫波時差；
Δtc—縱波時差；
ρ—密度。
工作原理：該頁巖儲層地應力預測建模方法，外部設備16通過GPRS網絡15可對單片機控制器1遠程控制，RAM存儲器4可對數據臨時存儲，ROM存儲器5可對數據永久存儲，便于數據調用，數據庫6內存儲有各種應力變化數據，計時器7可對數據采集時間以及預測建模的時間進行詳細的記錄，無線射頻收發器8可接收和發送無線信號，供電電源9可為系統供電，綜合數據采集器2可接收分別來自第一數據采集器11、第二數據采集器12、第三數據采集器13和第四數據采集器14的數據，距離傳感器17可檢測探測孔深度變化情況，孔隙壓力傳感器18可檢測頁巖儲層內孔隙壓力變化情況，應力傳感器19可檢測頁巖儲層的應力變化情況，應變式傳感器20可檢測頁巖儲層的構造應變變化情況，GPS定位器21可對探測孔位置進行定位，單片機控制器1接收綜合數據采集器2的數據，并反饋給數據處理器3，數據處理器3對數據綜合處理分析后，反饋給單片機控制器1，單片機控制器1將反饋數據反饋給應力模擬器10，應力模擬器10可對頁巖儲層的應力情況綜合模擬，直觀反映頁巖儲層的應力情況，整個系統智能化程度較高，對頁巖儲層的數據采集較為全面，保證了預測地應力的準確性，提高了預測效率。
本發明所建立的模型可更真實的反映頁巖氣儲層的復雜礦物以及復雜礦物造成的多孔隙類型，即以一定形態分散存在的干酪根的巖石物理特征，從而在使用該巖石物理模型進行儲層速度預測的過程中提高預測速度的精度；在此基礎上的應力預測使得在沒有測井橫波速度的情況下，根據常規測井曲線可以計算得到與測量值吻合的地下的最大、最小水平主應力，以及破裂壓力，是一種可靠的頁巖氣儲層地應力測井預測方法。
以上所述僅為本發明的較佳實施例而已，并不用以限制本發明，凡在本發明的精神和原則之內所作的任何修改、等同替換和改進等，均應包含在本發明的保護范圍之內。