本發(fā)明涉及一種基于壓汞法毛管壓力曲線的復(fù)雜巖石分形維數(shù)分布譜構(gòu)建方法，屬于石油工程及巖石物理領(lǐng)域。

背景技術(shù)：

分形理論自發(fā)展以來(lái)，已廣泛應(yīng)用于石油工程和巖石物理的各個(gè)領(lǐng)域，成為表征巖石孔隙結(jié)構(gòu)的重要手段之一，取得了較好的應(yīng)用效果。

毛管壓力曲線是分形理論用于表征巖石孔隙結(jié)構(gòu)的重要基礎(chǔ)資料，以往通常根據(jù)飽和度和壓力的冪指數(shù)關(guān)系構(gòu)建簡(jiǎn)單的單重分形維數(shù)模型并用最小二乘法求解。隨著研究的進(jìn)一步深入，所面臨的對(duì)象也日益復(fù)雜，從簡(jiǎn)單的高孔隙度、高滲透率砂巖向碳酸鹽巖、火山巖和致密砂巖等復(fù)雜巖石發(fā)展，分形理論在應(yīng)用中遇到一些瓶頸。傳統(tǒng)的分形理論和分形維數(shù)計(jì)算方法很難將復(fù)雜巖石的毛管壓力曲線進(jìn)行很好的擬合，特別是對(duì)于孔隙結(jié)構(gòu)極其復(fù)雜、非均質(zhì)性極強(qiáng)的巖石而言，它們的毛管壓力曲線呈明顯的分段特征。為克服傳統(tǒng)分形理論和模型在描述毛管壓力曲線中的缺陷，非常有必要對(duì)分形理論進(jìn)行深入分析，發(fā)展新的毛管壓力曲線表征方法。

技術(shù)實(shí)現(xiàn)要素：

基于上述技術(shù)問(wèn)題，本發(fā)明提供一種基于壓汞法毛管壓力曲線的復(fù)雜巖石分形維數(shù)分布譜構(gòu)建方法，為復(fù)雜巖石的孔隙結(jié)構(gòu)評(píng)價(jià)和分形特征研究提供參考。

本發(fā)明所采用的技術(shù)解決方案是：

一種基于壓汞法毛管壓力曲線的復(fù)雜巖石分形維數(shù)分布譜構(gòu)建方法，包括以下步驟：

a巖石預(yù)處理、基礎(chǔ)數(shù)據(jù)采集和壓汞實(shí)驗(yàn)

a1對(duì)巖石進(jìn)行鉆取、切割、打磨，得到呈柱塞樣的樣品，將樣品放入烘箱，烘干后取出，用游標(biāo)卡尺測(cè)量樣品的長(zhǎng)度、直徑，計(jì)算樣品的總體積；用電子天平測(cè)量樣品的質(zhì)量；

a2將樣品放入孔隙度儀，采用氦氣或氮?dú)庾鳛榻橘|(zhì)，應(yīng)用波義耳定律測(cè)量得到樣品的孔隙體積，計(jì)算樣品孔隙度；

a3將樣品放入壓汞儀，采用恒速或者恒壓法將汞注入樣品，記錄注入過(guò)程中進(jìn)汞壓力和汞體積的變化，進(jìn)汞壓力和汞體積采集的個(gè)數(shù)為m；

b將汞體積除以樣品的孔隙體積，得到進(jìn)汞飽和度；根據(jù)進(jìn)汞飽和度和進(jìn)汞壓力的冪指數(shù)關(guān)系，構(gòu)建多個(gè)分形維數(shù)的累加公式：

$S=\underset{j=1}{\overset{n}{\Σ}}{a}_j{P}_c^{-(3-D_j)}---\left(1\right)$

式(1)中：D_j為預(yù)設(shè)的第j個(gè)分形維數(shù)序列；a_j為D_j所對(duì)應(yīng)的比例；n是預(yù)設(shè)的分形維數(shù)的個(gè)數(shù)；P_c為進(jìn)汞壓力；S為進(jìn)汞壓力為P_c時(shí)的進(jìn)汞飽和度；

c將式(1)進(jìn)行變形，得到：

$S=\underset{j=1}{\overset{n}{\Σ}}{a}_j{e}^{-(3-D_j)ln\left(P_c\right)}---\left(2\right)$

d將式(2)表示為最優(yōu)化求解的矩陣形式，得到：

min||S-AX||² (3)

式(3)中：S是由進(jìn)汞飽和度形成的大小為m×1矩陣，已知量；A是由構(gòu)成的大小為m×n矩陣，已知量；X是分形維數(shù)的比例構(gòu)成的大小為n×1矩陣；m表示壓汞法毛管壓力曲線測(cè)量過(guò)程中進(jìn)汞壓力記錄的點(diǎn)數(shù)；

