本發(fā)明屬于地球動(dòng)力學(xué)研究技術(shù)領(lǐng)域，具體涉及一種接收函數(shù)和重力反演地殼厚度和波速比的方法。

背景技術(shù)：

地殼厚度和波速比等參數(shù)是研究地殼結(jié)構(gòu)和物質(zhì)組成的重要參數(shù),對(duì)研究地球動(dòng)力學(xué)過程具有重要意義。20世紀(jì)80年代以來，隨著寬頻帶地震臺(tái)網(wǎng)觀測(cè)技術(shù)的不斷完善和發(fā)展，遠(yuǎn)震體波接收函數(shù)方法的發(fā)展與應(yīng)用使地殼上地幔速度結(jié)構(gòu)的認(rèn)識(shí)得到了空前的提升。

接收函數(shù)是由遠(yuǎn)震徑向分量和垂直分量反褶積得到的時(shí)間序列，由于扣除了震源效應(yīng)和傳播路徑的影響，因此更適合用來研究地殼和上地幔間斷面的結(jié)構(gòu)信息。接收函數(shù)的H-κ疊加法,用于估計(jì)地殼厚度與波速比參數(shù)，通過波速比與泊松比的關(guān)系式可換算地殼泊松比參數(shù)。該方法不用人工挑選震相走時(shí),對(duì)間斷面上的速度差敏感，簡(jiǎn)單快速，近年來已廣泛用于研究全球多個(gè)地區(qū)的地殼厚度和波速比，是目前獲取地殼厚度(H)和波速比(κ)參數(shù)的常用方法。

然而實(shí)際應(yīng)用中，有些地區(qū)往往會(huì)由于臺(tái)基不穩(wěn)定、地震數(shù)據(jù)質(zhì)量差、臺(tái)站下方構(gòu)造復(fù)雜等問題，造成接收函數(shù)的多次反射震相不清晰、拾取困難，從而導(dǎo)致無法或很難從H-κ疊加圖中拾取出合理的H和κ參數(shù)值。如圖1(a)中最大值所在是一個(gè)區(qū)域，無法確定出相對(duì)應(yīng)的H和κ，圖1(b)中置信區(qū)間有三個(gè)，同樣無法確定出合理的H和κ值。其中，臺(tái)站下方構(gòu)造復(fù)雜包括較厚沉積層、殼內(nèi)低速層和傾斜地層等地質(zhì)情況。沉積層和基底之間存在明顯的速度界面，這在接收函數(shù)中會(huì)出現(xiàn)介于初至P波和轉(zhuǎn)換Ps波之間的一個(gè)波峰，從而對(duì)H-κ疊加結(jié)果產(chǎn)生較大影響，且沉積層越厚，這種影響就越大。

目前減少沉積層影響的方法有固定值域模擬搜索技術(shù)、添加沉積層速度結(jié)構(gòu)等。殼內(nèi)低速層引起的問題在接收函數(shù)的H-κ疊加中主要表現(xiàn)為出現(xiàn)雙Ps波峰，這種情況下就會(huì)出現(xiàn)疊加結(jié)果的嚴(yán)重偏差。

重力方法是研究地質(zhì)結(jié)構(gòu)和區(qū)域地質(zhì)構(gòu)造的重要方法之一，水平分辨率高。布格重力異常通常包含了地殼內(nèi)部物質(zhì)成分密度不均勻引起的重力異常(即地殼重力異常)，也包括了莫霍面起伏引起的重力異常(即莫霍面重力異常)。地殼厚度參數(shù)可描述莫霍面起伏變化，波速比參數(shù)則可描述地殼內(nèi)部物質(zhì)成分并與地殼密度參數(shù)相關(guān)，也即它可描述地殼內(nèi)部密度分布的不均勻性。因此,布格重力異常除了與地殼密度和厚度直接相關(guān)外，還與地殼波速比密切相關(guān)，通過布格重力異常反演理論上可以獲取這些參數(shù)。但是僅僅依靠重力反演一種方法不能獲得較為可信的地殼厚度與波速比。

