本公開涉及油氣地球物理勘探領域，更具體地，涉及一種反演裂縫密度的方法及系統(tǒng)。

背景技術：

在油氣地球物理勘探領域，對裂縫性儲層研究已成為日益重要的內容，而裂縫密度是裂縫性儲層定量表征中的重要參數。

在現有的預測裂縫性質的方法中，最常見的是橢圓擬合方法，對于HTI(atransversely isotropic medium with a horizontal axis of symmetry(具有對稱水平軸的橫向各向同性介質))介質，縱波的地震屬性(疊加速度和反射振幅)隨方位變化呈橢圓形特征，長軸方向反映裂縫發(fā)育主要方位，長、短軸之比反映了裂縫發(fā)育程度。但此方法只能估計裂縫的發(fā)育程度，并不能定量的表征裂縫密度。Isabel Varela對含有AVAZ(Amplitude variation with incident angle and azimuth(振幅隨入射角和方位角的變化))特征的數據運用SVD(sigular value decomposition奇異值分解(Singular Value Decomposition))反演裂縫密度。Morten Jakobsen等通過AVOZ分析估計裂縫儲層滲透率。朱培民等提出用2條正交測線上的P波含有AVO特征的觀測數據反演Thomsen參數，進而估計裂縫密度的方法。

發(fā)明人發(fā)現，上述幾種有關的反演裂縫密度的方法均是建立在可以從反射信息中取得準確的反射系數的基礎上，然而實際上地震勘探得到的是振幅信息，提取地震子波得到反射系數的過程不可避免也會對后續(xù)反演的精度產生影響，因此，有必要開發(fā)一種精確反演裂縫密度的方法。

公開于本公開背景技術部分的信息僅僅旨在加深對本公開的一般背景技術的理解，而不應當被視為承認或以任何形式暗示該信息構成已為本領域技術人 員所公知的現有技術。

技術實現要素：

本公開提出了一種反演裂縫密度的方法及系統(tǒng)，該方法可以基于相交測線上的反射振幅差值構造目標函數，基于目標函數進行反演得到與裂縫密度有關的Thomsen參數，以及基于Thomsen參數與裂縫密度的關系獲得裂縫密度，從而實現了對裂縫密度的精確反演。

根據本公開的一方面，提出了一種反演裂縫密度的方法，該方法可以包括以下步驟：基于地震資料獲取相交測線的縱波反射振幅差值，并基于擬合模型獲得縱波反射系數差值；基于縱波反射系數差值和縱波反射振幅差值，構造目標函數fitness；以及基于目標函數fitness進行反演計算以獲得裂縫密度，其中，目標函數fitness可以表示為：

其中，表示方位角；η_j表示相交測線的方位角差值；θ_k表示入射角；表示方位角分別為和入射角為θ_k時的縱波反射振幅差值；表示的最大值；表示反射系數差值或振幅差值最大時對應的入射角；Δγ、Δδ^(V)與Δε^(V)表示Thomsen參數，其中前置符號“Δ”代表上下兩層Thomsen參數的差值；表示擬合模型中兩條相交測線的縱波反射系數差值；表示的最大值；ω^ij表示權系數；I表示地震資料中方位角的測線總數量；J表示地震資料中方位角為時與其相交的測線的總數量；以及K 表示地震資料中入射角θ的總數量。

根據本公開的另一方面，提出了一種反演裂縫密度的系統(tǒng)，該系統(tǒng)可以包括以下單元：用于基于地震資料獲取相交測線的縱波反射振幅差值，并基于擬合模型獲得縱波反射系數差值的單元；用于基于縱波反射系數差值和縱波反射振幅差值，構造目標函數fitness的單元；以及用于基于目標函數fitness進行反演計算以獲得裂縫密度的單元，其中，目標函數fitness可以表示為：