e對(duì)式(3)中的A進(jìn)行奇異值分解，得到：

A＝UZV^T (4)

式(4)中：U是大小為m×m的正交矩陣，V是大小為n×n的正交矩陣；V^T是V的轉(zhuǎn)置；Z是大小為m×n的對(duì)角矩陣；

f設(shè)置閾值q將對(duì)角矩陣進(jìn)行截?cái)?，令?duì)角矩陣Z中所有小于q的元素為0，并將Z改寫(xiě)為Z_q，代入式(4)，可得：

A_q＝UZ_qV^T (5)

g采用矩陣變換法，將式(5)代入式(3)，得到X的初值：

$X_1={\mathrm{VZ}}_q^{-1}US---\left(6\right)$

式(6)中：是Z_q的逆矩陣；

h應(yīng)用迭代算法進(jìn)行求解，得到最優(yōu)的X

h1根據(jù)式(6)得到初值X₁，設(shè)置最小平均相對(duì)誤差為T(mén)和最大迭代次數(shù)為I，令i＝1；

h2計(jì)算誤差矩陣ΔS_i＝S-UZ_qV^TX_i和平均相對(duì)誤差ΔT_i＝||ΔS_i||/||S||；

h3如果ΔT_i＜T或者i＝I，停止迭代，輸出X_i；否則轉(zhuǎn)入步驟h4；

h4計(jì)算和X_i＝X_i+ΔX_i，令i＝i+1；并轉(zhuǎn)入步驟h2；

k完成優(yōu)化迭代步驟h后輸出X_i，繪制分形維數(shù)分布譜。

優(yōu)選的，預(yù)設(shè)的分形維數(shù)的個(gè)數(shù)n小于進(jìn)汞壓力記錄的點(diǎn)數(shù)m。

優(yōu)選的，預(yù)設(shè)的分形維數(shù)D_j最大值和最小值分別為3和2。

本發(fā)明的有益技術(shù)效果是：

本發(fā)明可實(shí)現(xiàn)對(duì)具有分段分形特征的毛管壓力曲線的細(xì)致刻畫(huà)，傳統(tǒng)的單重分形維數(shù)方法針對(duì)一塊樣品的毛管壓力曲線只能得到一個(gè)分形維數(shù)；而通過(guò)本發(fā)明，針對(duì)一塊樣品的毛管壓力曲線，可得到多個(gè)分形維數(shù)以及它們的比例，從而實(shí)現(xiàn)了分形維數(shù)分布譜的構(gòu)建，為復(fù)雜巖石的孔隙結(jié)構(gòu)分析和分形表征帶來(lái)新的思路。

附圖說(shuō)明

下面結(jié)合附圖與具體實(shí)施方式對(duì)本發(fā)明作進(jìn)一步說(shuō)明：

圖1為本發(fā)明提供的一種基于壓汞法毛管壓力曲線的復(fù)雜巖石分形維數(shù)分布譜構(gòu)建方法的技術(shù)流程圖；

圖2是某復(fù)雜碳酸鹽巖樣品的壓汞法毛管壓力曲線圖；

圖3是應(yīng)用本發(fā)明提供的方法對(duì)圖2所述樣品進(jìn)行處理得到的分形維數(shù)分布譜。

具體實(shí)施方式

結(jié)合附圖，一種基于壓汞法毛管壓力曲線的復(fù)雜巖石分形維數(shù)分布譜構(gòu)建方法，包括以下步驟：

a巖石預(yù)處理、基礎(chǔ)數(shù)據(jù)采集和壓汞實(shí)驗(yàn)

a1對(duì)巖石進(jìn)行鉆取、切割、打磨成直徑為2.54cm，長(zhǎng)度為約3-5cm的柱塞樣，將樣品放入烘箱，選擇溫度為95℃并恒溫24小時(shí)后，取出，用游標(biāo)卡尺測(cè)量樣品的長(zhǎng)度、直徑，計(jì)算樣品的總體積；用電子天平測(cè)量樣品的質(zhì)量。

a2將樣品放入氣測(cè)法孔隙度儀，采用氦氣或氮?dú)庾鳛榻橘|(zhì)，應(yīng)用波義耳定律測(cè)量得到樣品的孔隙體積，計(jì)算樣品孔隙度。

a3將樣品放入壓汞儀，采用恒速或者恒壓法將汞注入樣品，記錄注入過(guò)程中進(jìn)汞壓力和汞體積的變化，進(jìn)汞壓力和汞體積采集的個(gè)數(shù)為m。

b將汞體積除以樣品的孔隙體積，得到進(jìn)汞飽和度。根據(jù)進(jìn)汞飽和度和進(jìn)汞壓力的冪指數(shù)關(guān)系，構(gòu)建多個(gè)分形維數(shù)的累加公式：