總體來講，在地殼結(jié)構(gòu)復(fù)雜的情況下接收函數(shù)與重力反演任何一種方法單方面都不能獲得較為可信的地殼厚度與波速比。而重力可以為接收函數(shù)H-κ疊加法提供約束信息，降低反演的不確定性或多解性。

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

為了解決由于臺(tái)基不穩(wěn)定、地震數(shù)據(jù)質(zhì)量差、臺(tái)站下方構(gòu)造復(fù)雜等問題造成接收函數(shù)的多次反射震相不清晰、拾取困難而無法拾取出合理地殼厚度H和波速比κ的技術(shù)問題。本發(fā)明提供一種接收函數(shù)和重力聯(lián)合反演獲得地殼厚度及波速比的方法；本發(fā)明通過利用布格重力異常反演莫霍面深度并轉(zhuǎn)換成地殼厚度，在接收函數(shù)H-k疊加圖中拾取地殼厚度相對(duì)應(yīng)的波速比，以區(qū)域內(nèi)所有臺(tái)站的H和κ作為初值，計(jì)算密度參數(shù)與噪聲參數(shù)，然后在H和κ的參數(shù)空間內(nèi)掃面H和κ，并用該H和κ替換窗口中心臺(tái)站的H和κ，計(jì)算相應(yīng)的重力似然圖，然后將該重力似然圖與接收函數(shù)H-κ疊加圖歸一化并點(diǎn)乘即可得到聯(lián)合反演的最優(yōu)H和κ。

本發(fā)明是通過以下技術(shù)方案實(shí)現(xiàn)的：

一種準(zhǔn)確獲取地殼厚度和波速比的方法，所述方法用于降低接收函數(shù)H-κ疊加法的不確定性，提高獲取地殼厚度H和波速比κ的精度和效率，所述方法基于接收函數(shù)和重力進(jìn)行聯(lián)合反演，所述方法首先利用布格重力異常反演獲取地殼厚度H初值，進(jìn)而獲得波速比κ初值；然后，對(duì)選擇的任意一個(gè)待反演地震臺(tái)站為中心進(jìn)行單一臺(tái)站聯(lián)合反演，得到單一臺(tái)站的聯(lián)合反演的H-κ圖，從而獲取該地震臺(tái)站的最佳估計(jì)的地殼厚度H和波速比κ值；最后，對(duì)研究區(qū)域內(nèi)所有待反演站臺(tái)進(jìn)行上述單一臺(tái)站聯(lián)合反演過程，以獲得研究區(qū)域內(nèi)所有待反演地震臺(tái)站的地殼厚度H和波速比κ；所述聯(lián)合反演的H-κ圖中地殼厚度H和波速比κ拾取簡(jiǎn)單；所述方法能夠?yàn)榈貧ぴ龊窈偷貧ぷ冃螜C(jī)制模式的研究提供高精度的地殼厚度和波速比數(shù)值。

進(jìn)一步地，所述方法包括如下步驟：

利用布格重力異常反演研究區(qū)域內(nèi)的莫霍面深度并轉(zhuǎn)換成地殼厚度H，提取待反演臺(tái)站位置所對(duì)應(yīng)的H作為聯(lián)合反演的地殼厚度H初值；

以所述地殼厚度H為引導(dǎo)，在接收函數(shù)H-κ疊加圖中拾取與所有待反演地震臺(tái)站的地殼厚度相對(duì)應(yīng)的波速比κ，作為聯(lián)合反演的波速比κ初值；

選擇任意一個(gè)待反演的地震臺(tái)站為中心，搜索一定窗口內(nèi)所有地震臺(tái)站的地殼厚度H、波速比κ和布格重力異常數(shù)據(jù)；