其中，表示方位角；η_j表示相交測線的方位角差值；θ_k表示入射角；表示方位角分別為和入射角為θ_k時的縱波反射振幅差值；表示的最大值；θ_km表示反射系數差值或振幅差值最大時對應的入射角；Δγ、Δδ^(V)與Δε^(V)表示Thomsen參數；表示擬合模型中兩條相交測線的縱波反射系數差值；表示的最大值；ω^ij表示權系數；I表示地震資料中方位角的測線總數量；J表示地震資料中方位角為時與其相交的測線的總數量；以及K表示地震資料中入射角θ的總數量。

本公開的方法和裝置具有其它的特性和優(yōu)點，這些特性和優(yōu)點從并入本文中的附圖和隨后的具體實施方案中將是顯而易見的，或者將在并入本文中的附圖和隨后的具體實施方案中進行詳細陳述，這些附圖和具體實施方案共同用于解釋本公開的特定原理。

附圖說明

通過結合附圖對本公開示例性實施方案進行更詳細的描述，本公開的上述 以及其它目的、特征和優(yōu)勢將變得更加明顯，其中，在本公開示例性實施方案中，相同的參考標號通常代表相同部件。

圖1示出了根據本公開的一個實施方案的反演裂縫密度的方法的步驟的流程圖。

圖2a-2b為示出了根據本公開的一個示例的Thomsen參數與裂縫密度關系的示意圖，其中，圖2a為裂縫含水，圖2b為裂縫含氣。

圖3a-3e為示出了根據本公開的一個示例的Z分量地震記錄的示意圖，其中，圖3a為方位角0°；圖3b為方位角30°；圖3c為方位角45°；圖3d為方位角60°；以及圖3e為方位角90°。

圖4a-4h為示出了根據本公開的一個示例的相交測線上HTI介質反射振幅差值的示意圖，其中，圖4a為方位角0°和90°；圖4b為方位角30°和60°；圖4c為方位角30°和90°；圖4d為方位角0°和30°；圖4e為方位角0°和45°；圖4f為方位角45°和90°；圖4g為方位角45°和60°；以及圖4h為方位角30°和45°。

具體實施方式

下面將參照附圖更詳細地描述本公開的優(yōu)選實施方案。雖然附圖中顯示了本公開的優(yōu)選實施方案，然而應該理解，可以以各種形式實現本公開而不應被這里闡述的實施方案所限制。相反，提供這些實施方案是為了使本公開更加透徹和完整，并且能夠將本公開的范圍完整地傳達給本領域的技術人員。

實施方案1

圖1示出了根據本公開的一個實施方案的反演裂縫密度的方法的流程圖。根據本公開的實施方案的反演裂縫密度的方法可以包括以下步驟：步驟101，基于地震資料獲取相交測線的縱波反射振幅差值和擬合模型中縱波反射系數差值；步驟102，基于縱波反射系數差值和縱波反射振幅差值，構造目標函數fitness； 以及步驟103，基于目標函數fitness進行反演計算以獲得裂縫密度，其中，目標函數fitness可以表示為：

其中，表示方位角；η_j表示相交測線的方位角差值；θ_k表示入射角；表示方位角分別為和入射角為θ_k時的縱波反射振幅差值；表示的最大值；θ_km表示反射系數差值或振幅差值最大時對應的入射角；Δγ、Δδ^(V)與Δε^(V)表示Thomsen參數，其中前置符號“Δ”代表上下兩層Thomsen參數的差值；表示擬合模型中兩條相交測線的縱波反射系數差值；表示的最大值；ω^ij表示權系數；I表示地震資料中方位角的測線總數量；J表示地震資料中方位角為時與其相交的測線的總數量；以及K表示地震資料中入射角θ的總數量。

本實施方案可以通過利用縱波反射系數差值和縱波反射振幅差值來構造目標函數，并基于目標函數進行計算以獲得裂縫密度，實現了對裂縫密度的精確反演。

獲取縱波反射系數差值和縱波反射振幅差值

在一個示例中，可以基于擬合模型獲得縱波反射系數差值。其中，可以使用Rüger反射系數近似公式來計算縱波反射系數差值。

根據Rüger反射系數近似公式，當方位角為和時，縱波反射系數差值可以用如下公式(1)表示：

