基于交叉梯度約束的地震初至縱波與瑞雷面波聯(lián)合反演

傅磊， 劉四新

吉林大學(xué) 地球探測科學(xué)與技術(shù)學(xué)院， 長春 130026

?

基于交叉梯度約束的地震初至縱波與瑞雷面波聯(lián)合反演

傅磊， 劉四新

吉林大學(xué) 地球探測科學(xué)與技術(shù)學(xué)院， 長春 130026

本文提出了一種初至縱波(P波)與瑞雷面波的交叉梯度聯(lián)合反演策略.通過對初至P波進(jìn)行全波形反演可以獲得近地表P波速度結(jié)構(gòu)；通過對僅含瑞雷面波信息的地震數(shù)據(jù)轉(zhuǎn)換到頻率-波數(shù)域進(jìn)行加窗振幅波形反演(Windowed-Amplitude Waveform Inversion，w-AWI)可獲得近地表橫波(S波)速度結(jié)構(gòu).在二者反演的目標(biāo)函數(shù)中均加入P波速度和S波速度的交叉梯度作為正則化約束項，使得在反演過程中P波速度和S波速度相互制約，相互約束，從而實(shí)現(xiàn)對地震初至P波與瑞雷面波的聯(lián)合反演.數(shù)值模擬結(jié)果表明交叉梯度聯(lián)合反演可以提高S波速度反演分辨率，而P波速度反演結(jié)果并沒有得到提高.實(shí)際資料的反演結(jié)果表明，交叉梯度聯(lián)合反演能夠獲得更加可信的近地表速度結(jié)構(gòu).

初至縱波； 瑞雷面波； 交叉梯度； 聯(lián)合反演

1 引言

在勘探地震的偏移成像中，近地表速度(P波速度和S波速度)結(jié)構(gòu)(劉伊克等, 2001)的精確度制約著深層地質(zhì)結(jié)構(gòu)的構(gòu)建.目前為止，大多數(shù)面波反演(Xia et al., 1999; 石耀霖和金文, 1995; 胡家富等, 1999; 李明明和何玉梅, 2011)基于層狀模型假設(shè)，這一假設(shè)在實(shí)際復(fù)雜的地質(zhì)環(huán)境中并不適用.為了克服這個層狀模型假設(shè)，許多研究人員嘗試采用波形反演的思路對面波進(jìn)行反演.Sch?fer等(2013)采用全波形反演方法對兩個實(shí)際數(shù)據(jù)進(jìn)行測試，結(jié)果表明反演過程收斂速度極慢以及反演結(jié)果往往陷入局部最小值.Solano等(2013)提出了一種振幅波形反演(Amplitude Waveform Inversion)的方法，其目標(biāo)函數(shù)為觀測數(shù)據(jù)振幅譜與預(yù)測數(shù)據(jù)振幅譜的最小二范數(shù)，而傳統(tǒng)的全波形反演(Tarantola,1984; Pratt et al., 1998; Virieux and Operto, 2009)的目標(biāo)函數(shù)為觀測數(shù)據(jù)與預(yù)測數(shù)據(jù)的最小二范數(shù)，結(jié)果表明振幅波形反演具有更好的收斂性質(zhì).這是因為傳統(tǒng)的全波形反演匹配的是所有波形，而振幅波形反演只匹配振幅譜，振幅譜僅僅是所有波形信息中的一部分，這種簡化的數(shù)據(jù)(Luo and Schuster, 1991a, 1991b)使得振幅波形反演具有更快的收斂速度和更好的收斂效果.Solano等(2014)提出了一種加窗振幅波形反演方法(Windowed-Amplitude Waveform Inversion)，與振幅波形反演相比，該方法更加穩(wěn)定.然而與傳統(tǒng)全波形反演相比，振幅波形反演是以犧牲分辨率為代價的.為了部分恢復(fù)振幅波形反演丟失的分辨率，本文提出了一種初至P波與瑞雷面波的交叉梯度聯(lián)合反演(于鵬等, 2009; 彭淼等, 2013; 李兆祥等, 2015)方法，采用交叉梯度同時作為初至P波反演和瑞雷面波加窗振幅波形反演的約束項，將分辨率較低的瑞雷面波振幅波形反演與高分辨率的P波波形反演結(jié)合起來，可以獲得更高分辨率的S波速度結(jié)構(gòu).