對(duì)上述選擇的地震臺(tái)站進(jìn)行單一臺(tái)站聯(lián)合反演，使用單一臺(tái)站聯(lián)合反演中獲得的最佳估計(jì)的H和κ值替換中心地震臺(tái)站的原H和κ值；其中，所述單一臺(tái)站聯(lián)合反演基于重力反演的H-κ似然圖與接收函數(shù)反演的H-κ疊加圖；

然后選擇下一個(gè)待聯(lián)合反演的臺(tái)站，重復(fù)所述單一臺(tái)站聯(lián)合反演的過程，直到研究區(qū)域內(nèi)所有待反演地震臺(tái)站的H和κ被替換為最佳估計(jì)的H和κ值，從而獲得研究區(qū)域內(nèi)所有待反演地震臺(tái)站的地殼厚度H和波速比κ。

進(jìn)一步地，所述單一臺(tái)站聯(lián)合反演的過程具體為：

(1)將以所述待反演地震臺(tái)站為中心的一定窗口內(nèi)的所有地震臺(tái)站的地殼厚度H、波速比κ和布格重力異常數(shù)據(jù)，分別網(wǎng)格化成相同大小的網(wǎng)格數(shù)據(jù)；

(2)用步驟(1)獲得的網(wǎng)格化后的地殼厚度H、波速比κ和布格重力異常，結(jié)合線性迭代算法計(jì)算出莫霍面上下密度差Δρ_Moho和地殼密度隨波速比的變化率

(3)利用求得的所述莫霍面上下密度差Δρ_Moho和地殼密度隨波速比的變化率計(jì)算出模型理論重力異常；

(4)將布格重力異常與所述模型理論重力異常作差，求得異常殘差，采用似然估計(jì)法分別計(jì)算出異常殘差的均值μ與方差σ²，再計(jì)算似然函數(shù)值；

(5)以接收函數(shù)H-κ疊加圖相同的范圍和步長(zhǎng)選取H和κ值，重復(fù)步驟(3)和(4)，形成重力反演的H-κ似然圖；

(6)將接收函數(shù)反演的H-κ疊加圖與重力反演的H-κ似然圖相乘，得到聯(lián)合反演的H-κ圖，從聯(lián)合反演的H-κ圖中拾取最佳估計(jì)的H和κ值。

進(jìn)一步地，利用布格重力異常反演研究區(qū)域內(nèi)的莫霍面深度并轉(zhuǎn)換成地殼厚度H，具體為：

用頻率域?yàn)V波法對(duì)研究區(qū)域內(nèi)的布格重力異常作異常分離獲得莫霍面重力異常，再應(yīng)用密度界面反演方法獲得研究區(qū)域內(nèi)的莫霍面深度分布，并加上地形換算為地殼厚度分布，提取研究區(qū)內(nèi)所有待反演臺(tái)站的位置所對(duì)應(yīng)的地殼厚度值作為聯(lián)合反演的地殼厚度H初值。

進(jìn)一步地，步驟(2)中，所述布格重力異常由莫霍面重力異常與地殼重力異常組成，所述莫霍面上下密度差Δρ_Moho與地殼密度隨波速比的變化率的計(jì)算使用以下公式：

Δg_B＝Δg_Moho+Δg_Crust (1)

其中,Δg_B、Δg_Moho及Δg_Crust分別代表布格重力異常、莫霍面重力異常及地殼重力異常；Δρ_Moho為莫霍面上下密度差，為地殼密度隨波速比的變化率，為研究區(qū)地殼厚度的平均值，為研究區(qū)波速比的平均值，F(xiàn){·}和F^-1{·}分別表示傅里葉變換和反傅里葉變換，G為萬有引力常數(shù)，f為波數(shù)，c為地殼厚度變化的校正因子；H為地殼厚度，κ為波速比。

進(jìn)一步地，步驟(3)中計(jì)算出模型理論重力異常具體為：將步驟(2)中計(jì)算獲得的密度參數(shù)Δρ_Moho與反代入步驟(2)的公式(2)、(3)中計(jì)算出Δg_Moho和Δg_Crust，計(jì)算所述模型理論重力異常Δg_Model的公式為：