$S=\underset{j=1}{\overset{n}{\Σ}}{a}_j{P}_c^{-(3-D_j)}---\left(1\right)$

式(1)中：D_j為預(yù)設(shè)的第j個(gè)分形維數(shù)序列；a_j為D_j所對(duì)應(yīng)的比例；n是預(yù)設(shè)的分形維數(shù)的個(gè)數(shù)；P_c為進(jìn)汞壓力；S為進(jìn)汞壓力為P_c時(shí)的進(jìn)汞飽和度。

c將式(1)進(jìn)行變形，得到：

$S=\underset{j=1}{\overset{n}{\Σ}}{a}_j{e}^{-(3-D_j)ln\left(P_c\right)}---\left(2\right)$

d將式(2)表示為最優(yōu)化求解的矩陣形式，得到：

min||S-AX||² (3)

式(3)中：S是由進(jìn)汞飽和度形成的大小為m×1矩陣，已知量；A是由構(gòu)成的大小為m×n矩陣，已知量；X是分形維數(shù)的比例構(gòu)成的大小為n×1矩陣；m表示壓汞法毛管壓力曲線測(cè)量過(guò)程中進(jìn)汞壓力記錄的點(diǎn)數(shù)。

e對(duì)式(3)中的A進(jìn)行奇異值分解，得到：

A＝UZV^T (4)

式中：U是大小為m×m的正交矩陣，V是大小為n×n的正交矩陣；V^T是V的轉(zhuǎn)置；Z是大小為m×n的對(duì)角矩陣。

f設(shè)置閾值q將對(duì)角矩陣進(jìn)行截?cái)?，令?duì)角矩陣Z中所有小于q的元素為0，并將Z改寫(xiě)為Z_q，代入式(4)，可得：

A_q＝UZ_qV^T (5)

g采用矩陣變換法，將式(5)代入式(3)，得到X的初值為：

式中：是Z_q的逆矩陣。

h應(yīng)用迭代算法進(jìn)行求解，得到最優(yōu)的X

h1根據(jù)式(6)得到初值X₁，設(shè)置最小平均相對(duì)誤差為T(mén)和最大迭代次數(shù)為I，令i＝1。

h2計(jì)算誤差矩陣ΔS_i＝S-UZ_qV^TX_i和平均相對(duì)誤差ΔT_i＝||ΔS_i||/||S||。

h3如果ΔT_i＜T或者i＝I，停止迭代，輸出X_i；否則轉(zhuǎn)入步驟h4。

h4計(jì)算和X_i＝X_i+ΔX_i，令i＝i+1；并轉(zhuǎn)入步驟h2。

k完成優(yōu)化迭代步驟h后輸出X_i，繪制分形維數(shù)分布譜。

上述步驟中：預(yù)設(shè)的分形維數(shù)的個(gè)數(shù)n小于進(jìn)汞壓力記錄的點(diǎn)數(shù)m。

上述步驟中：預(yù)設(shè)的分形維數(shù)D_j最大值和最小值分別為3和2。

上述方法不僅僅限于壓汞法毛管壓力曲線的處理，還適用于半滲透率隔板法毛管壓力曲線和離心法毛管壓力曲線的處理。

圖1是一種基于壓汞法毛管壓力曲線的巖石分形維數(shù)分布譜構(gòu)建方法的技術(shù)流程圖，主要包括巖石鉆取、切割、打磨、烘干等預(yù)處理，巖石直徑、長(zhǎng)度、質(zhì)量、孔隙度等基本參數(shù)測(cè)量，巖石壓汞法毛管壓力測(cè)量與數(shù)據(jù)整理，巖石分形維數(shù)的分布譜求解等四個(gè)部分，這四個(gè)部分缺一不可，且順序不可顛倒。

圖2是某碳酸鹽巖樣品實(shí)測(cè)的壓汞法毛管壓力曲線，壓汞實(shí)驗(yàn)采用恒速法。需要說(shuō)明的是，采用恒壓法得到的毛管壓力曲線，也可用本發(fā)明提供的方法進(jìn)行處理。

圖3是應(yīng)用本發(fā)明所提供的方法對(duì)圖2所述的某碳酸鹽巖樣品進(jìn)行處理得到的分形維數(shù)的分布譜，從圖中可知，該樣品的分形維數(shù)主要在2到2.45之間，呈單峰分布。

上述方式中未述及的部分采取或借鑒已有技術(shù)即可實(shí)現(xiàn)。

需要說(shuō)明的是，在本說(shuō)明書(shū)的教導(dǎo)下，本領(lǐng)域技術(shù)人員所作出的任何等同替代方式，或明顯變型方式，均應(yīng)在本發(fā)明的保護(hù)范圍之內(nèi)。