本公開涉及石油地震勘探速度建模技術(shù)，尤其涉及一種高精度的基于地面地震構(gòu)造約束的vsp(verticalseismicprofilling，垂直地震剖面)數(shù)據(jù)全波形反演建模方法。

背景技術(shù)：

隨著油氣勘探的日益精細(xì)化、復(fù)雜化，高精度的速度建模在地球物理處理中扮演著越來越重要的角色。由于地震波速度不僅決定偏移成像的質(zhì)量，還與巖石性質(zhì)關(guān)系密切，能夠反映巖石類別和富含流體情況，最終還可以影響地震解釋結(jié)果的可信性，因而地震波速度是一個非常重要的地層物性參數(shù)。

近年來基于地震數(shù)據(jù)振幅、相位、波形匹配的全波形反演方法由于是基于有限頻近似，未采用高頻近似假設(shè)，因而理論上具有高精度刻畫地下構(gòu)造的能力，被業(yè)界認(rèn)為是速度建模精度最高的一種方法。目前該方法在常規(guī)地面地震中已經(jīng)得到了較好的應(yīng)用，特別是海上地震資料的全波形反演已有大量的成功應(yīng)用案例。但是采用vsp數(shù)據(jù)進(jìn)行全波形反演以得到高精度速度模型目前還鮮有研究案例。

目前業(yè)內(nèi)常規(guī)利用vsp數(shù)據(jù)對地下介質(zhì)構(gòu)造進(jìn)行建模的方式有很多。在地震數(shù)據(jù)解釋中，業(yè)內(nèi)專家學(xué)者陳信平(1992)、周熙襄等(1992)采用零偏vsp的初至旅行時來計算不同地層間的層速度。當(dāng)?shù)叵聵?gòu)造復(fù)雜時，采用這種方法求解的層速度誤差較大。另外一種利用vsp數(shù)據(jù)進(jìn)行速度建模的方法是旅行時層析。國外學(xué)者stewart(1984)、salo等(1989)、lee(1990)詳實(shí)系統(tǒng)地研究了如何應(yīng)用vsp數(shù)據(jù)進(jìn)行旅行時層析以獲取地層層速度，此后國內(nèi)學(xué)者劉清林(1996)采用多偏移距vsp數(shù)據(jù)完成了井間層析成像，朱健等(1999)利用vsp數(shù)據(jù)上行波的旅行時同時對層速度和界面形態(tài)參數(shù)進(jìn)行了反演。黃光南等(2012)發(fā)展了變阻尼約束層析成像方法，并把該方法應(yīng)用在了vsp數(shù)據(jù)速度建模當(dāng)中，在一定程度上改善了射線不均勻覆蓋帶來的影響。由于vsp數(shù)據(jù)的角度覆蓋有限，且射線層析類速度建模方案往往基于高頻近似假設(shè)，因而采用vsp地震數(shù)據(jù)層析建模往往只能獲取有限角度范圍內(nèi)的低波數(shù)背景速度，而無法精細(xì)刻畫異常速度體的準(zhǔn)確位置和構(gòu)造形態(tài)。雖然全波形反演可以有效刻畫模型的構(gòu)造細(xì)節(jié)，但是當(dāng)初始模型不準(zhǔn)或數(shù)據(jù)不完備時，梯度會產(chǎn)生噪音，造成反演結(jié)果陷入局部極小值問題。

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

本公開的目的是通過在vsp數(shù)據(jù)全波形反演的梯度中引入地面地震構(gòu)造約束，從而有效地解決常規(guī)地面地震數(shù)據(jù)全波形反演面臨的局部極小值問題，并在一定程度上突破了vsp數(shù)據(jù)角度覆蓋局限問題。

本公開提供一種基于地面地震構(gòu)造約束的vsp數(shù)據(jù)全波形反演建模方法，包括以下步驟：

步驟s1：基于地面地震構(gòu)造約束進(jìn)行vsp數(shù)據(jù)全波形反演，獲得中間速度模型m1；

其中所述步驟s1包括以下子步驟s11至s16：

子步驟s11：采用初始速度模型m0以及地面地震數(shù)據(jù)進(jìn)行疊前深度偏移，獲得地下地質(zhì)構(gòu)造模式f(x)，其中x表示空間向量；

子步驟s12：構(gòu)建所述地下地質(zhì)構(gòu)造模式f(x)的結(jié)構(gòu)特征張量算子；