其中，R可以表示縱波反射系數，可以表示方位角為時的測線的縱波反射系數，可以表示方位角為時的測線的縱波反射系數，ΔR可以表示方位角為和時的相交測線的縱波反射系數差值，α可以表示HTI介質各向同性面的縱波速度，β可以表示HTI介質各向同性面的SH波速度，可以表示HTI介質各向同性面的縱波速度的平均值，可以表示HTI介質各向同性面的SH波速度的平均值，θ可以表示入射角。其中，γ、δ^(V)與ε^(V)可以表示Thomsen參數，γ可以表示在水平方向傳播的SH波與SV波速度的差別，，δ^(V)可以表示P波相速度在垂直入射時的二階導數，ε^(V)可以表示縱波速度在垂直和水平方向的差別，Thomsen參數上角標V代表區(qū)分針對VTI介質的Thomsen參數，前置符號“Δ”代表上下兩層Thomsen參數的差值。

本領域技術人員應該理解，獲得縱波反射系數差值的方法并不局限于此，而是可以采用本領域技術人員所公知的各種反射系數近似方式來獲得縱波反射系數的差值。

從公式(1)可以獲得縱波反射系數差值，并可以使用表示兩條相交測線的縱波反射系數差值，其中(i＝1,2,...,I；j＝1,2,...,J；k＝1,2,...,K)，I可以表示地震資料中方位角的測線總數量，J可以表示地震資料中方位角為時與其相交的測線的總數量，K可以表示地震資料中入射角θ的總數量。

在一個示例中，可以基于地震資料獲取相交測線的縱波反射振幅差值。其中，可以使用觀測地震數據中目的層位上的相交測線反射振幅來計算縱波反射 振幅差值，其可以使用如下公式(2)表示：

其中，A可以表示從反射信息中提取到的縱波反射振幅，可以表示方位角分別為和入射角為θ_k時的縱波反射振幅差值。本領域技術人員應當理解，本公開并不限于此，而可以使用本領域公知的任何地震數據處理與解釋手段來計算縱波反射振幅差值。

構造目標函數

在一個示例中，可以基于縱波反射系數差值和縱波反射振幅差值，構造目標函數。所構建的目標函數可以用如下公式(3)表示：

其中，ω^ij可以表示權系數，可以表示在K個入射角中反射系數差值或振幅差值最大時對應的入射角，可以表示方位角分別為和入射角為θ_k時的縱波反射振幅差值的最大值，可以表示兩條相交測線的縱波反射系數差值的最大值。

可以使用如下公式(4)表示：

可以使用如下公式(5)表示：

其中，從擬合模型中獲得的縱波反射系數差值和從地震資料中獲取的縱波反射振幅差值兩者均可以是無量綱的物理量。由于反射系數與反射振幅呈線性 關系，因此可以利用反射振幅差值中的最大值點為參考，對應擬合模型中反射系數差值的最大值，可以計算出反射振幅差值與反射系數差值的比值關系，再根據此比值關系將其它反射系數差值轉換為擬合模型的反射振幅差值。因而公式(3)中的后一項可以為擬合模型的反射振幅差值，與前一項(實際的縱波反射振幅差值)可以是同一個物理量。因此，可以使用與的差值最小的情況(也即目標函數fitness取最小值)下的Thomsen參數計算裂縫密度。

在一個示例中，權系數ω^ij可以用于控制一組不同方位角的縱波反射振幅差值的可信程度。其中，權系數ω^ij可以由用戶根據縱波反射振幅差值的可信程度自行設定，例如ω^ij可以設定為0-1之間的值。某一入射角或方位角的觀測數據的質量越好，也即可信程度越高，則權系數可以越大。具體地，可以通過地震數據的信號與噪音情況、某一地震面元上各個方位上的地震數據是否足夠多來判斷數據是否可信，權系數可以是主觀參數。在實際情況下，如果該地區(qū)數據采集是規(guī)律的，即每個面元的炮點和檢波點分布一致，則某一個i或者某一個j對應的權系數是可以一樣的。