Δg_Model＝Δg_Moho+Δg_Crust (4)

其中：Δg_Moho是莫霍面重力異常，Δg_Crust為地殼重力異常。

進(jìn)一步地，步驟(4)中，假設(shè)異常殘差服從高斯分布，因而采用下式計(jì)算μ與σ²：

其中：是待反演地震臺(tái)站所在點(diǎn)上的重力異常殘差，n代表計(jì)算窗口內(nèi)的網(wǎng)格點(diǎn)數(shù)。

進(jìn)一步地，步驟(5)中，在與接收函數(shù)H-κ疊加圖相同的H和κ參數(shù)范圍內(nèi)，變動(dòng)所述待反演地震臺(tái)站的H和κ，當(dāng)其似然值最大時(shí)，擬合出的模型重力異常與布格重力異常最接近，故此時(shí)的H和κ為最佳估計(jì)值，重力反演的H-κ似然圖采用如下算法：

其中：L(μ,σ²)表示最大似然值；μ與σ²分別表示異常殘差的均值和方差，是待反演地震臺(tái)站所在點(diǎn)上的重力異常殘差；n代表計(jì)算窗口內(nèi)的網(wǎng)格點(diǎn)數(shù)。

進(jìn)一步地，步驟(6)具體為：將重力反演的H-κ似然圖與接收函數(shù)H-k疊加圖均歸一化，通過兩矩陣點(diǎn)乘實(shí)現(xiàn)相互加權(quán)，得到聯(lián)合反演的H-κ圖，從聯(lián)合反演的H-κ圖中拾取最佳估計(jì)的H和κ值。

進(jìn)一步地，不斷更新并迭代研究區(qū)域內(nèi)所有待反演地震臺(tái)站的H和κ，以避免由于地殼厚度H和波速比κ初值誤差太大而引起的聯(lián)合反演結(jié)果不準(zhǔn)確的情況。

本發(fā)明的有益技術(shù)效果：

(1)重力似然反演采用滑動(dòng)窗口，每次反演只對(duì)窗口內(nèi)的小面積進(jìn)行，提高計(jì)算精度；

(2)考慮布格重力異常由莫霍面重力異常和地殼重力異常組成，簡(jiǎn)化計(jì)算算法；

(3)針對(duì)接收函數(shù)計(jì)算中多次波不明顯而無法估計(jì)H、κ的問題，采用重力異常界面反演得到莫霍面深度，進(jìn)一步換算出地殼厚度初值，再引導(dǎo)從接收函數(shù)H-κ疊加圖中拾取出波速比參數(shù)，減少計(jì)算時(shí)間；

(4)布格重力異常反演地殼厚度和波速比過程中先用似然估計(jì)估計(jì)出噪聲的均值和方差，然后在H和κ掃描過程中計(jì)算相應(yīng)的最大似然值；

(5)區(qū)域內(nèi)所有待反演地震臺(tái)站的多次迭代計(jì)算有效提高了地殼厚度與波速比的反演精度；解決了在天然地震接收函數(shù)數(shù)據(jù)量少、數(shù)據(jù)質(zhì)量差、臺(tái)站下方存在較厚沉積層或者地殼內(nèi)部存在低速層時(shí)，計(jì)算地殼厚度和波速比不準(zhǔn)確的問題，能夠?yàn)榈貧ぴ龊窈偷貧ぷ冃螜C(jī)制模式的研究提供相對(duì)可靠的地殼厚度和波速比數(shù)值，為檢驗(yàn)地殼增厚和變形機(jī)制模式提供重要依據(jù)。

(6)本發(fā)明通過重力似然圖與接收函數(shù)疊加圖相互加權(quán)實(shí)現(xiàn)最優(yōu)H和κ的拾取，為解決復(fù)雜情況下接收函數(shù)反演的多解性問題提供有效途徑。