地球物理數(shù)據(jù)的聯(lián)合反演可以分為兩大類.第一類是基于不同參數(shù)之間的巖石物理關(guān)系，例如飽和度和孔隙度(Hoversten et al., 2006)，電導(dǎo)率和介電常數(shù)(Kowalsky et al., 2005).例如，P波和S波速度之間的關(guān)系是依賴于介質(zhì)屬性的.Castagna等(1985)指出，在砂巖中P波與S波比約為1.5～1.7，在頁巖中該比值大于2.這就意味著在聯(lián)合反演中，采用不同參數(shù)之間錯誤的經(jīng)驗關(guān)系可能會帶來嚴(yán)重的耦合誤差.

第二類聯(lián)合反演利用不同參數(shù)之間的結(jié)構(gòu)相似性，同一地質(zhì)模型可以用不同的物性參數(shù)來描述，例如密度、P波速度、電導(dǎo)率等，這些參數(shù)存在著結(jié)構(gòu)相似性.Gallardo和Meju(2003)采用交叉梯度聯(lián)合反演方法成功地對直流電法和地震走時數(shù)據(jù)進(jìn)行了反演.Colombo和De Stefano(2007)等利用交叉梯度聯(lián)合反演方法同時對地震、重力和大地電磁數(shù)據(jù)進(jìn)行了反演.本文采用第二類聯(lián)合反演策略對初至地震P波和瑞雷面波進(jìn)行交叉梯度聯(lián)合反演.

對于初至縱波，本文采用全波形反演獲得地下P波速度結(jié)構(gòu)，采用加窗振幅波形反演方法對瑞雷面波進(jìn)行反演，同時采用P波速度和S波速度的交叉梯度作為二者的約束項.結(jié)果顯示，與單獨(dú)的瑞雷面波振幅波形反演相比，交叉梯度聯(lián)合反演的S波速度具有更高的精度.本文第2節(jié)和第3節(jié)將分別給出交叉梯度聯(lián)合反演的方法原理以及數(shù)值模擬結(jié)果.之后，將該聯(lián)合反演方法應(yīng)用于實(shí)際數(shù)據(jù).最后是結(jié)論.

2 方法原理

對于地震數(shù)據(jù)，初至P波可以通過求解二維常密度聲波方程來進(jìn)行模擬，

(1)

其中v是縱波速度，p(x,t)是壓力場，f(x,t)是震源項.

另一方面，瑞雷面波可以通過求解彈性波方程進(jìn)行模擬獲得，

(2)

其中ρ是密度，ui是第i個位移分量，cijkl是剛度矩陣，si是震源項.

2.1 交叉梯度

交叉梯度的定義為兩個模型參數(shù)梯度的叉積：

(3)

其中mP代表P波速度，而mS代表S波速度.交叉梯度正則項是基于這樣的事實(shí)，當(dāng)兩個模型參數(shù)的交叉梯度達(dá)到最小化時，它們的結(jié)構(gòu)相似度將達(dá)到最大化.如果在P波和S波速度中存在相同的邊界，且該邊界的位置和方向一致，那么它們的交叉梯度將等于零.

對于二維問題，方程(3)可以改寫為：

(4)

如圖1a為P波速度模型，波速從左到右逐漸增大；圖1b為S波速度模型，一個低速方形異常體、一個高速圓形異常體及一個高速方形異常體嵌于一均勻背景速度模型；圖1c為二者的交叉梯度值.以左側(cè)低速異常體為例，雖然在S波速度中四條邊界處的梯度值均不為零；然而在P波速度中低速方形所在位置的左、右邊界處的梯度與S波速度中對應(yīng)的梯度平行，故而交叉梯度值為零，僅在上下邊界處其交叉梯度值不為零.

圖1 (a)P波速度模型，(b)S波速度模型，(c)相應(yīng)P波速度與S波速度的交叉梯度Fig.1 (a) P-wave velocity model, (b) S-wave velocity model, (c) Corresponding cross-gradients of P-wave velocity and S-wave velocity models