獲得裂縫密度

在一個示例中，可以基于目標函數fitness進行反演計算以獲得裂縫密度。通過非線性反演算法求得公式(3)所示的目標函數fitness的最小值，基于目標函數fitness的最小值可以獲得與該最小值對應的一組Thomsen參數Δγ，Δδ^(V)和Δε^(V)的值。本領域技術人員應當理解，求得目標函數fitness的最小值的方法并不限于此，而可以使用本領域公知的任何方法求得目標函數fitness的最小值。

如果上層介質為均勻各向同性介質，下層為HTI介質，則Δγ，Δδ^(V)和Δε^(V) 均等于下層介質的Thomsen參數γ，δ^(V)和ε^(V)。因此，根據所獲得的Thomsen參數Δγ，Δδ^(V)和Δε^(V)的值，可以基于如下公式(6)和(7)中的HTI介質Thomsen參數與裂縫密度e的關系，獲得裂縫密度e：

$e=\frac{3(3-2g)\mathrm{\γ}}{8(1+2\mathrm{\γ})}---\left(6\right)$

$g=\frac{\mathrm{\μ}}{\mathrm{\λ}+2\mathrm{\μ}}=\frac{V_S^2}{V_P^2}---\left(7\right)$

其中，V_P和V_S可以為背景介質的縱波和橫波速度，e可以表示裂縫密度，λ和μ是拉梅常數，g可以表示橫波與縱波速度比值的平方。

應用示例

為便于理解本公開實施方案的方案及其效果，以下給出一個具體應用示例。本領域技術人員應理解，該示例僅為了便于理解本公開，其任何具體細節(jié)并非意在以任何方式限制本公開。

圖2a和圖2b為示出了根據本公開的一個示例的Thomsen參數與裂縫密度關系的示意圖。其中，圖2a為裂縫含水，圖2b為裂縫含氣。在圖2a和圖2b中，縱波與橫波背景速度分別為3800m/s和2000m/s，密度為2.4g/cm3，裂縫縱橫比為10-3，縱軸可以表示Thomsen參數，圖2a的橫軸可以表示含水裂縫密度，圖2b的橫軸可以表示含氣裂縫密度。

圖3a-3e為示出了根據本公開的一個示例的Z分量地震記錄的示意圖，其中，圖3a為方位角0°；圖3b為方位角30°；圖3c為方位角45°；圖3d為方位角60°；以及圖3e為方位角90°。其中，模型上層為均勻各向同性介質，下層為HTI介質。

圖4a-4h為示出了根據本公開的一個示例的相交測線上HTI介質反射振幅差值的示意圖，其中，圖4a為方位角0°和90°；圖4b為方位角30°和60°；圖4c為方位角30°和90°；圖4d為方位角0°和30°；圖4e為方位角0°和 45°；圖4f為方位角45°和90°；圖4g為方位角45°和60°；以及圖4h為方位角30°和45°。

如圖4a-4h所示，反演得到的Thomsen參數結果：γ相對誤差為1.3％；δ^(V)相對誤差為28％；ε^(V)相對誤差為26.6％。其中相對誤差公式為：|(真實值-反演值)/真實值|×100％。顯然，三個Thomsen參數反演值中，γ反演結果最接近真實值。根據公式(4)中裂縫密度e與γ的關系，可估計裂縫密度為相對誤差為9％。數值反演結果表明：相交測線上反射振幅差值正演值與反演值擬合情況較好，基于縱波反射系數差值和縱波反射振幅差值來反演裂縫密度的方法是可行且穩(wěn)定的。

本領域技術人員應理解，上面對本公開的實施方案的描述的目的僅為了示例性地說明本公開的實施方案的有益效果，并不意在將本公開的實施方案限制于所給出的任何示例。

實施方案2