附圖說明

圖1：a圖為接收函數(shù)數(shù)據(jù)少或者質(zhì)量較差時(shí)的接收函數(shù)H-κ疊加圖；b圖為存在殼內(nèi)低速層時(shí)的接收函數(shù)H-κ疊加圖；

圖2為重力反演的H-κ似然圖；

圖3分別為簡(jiǎn)單模型和復(fù)雜模型的接收函數(shù)與重力聯(lián)合反演圖；

圖4為單一臺(tái)站處理流程圖；

圖5為研究區(qū)域內(nèi)所有待反演地震臺(tái)站處理流程圖。

具體實(shí)施方式

為了使本發(fā)明的目的、技術(shù)方案及優(yōu)點(diǎn)更加清楚明白，以下結(jié)合附圖及實(shí)施例，對(duì)本發(fā)明進(jìn)行進(jìn)一步詳細(xì)描述。應(yīng)當(dāng)理解，此處所描述的具體實(shí)施例僅僅用于解釋本發(fā)明，并不用于限定本發(fā)明。

相反，本發(fā)明涵蓋任何由權(quán)利要求定義的在本發(fā)明的精髓和范圍上做的替代、修改、等效方法以及方案。進(jìn)一步，為了使公眾對(duì)本發(fā)明有更好的了解，在下文對(duì)本發(fā)明的細(xì)節(jié)描述中，詳盡描述了一些特定的細(xì)節(jié)部分。對(duì)本領(lǐng)域技術(shù)人員來說沒有這些細(xì)節(jié)部分的描述也可以完全理解本發(fā)明。

實(shí)施例1

地殼厚度和波速比是研究地殼結(jié)構(gòu)和物質(zhì)組成的重要參數(shù)，對(duì)研究地球動(dòng)力學(xué)機(jī)制有著重要的意義。在檢驗(yàn)地殼增厚和地殼變形機(jī)制模式的研究中，需要以地殼厚度和波速比作為研究基礎(chǔ)。以地殼厚度和波速比在青藏高原東北緣的研究中的應(yīng)用為例，青藏高原東北緣是青藏高原向東北方向生長(zhǎng)的最新的和正在變形的重要組成部分，廣泛發(fā)育的第四系褶皺和逆沖斷裂等表明該區(qū)正遭受著地殼縮短，并且伴隨著垂直差異性隆升作用，地殼明顯增厚。國(guó)內(nèi)外學(xué)者通過研究利用地殼厚度和波速比對(duì)青藏高原東北緣地殼增厚和變形機(jī)制提出了很多模式，具有代表性的有：垂直連續(xù)變形說、上地殼增厚說和下地殼流說。其中，垂直連續(xù)變形說對(duì)應(yīng)的耦合地殼增厚模型假設(shè)上地殼和下地殼的增厚速率一致，波速比和密度不會(huì)隨地殼的增厚有明顯變化。上地殼增厚說對(duì)應(yīng)的上地殼增厚模型假設(shè)上地殼增厚速率明顯大于下地殼，平均波速比和密度隨地殼增厚而降低。下地殼流說對(duì)應(yīng)的下地殼流模型假設(shè)地殼增厚主要是因?yàn)橄噜彴鍓K侵入到下地殼而引起，地殼平均波速比和密度隨地殼增厚而增加。

由以上青藏高原東北緣的例子可知，地殼厚度和波速比是地殼增厚和地殼變形機(jī)制模式的研究的基礎(chǔ)。并且，地殼厚度和波速比在地殼增厚和地殼變形機(jī)制模式的研究中的應(yīng)用并不以地區(qū)為限制，可以被廣泛應(yīng)用于不同地區(qū)的研究。

另外，地殼厚度和波速比還能夠用于研究地殼組分，例如，有學(xué)者對(duì)華北克拉通中西部地區(qū)的地殼厚度與波速比進(jìn)行研究，結(jié)果顯示該地區(qū)山區(qū)平均波速比主要集中分布于1.70～1.77之間,暗示山區(qū)塊體較東部盆地地區(qū)地殼組成更富長(zhǎng)英質(zhì),而缺少鐵鎂質(zhì)成分。