子步驟s13：基于時間域數(shù)據(jù)匹配的全波形反演構(gòu)建vsp數(shù)據(jù)殘差，基于所述vsp數(shù)據(jù)殘差計算所述速度模型m0的梯度gv(x)，其中所述vsp數(shù)據(jù)殘差是模擬的vsp數(shù)據(jù)與觀測的vsp數(shù)據(jù)之差；

子步驟s14：基于各向異性擴(kuò)散方程構(gòu)建沿所述地下地質(zhì)構(gòu)造模式f(x)的結(jié)構(gòu)特征張量算子平滑的地面地震構(gòu)造約束濾波算子，對所述梯度gv(x)沿所述地面地震構(gòu)造約束方向進(jìn)行平滑濾波，得到基于地面地震構(gòu)造約束的速度模型梯度gv(x)；

子步驟s15：通過p-l-bfgs算法更新所述速度模型m0，獲得更新的速度模型m0’；

子步驟s16：判斷是否完成給定次數(shù)的迭代，若是則將更新的速度模型作為中間速度模型m1，并繼續(xù)到步驟s2，否則返回子步驟s11，以更新的速度模型m0’作為子步驟s11中的初始速度模型；

步驟s2：針對所述中間速度模型m1，進(jìn)行沒有地面地震構(gòu)造約束的vsp數(shù)據(jù)全波形反演，得到最終速度模型m。

優(yōu)選地，采用克?；舴蚱扑阕?、單程波偏移算子或者逆時偏移算子進(jìn)行疊前深度偏移。

優(yōu)選地，所述結(jié)構(gòu)特征張量算子包括所述地下地質(zhì)構(gòu)造模式f(x)在每個點(diǎn)的一階對稱半正定結(jié)構(gòu)張量矩陣的特征向量和特征值。

優(yōu)選地，采用聲波正演模擬方程模擬vsp數(shù)據(jù)和正向傳播波場，其中所述聲波正演模擬方程表示為公式(3)：

其中，

p(x,t)＝[vx(x,t),vy(x,t),vz(x,t),p(x,t)]^t代表正向傳播波場，

s(x,t)＝[0,0,0,s(xs,t)]^t代表震源矢量，xs為震源位置坐標(biāo)，

b代表地下介質(zhì)密度的倒數(shù)，

κ＝v²/b代表體積模量，

v代表地下介質(zhì)的速度參數(shù)。

優(yōu)選地，根據(jù)所述vsp數(shù)據(jù)殘差求解逆時傳播的伴隨波場，以及基于所述伴隨波場根據(jù)公式(4)求解所述速度模型m0的梯度gv(x)：

其中，ns表示震源個數(shù)，t代表波場傳播的終止時刻，q(x)是所述逆時傳播的伴隨波場的分量。

優(yōu)選地，所述地面地震構(gòu)造約束濾波算子表述為公式(5)：

其中，d(x)表示擴(kuò)散張量，且擴(kuò)散張量d(x)與所述結(jié)構(gòu)特征張量算子的特征值一致，α代表光滑因子。

優(yōu)選地，采用非精確線性搜索法計算對所述速度模型m0進(jìn)行更新的步長。

優(yōu)選地，采用p-l-bfgs算法求解海塞矩陣的逆矩陣，并基于所述逆矩陣更新所述速度模型m0。

優(yōu)選地，所述步驟2包括以下子步驟s21至s23：

子步驟s21：基于所述時間域數(shù)據(jù)匹配的全波形反演算子構(gòu)建vsp數(shù)據(jù)殘差，基于所述vsp數(shù)據(jù)殘差計算所述中間速度模型m1的梯度；

子步驟s22：通過p-l-bfgs算法更新所述速度模型m1，獲得更新的速度模型m1’；

子步驟s23：判斷是否完成給定次數(shù)的迭代，若是則將所述更新的速度模型m1’作為最終速度模型m并輸出；否則返回子步驟s21，以更新的速度模型m1’作為子步驟s21中的中間速度模型。

與現(xiàn)有技術(shù)相比，本發(fā)明的有益效果是：采用了基于波動理論的全波形反演算子，相對于傳統(tǒng)的vsp數(shù)據(jù)建模方式，建模精度更高，能夠更加清晰的刻畫模型的構(gòu)造細(xì)節(jié)；引入地面地震構(gòu)造約束濾波算子，不僅克服了傳統(tǒng)全波形反演的數(shù)據(jù)不完備引起的建模噪音問題，以及低頻缺失造成的反演陷入局部極值問題，而且在一定程度上解決了vsp數(shù)據(jù)角度覆蓋有限的問題，提高了建模的角度范圍；地面地震構(gòu)造約束的vsp數(shù)據(jù)全波形反演與不加地面地震構(gòu)造約束的vsp數(shù)據(jù)全波形反演串聯(lián)策略的實(shí)施，可以進(jìn)一步提高vsp數(shù)據(jù)建模的精度。

