彭恒初,胡家富,楊海燕,文麗敏
云南大學地球物理系,昆明 650091
接收函數(shù)反演地殼S波速度結(jié)構(gòu)的有效約束方法
彭恒初,胡家富*,楊海燕,文麗敏
云南大學地球物理系,昆明 650091
本文通過對徑向接收函數(shù)和垂直向接收函數(shù)進行低通濾波,獲取了S波視速度隨低通濾波參數(shù)的變化曲線,然后利用經(jīng)驗關(guān)系將它轉(zhuǎn)換成了臺站下方的S波速度結(jié)構(gòu),并以此作為接收函數(shù)反演的初始模型.理論數(shù)值實驗表明:由于初始模型的S波速度值提供了有效的約束,即使Moho面深度并不準確,但反演迭代過程還是快速地向真解逼近.另外,通過給觀測波形加入10%的噪聲,在保持S波速度不變的情況下,分別對波速比進行5%的正負擾動(即泊松比分別擾動為0.23和0.27),反演結(jié)果仍然快速向真解收斂.對保山臺記錄的遠震接收函數(shù)反演結(jié)果表明:用本文方法反演所得結(jié)果與測深結(jié)果較為一致.這充分說明只要S波速度值(而非泊松比)能夠提供有效的約束,接收函數(shù)的反演過程對P波速度的選取并不敏感.
接收函數(shù),反演,初始模型
遠震P波接收函數(shù)在很大程度上消除了震源時間函數(shù)和傳播路徑的影響,主要由觀測臺站下方的地殼、上地幔速度界面所產(chǎn)生的Ps轉(zhuǎn)換波及多次反射震相組成,并且對近臺站的剪切波速度結(jié)構(gòu)非常敏感[1].因而,在過去十幾年里,用接收函數(shù)反演臺站下方的S波速度結(jié)構(gòu)是探測殼幔結(jié)構(gòu)的有效方法之一.不過,接收函數(shù)反演是一個強非線性問題,反演結(jié)果有很強的非惟一性[1].對于非線性問題,目前主要采用以下兩種處理辦法,一種是對反演方程進行線性化,將其轉(zhuǎn)變?yōu)橐粋€線性反演問題;而另一種就是直接進行非線性反演.
為了有效抑制接收函數(shù)反演問題的非惟一性,近年來發(fā)展了多種非線性反演方法.其中,一種最常見的方法是利用半空間上覆蓋一層的最簡單模型,然后基于二維網(wǎng)格進行搜索[2].這種方法由于比較簡單,被廣泛用于確定殼幔邊界.但是這種最簡化的地殼模型忽略了地殼結(jié)構(gòu)本身的復雜性,因此這種近似有很大的局限性.為了避免反演問題的解局限于目標函數(shù)的局部最小域內(nèi),很多全局優(yōu)化方法,如蒙特卡洛方法,包括屬于經(jīng)典蒙特卡洛方法范疇的模擬退火算法[3-4]、遺傳算法[5-9],以及一些新方法,如相鄰算法[10-11]等也被應用于接收函數(shù)反演問題.劉啟元等[12]根據(jù)Tarantola[13]的非線性反演理論,提出了根據(jù)接收函數(shù)徑向與垂向分量復譜比的非線性反演.然而,在采用這些非線性反演方法的過程中,對于反演結(jié)果的非惟一性,一般采用兩種辦法進行處理.一種是基于誤差相對于研究參數(shù)變化統(tǒng)計關(guān)系的一些假定,利用經(jīng)典的線性化理論進行分析計算[2];另一種是在網(wǎng)格搜索過程中引入bootstrap重新采樣[14].
對于線性化反演方法,其中以Ammon[1]的方法最具代表性,由于采用了Randall的高效算法[15]計算微分地震圖,以及Shaw和Orcutt的跳躍反演技術(shù)[16],因此反演效率極高,并且被地球物理學家廣泛采用.Ammon的研究[1]表明,接收函數(shù)反演結(jié)果在很大程度上取決于初始模型的選取.如果初始模型選取合理,則只需要有限次數(shù)的迭代即可收斂到非常接近于真實結(jié)構(gòu),因此如何獲取合適的初始模型是反演計算的關(guān)鍵.為了抑制反演的非惟一性,Julia[17]充分利用接收函數(shù)對深度變化敏感、面波頻散對S波速度值變化敏感的特點,采用了接收函數(shù)與面波頻散聯(lián)合反演[17].雖然聯(lián)合反演在一定程度上抑制了反演的非惟一性,但反演結(jié)果還是很依賴于初始模型.另外,聯(lián)合反演也受到觀測條件的限制.針對這一反演問題,吳慶舉等[18-19]基于小波變換多分辨率分析的特點,提出了一種用二進離散小波變換反演寬頻帶接收函數(shù)的方法,將寬頻帶接收函數(shù)在5階分辨尺度上展開,分解到不同頻段,從而將寬頻帶地震波形反演問題分解成不同分辨尺度的反問題,由高階到低階分別對不同分辨尺度的地震波形進行反演,并以n階尺度的反演結(jié)果作為n-1階尺度的初始模型.二進離散小波變換反演是一種廣義的線性化反演方法,對反演結(jié)果的非惟一性有一定的約束,但對第5階尺度進行反演時,還是需要提供初始模型,因此,最終的反演結(jié)果還是依賴于初始模型的選取.
在實際反演中,選取反演的初始模型主要參考已有的研究結(jié)果,為了最大限度地消除接收函數(shù)反演的不惟一性,選取充分接近真實模型的初始模型是關(guān)鍵.本文用不同的低通濾波參數(shù)對接收函數(shù)進行濾波,并獲取相應的視入射角,基于半空間的入射波理論得到介質(zhì)的S波視速度[20],從接收函數(shù)中提取模型信息.通過數(shù)值實驗總結(jié)出低通濾波參數(shù)與影響深度的關(guān)系,利用Ps震相的延時與深度的關(guān)系[2],由淺到深剝離出各層的S波速度,并以此作為反演的初始模型.
接收函數(shù)反演就是將波形中包含的信息,轉(zhuǎn)化為參數(shù)化的地殼模型結(jié)構(gòu),其對應的正演問題可表達為[1]
式中,d是接收函數(shù)波形構(gòu)成的矢量,F(xiàn)代表作用于模型m以生成接收函數(shù)d的泛函,m代表參數(shù)化后的模型空間.(1)式所表達的泛函關(guān)系是非線性的,但在初始模型m0非常靠近真實模型m的情況下,可將觀測接收函數(shù)在初始模型m0附近進行Taylor展開,將問題線性化[1]
這里,δm是模型修正矢量,D是泛函F在m0處的偏導數(shù)矩陣,即Jacobi矩陣.線性化反演就是利用(2)式逐步迭代,不斷修正參數(shù)化的模型矢量,最終逼近真解.因此,如何獲取一個充分接近真實模型的反演初始模型是線性化反演問題的關(guān)鍵.
一束平面P波從均勻半空間入射到自由表面,真入射角iP、視入射角iP、水平慢度(射線參數(shù))p、半空間的S波速度VS以及P波速度VP之間有如下關(guān)系[20]:
變換(3)式可以得到
這里,視入射角iP可表達為[20]
其中,RRF代表徑向接收函數(shù),ZRF代表垂直向接收函數(shù).利用接收函數(shù)計算視入射角要比直接利用原始記錄來得方便,免去了直接用原始記錄讀取震相的麻煩.給接收函數(shù)加上低通濾波平滑因子后,可將視入射角和半空間S波速度表達為低通濾波參數(shù)T的函數(shù),L.Svenningsen和B.H.Jacobsen[20]將這樣得到的半空間S波速度稱為S波視速度VS,app.
為了檢驗上述方法,我們以一層地殼覆蓋在半空間之上的模型為例,圖1a中給出的是S波速度隨深度的變化規(guī)律,地殼的S波速度為3.50km/s,上地幔的S波速度為4.50km/s,若泊松比σ取0.25,可利用經(jīng)驗關(guān)系:
給出密度ρ,以此確定地殼模型.我們利用Kennett提出的傳播矩陣法[21]來計算地震響應地震圖,分別對徑向分量和垂直分量進行反褶積得到相應的接收函數(shù),利用方程(6)對相應的接收函數(shù)進行平滑,得到如圖1b所示的視速度隨濾波參數(shù)的變化曲線.從圖1中可看出,濾波參數(shù)很小時,視速度值與模型值非常接近,隨著濾波參數(shù)的增大視速度也呈增加的趨勢.視速度VS,app關(guān)于濾波參數(shù)T的一階導數(shù)在一定程度上體現(xiàn)了S波速度隨深度的變化率,在導數(shù)最大的地方就是S波速度相對于深度變化最劇烈的地方,這意味著該影響深度可能就是界面的深度.通過計算,我們發(fā)現(xiàn)在4s附近視速度隨濾波參數(shù)的變化率最大,而理論模型的地殼厚度是20km,若以濾波參數(shù)乘以倍數(shù)因子5作為其對應的影響深度,則視速度的最大變化率發(fā)生在20km附近(如圖1),這與理論模型具有一定的相關(guān)性.由(6)式得到的VS,app,在濾波參數(shù)小于Ps震相走時的情況下,只有直達P波對徑向接收函數(shù)RRF(t=0)有貢獻,此時得到的VS,app就是地表附近的絕對S波速度;而隨著濾波參數(shù)的增大,影響深度越來越大,得到的VS,app是等效半空間所包含的各層的平均效果[20].因此,在第一個分層點得到的S波速度可認為是第一層的絕對S波速度,而在第二個分層點得到的是將一、二兩層等效為一層時的結(jié)果,對于多層的情況,依次類推.考慮到Ps震相的延時與層厚和S波速度有如下關(guān)系[2]:
一般情況下,波速比κ可取值1.732,我們可以先計算出第一層的,以及一、二兩層等效為一層的,由(8)式可知即為第二層的.在已知厚度和波速比的情況下,由(8)式即可反推第二層的S波速度.依次類推,即可由淺入深剝離出各層的絕對S波速度.
另外,隨著濾波參數(shù)T的增加,影響深度不斷增大,越來越多的震相對徑向接收函數(shù)RRF(t=0)產(chǎn)生影響,淺層結(jié)構(gòu)的貢獻越來越小,因此可取濾波參數(shù)較大時VS,app隨濾波參數(shù)變化的漸近趨勢作為Moho面以下的S波速度,由此得到速度隨深度變化的模型如圖1c所示.對比圖1a與1c,從接收函數(shù)中所提取的地殼S波速為3.55km/s,上地幔S波速度為4.40km/s,這一結(jié)果與理論模型非常接近,然而,厚度卻有一定的差異.
一般而言,在地球內(nèi)部速度隨深度增加而增加.如圖2所示,若地殼由厚度分別為10km、10km、20km,S波速度分別為3.00km/s、3.50km/s、3.80km/s的三層組成,波速比取為1.732,按上述的經(jīng)驗公式確定各層的密度.利用該模型合成的理論地震圖經(jīng)反褶積后得到相應徑向接收函數(shù)如圖2所示.反演的任務就是通過擬合圖2中的接收函數(shù),以提取相應的地殼S波速度結(jié)構(gòu).由于地殼結(jié)構(gòu)比較簡單,各界面上產(chǎn)生的Ps震相很容易分辨,Moho面的Ps震相的延時為5s,經(jīng)Ps相與直達P的延時估算地殼厚度為40km.若沒有其他約束信息,可考慮到地殼S波平均速度在3.5km/s左右,我們給出一個由20層厚度為2km、S波速度為3.5km/s的模型作為反演的初始模型,利用Ammon[1]的線性化反演程序進行反演計算,整個迭代過程快速地向真解逼近,第5次反演結(jié)果在Moho面以上已基本與真解重合,接收函數(shù)波形也達到99%的擬合度(如圖3).然而,若將地殼的平均速度取為3.6km/s,再次進行反演時,反演結(jié)果卻與真解出現(xiàn)了較大的偏差.如圖4所示,雖然接收函數(shù)波形達到了99%的擬合度,但第2、3層的速度值有明顯差異,并且Moho面的深度也不準確.兩次反演使用的初始模型Moho面深度都與真實模型一致,惟一的區(qū)別就在于S波速度值,取3.5km/s則剛好與真實模型的第二層一致,而取3.6km/s則與真實模型每一層的速度值都有差異.這個試驗充分說明接收函數(shù)反演嚴重依賴于初始模型的速度值.
圖1 (a)半空間上覆蓋一層的理論模型;(b)S波視速度與濾波參數(shù)的關(guān)系曲線;(c)從S波視速度與濾波參數(shù)關(guān)系曲線獲取的S波速度結(jié)構(gòu)Fig.1 (a)The theoretical model of the single layer over a half-space;(b)The relation curve between the apparent velocity of S wave versus to filtering parameter;(c)S-velocity structure acquired from the curve illustrated in(b)
圖2 (a)由(b)所示地殼模型合成的徑向接收函數(shù)(高斯參數(shù)α=2.0);(b)半空間上覆蓋三個水平層的地殼模型Fig.2 (a)Synthetic radial receiver function calculated from the model(b)(Gaussian parameterα=2.0);(b)The crustal model of three-layer over a half-space
圖3 (a)最終反演模型(第5次迭代)的合成波形與觀測波形的擬合情況(高斯參數(shù)α=2.0);(b)反演所得到的S波速度結(jié)構(gòu).黑色粗線表示圖2b所示的真實模型、黑色虛線表示由經(jīng)驗方法獲取的初始模型(泊松比σ=0.25,地殼平均S波速度為3.5km/s)、黑色細線表示1~4次反演迭代的結(jié)果,灰色粗線表示最終反演結(jié)果Fig.3 (a)The observed waveform corresponding to the model illustrated in Fig.2band the synthetic waveform calculated from the inversion result(Gaussian parameterα=2.0);(b)The S-wave velocity structure from the inversion iteration.Black thick line represents true model illustrated in Fig.2b,black dash line represents the initial model(Poisson′s ratioσ=0.25,the average crustal velocity of S wave is 3.5km/s)obtained from the empirical method,black thin lines represent results of the 1~4th inversion iterations,gray thick line represents the final result
圖4 同圖3,但初始模型地殼S波速度變?yōu)?.6km/sFig.4 Same as Fig.3,but average S-wave velocity of the crust is perturbed to 3.6km/s
圖5 (a)真實模型(如圖2b所示);(b)S波視速度與濾波參數(shù)的關(guān)系曲線;(c)從S波視速度與濾波參數(shù)關(guān)系曲線(b)中獲取的S波速度結(jié)構(gòu)Fig.5 (a)The true model illustrated in Fig.2b;(b)The relation curve between the apparent velocity of S wave versus to filtering parameter;(c)S-velocity structure acquired from the curve illustrated in(b)
圖6 同圖3,但采用如圖5c所示初始模型Fig.6 Same as Fig.3,but the initial model is replaced by the one illustrated in Fig.5c
對于同一個例子,作為對比分析,我們利用等效半空間理論得到S波視速度VS,app和濾波參數(shù)T的關(guān)系如圖5b所示.從圖5中可看出,濾波參數(shù)T<2s時的視速度值與真實模型的第一層速度值很接近,濾波參數(shù)T>19s時的視速度值與真實模型的半空間的速度值很接近.利用分層點選取原則建立的初始模型與真實模型值已很接近,仍有三層結(jié)構(gòu),但是,初始模型不僅層厚與真實模型有一定的差異,而且Moho面的深度僅為32km,這與真實模型存在較大的偏差.以圖5c所示的S波速度結(jié)構(gòu)為基礎(chǔ),將每一層分為厚度為2km的薄層,利用前述方法確定P波速度和密度以建立反演的初始模型.如圖6所示,經(jīng)過五次反演迭代計算,反演過程快速地向真解逼近,最終反演結(jié)果(第5次迭代)Moho面以上各層的S波速度值已基本與真實模型重合,僅在Moho面以下稍有差異.
上面的例子正如Ammon[1]所論述的那樣,接收函數(shù)反演結(jié)果嚴重依賴于初始模型的選取.為了對比分析,現(xiàn)在選取Ammon[1]所采用的包含低速層的復雜地殼模型進行數(shù)值實驗.如圖7所示,S波速度隨著深度增加而增加,但深度在30km附近存在一厚度為4km的低速層,從模型中計算出的理論接收函數(shù)也表明,Ps相較為復雜,存在多個Ps相,其中,來自Moho面的Ps出現(xiàn)在5.5s處.下面就嘗試用不同的初始模型來擬合該理論接收函數(shù).
由于本例中的地殼厚度為44km,與上一節(jié)中的很接近,我們同樣給出一個由22層厚度為2km、S波速度為3.5km/s的地殼模型作為反演的初始模型,進行迭代.雖然初始模型的Moho面深度與真實模型一致,如圖8所示,經(jīng)過五次反演迭代后,波形擬合度達到了98%,反演結(jié)果中也捕捉到了低速層,但低速層以及Moho面的深度是不準確的.說明反演過程中,接收函數(shù)對深度變化敏感的特點都沒有得到發(fā)揮,另外,各層的S波速度值和真實模型相差較大.對于這樣一個復雜的模型,若沒有可靠的速度值作為約束條件,即使Moho面深度精確給出,也不能有效抑制反演的不惟一性.
圖7 (a)由(b)所示地殼模型合成的徑向接收函數(shù)(高斯參數(shù)α=2.0);(b)Ammon[1]采用的含低速層的地殼模型Fig.7 (a)Synthetic radial receiver function calculated from the model(b)(Gaussian parameterα=2.0);(b)The crustal model with a low velocity used by Ammon[1]
圖8 同圖3,但采用如圖7所示地殼結(jié)構(gòu)和接收函數(shù)作為真實模型和觀測波形Fig.8 Same as Fig.3,but the true model and observed waveform are replaced by the ones illustrated in Fig.7
對于同一個模型,我們從低通濾波后的接收函數(shù)中獲取S波視速度VS,app和濾波參數(shù)T的關(guān)系曲線如圖9所示,并利用前述分層點選取原則得到相應的初始模型.初始模型(圖9c)由4層組成,層厚與真實模型有較大差異,而且不含低速層.不過,在淺部,初始模型的速度值與真實模型很接近.為了更好地擬合接收函數(shù),我們將初始模型的每一層均拆分為厚度為2km的薄層,經(jīng)過五次反演迭代,反演過程快速地逼近“真實”模型(圖10),并準確地捕捉到了低速層.
為了檢驗反演過程對模型參數(shù)的敏感性,我們在理論接收函數(shù)濾波頻段上增加了10%的隨機噪聲,并且對初始模型的波速比分別進行5%的正負擾動,即泊松比分別取為0.23和0.27,再次進行反演計算.反演迭代過程如圖11和圖12所示,經(jīng)過5次迭代,雖然波形擬合度不高,但反演結(jié)果仍然快速地向真實模型逼近.另外,低速層和Moho面的捕捉仍然非常準確,殼內(nèi)各層的S波速度值也與真實解吻合得非常好.
圖9 (a)“真實”模型(如圖7b所示);(b)S波視速度與濾波參數(shù)的關(guān)系曲線;(c)從S波視速度與濾波參數(shù)關(guān)系曲線(b)中獲取的S波速度結(jié)構(gòu)Fig.9 (a)The“ture”model illustrated in Fig.7b;(b)The relation curve between the apparent velocity of S wave versus to filtering parameter;(c)S-velocity structure acquired from the curve illustrated in(b)
圖10 (a)最終反演模型(第5次迭代)的合成波形與觀測波形的擬合情況(高斯參數(shù)α=2.0);(b)反演所得到的S波速度結(jié)構(gòu).黑色粗線表示圖7b所示的真實模型、黑色虛線表示用本文方法獲取的初始模型(即圖9c,泊松比σ=0.25)、黑色細線表示1~4次反演迭代的結(jié)果,灰色粗線表示最終反演結(jié)果Fig.10 (a)The observed waveform corresponding to the model illustrated in Fig.2band the synthetic waveform calculated from the inversion result(Gaussian parameterα=2.0);(b)The S-wave velocity structure from the inversion iteration.Black thick line represents true model illustrated in Fig.2b,black dash line represents the initial model obtained from the study(showed in Fig.9c,Poisson′s ratioσ=0.25),black thin lines represent results of the 1~4th inversion iterations,gray thick line represents the final result
云南地區(qū)有長達幾百公里的測深剖面,為檢驗反演結(jié)果提供了條件.雖然數(shù)值實驗取得了較好的反演結(jié)果,為了檢驗實際觀測數(shù)據(jù)的反演效果,我們從云南保山臺記錄的遠震事件提取了接收函數(shù)波形(如圖13所示),事件震中距為6520km.由于震相比較清晰,可以看到Ps震相延時約為5s左右,由此可估計Moho面深度約為40km.我們先采用本文的方法,用不同的低通濾波參數(shù)對接收函數(shù)進行濾波,得到S波視速度與濾波參數(shù)的關(guān)系曲線,再利用剝層技術(shù)獲取用于反演的初始模型.最后,反演得到了保山臺下方的S波速度結(jié)構(gòu)(如圖13c).為了與前人的結(jié)果作對比,我們收集了胡家富等[22]的接收函數(shù)反演結(jié)果(圖13c細實線),以及人工測深的結(jié)果[23].由于人工測深給出的是P波速度,這里利用保山臺的泊松比[22]換算為S波速度(圖13c粗虛線和細虛線).與已有的接收函數(shù)反演結(jié)果相比,雖然Moho面的深度較為一致,但本文所得結(jié)果更接近人工測深結(jié)果,尤其是在上地幔里更接近測深結(jié)果.另外,胡家富等的反演結(jié)果在22km附近有一個明顯的低速層,而本文的反演結(jié)果在這個位置沒有明顯的低速結(jié)構(gòu),這與爆破地震研究結(jié)果較為一致.
圖11 同圖10,但在觀測波形的濾波頻段上加入了10%的隨機噪聲,初始模型泊松比取0.23Fig.11 Same as Fig.10,but the observed waveform is mixed with 10%random noise,and the Poisson's ratio of the initial model is perturbed to 0.23
圖12 同圖10,但在觀測波形的濾波頻段上加入了10%的隨機噪,初始模型泊松比取0.27Fig.12 Same as Fig.10,but the observed waveform is mixed with 10%random noise,and the Poisson's ratio of the initial model is perturbed to 0.27
Ammon[1]的線性化反演算法在體現(xiàn)高效的同時,由于引入了平滑約束,在初始模型的選取上給予了一定程度的靈活性.不同的初始模型只要提供了足夠相容的約束信息,都會向同一個解收斂.因此當初始模型發(fā)生擾動時,不會使最終結(jié)果過于發(fā)散,從而保障了反演結(jié)果的穩(wěn)定性和可靠性.數(shù)值實驗表明,當?shù)貧そY(jié)構(gòu)簡單時,通過Ps震相的延時判斷出Moho面深度,并取地殼平均S波速度為3.5km/s作為初始模型,也可能使反演過程向真解收斂,但前提是地殼平均S波速度值能為反演提供足夠準確的約束.但在地殼結(jié)構(gòu)比較復雜的情況下,選取一個簡單的初始模型將導致反演嚴重偏離真實解.其原因主要是各層真實速度值偏離地殼平均速度值較大所致,此時僅依靠Moho面深度和地殼平均S波速度的約束肯定是不夠的,還必須提供更多的分層信息.
圖13 (a)保山臺記錄的遠震事件三分量波形(震中距6520km);(b)觀測徑向接收函數(shù)(實線)和反演模型的合成波形(虛線)的擬合情況(高斯參數(shù)α=2.0);(c)細實線表示參考模型[22],粗虛線為爆破地震保山(南)炮點的速度結(jié)構(gòu)[23],細虛線為爆破地震保山(北)炮點的速度結(jié)構(gòu)[23],粗實線為本文的反演結(jié)果Fig.13 (a)Teleseismic 3-component waveforms recorded by Station BS(dist=6520km);(b)Radial receiver function(solid line)calculated from(a)and synthetic radial receiver function(dash line)calculated from the final inversion result(α=2.0);(c)Thin solid line represents the reference model[22],thick dash line represents S-velocity structure of the south explosion point of Baoshan[23],thin dash line represents S-velocity structure of the north explosion point of Baoshan[23],thick solid line represents the inversion result calculated from the initial model acquired with our method.
由于接收函數(shù)具有對深度變化敏感、而對速度值變化不敏感的特點,如果初始模型能在幾個地層的S波速度值上提供準確的約束,即使深度不準確,反演結(jié)果仍然可以向真解收斂.對于包含低速層的復雜地殼模型,僅僅依靠Moho面深度和地殼平均S波速度值的約束是不夠的,Moho面以上需要區(qū)分出更多的層次才能提供足夠的約束.本文根據(jù)不同濾波參數(shù)下,從低通濾波后的接收函數(shù)獲取的視入射角不同,影響深度也不同,利用經(jīng)驗關(guān)系從視速度隨濾波參數(shù)的變化曲線中提取了S波速度結(jié)構(gòu),與真實模型相比,速度值整體變化趨勢非常接近,但每一層的深度可能存在較大差異.這樣選取有利于充分發(fā)揮接收函數(shù)對深度變化敏感的優(yōu)勢,由于速度值提供了有效的約束,雖然Moho面深度并不準確,但反演迭代過程還是快速地向真解逼近.另外,通過給真實波形加入10%的噪聲,在保待S波速度值不變的情況下,分別對波速比進行5%的正負擾動(即泊松比分別擾動為0.23和0.27),反演結(jié)果仍然是快速向真解收斂.最后我們用本文方法對保山臺實測遠震事件數(shù)據(jù)進行實驗,并與前人的接收函數(shù)反演結(jié)果和測深結(jié)果進行對比,可以發(fā)現(xiàn)用本文方法反演得到的結(jié)果與測深結(jié)果更加吻合.這充分說明只要S波速度值能夠提供準確的約束,接收函數(shù)的反演過程對P波速度的選取并不敏感.
接收函數(shù)反演具有很強的非惟一性,正演合成波形的擬合度并不足以說明解的好壞,反演結(jié)果是否接近于真實情況在很大程度上取決于初始模型的選取,這是由于線性化過程的特點所決定的.當前很多學者一般都傾向于采用其他研究者的結(jié)果作為反演計算的初始模型,但這樣引用的初始模型是否就一定會向真解收斂是很難論證的,甚至會陷入一個死循環(huán),從而導致結(jié)果更加偏離真實解.劉啟元等[24]發(fā)展了基于貝葉斯反演理論的接收函數(shù)與環(huán)境噪聲的非線性聯(lián)合反演方法,在反演過程中利用環(huán)境噪聲的相速度頻散數(shù)據(jù)對地殼S波絕對速度提供了有力的約束,使得反演結(jié)果對初始模型的依賴程度大大降低.然而,該方法在實施過程中需要密集的地震臺網(wǎng)覆蓋,對觀測條件的要求較高.本文嘗試從接收函數(shù)所包含的信息中提取初始模型,使初始模型就局限在最靠近真解的解空間內(nèi),通過S波速度值提供有效的約束,有效地抑制了反演的非惟一性.
(References)
[1] Ammon C J,Randall G E,Zandt G.On the nonuniqueness of receiver function inversions.J.Geophys.Res.,1990,95(B10):15303-15318.
[2] Zhu L P,Kanamori H.Moho depth variation in southern California from teleseismic receiver functions.J.Geophys.Res.,2000,105(B2):2969-2980.
[3] Sen M K,Stoffa P L.Nonlinear one-dimensional seismic waveform inversion using simulated annealing.Geophysics,1991,56:1624-1638.
[4] 劉鵬程,紀晨,Hartzell S H.改進的模擬退火-單純形綜合反演方法.地球物理學報,1995,38(2):199-205.Liu P C,Ji C,Hartzell S H.An improved simulated annealing-downhill simplex hybrid global inverse algorithm.Chinse J.Geophys.(in Chinese),1995,38(2):199-205.
[5] Sambridge M,Drijkoningen G.Genetic algorithms in seismic waveform inversion.Geophys.J.Int.,1992,109(2):323-342.
[6] Sen M K,Stoffa P L.Rapid sampling of model space using genetic algorithms:examples from seismic waveform inversion.Geophys.J.Int.,1992,108(1):281-292.
[7] 石耀霖,金文.面波頻散反演地球內(nèi)部構(gòu)造的遺傳算法.地球物理學報,1995,38(2):189-198.Shi Y L,Jin W.Genetic algorithms inversion of lithospheric structure from surface wave dispersion.Chinese J.Geophys.(in Chinse),1995,38(2):189-198.
[8] 吳建平,明躍紅,湯毅.遺傳算法在上地幔速度結(jié)構(gòu)研究中的應用.地震地磁觀測與研究,1997,18(2):11-29.Wu J P,Ming Y H,Tang Y.Application of genetic algorithm in study of upper mantle velocity structure.Seismological and Geomagnetic Observation and Research(in Chinese),1997,18(2):11-29.
[9] 吳建平,明躍紅,曾融生.遺傳算法中的光滑約束反演及其在青藏高原面波研究中的應用.地震學報,2001,23(1):45-53.Wu J P,Ming Y H,Zeng R S.Smooth constraint inversion technique in genetic algorithms and its application to surface wave study in the Tibetan plateau.Acta Seismologica Sinica(in Chinese),2001,23(1):45-53.
[10] Sambridge M.Geophysical inversion with a neighbourhood algorithm—Ⅰ.Searching aparameter space.Geophys.J.Int.,1999,138(2):479-494.
[11] Sambridge M.Geophysical inversion with a neighbourhood algorithm-II.Appraising the ensemble.Geophys.J.Int.,1999,138:727-746.
[12] 劉啟元,Kind R,李順成.接收函數(shù)復譜比的最大或然性估計及非線性反演.地球物理學報,1996,39(4):502-513.Liu Q Y,Kind R,Li S C.Maximal likelihood estimation and nonlinear inversion of the complex receiver function spectrum ratio.Chinese J.Geophys.(in Chinese),1996,39(4):502-513.
[13] Tarantola A.Inverse Problem Theory,Method for Data Fitting and Model Parameter Estimation.Amsterdam:Elsevier,1987.
[14] Sandvol E,Seber D,Calvert A,et al.Grid search modeling of receiver functions:implications for crustal structure in the Middle East and North Africa.J.Geophys.Res.,1998,103(B11):26899-26917.
[15] Randall G E.Efficient calculation of differential seismograms for lithospheric receiver functions.Geophys.J.Int.,1989,99(3):469-481.
[16] Shaw P R,Orcutt J A.Waveform inversion of seismic refraction data and applications to young Pacific crust.Geophys.J.R.Astron.Soc.,1985,82(3):375-414.
[17] JuliàJ,Ammon C J,Herrmann R B,et al.Joint inversion of receiver function and surface wave dispersion observations.Geophys.J.Int.,2000,143(1):99-112.
[18] Wu Q J,Li Y H,Zhang R Q,et al.Wavelet modeling of broad-band receiver functions.Geophys.J.Int.,2007,170(2):534-544.
[19] 吳慶舉,田小波,張乃鈴等.用小波變換方法反演接收函數(shù).地震學報,2003,25(6):601-607.Wu Q J,Tian X B,Zhang N L,et al.Inversion of receiver function by wavelet transformation.Acta Seismologica Sinica(in Chinese),2003,25(6):601-607.
[20] Svenningsen L,Jacobsen B H.Absolute S-velocity estimation from receiver functions.Geophys.J.Int.,2007,170(3):1089-1094.
[21] Kennett B L N.Seismic Wave Propagation in Stratified Media.Cambridge:Cambridge University Press,1983.
[22] 胡家富,蘇有錦,朱雄關(guān)等.云南的地殼S波速度與泊松比結(jié)構(gòu)及其意義.中國科學(D輯),2003,33(8):714-722.Hu J F,Su Y J,Zhu X G,et al.S-wave velocity and Poisson′s ratio structure of crust in Yunnan and its implication.Science in China(Series D)(in Chinese),2002,33(8):714-722.
[23] 胡鴻翔,陸涵行,王椿鏞等.滇西地區(qū)地殼結(jié)構(gòu)的爆破地震研究.地球物理學報,1986,29(2):133-144.Hu H X,Lu H X,Wang C Y,et al.Explosion investigation of the crustal structure in western Yunnan Province.Chinese J.Geophys.(in Chinese),1986,29(2):133-144.
[24] 劉啟元,李昱,陳九輝等.基于貝葉斯理論的接收函數(shù)與環(huán)境噪聲聯(lián)合反演.地球物理學報,2010,53(11):2603-2612.Liu Q Y,Li Y,Chen J H,et al.Joint inversion of receiver function and ambient noise based on Bayesian theory.Chinese J.Geophys.(in Chinese),2010,53(11):2603-2612.
An effective technique to constrain the non-uniqueness of receiver function inversion
PENG Heng-Chu,HU Jia-Fu*,YANG Hai-Yan,WEN Li-Min
Department of Geophysics,Yunnan University,Kunming650091,China
The relation curve between the apparent velocity of S wave and the filter-parameter is obtained by low-pass filtering the radial and vertical receiver function,then it is transformed into the S velocity structure beneath the station with empirical formula;moreover,the S velocity structure will be regarded as initial model for receiver function inversion.The numerical simulation indicated that the inversion processes would approach quickly to the true model once the initial S velocity can supply the effective constraint to inversion,even if the depth of Moho is not accurate.In addition,for the same initial S structure,the inversion results also approach to the true model quickly,even adding 10%of random noises to the“observed”receiver function,and 5%of positive and negative perturbation to the velocity ratio(corresponding to the Poisson′s ratio of 0.27and 0.23)respectively.We processed the teleseismic receiver function recorded by Baoshan station.The result from our method is similar to that of sounding section.These confirmed that the inversion processes of receiver function are not sensitive to the initial P velocity if the initial S wave velocity can provide the effective constraint information.
Receiver function,Inversion,Initial model
10.6038/j.issn.0001-5733.2012.04.013
P315
2011-08-08,2012-03-09收修定稿
國家自然科學基金(U0933602)和云南省基金資助項目(2007D013M)資助.
彭恒初,男,講師,2004年畢業(yè)于北京大學,獲碩士學位,研究方向為深部地球物理.E-mail:terrypku@163.com
*通訊作者胡家富,E-mail:jfhu@ynu.edu.cn
彭恒初,胡家富,楊海燕等.接收函數(shù)反演地殼S波速度結(jié)構(gòu)的有效約束方法.地球物理學報,2012,55(4):1168-1178,
10.6038/j.issn.0001-5733.2012.04.013.
Peng H C,Hu J F,Yang H Y,et al.An effective technique to constrain the non-uniqueness of receiver function inversion.Chinese J.Geophys.(in Chinese),2012,55(4):1168-1178,doi:10.6038/j.issn.0001-5733.2012.04.013.
(本文編輯 胡素芳)