而目前使用最廣泛的天然地震接收函數(shù)H-κ疊加法在實(shí)際應(yīng)用中，往往會(huì)由于臺(tái)基不穩(wěn)定、地震數(shù)據(jù)質(zhì)量差、臺(tái)站下方構(gòu)造復(fù)雜等問題，造成接收函數(shù)的多次反射震相不清晰、拾取困難，從而導(dǎo)致無法或很難從H-κ疊加圖中拾取出合理的地殼厚度H和波速比κ參數(shù)值。錯(cuò)誤或不準(zhǔn)確的地殼厚度與波速比對(duì)地殼增厚和地殼變形機(jī)制模式帶來了很大的阻礙，會(huì)造成不能夠獲得正確的地殼增厚和地殼變形機(jī)制模式。

綜上所述，在地殼增厚和地殼變形機(jī)制模式的研究中，為了獲得正確的地殼增厚和地殼變形機(jī)制模式，亟需一種方法能夠準(zhǔn)確測(cè)定地殼厚度和波速比。

本實(shí)施例提供一種一種接收函數(shù)和重力反演地殼厚度和波速比的方法，所述方法能夠更加準(zhǔn)確地獲取地殼厚度和波速比，能夠確保為地殼增厚和地殼變形機(jī)制模式的研究提供高精度的地殼厚度和波速比數(shù)值，以獲得正確的地殼增厚和地殼變形機(jī)制模式。所述方法包括如下步驟：

利用布格重力異常反演研究區(qū)域內(nèi)的莫霍面深度并轉(zhuǎn)換成地殼厚度H，提取待反演臺(tái)站位置所對(duì)應(yīng)的H作為聯(lián)合反演的地殼厚度H初值；具體為：用頻率域?yàn)V波法對(duì)研究區(qū)域內(nèi)的布格重力異常作異常分離獲得莫霍面重力異常，再應(yīng)用密度界面反演方法獲得研究區(qū)域內(nèi)的莫霍面深度分布，提取待反演臺(tái)站位置的莫霍面深度并加上地形換算為地殼厚度分布，作為聯(lián)合反演的地殼厚度H初值。

以所述地殼厚度H為引導(dǎo)，在接收函數(shù)H-κ疊加圖中拾取與所有待反演地震臺(tái)站的地殼厚度相對(duì)應(yīng)的波速比κ，作為聯(lián)合反演的波速比κ初值；

選擇任意一個(gè)待反演的地震臺(tái)站為中心，搜索一定窗口內(nèi)所有地震臺(tái)站的地殼厚度H、波速比κ和布格重力異常數(shù)據(jù)；

對(duì)上述選擇的地震臺(tái)站進(jìn)行單一臺(tái)站聯(lián)合反演，使用單一臺(tái)站聯(lián)合反演中獲得的最佳估計(jì)的H和κ值替換中心地震臺(tái)站的原H和κ值；其中，所述單一臺(tái)站聯(lián)合反演基于重力反演的H-κ似然圖與接收函數(shù)反演的H-κ疊加圖；

然后選擇下一個(gè)待聯(lián)合反演的臺(tái)站，重復(fù)所述單一臺(tái)站聯(lián)合反演的過程，直到研究區(qū)域內(nèi)所有待反演地震臺(tái)站的H和κ被替換為最佳估計(jì)的H和κ值，從而獲得研究區(qū)域內(nèi)所有待反演地震臺(tái)站的地殼厚度H和波速比κ。

其中，所述單一臺(tái)站聯(lián)合反演的過程具體為：

(1)將以所述待反演地震臺(tái)站為中心的一定窗口內(nèi)的所有地震臺(tái)站的地殼厚度H、波速比κ和布格重力異常數(shù)據(jù)，分別網(wǎng)格化成相同大小的網(wǎng)格數(shù)據(jù)；