2.2 初至P波全波形反演

波動方程波形反演的目標(biāo)函數(shù)為：

‖Δdrs‖2,

(5)

其中Δdrs=dcal(xr,xs)-dobs(xr,xs)是初至P波殘差，它表示為預(yù)測數(shù)據(jù)與觀測數(shù)據(jù)之差.對于P波速度，它可以通過任何一種梯度類方法進(jìn)行迭代更新，例如通過共軛梯度(Hestenes and Stiefel, 1952; 胡祖志, 2006)方法對P波速度進(jìn)行更新：

(6)

(7)

(8)

(9)

其中，

(10)

2.3 瑞雷面波加窗振幅波形反演

加窗振幅波形反演(Solano et al., 2014)是在全波形反演的基礎(chǔ)上衍生出來的，它的目標(biāo)函數(shù)為：

(11)

(12)

(13)

(14)

(15)

其中

(16)

δd(xr,t|xs)=

(17)

2.4 交叉梯度聯(lián)合反演目標(biāo)函數(shù)及梯度

為了提高反演結(jié)果的精度，本文將初至P波與S波速度進(jìn)行聯(lián)合同步反演，對于初至P波，其全波形聯(lián)合反演目標(biāo)函數(shù)為：

(18)

同理對于瑞雷面波，其加窗振幅波形聯(lián)合同步反演的目標(biāo)函數(shù)則為：

ΦS(vS,vP)=

(19)

其中λ1和λ2分別為各自的正則化系數(shù)，當(dāng)這兩個正則化系數(shù)同時置為零時，方程(18)和(19)退化為無耦合的單獨(dú)反演方程.對于交叉梯度聯(lián)合反演，關(guān)鍵問題是計算目標(biāo)函數(shù)(18)和(19)分別對P波速度和S波速度的偏導(dǎo)，即獲得各自的梯度.注意到交叉梯度聯(lián)合反演目標(biāo)函數(shù)均包含兩項，原始目標(biāo)函數(shù)項和正則約束項，而原始目標(biāo)函數(shù)對應(yīng)的梯度已知，故關(guān)鍵問題是求取正則約束項分別對P波速度和S波速度的梯度，根據(jù)Hu等(2009)，經(jīng)過數(shù)學(xué)推導(dǎo)，正則項的梯度分別表示為：

(20)

(21)

3 數(shù)值實(shí)驗

3.1 目標(biāo)函數(shù)形態(tài)

全波形反演是病態(tài)的非線性反演問題，非線性、跳周、子波估計不準(zhǔn)等現(xiàn)象往往使得波形反演問題陷入局部最小值.本文選取加窗振幅波形反演對瑞雷面波進(jìn)行反演，是因為與全波形反演相比，它更容易收斂到全局最小值.為了比較全波形反演、振幅波形反演及加窗振幅波形反演三種目標(biāo)函數(shù)的形態(tài)，首先構(gòu)建一簡單模型如圖2所示，圖2為S波速度模型，一條帶狀異常體速度vS1嵌于均勻背景速度vS0中，真實(shí)速度模型和初始速度模型同時取決于vS1和vS0.假設(shè)真實(shí)背景速度大小為800 m·s-1，真實(shí)異常體速度大小為960 m·s-1；通過改變vS1和vS0獲得不同的初始模型，其中vS1的變化范圍為400～1400 m·s-1，vS0的變化范圍為600～1400 m·s-1.同時P波速度設(shè)為2000 m·s-1，密度設(shè)為1000 kg·m-3.圖中紅色星號代表震源所在位置，采用中心頻率為20 Hz的雷克子波作為震源；檢波器均勻分布于地表.圖3為三種不同目標(biāo)函數(shù)的形態(tài)，十字虛白線交叉點(diǎn)為真模型；圖3a為全波形反演目標(biāo)函數(shù)形態(tài)；圖3b為振幅波形反演目標(biāo)函數(shù)形態(tài)；圖3c為加窗振幅波形反演目標(biāo)函數(shù)形態(tài)，可見加窗目標(biāo)函數(shù)具有最寬的全局最小范圍，最容易收斂到全局最小.相比于初至P波，瑞雷面波具有更強(qiáng)的非線性性質(zhì)，故而本文對瑞雷面波采用加窗振幅波形反演方法；而對初至P波采用全波形反演方法.