在該實施方案中，提供了一種反演裂縫密度的系統(tǒng)，該系統(tǒng)可以包括以下單元：用于基于地震資料獲取相交測線的縱波反射振幅差值，并基于擬合模型獲得縱波反射系數差值的單元；用于基于縱波反射系數差值和縱波反射振幅差值，構造目標函數fitness的單元；以及用于基于目標函數fitness進行反演計算以獲得裂縫密度的單元。其中，目標函數fitness可以表示為：

其中，表示方位角；η_j表示相交測線的方位角差值；θ_k表示入射角；表示方位角分別為和入射角為θ_k時的縱波反射振幅差值；表示的最大值；表示反射系數差值或振幅差值最大時對 應的入射角；Δγ、Δδ^(V)與Δε^(V)表示Thomsen參數；表示擬合模型中兩條相交測線的縱波反射系數差值；表示的最大值；ω^ij表示權系數；I表示地震資料中方位角的測線總數量；J表示地震資料中方位角為時與其相交的測線的總數量；以及K表示地震資料中入射角θ的總數量。

本實施方案可以通過利用縱波反射系數差值和縱波反射振幅差值來構造目標函數，并基于目標函數進行計算以獲得裂縫密度，實現了對裂縫密度的精確反演。

在一個示例中，可以使用Rüger反射系數近似公式來計算縱波反射系數差值。

在一個示例中，可以使用觀測地震數據中目的層位上的相交測線的反射振幅來計算縱波反射振幅差值。

在一個示例中，目標函數fitness中的權系數可以用于控制一組不同方位角的縱波反射振幅差值的可信程度。

在一個示例中，基于所述目標函數fitness進行反演計算以獲得裂縫密度可以包括：通過非線性反演算法獲得目標函數fitness的最小值；基于目標函數fitness的最小值獲得Thomsen參數；以及基于Thomsen參數進行計算以獲得裂縫密度。

本領域技術人員應理解，上面對本公開的實施方案的描述的目的僅為了示例性地說明本公開的實施方案的有益效果，并不意在將本公開的實施方案限制于所給出的任何示例。

本公開可以是系統(tǒng)、方法和/或計算機程序產品。計算機程序產品可以包括計算機可讀存儲介質，其上載有用于使處理器實現本公開的各個方面的計算機可讀程序指令。

計算機可讀存儲介質可以是可以保持和存儲由指令執(zhí)行設備使用的指令 的有形設備。計算機可讀存儲介質例如可以是――但不限于――電存儲設備、磁存儲設備、光存儲設備、電磁存儲設備、半導體存儲設備或者上述的任意合適的組合。計算機可讀存儲介質的更具體的例子(非窮舉的列表)包括：便攜式計算機盤、硬盤、隨機存取存儲器(RAM)、只讀存儲器(ROM)、可擦式可編程只讀存儲器(EPROM或閃存)、靜態(tài)隨機存取存儲器(SRAM)、便攜式壓縮盤只讀存儲器(CD-ROM)、數字多功能盤(DVD)、記憶棒、軟盤、機械編碼設備、例如其上存儲有指令的打孔卡或凹槽內凸起結構、以及上述的任意合適的組合。這里所使用的計算機可讀存儲介質不被解釋為瞬時信號本身，諸如無線電波或者其他自由傳播的電磁波、通過波導或其他傳輸媒介傳播的電磁波(例如，通過光纖電纜的光脈沖)、或者通過電線傳輸的電信號。

這里所描述的計算機可讀程序指令可以從計算機可讀存儲介質下載到各個計算/處理設備，或者通過網絡、例如因特網、局域網、廣域網和/或無線網下載到外部計算機或外部存儲設備。網絡可以包括銅傳輸電纜、光纖傳輸、無線傳輸、路由器、防火墻、交換機、網關計算機和/或邊緣服務器。每個計算/處理設備中的網絡適配卡或者網絡接口從網絡接收計算機可讀程序指令，并轉發(fā)該計算機可讀程序指令，以供存儲在各個計算/處理設備中的計算機可讀存儲介質中。