窗口大小設(shè)置一般遵循如下規(guī)律：網(wǎng)格數(shù)大于11×11網(wǎng)格，網(wǎng)格間距大于1/4的臺(tái)站平均間距且小于臺(tái)站平均間距；

(2)用步驟(1)獲得的網(wǎng)格化后的地殼厚度H、波速比κ和布格重力異常，結(jié)合線性迭代算法計(jì)算出莫霍面密度差Δρ_Moho和地殼密度隨波速比變化率

所述地殼密度隨波速比變化率的計(jì)算使用以下公式：

Δg_B＝Δg_Moho+Δg_Crust (1)

其中,Δg_B、Δg_Moho及Δg_Crust分別代表布格重力異常、莫霍面重力異常及地殼重力異常；Δρ_Moho為莫霍面上下密度差,為地殼密度隨波速比的變化率,為研究區(qū)地殼厚度的平均值,為研究區(qū)波速比的平均值,F{·}和F^-1{·}分別表示傅里葉變換和反傅里葉變換,G為萬有引力常數(shù),f為波數(shù),c為地殼厚度變化的校正因子；H為地殼厚度，κ為波速比。附圖4中

(3)利用求得的莫霍面上下密度差Δρ_Moho和地殼密度隨波速比的變化率計(jì)算出模型理論重力異常；模型理論重力異常的計(jì)算具體為：將步驟(2)中計(jì)算獲得的密度參數(shù)反代入步驟(2)的公式(2)、(3)中計(jì)算出Δg_Moho和Δg_Crust，計(jì)算所述模型理論重力異常Δg_Model公式為：

Δg_Model＝Δg_Moho+Δg_Crust (4)

其中：Δg_Moho是莫霍面重力異常，Δg_Crust為地殼重力異常。

(4)將布格重力異常與所述模型理論重力異常作差，求得異常殘差，即異常噪音，采用似然估計(jì)法分別計(jì)算出異常殘差的均值μ與方差σ²，再計(jì)算重力似然函數(shù)值；

假設(shè)異常殘差服從高斯分布，因而采用下式計(jì)算μ與σ²：

其中：是待反演地震臺(tái)站所在點(diǎn)上的重力異常殘差，n代表計(jì)算窗口內(nèi)的網(wǎng)格點(diǎn)數(shù)。

(5)以接收函數(shù)H-κ疊加圖相同的范圍和步長(zhǎng)選取H和κ值，重復(fù)步驟(3)和(4)，形成重力反演的H-κ似然圖(即附圖4中的MLE(maximum likelihood estimation))，重復(fù)步驟(3)時(shí)繼續(xù)使用步驟(2)中第一次計(jì)算的密度參數(shù)Δρ_Moho和重復(fù)步驟(4)時(shí)繼續(xù)使用步驟(4)中第一次計(jì)算的噪聲參數(shù)μ與σ²；

在與接收函數(shù)H-κ疊加圖中H和κ相同的參數(shù)范圍內(nèi)，變動(dòng)所述待反演地震臺(tái)站的H和κ，當(dāng)似然值最大時(shí)，擬合出的模型重力異常與布格重力異常最接近，故此時(shí)的H和κ為最佳估計(jì)值，重力反演的H-κ似然圖采用如下算法：

其中：L(μ,σ²)表示最大似然值；μ與σ²分別表示異常殘差的均值和方差，是待反演地震臺(tái)站所在點(diǎn)上的重力異常殘差；n代表計(jì)算窗口內(nèi)的網(wǎng)格點(diǎn)數(shù)。

(6)將重力反演的H-κ似然圖與接收函數(shù)H-k疊加圖均歸一化，通過兩矩陣點(diǎn)乘實(shí)現(xiàn)相互加權(quán)，得到聯(lián)合反演的H-κ圖，從聯(lián)合反演的H-κ圖中拾取最佳估計(jì)的H和κ值。