3.2 數(shù)值算例

首先將交叉梯度聯(lián)合反演方法應(yīng)用于簡單數(shù)值模型圖4.其中圖4a是P波速度模型，圖4b為S波速度模型.該模型包括三層巖層，從地表到深部其速度逐漸遞增，且不同巖層之間界面為不規(guī)則界面.模型被離散為30×120的等體積立方體網(wǎng)格，網(wǎng)格間距為3 m.圖中紅色X表示震源所在位置，其間距為6 m；地表均勻布設(shè)120個檢波器，檢波器間距為3 m.源子波為中心頻率為20 Hz的雷克子波.初始速度模型見圖5，它是將真實(shí)速度模型進(jìn)行50次高斯濾波圓滑獲得.

圖2 用于計算目標(biāo)函數(shù)形態(tài)的S波速度模型，v S0代表背景速度，v S1代表異常體速度Fig.2 S velocity models used to compute different target function shapes. vS0 represents background S velocity. vS1 represents abnormal S velocity

圖3 不同方法歸一化目標(biāo)函數(shù)(a) 全波形反演； (b) 振幅波形反演； (c) 加窗振幅波形反演.Fig.3 Normalized objective functions by different methods

圖4 真實(shí)速度模型其中紅色X代表震源所在位置.(a) P波速度模型；(b) S波速度模型.Fig.4 Real velocity modelsRed X indicates source location. (a) P-wave velocity model; (b) S-wave velocity model.

圖5 (a) 初始P波速度模型， (b) 初始S波速度模型Fig.5 Initial P-wave velocity model (a) and S-wave velocity model (b)

圖6 初至P波全波形反演結(jié)果(a) 單獨(dú)反演第10次迭代P波速度結(jié)果； (b) 聯(lián)合反演第10次迭代P波速度結(jié)果.Fig.6 Full-waveform inversion results of first arrival P waves(a) Single inversion of P-wave velocity after 10 iterations; (b) Joint inversion of P-wave velocity after 10 iterations.

圖7 瑞雷面波反演結(jié)果(a) 單獨(dú)反演第10次迭代S波速度結(jié)果； (b) 聯(lián)合反演第10次迭代S波速度結(jié)果.Fig.7 Inversion results for Rayleigh waves(a) Single inversion of S-wave velocity after 10 iterations; (b) Joint inversion of S-wave velocity after 10 iterations.

圖6a是經(jīng)過10次迭代后單獨(dú)反演的P波速度結(jié)果，由圖可知，P波反演結(jié)果中基本可以識別出三層巖層，其中淺層形態(tài)與真實(shí)模型吻合較好；底層高速層的形態(tài)及反演結(jié)果與真實(shí)模型基本一致；中間層的形態(tài)與真實(shí)模型基本一致，然而其反演結(jié)果與真實(shí)模型存在一定誤差.圖6b是經(jīng)過10次迭代后交叉梯度聯(lián)合反演P波速度結(jié)果，其結(jié)果同樣并沒有得到提高.

圖7所示為經(jīng)過10次迭代反演S波速度結(jié)果.圖7a為單獨(dú)反演S波速度結(jié)果，由圖可知，淺層速度結(jié)構(gòu)刻畫得較好，淺層與中間層的界面與真實(shí)模型基本吻合，但是中間層及底層的反演結(jié)果與真實(shí)模型存在較大誤差.這是因為瑞雷面波的穿透深度僅為其波長的1/3至1/2，而這里淺層的速度為800 m·s-1，對應(yīng)的波長為40 m，故而其探測深度應(yīng)在20 m范圍內(nèi).圖7b為交叉梯度聯(lián)合反演S波速度結(jié)果，可以看到：交叉梯度聯(lián)合反演S波速度結(jié)果中，淺層界面的連續(xù)性好于單獨(dú)反演結(jié)果.