用于執(zhí)行本公開操作的計算機程序指令可以是匯編指令、指令集架構(ISA)指令、機器指令、機器相關指令、微代碼、固件指令、狀態(tài)設置數據、或者以一種或多種編程語言的任意組合編寫的源代碼或目標代碼，所述編程語言包括面向對象的編程語言—諸如Smalltalk、C++等，以及常規(guī)的過程式編程語言—諸如“C”語言或類似的編程語言。計算機可讀程序指令可以完全地在用戶計算機上執(zhí)行、部分地在用戶計算機上執(zhí)行、作為一個獨立的軟件包執(zhí)行、部分在用戶計算機上部分在遠程計算機上執(zhí)行、或者完全在遠程計算機或服務器上執(zhí)行。在涉及遠程計算機的情形中，遠程計算機可以通過任意種類的網絡—包括局域網(LAN)或廣域網(WAN)—連接到用戶計算機，或者，可以連接到外部計算機(例如利用因特網服務提供商來通過因特網連接)。在一些實施方案中，通 過利用計算機可讀程序指令的狀態(tài)信息來個性化定制電子電路，例如可編程邏輯電路、現場可編程門陣列(FPGA)或可編程邏輯陣列(PLA)，該電子電路可以執(zhí)行計算機可讀程序指令，從而實現本公開的各個方面。

這里參照根據本公開實施方案的方法、裝置(系統(tǒng))和計算機程序產品的流程圖和/或框圖描述了本公開的各個方面。應當理解，流程圖和/或框圖的每個方框以及流程圖和/或框圖中各方框的組合，都可以由計算機可讀程序指令實現。

這些計算機可讀程序指令可以提供給通用計算機、專用計算機或其它可編程數據處理裝置的處理器，從而生產出一種機器，使得這些指令在通過計算機或其它可編程數據處理裝置的處理器執(zhí)行時，產生了實現流程圖和/或框圖中的一個或多個方框中規(guī)定的功能/動作的裝置。也可以把這些計算機可讀程序指令存儲在計算機可讀存儲介質中，這些指令使得計算機、可編程數據處理裝置和/或其他設備以特定方式工作，從而，存儲有指令的計算機可讀介質則包括一個制造品，其包括實現流程圖和/或框圖中的一個或多個方框中規(guī)定的功能/動作的各個方面的指令。

也可以把計算機可讀程序指令加載到計算機、其它可編程數據處理裝置、或其它設備上，使得在計算機、其它可編程數據處理裝置或其它設備上執(zhí)行一系列操作步驟，以產生計算機實現的過程，從而使得在計算機、其它可編程數據處理裝置、或其它設備上執(zhí)行的指令實現流程圖和/或框圖中的一個或多個方框中規(guī)定的功能/動作。

附圖中的流程圖和框圖顯示了根據本公開的多個實施方案的系統(tǒng)、方法和計算機程序產品的可能實現的體系架構、功能和操作。在這點上，流程圖或框圖中的每個方框可以代表一個模塊、程序段或指令的一部分，所述模塊、程序段或指令的一部分包含一個或多個用于實現規(guī)定的邏輯功能的可執(zhí)行指令。在有些作為替換的實現中，方框中所標注的功能也可以以不同于附圖中所標注的順序發(fā)生。例如，兩個連續(xù)的方框實際上可以基本并行地執(zhí)行，它們有時也可以按相反的順序執(zhí)行，這依所涉及的功能而定。也要注意的是，框圖和/或流程 圖中的每個方框、以及框圖和/或流程圖中的方框的組合，可以用執(zhí)行規(guī)定的功能或動作的專用的基于硬件的系統(tǒng)來實現，或者可以用專用硬件與計算機指令的組合來實現。

以上已經描述了本公開的各實施方案，上述說明是示例性的，并非窮盡性的，并且也不限于所披露的各實施方案。在不偏離所說明的各實施方案的范圍和精神的情況下，對于本技術領域的普通技術人員來說許多修改和變更都是顯而易見的。本文中所用術語的選擇，旨在最好地解釋各實施方案的原理、實際應用或對市場中的技術的技術改進，或者使本技術領域的其它普通技術人員能理解本文披露的各實施方案。