在本實(shí)施例所述方法中，不斷更新并迭代研究區(qū)域內(nèi)所有待反演地震臺(tái)站的H和κ，以避免由于地殼厚度H和波速比κ初值誤差而引起的聯(lián)合反演結(jié)果不準(zhǔn)確的情況。

對(duì)上述方法進(jìn)行模型測(cè)試：

步驟一、建立一個(gè)平均深度為35km的地殼厚度模型，一個(gè)平均波速比為1.8的波速比模型，網(wǎng)格大小均為21×21，網(wǎng)格間距為15km。設(shè)置Δρ_Moho＝0.5g/cm³，利用如下公式正演重力異常：

Δg_B＝Δg_Moho+Δg_Crust

其中,Δg_B、Δg_Moho及Δg_Crust分別代表布格重力異常、莫霍面重力異常及地殼重力異常；Δρ_Moho為莫霍面上下密度差,為地殼密度隨波速比的變化率,為研究區(qū)地殼厚度的平均值,為研究區(qū)波速比的平均值,F{·}和F^-1{·}分別表示傅里葉變換和反傅里葉變換,G為萬有引力常數(shù),f為波數(shù),c為地殼厚度變化的校正因子；H為地殼厚度，κ為波速比。

步驟二、以區(qū)域中心的地殼厚度和波速比為參考建立兩種地殼模型，分別為簡(jiǎn)單模型和復(fù)雜模型，其地殼厚度均為35km，然后正演接收函數(shù)；

步驟三、模型區(qū)域內(nèi)隨機(jī)產(chǎn)生65個(gè)點(diǎn)位作為臺(tái)站位置，其中有一個(gè)位于區(qū)域中心，假設(shè)為需要聯(lián)合反演的臺(tái)站，提取這些點(diǎn)上的地殼厚度與波速比，以此作為初始值；

步驟四、將65個(gè)臺(tái)站的地殼厚度與波速比網(wǎng)格化成與重力一樣的網(wǎng)格大小，即網(wǎng)格大小為21×21，網(wǎng)格間距為15km；

步驟五、使用上述網(wǎng)格化后布格重力異常及H、κ,用線性迭代算法計(jì)算出Δρ_Moho與計(jì)算出Δρ_Moho＝0.52g/cm³，

步驟六、利用求得的密度參數(shù)Δρ_Moho與計(jì)算出模型理論重力異常；

步驟七、將布格重力異常與模型重力異常作差求得異常噪音，采用似然估計(jì)法分別計(jì)算出μ＝0與σ²＝4.85，將這兩個(gè)參數(shù)帶入以下步驟中計(jì)算重力似然值；

步驟八、以接收函數(shù)H-κ疊加相同的范圍和步長(zhǎng)選取H和κ值，重復(fù)步驟五和六，形成重力反演的H-κ似然圖，如圖2，重復(fù)步驟五時(shí)繼續(xù)使用第一次計(jì)算的密度參數(shù)，重復(fù)步驟六時(shí)繼續(xù)使用第一次計(jì)算的噪聲參數(shù)；

步驟九、分別將兩種地殼模型的接收函數(shù)H-κ疊加圖與重力反演的H-κ似然圖歸一化并點(diǎn)乘，得到聯(lián)合反演的H-κ圖，如圖3，可見無論是對(duì)應(yīng)于接收函數(shù)數(shù)據(jù)數(shù)量和質(zhì)量問題的簡(jiǎn)單模型，還是對(duì)應(yīng)于存在沉積層和殼內(nèi)低速層的復(fù)雜模型，與重力聯(lián)合后均得到了一個(gè)相對(duì)穩(wěn)定的值，從中拾取最佳估計(jì)的H和κ值分別為35km/1.85與36km/1.84，與模型設(shè)定值的誤差在可接受范圍內(nèi)，因此，模型試驗(yàn)結(jié)果證明該方法有效。