4 實(shí)際資料

4.1 數(shù)據(jù)采集及預(yù)處理

此二維地震勘探的目的是為了探明地下某隱伏斷層在地表的位置.該二維地震剖面包括240個垂直分量檢波器，檢波器間距為5 m.時間采樣率為1 ms，記錄長度為3 s，每個地震記錄進(jìn)行32次疊加，總共240炮地震記錄.圖8a，8b，8c，8d分別為偏移距為10 m，30 m，50 m，70 m的共偏移距地震剖面.通過對比不同共偏移距剖面，我們認(rèn)為該隱伏斷層可能位于第55炮所在位置，即位于水平位置275 m處，即如圖中紅色虛線所示.故本文選取第21炮至第81炮地震記錄進(jìn)行反演.

圖9a為第21炮原始地震記錄，圖中紅色虛線所示為初至P波，而藍(lán)色虛線所示為瑞雷面波.通過對兩條虛線的斜率進(jìn)行估算，可以得到近地表的P波速度和S波速度分別約為2900 m·s-1和460 m·s-1.在數(shù)據(jù)預(yù)處理之前，首先手動提取初至P波和瑞雷面波走時，提取的走時將用于后續(xù)提取初至P波和瑞雷面波波形.

數(shù)據(jù)預(yù)處理的第一步是選取3～50 Hz濾波器對原始數(shù)據(jù)進(jìn)行帶通濾波；其次對數(shù)據(jù)進(jìn)行三維到二維轉(zhuǎn)換及幾何擴(kuò)散校正；第三步去除壞道并采用POCS(Projection onto Convex Sets)(Abma and Kabir, 2006)對缺失道進(jìn)行數(shù)據(jù)插值；第四步在時間域?qū)?shù)據(jù)進(jìn)行插值，使得時間采樣率滿足正演模擬需求；最后切除近偏移距地震道以及根據(jù)提取的走時將初至P波和瑞雷面波提取出.經(jīng)過預(yù)處理后的數(shù)據(jù)如圖9b所示.

圖8 經(jīng)過帶通濾波及道內(nèi)歸一化后的共偏移距剖面紅色虛線指示斷層的可能位置. (a) 偏移距為10 m； (b) 偏移距為30 m； (c) 偏移距為50 m； (d) 偏移距為70 m.Fig.8 Common offset sections from processed data after bandpass filtering and trace normalizationRed dash lines indicate the possible fault location.

圖9 (a) 第21炮原始數(shù)據(jù)記錄，紅色虛線為提取的P波走時，藍(lán)色虛線為提取的瑞雷面波走時； (b) 經(jīng)過數(shù)據(jù)預(yù)處理后的數(shù)據(jù)記錄Fig.9 (a) Raw record for 21th shot with picked traveltimes for first-arrival P waves (red dashed line) and Rayleigh waves (blue dashed line); (b) Data after pre-processing

4.2 反演結(jié)果

通過對瑞雷面波頻散曲線進(jìn)行反演(Zhang et al., 2015)獲得的S波速度可以作為聯(lián)合反演的初始S波速度模型，如圖10b所示，該初始S波速度模型被離散為30×121等體積立方體網(wǎng)格，網(wǎng)格間距為2.5 m.從圖可知，在30 m深度范圍內(nèi)，從左到右，S波速度從420 m·s-1遞增到530 m·s-1，低速異常體位于水平方向100～150 m，這意味著斷層地表水平位置可能位于150～170 m；深層為一高速帶，故在45 m以下的區(qū)域設(shè)置為650 m·s-1常速度.基于此S波速度模型，以及從圖9a估算的P波速度值(2900 m·s-1)，我們可以構(gòu)建初始P波速度模型，如圖10a所示.

圖11a是經(jīng)過15次迭代后單獨(dú)反演P波速度結(jié)果，圖11b是經(jīng)過15次迭代后交叉梯度聯(lián)合反演P波速度結(jié)果.由圖可知，相比于單獨(dú)反演，交叉梯度聯(lián)合反演P波速度結(jié)果并沒有得到提高.然而與初始P波速度相比，反演結(jié)果確有很大提高，其中低速異常體面積擴(kuò)大，從中可以推斷隱伏斷層水平位置可能位于230～280 m，這與共偏移距結(jié)果指示的隱伏斷層位置更加吻合.