附圖說明

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

圖1示出了根據(jù)示例性實(shí)施例的基于地面地震構(gòu)造約束的vsp數(shù)據(jù)全波形反演建模方法的流程圖；

圖2示出了根據(jù)示例性實(shí)施例的基于地面地震構(gòu)造約束的vsp數(shù)據(jù)全波形反演建模方法的步驟1的詳細(xì)流程圖；

圖3示出了根據(jù)示例性實(shí)施例的基于地面地震構(gòu)造約束的vsp數(shù)據(jù)全波形反演建模方法的步驟2的詳細(xì)流程圖；

圖4示出了根據(jù)示例性實(shí)施例的復(fù)雜模型的逆時偏移結(jié)果；

圖5示出了根據(jù)示例性實(shí)施例的復(fù)雜構(gòu)造模型偏移結(jié)果及其結(jié)構(gòu)特征張量算子的方向性示意圖；

圖6示出了根據(jù)現(xiàn)有技術(shù)的復(fù)雜構(gòu)造的真實(shí)速度模型；

圖7示出了根據(jù)現(xiàn)有技術(shù)的常梯度初始速度模型；

圖8和圖9分別顯示根據(jù)現(xiàn)有技術(shù)，采用vsp數(shù)據(jù)以常梯度初始速度模型進(jìn)行傳統(tǒng)全波形反演所得到的梯度結(jié)果和最終速度模型結(jié)果；

圖10示出了根據(jù)示例性實(shí)施例的采用地面地震構(gòu)造約束濾波算子濾波后的速度模型梯度；

圖11示出了根據(jù)示例性實(shí)施例的基于地面地震構(gòu)造約束的vsp數(shù)據(jù)全波形反演速度建模結(jié)果；以及

圖12示出了根據(jù)示例性實(shí)施例的基于地面地震構(gòu)造約束的vsp數(shù)據(jù)全波形反演聯(lián)合無構(gòu)造約束的vsp數(shù)據(jù)全波形反演的速度建模結(jié)果。

具體實(shí)施方式

下面將參照附圖更詳細(xì)地描述本公開的優(yōu)選實(shí)施例。雖然附圖中顯示了本公開的優(yōu)選實(shí)施例，然而應(yīng)該理解，可以以各種形式實(shí)現(xiàn)本公開而不應(yīng)被這里闡述的實(shí)施例所限制。相反，提供這些實(shí)施例是為了使本公開更加透徹和完整，并且能夠?qū)⒈竟_的范圍完整地傳達(dá)給本領(lǐng)域的技術(shù)人員。

圖1-3示出了根據(jù)示例性實(shí)施例的基于地面地震構(gòu)造約束的vsp數(shù)據(jù)全波形反演建模方法的流程圖。如圖1-3所示，根據(jù)示例性實(shí)施例的基于地面地震構(gòu)造約束的vsp數(shù)據(jù)全波形反演建模方法包括以下步驟：

步驟s1：基于地面地震構(gòu)造約束進(jìn)行vsp數(shù)據(jù)全波形反演，獲得中間速度模型m1

步驟s1包括以下子步驟s11至s16：

子步驟s11：采用初始速度模型m0以及地面地震數(shù)據(jù)進(jìn)行疊前深度偏移，獲得地下地質(zhì)構(gòu)造模式f(x)

采用初始速度模型m0以及地面地震數(shù)據(jù)進(jìn)行疊前深度域偏移，可以得到地下地質(zhì)構(gòu)造模式f(x)，其中x表示空間向量。在進(jìn)行第一次運(yùn)算時，初始速度模型m0可以選擇本領(lǐng)域現(xiàn)有的速度模型，例如常梯度速度模型，在經(jīng)過后續(xù)的迭代更新后，m0是經(jīng)過更新的速度模型。在進(jìn)行疊前深度偏移時采用的成像算子可以是克?；舴?kirchhoff)偏移算子，也可以是單程波偏移算子或逆時偏移算子。圖4顯示根據(jù)示例性實(shí)施例的復(fù)雜模型的逆時偏移結(jié)果，經(jīng)過逆時偏移可以清楚獲取地下地質(zhì)構(gòu)造模式。

子步驟s12：構(gòu)建地下地質(zhì)構(gòu)造模式f(x)的結(jié)構(gòu)特征張量算子

