基于地面地震構造約束的VSP數據全波形反演建模方法與流程

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

背景技術：

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

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

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

技術實現要素：

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

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

步驟s1：基于地面地震構造約束進行vsp數據全波形反演，獲得中間速度模型m1；

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

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

子步驟s12：構建所述地下地質構造模式f(x)的結構特征張量算子；

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

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

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

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

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

優選地，采用克希霍夫偏移算子、單程波偏移算子或者逆時偏移算子進行疊前深度偏移。

優選地，所述結構特征張量算子包括所述地下地質構造模式f(x)在每個點的一階對稱半正定結構張量矩陣的特征向量和特征值。

優選地，采用聲波正演模擬方程模擬vsp數據和正向傳播波場，其中所述聲波正演模擬方程表示為公式(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為震源位置坐標，

b代表地下介質密度的倒數，

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

v代表地下介質的速度參數。

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

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

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

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

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

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

優選地，所述步驟2包括以下子步驟s21至s23：

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

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

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

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

附圖說明

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

圖1示出了根據示例性實施例的基于地面地震構造約束的vsp數據全波形反演建模方法的流程圖；

圖2示出了根據示例性實施例的基于地面地震構造約束的vsp數據全波形反演建模方法的步驟1的詳細流程圖；

圖3示出了根據示例性實施例的基于地面地震構造約束的vsp數據全波形反演建模方法的步驟2的詳細流程圖；

圖4示出了根據示例性實施例的復雜模型的逆時偏移結果；

圖5示出了根據示例性實施例的復雜構造模型偏移結果及其結構特征張量算子的方向性示意圖；

圖6示出了根據現有技術的復雜構造的真實速度模型；

圖7示出了根據現有技術的常梯度初始速度模型；

圖8和圖9分別顯示根據現有技術，采用vsp數據以常梯度初始速度模型進行傳統全波形反演所得到的梯度結果和最終速度模型結果；

圖10示出了根據示例性實施例的采用地面地震構造約束濾波算子濾波后的速度模型梯度；

圖11示出了根據示例性實施例的基于地面地震構造約束的vsp數據全波形反演速度建模結果；以及

圖12示出了根據示例性實施例的基于地面地震構造約束的vsp數據全波形反演聯合無構造約束的vsp數據全波形反演的速度建模結果。

具體實施方式

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

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

步驟s1：基于地面地震構造約束進行vsp數據全波形反演，獲得中間速度模型m1

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

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

采用初始速度模型m0以及地面地震數據進行疊前深度域偏移，可以得到地下地質構造模式f(x)，其中x表示空間向量。在進行第一次運算時，初始速度模型m0可以選擇本領域現有的速度模型，例如常梯度速度模型，在經過后續的迭代更新后，m0是經過更新的速度模型。在進行疊前深度偏移時采用的成像算子可以是克希霍夫(kirchhoff)偏移算子，也可以是單程波偏移算子或逆時偏移算子。圖4顯示根據示例性實施例的復雜模型的逆時偏移結果，經過逆時偏移可以清楚獲取地下地質構造模式。

子步驟s12：構建地下地質構造模式f(x)的結構特征張量算子

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

對于二維偏移成像數據體，它在每個點的一階對稱半正定結構張量矩陣表述為公式(1)：

根據矩陣的特征分解定理把結構張量矩陣進行特征值分解，得到公式(2)：

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

其中，λu≥λv≥0表示特征值，u，v分別為兩個特征值對應的特征向量。特征值及特征向量刻畫了局部圖像的特征。當兩個特征值都為0時，圖像區域為常數；當兩個特征值相等并都大于零時，表示圖像為各向同性；當兩個特征值不相等并均大于零時，表示圖像為各向異性，且局部圖像區域內存在主方向，較大特征值對應的特征向量表示圖像梯度變化最快的方向，其與圖像局部的線性特征方向相垂直，相對而言較小特征值對應的特征向量代表了圖像局部的線性方向。圖5顯示了復雜構造模型偏移成像結果及其結構特征張量算子的方向性示意圖。

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

采用聲波正演模擬方程模擬vsp數據和正向傳播波場，然后計算vsp觀測數據和模擬的vsp數據在時間空間域中的數據殘差，即vsp數據殘差。其中聲波正演模擬方程表述為公式(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為震源位置坐標，

b為地下介質密度的倒數，

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

v是地下介質的速度參數。

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

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

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

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

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

為了提升vsp數據反演結果的質量，改善有限的角度覆蓋范圍，需要對速度模型m0的梯度gv(x)進行預處理。在此預處理方式采用沿地面地震構造約束方向對梯度gv(x)進行平滑濾波，地震構造約束方向即子步驟12中求解出的較小特征值對應的特征向量代表的局部線性方向。

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

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

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

ηgv(x)＝gv(x)

其中η是與擴散張量有關的稀疏矩陣。

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

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

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

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

圖11為基于地面地震構造約束的vsp數據全波形反演速度建模結果。

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

步驟2包括以下子步驟：

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

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

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

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

圖12為聯合反演得到的最終速度建模結果，其速度模型的細節刻畫更為清晰、更加準確。

上述技術方案只是本發明的一種實施例，對于本領域內的技術人員而言，在本發明公開的原理的基礎上，很容易做出各種類型的改進或變形，而不僅限于本發明上述具體實施例的描述，因此前面的描述只是優選的，而并不具有限制性的意義。