圖12所示為經(jīng)過15次迭代反演S波速度結(jié)果.圖12a為單獨(dú)反演S波速度結(jié)果，可見在水平位置100～200 m處仍然存在一低速異常體，但是該異常體的形態(tài)與P波反演結(jié)果形態(tài)略有不同.并且在淺部水平位置250～300 m處存在一高速異常體，這與P波反演結(jié)果不同.圖12b為交叉梯度聯(lián)合反演S波速度結(jié)果，該結(jié)果中低速異常體的位置、形態(tài)與P波速度反演結(jié)果非常相似，而且淺部也未見局部高速異常體.單獨(dú)反演中，S波速度反演結(jié)果與P波速度反演結(jié)果相似度不高；而在交叉梯度聯(lián)合反演中，S波速度與P波速度反演結(jié)果相似度極高，其反演結(jié)果更真實(shí)可靠.

5 結(jié)論與討論

本文提出了一種初至P波和瑞雷面波聯(lián)合反演的方法，對初至P波進(jìn)行全波形反演，以及對瑞雷面波采用加窗振幅波形反演，同時將P波速度和S波速度的交叉梯度作為二者反演的正則化約束項.數(shù)值模擬結(jié)果表明，與單獨(dú)反演相比，聯(lián)合反演可以消除S波反演結(jié)果中的偽異常，提高S波反演精度；對于P波速度，聯(lián)合反演并不能獲得更好的反演結(jié)果.實(shí)際資料反演結(jié)果表明，交叉梯度聯(lián)合反演的S波速度與P波速度反演結(jié)果相似度更高，其反演結(jié)果更可信.

交叉梯度聯(lián)合反演的缺點(diǎn)是，當(dāng)真實(shí)的P波速度模型與S波速度模型在結(jié)構(gòu)上相似度不高時，采用交叉梯度作為反演約束項會得到錯誤的反演結(jié)果.

圖10 初始速度模型(a) 基于初始S波速度模型估算出的初始P波速度模型； (b) 通過面波頻散曲線反演獲得的S波速度作為初始S波速度.Fig.10 Initial velocity models(a) Initial P-wave velocity model estimated based on initial S-wave velocity model; (b) Initial S-wave velocity model computed from surface-wave dispersion inversion.

圖11 初至P波波動方程走時反演結(jié)果(a) 單獨(dú)反演第15次迭代P波速度結(jié)果； (b) 聯(lián)合反演第15次迭代P波速度結(jié)果.Fig.11 Inversion results for first-arrival P waves(a) Single inversion of P-wave velocity after 15 iterations; (b) Joint inversion of P-wave velocity after 15 iterations.

圖12 S波波動方程走時反演結(jié)果(a) 單獨(dú)反演第15次迭代S波速度結(jié)果； (b) 聯(lián)合反演第15次迭代S波速度結(jié)果.Fig.12 Inversion results for S-wave velocity(a) Single inversion of S-wave velocity after 15 iterations; (b) Joint inversion of S-wave velocity after 15 iterations.

Abma R, Kabir N. 2006. 3D interpolation of irregular data with a POCS algorithm.Geophysics, 71(6): E91-E97.

Castagna J P, Batzle M L, Eastwood R L. 1985. Relationships between compressional-wave and shear-wave velocities in clastic silicate rocks.Geophysics, 50(4): 571-581. Colombo D, De Stefano M. 2007. Geophysical modeling via simultaneous joint inversion of seismic, gravity, and electromagnetic data: Application to prestack depth imaging.TheLeadingEdge, 26(3): 326-331. Gallardo L A, Meju M A. 2003. Characterization of heterogeneous near-surface materials by joint 2D inversion of DC resistivity and seismic data.GeophysicalResearchLetters, 30(13): 1658. Hestenes M R, Stiefel E. 1952. Methods of conjugate gradients for solving linear systems.JournalofResearchoftheNationalBureauofStandards, 49(6): 409-436.