基于圖像結(jié)構(gòu)特征理論，在圖像學(xué)領(lǐng)域中，可以把子步驟s12中通過偏移成像得到的地下地質(zhì)構(gòu)造模式f(x)看作一幅圖像，該圖像具有方向性紋理結(jié)構(gòu)，通過求解紋理結(jié)構(gòu)可以得到結(jié)構(gòu)特征張量算子，其包括地下地質(zhì)構(gòu)造模式f(x)的在每個點(diǎn)的一階對稱半正定結(jié)構(gòu)張量矩陣的特征向量和特征值，用于表征圖像的局部結(jié)構(gòu)信息，即圖像特征的主要變化方向及變化速率。

對于二維偏移成像數(shù)據(jù)體，它在每個點(diǎn)的一階對稱半正定結(jié)構(gòu)張量矩陣表述為公式(1)：

根據(jù)矩陣的特征分解定理把結(jié)構(gòu)張量矩陣進(jìn)行特征值分解，得到公式(2)：

t＝λuuu^t+λvvv^t(2)

其中，λu≥λv≥0表示特征值，u，v分別為兩個特征值對應(yīng)的特征向量。特征值及特征向量刻畫了局部圖像的特征。當(dāng)兩個特征值都為0時，圖像區(qū)域?yàn)槌?shù)；當(dāng)兩個特征值相等并都大于零時，表示圖像為各向同性；當(dāng)兩個特征值不相等并均大于零時，表示圖像為各向異性，且局部圖像區(qū)域內(nèi)存在主方向，較大特征值對應(yīng)的特征向量表示圖像梯度變化最快的方向，其與圖像局部的線性特征方向相垂直，相對而言較小特征值對應(yīng)的特征向量代表了圖像局部的線性方向。圖5顯示了復(fù)雜構(gòu)造模型偏移成像結(jié)果及其結(jié)構(gòu)特征張量算子的方向性示意圖。

子步驟s13：基于時間域數(shù)據(jù)匹配的全波形反演構(gòu)建vsp數(shù)據(jù)殘差，基于vsp數(shù)據(jù)殘差計算速度模型m0的梯度gv(x)

采用聲波正演模擬方程模擬vsp數(shù)據(jù)和正向傳播波場，然后計算vsp觀測數(shù)據(jù)和模擬的vsp數(shù)據(jù)在時間空間域中的數(shù)據(jù)殘差，即vsp數(shù)據(jù)殘差。其中聲波正演模擬方程表述為公式(3)：

其中：

p(x,t)＝[vx(x,t),vy(x,t),vz(x,t),p(x,t)]^t代表正向傳播波場，

s(x,t)＝[0,0,0,s(xs,t)]^t代表震源矢量，xs為震源位置坐標(biāo)，

b為地下介質(zhì)密度的倒數(shù)，

κ＝v²/b為體積模量，

v是地下介質(zhì)的速度參數(shù)。

采用聲波正演模擬方程模擬vsp數(shù)據(jù)后，可以計算vsp觀測數(shù)據(jù)和模擬的vsp數(shù)據(jù)在時間空間域中的數(shù)據(jù)殘差，即vsp數(shù)據(jù)殘差。然后，根據(jù)vsp數(shù)據(jù)殘差求解逆時傳播的伴隨波場，以及基于所述伴隨波場求解速度模型m0的梯度gv(x)，速度模型m0的梯度gv(x)的表達(dá)式為公式(4)：

其中，ns表示震源個數(shù)，t代表波場傳播的終止時刻，q(x)是逆時傳播的伴隨波場的分量。

子步驟s14：基于各向異性擴(kuò)散方程構(gòu)建沿地下地質(zhì)構(gòu)造模式f(x)的結(jié)構(gòu)特征張量算子平滑的地面地震構(gòu)造約束濾波算子，并對子步驟s13中得到的速度模型m0的梯度gv(x)沿地面地震構(gòu)造約束方向進(jìn)行平滑濾波，得到基于地面地震構(gòu)造約束的速度模型梯度gv(x)

實(shí)際地震資料往往缺少道集或存在較強(qiáng)的噪音干擾，而且其低頻成分信息往往不可靠，因而采用這種數(shù)據(jù)進(jìn)行傳統(tǒng)全波形反演時會使反演陷入局部極值，反演結(jié)果中出現(xiàn)較強(qiáng)的噪音，而這種噪音是每次迭代的梯度中的噪音引起的。