Hoversten G M, Cassassuce F, Gasperikova E, et al. 2006. Direct reservoir parameter estimation using joint inversion of marine seismic AVA and CSEM data.Geophysics, 71(3): C1-C13.

Hu J F, Duan Y K, Hu Y L, et al. 1999. Inversion of Shear-wave velocity structure in shallow soil from Rayleigh waves.ChineseJ.Geophys. (in Chinese), 42(3): 393-400.

Hu W Y, Abubakar A, Habashy T M. 2009. Joint electromagnetic and seismic inversion using structural constraints.Geophysics, 74(6): R99-R109.

Hu Z Z, Hu X Y, He Z X. 2006. Pseudo-Three-Dimensional magnetotelluric inversion using nonlinear conjugate gradients.ChineseJ.Geophys. (in Chinese), 49(4): 1226-1234.

Kowalsky M B, Finsterle S, Peterson J, et al. 2005. Estimation of field-scale soil hydraulic and dielectric parameters through joint inversion of GPR and hydrological data.WaterResourcesResearch, 41(11): W11425, doi:10.1029/2005WR004237.

Li M M, He Y M. 2011. Lithospheric structure beneath northeastern boundary region of the north China craton from Rayleigh wave dispersion inversion.ActaSeismologicaSinica(in Chinese), 33(2): 143-155.

Li Z X, Tan H D, Fu S S, et al. 2015. Two-dimensional synchronous inversion of TDIP with cross-gradient constraint.ChineseJ.Geophys. (in Chinese), 58(12): 4718-4726, doi: 10.6038/cjg20151232.

Liu Y K, Chang X, Wang H, et al. 2001. Estimation of near-surface velocity and seismic tomographic static corrections.ChineseJ.Geophys. (in Chinese), 44(2): 272-278. Luo Y, Schuster G T. 1991a. Wave equation inversion of skeletalized geophysical data.GeophysicalJournalInternational, 105(2): 289-294. Luo Y, Schuster G T. 1991b. Wave-equation traveltime inversion.Geophysics, 56(5): 645-653.

Peng M, Tan H D, Jiang M, et al. 2013. Three-dimensional joint

inversion of magnetotelluric and seismic travel time data with cross-gradient constraints.ChineseJ.Geophys. (in Chinese), 56(8): 2728-2738, doi: 10.6038/cjg20130821.

Pratt R G, Shin C, Hick G J. 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion.GeophysicalJournalInternational, 133(2): 341-362.

Sch?fer M, Groos L, Forbriger T, et al. 2013. 2D Full waveform inversion of recorded shallow seismic Rayleigh waves on a significantly 2D structure. ∥ Near Surface Geoscience 2013-19th EAGE European Meeting of Environmental and Engineering Geophysics. EAGE.Shi Y L, Jin W. 1995. Genetic algorithms inversion of lithospheric structure from surface wave dispersion.ChineseJ.Geophys. (ActaGeophysicaSinica) (in Chinese), 38(2): 189-198.

Solano C A P, Donno D, Chauris H. 2013. 2D surface wave inversion in the F-K domain. ∥ 75th EAGE Conference & Exhibition incorporating SPE EUROPEC. EAGE.

Solano C A P, Donno D, Chauris H. 2014. Alternative waveform inversion for surface wave analysis in 2-D media.GeophysicalJournalInternational, 198(3): 1359-1372.

Tarantola A. 1984. Inversion of seismic reflection data in the acoustic approximation.Geophysics, 49(8): 1259-1266.

Virieux J, Operto S. 2009. An overview of full-waveform inversion in exploration geophysics.Geophysics, 74(6): WCC1-WCC26.

Xia J H, Miller R D, Park C B. 1999. Estimation of near-surface shear-wave velocity by inversion of Rayleigh waves.Geophysics, 64(3): 691-700.Yu P, Dai M G, Wang J L, et al. 2009. Joint inversion of

magnetotelluric and seismic data based on random resistivity and velocity distributions.ChineseJ.Geophys. (in Chinese), 52(4): 1089-1097.Zhang Z D, Liu Y K, Schuster G. 2015. Wave equation inversion of skeletonized surface waves. ∥ 2015 SEG Annual Meeting. New Orleans, Louisiana: Society of Exploration Geophysicists.

附中文參考文獻(xiàn)

胡家富, 段永康, 胡毅力等. 1999. 利用Rayleigh波反演淺土層的剪切波速度結(jié)構(gòu). 地球物理學(xué)報, 42(3): 393-400.

胡祖志, 胡祥云, 何展翔. 2006. 大地電磁非線性共軛梯度擬三維反演. 地球物理學(xué)報, 49(4): 1226-1234.

李明明, 何玉梅. 2011. 利用瑞雷面波反演華北克拉通東北部邊界的巖石圈結(jié)構(gòu). 地震學(xué)報, 33(2): 143-155.

李兆祥, 譚捍東, 付少帥等. 2015. 基于交叉梯度約束的時間域激發(fā)極化法二維同步反演. 地球物理學(xué)報, 58(12): 4718-4726, doi: 10.6038/cjg20151232.

劉伊克, 常旭, 王輝等. 2001. 三維復(fù)雜地形近地表速度估算及地震層析靜校正. 地球物理學(xué)報, 44(2): 272-278.

彭淼, 譚捍東, 姜枚等. 2013. 基于交叉梯度耦合的大地電磁與地震走時資料三維聯(lián)合反演. 地球物理學(xué)報, 56(8): 2728-2738, doi: 10.6038/cjg20130821.

石耀霖, 金文. 1995. 面波頻散反演地球內(nèi)部構(gòu)造的遺傳算法. 地球物理學(xué)報, 38(2): 189-198.

于鵬, 戴明剛, 王家林等. 2009. 電阻率和速度隨機(jī)分布的MT與地震聯(lián)合反演. 地球物理學(xué)報, 52(4): 1089-1097.

(本文編輯 何燕)

Joint inversion of first arrival P waves and Rayleigh waves based on cross-gradient constraint

FU Lei, LIU Si-Xin

CollegeofGeo-explorationScienceandTechnology,JilinUniversity,Changchun130026,China

In seismic exploration, knowledge of the near-surface velocity can improve the construction of an accurate geological model. Until recently, most surface-wave inversion approaches assumed a layered model, which is not justified in a complex geological environment.To overcome this 1D limitation, researchers have begun to explore the benefits of 2D waveform inversion of surface waves. Windowed-Amplitude waveform inversion is a replacement to classical Full Wave form Inversion (FWI), which has faster convergence, however, lower spatial resolution in the tomogram is deduced. To partly recover this lost resolution, we propose a joint inversion of first arrival P-wave data and Rayleigh wave data that uses the cross-product between the P-velocity and S-velocity gradients as a regularization term. First, the P-velocity model is inverted by the full waveform inversion method, in which the first arrivals are those of the refraction and direct arrivals. Then, the amplitude spectra of the surface waves in the frequency-wavenumber domain are inverted for the S-velocity model with a cross-product constraint. This constraint regularization insists that the S-velocity gradient is closely parallel to that of the P-velocity gradient. Results show that cross-product regularization provides a significant reduction in artifacts in the S-velocity tomogram. In this case the less robust inversion of surface waves is combined with high resolution P-wave tomography to achieve better spatial resolution of the S-wave velocity tomogram.

First-arrival P wave; Rayleigh wave; Cross-gradient; Joint inversion

10.6038/cjg20161209.

國家自然科學(xué)基金項目(3A515AV44423),中國博士后科學(xué)基金(801151010423)共同資助.

傅磊，男，1985年生，講師，主要從事地球物理正反演及成像研究.E-mail: leifu@jlu.edu.cn

10.6038/cjg20161209

P631

2016-02-01，2016-11-14收修定稿

傅磊，劉四新. 2016. 基于交叉梯度約束的地震初至縱波與瑞雷面波聯(lián)合反演. 地球物理學(xué)報，59(12):4464-4472,

Fu L, Liu S X. 2016. Joint inversion of first arrival P waves and Rayleigh waves based on cross-gradient constraint.ChineseJ.Geophys. (in Chinese),59(12):4464-4472,doi:10.6038/cjg20161209.