圖8和圖9分別顯示了采用vsp數(shù)據(jù)以常梯度初始速度模型(圖7)進(jìn)行傳統(tǒng)全波形反演所得到的梯度結(jié)果和最終速度模型結(jié)果，從圖中可以看出梯度圖像中存在較強(qiáng)的噪音干擾，而反演結(jié)果陷入了局部極值，與真實(shí)速度模型(圖6)差距較遠(yuǎn)。

為了提升vsp數(shù)據(jù)反演結(jié)果的質(zhì)量，改善有限的角度覆蓋范圍，需要對速度模型m0的梯度gv(x)進(jìn)行預(yù)處理。在此預(yù)處理方式采用沿地面地震構(gòu)造約束方向?qū)μ荻萭v(x)進(jìn)行平滑濾波，地震構(gòu)造約束方向即子步驟12中求解出的較小特征值對應(yīng)的特征向量代表的局部線性方向。

地面地震構(gòu)造約束濾波算子采用近似的各向異性擴(kuò)散方程來構(gòu)建，如公式(5)所示：

其中，d(x)表示擴(kuò)散張量，為保證擴(kuò)散方向沿著地面地震構(gòu)造方向進(jìn)行，擴(kuò)散張量d(x)應(yīng)與子步驟s12求解出的結(jié)構(gòu)特征張量算子的特征值一致，gv(x)表示經(jīng)過濾波后的速度模型梯度，α表示光滑因子，α越大，原始圖像沿構(gòu)造方向擴(kuò)散程度越大，濾波后的圖像越平滑，當(dāng)α＝0時代表擴(kuò)散程度為0，此時并未對原始梯度做任何預(yù)處理操作。

經(jīng)過雙線性變換之后可得到近似各向異性方程的緊致矩陣算子表達(dá)式：

ηgv(x)＝gv(x)

其中η是與擴(kuò)散張量有關(guān)的稀疏矩陣。

采用共軛梯度法求解上述線性方程組，可得到沿地面地震構(gòu)造約束方向經(jīng)過濾波的速度模型梯度gv(x)。圖10顯示了采用地面地震構(gòu)造約束濾波算子濾波后的速度模型梯度，其不僅消除了原始梯度的噪音，而且保留了地面地震的構(gòu)造形態(tài)。

子步驟s15：通過p-l-bfgs算法更新速度模型m0，獲得更新的速度模型m0’

在對梯度進(jìn)行濾波之后，采用非精確線性搜索法求解對速度模型進(jìn)行更新的步長，采用最優(yōu)化理論的p-l-bfgs優(yōu)化迭代算法近似求解海塞矩陣的逆矩陣，并基于該逆矩陣更新速度模型m0，獲得更新的速度模型m0’。

子步驟s16：判斷是否完成給定次數(shù)的迭代，若是則繼續(xù)到步驟s2；若未完成給定次數(shù)的迭代，則返回子步驟s11，以更新的速度模型作為子步驟s11中的初始的速度模型。

圖11為基于地面地震構(gòu)造約束的vsp數(shù)據(jù)全波形反演速度建模結(jié)果。

步驟2：針對步驟1中獲得的中間速度模型m1，進(jìn)行沒有地面地震構(gòu)造約束的vsp數(shù)據(jù)全波形反演，得到最終速度模型m

步驟2包括以下子步驟：

子步驟s21：基于時間域數(shù)據(jù)匹配的全波形反演算子構(gòu)建vsp數(shù)據(jù)殘差，基于vsp數(shù)據(jù)殘差計算中間速度模型m1的梯度；

子步驟s22：基于最優(yōu)化理論，通過p-l-bfgs算法更新速度模型m1，獲得更新的速度模型m1’；

子步驟s23：判斷是否完成給定次數(shù)的迭代，將更新的速度模型m1’作為最終速度模型m并輸出；若未完成給定次數(shù)的迭代，則返回子步驟s21，以更新的速度模型m1’作為子步驟s21中的中間速度模型。

子步驟s21至s23的細(xì)節(jié)與前述子步驟s13、s15和s16相同，只是基于中間速度模型m1來執(zhí)行，在此不再贅述。

圖12為聯(lián)合反演得到的最終速度建模結(jié)果，其速度模型的細(xì)節(jié)刻畫更為清晰、更加準(zhǔn)確。

上述技術(shù)方案只是本發(fā)明的一種實(shí)施例，對于本領(lǐng)域內(nèi)的技術(shù)人員而言，在本發(fā)明公開的原理的基礎(chǔ)上，很容易做出各種類型的改進(jìn)或變形，而不僅限于本發(fā)明上述具體實(shí)施例的描述，因此前面的描述只是優(yōu)選的，而并不具有限制性的意義。