伍 晗，寧杰遠(yuǎn),2*

（1.北京大學(xué)地球與空間科學(xué)學(xué)院，北京 100871；2.河北紅山地球物理國家野外科學(xué)觀測研究站，河北 邢臺 054000）

0 引言

背景噪聲成像作為一種不依賴于天然地震分布、對淺部結(jié)構(gòu)更敏感的成像技術(shù)，近些年來得到了廣泛關(guān)注[1-18]。然而背景噪聲成像的準(zhǔn)確度，受到各種因素的影響。在許多研究中[3-5]，假定背景噪聲源是均勻分布的。然而，噪聲源的實(shí)際分布既不是各向同性的，也不是始終穩(wěn)定分布的，而是在方位和時(shí)間上都有變化的[15]；而且周期不同，噪聲源的分布也會有所差異[19]。

研究噪聲源的時(shí)空變化，可以用于監(jiān)測地殼運(yùn)動、海洋運(yùn)動和人為活動的變化。此外，對結(jié)構(gòu)的精細(xì)成像也需要考慮噪聲源的不均勻分布帶來的測量誤差。有些學(xué)者定量地從理論和實(shí)際數(shù)據(jù)出發(fā)，對源的不均勻分布對背景噪聲成像結(jié)果的影響進(jìn)行了研究[11,15,20-21]，多數(shù)人認(rèn)為噪聲源的不均勻分布所導(dǎo)致相速度測量偏差不大，但需進(jìn)一步確認(rèn)[15,21]。噪聲源不均勻也會影響從背景噪聲互相關(guān)中提取用于研究介質(zhì)衰減的振幅信息。在定量描述噪聲源的不均勻分布造成的上述偏差時(shí)，首先要對噪聲源分布有定量的估計(jì)。

本文將在速度模型和臺陣內(nèi)部臺站對地震波形互相關(guān)函數(shù)均給定的條件下，嘗試基于Tikhonov正則化[22]方法，利用互相關(guān)函數(shù)和噪聲源分布的線性關(guān)系，反演背景噪聲源的分布。該方法相較于前人的方法[11,15]，在噪聲源分布反演時(shí)可以得到更為穩(wěn)定的解，有利于通過迭代得到更加正確的噪聲源分布及波速反演結(jié)果。

1 臺站對互相關(guān)函數(shù)的推導(dǎo)及合成

1.1 平面波假設(shè)

當(dāng)主要噪聲源和臺陣之間的距離遠(yuǎn)大于臺陣的規(guī)模時(shí)，如果不考慮臺陣內(nèi)或附近的不均勻結(jié)構(gòu)引起的局部散射，也不考慮臺陣附近的噪聲源對背景噪聲干涉的貢獻(xiàn)，可以假設(shè)彈性波能量是以平面波的形式傳播。在區(qū)域臺陣面波分析和背景噪聲干涉理論推導(dǎo)中，為了簡化，平面波近似被廣泛使用[8,11,15,21]。在平面波假設(shè)下，噪聲源的能量分布僅僅與周期和入射角度有關(guān)，用E(ω,θ)表示。其中：ω表示頻率，θ表示平面波入射角度和正東方向的夾角（圖1）。

圖1 平面波入射模型示意圖

1.2 互相關(guān)函數(shù)的推導(dǎo)

對于某一個(gè)固定頻率 ω，入射角固定的二維平面波可以表示為：

式中：A表 示振幅；ω表示頻率；表示平面波的傳播方向；表示觀測點(diǎn)的位置。

對于臺站對AB，從某一方向入射平面波的互相關(guān)函數(shù)為

將（1）式代入（2）式，可以得到

式中：δt是入射角為 θ的平面波到臺站對AB的走時(shí)差。由圖1 可知，在射線近似下兩者的時(shí)間差 δt滿足

式（3）是對于某一個(gè)角度入射平面波對于互相關(guān)函數(shù)的貢獻(xiàn)，那么對 θ進(jìn)行積分就能得到AB臺站對的互相關(guān)函數(shù)：

根據(jù)Yao et al.[15]，對臺站對的互相關(guān)函數(shù)離散化，可以表示為：

式中：Cn(ω,t)表示互相關(guān)函數(shù)；Δn表示第n個(gè)臺站對的臺間距；Wn(t,Δn)表示窗函數(shù)；Em(θm)是角度為θm頻率為ω的噪聲源能量；δtnm是第m個(gè)入射角對應(yīng)的平面波在第n個(gè)臺站對的走時(shí)差；Hnm(t,δtnm)是taper函數(shù)。

相應(yīng)的taper 系數(shù)和窗函數(shù)定義為：

式中：T為平面波的周期；T*是Hnm(t,δtmn)的寬度，這里取T*=5T；vminvmax表示最大和最小的參考速度，用于確定主要的面波時(shí)窗，這里vmin=2 km/s,vmax=5 km/s。

基于公式（6），在已知速度模型和噪聲源分布的情況下，可以理論地合成固定頻率T=30 s時(shí)臺站對的互相關(guān)函數(shù)。當(dāng)噪聲源隨方位角分布的輸入模型如圖2a 所示，速度模型和臺站分布如圖2b 所示，可以計(jì)算出臺站對所對應(yīng)的互相關(guān)函數(shù)結(jié)果（圖2c）。從圖2 中可以知道，即使速度結(jié)構(gòu)是均勻的，噪聲源的不均勻分布也會導(dǎo)致互相關(guān)函數(shù)有明顯的不同，會對到時(shí)和振幅的測量產(chǎn)生不可忽略的影響。

圖2 互相關(guān)函數(shù)合成示意圖

2 利用Tikhonov 正則化對噪聲源進(jìn)行反演

2.1 利用Tikhonov 正則化對噪聲源分布進(jìn)行估計(jì)

在已知臺站對的互相關(guān)函數(shù)和參考速度模型時(shí)，可以線性反演出噪聲源隨方位角的分布。為了方便計(jì)算，將時(shí)間域的記錄轉(zhuǎn)換到頻率域。

利用傅里葉變換，將（6）式變換到頻率域：

取傅里葉變換后，頻率 ω對應(yīng)的實(shí)部和虛部，式（9）即可用矩陣的方式表達(dá)為

式中：

DR=Re{F{Cn(ω,t)Wn(t,Δn)}}，是N×1的向量；

根據(jù)以上描述，可以將互相關(guān)函數(shù)、噪聲源分布、速度模型參數(shù)用線性方程組聯(lián)系起來。在已知互相關(guān)函數(shù)和速度模型的情況下，對噪聲源分布進(jìn)行線性反演。該問題由于矩陣非零的最大和最小奇異值之比很大，是離散病態(tài)問題，不能通過標(biāo)準(zhǔn)的方法直接求出有意義的解，可以利用更復(fù)雜的正則化方法，得到有意義的解。本文將采用Tikhonov正則化方法[22]，對噪聲源的能量分布進(jìn)行反演。

Tikhonov 正則化針對離散病態(tài)問題，能夠估計(jì)出穩(wěn)定的解。定義目標(biāo)函數(shù)χ(E)：

式中：λ是約束系數(shù)，由L-曲線準(zhǔn)則[23]確定。

對目標(biāo)函數(shù)求極小值，可以得到反演估計(jì)的解E*的表達(dá)式：

式中：I是單位矩陣。

至此，本文利用Tikhonov 正則化，通過（12）式可以對噪聲源方位角分布進(jìn)行估計(jì)。

下面展示一個(gè)例子，來看噪聲源分布的反演結(jié)果與輸入模型的比較。噪聲源的輸入模型和速度模型都如圖2 所示，均勻模型、噪聲源能量隨角度不均勻分布，反演時(shí)用到的速度模型與輸入的速度模型一致。利用（12）式進(jìn)行反演，得到的反演結(jié)果如圖3 所示。噪聲源分布的反演結(jié)果與輸入模型基本完全吻合，說明在利用Tikhonov 正則化能夠得到一個(gè)穩(wěn)定的，與真解相近的解。

圖3 噪聲源分布反演結(jié)果與輸入模型對比圖

2.2 對背景噪聲源反演結(jié)果的評價(jià)

由于噪聲源分布反演問題的解具有非唯一性，在完成解的估計(jì)之后，還需要對噪聲源反演的結(jié)果進(jìn)行評價(jià)。

定義正演算子A：

定義反演算子A*：

那么相應(yīng)的模型分辨率矩陣Rmodel，數(shù)據(jù)分辨率矩陣Rdata，單位協(xié)方差矩陣Cov(E*)為：

在圖3 的例子中，可以計(jì)算出相應(yīng)的模型分辨率矩陣、數(shù)據(jù)分辨率矩陣、單位協(xié)方差矩陣（圖4）。

模型分辨率矩陣展現(xiàn)了模型的估計(jì)值對模型真實(shí)值的逼近程度，如果Rmodel等于單位矩陣I，說明模型可以完全分辨，否則模型參數(shù)的估計(jì)值就是真實(shí)模型參數(shù)的加權(quán)平均值；如圖4a 所示，模型分辨率矩陣的每行有中心位于對角線上的窄峰，接近于I，說明在反演中可以很好地分辨模型。數(shù)據(jù)分辨率矩陣描述了預(yù)測數(shù)據(jù)與觀測數(shù)據(jù)的擬合程度，如果Rdata等于單位矩陣I，說明數(shù)據(jù)預(yù)測的誤差為0；如圖4b 所示，數(shù)據(jù)分辨率矩陣的每行有一個(gè)尖銳的極大值，其中心位于對角線上，說明數(shù)據(jù)在反演中也得到了很好的擬合。數(shù)據(jù)的單位協(xié)方差矩陣反映了數(shù)據(jù)中誤差的放大程度，如果數(shù)據(jù)中各分量是獨(dú)立互不相關(guān)的，那么協(xié)方差矩陣就是一個(gè)單位矩陣；如圖4c 所示，單位協(xié)方差矩陣不是單位矩陣，說明數(shù)據(jù)中的分量是相關(guān)的，單位協(xié)方差矩陣的值主要集中在對角線附近，說明數(shù)據(jù)中距離相近的分量相關(guān)性比較強(qiáng)，單位協(xié)方差矩陣的幅值表明了誤差的放大程度；圖4c 所示的幅值最大為8 左右，說明反演方法對數(shù)據(jù)誤差的響應(yīng)并不大，反演結(jié)果比較穩(wěn)定。

3 對反演方法的穩(wěn)定性測試

為了測試上述反演方法的穩(wěn)定性，設(shè)計(jì)了一系列不同的模型來展示在不同情況下解的估計(jì)。主要是從速度模型、不同噪聲源分布模型、周期和增加噪聲的角度去測試。

3.1 不同速度模型

為了測試反演方法在不同的速度模型下的穩(wěn)定性，在保證噪聲源分布如圖2a 所示，周期T=30 s等其他條件不變的情況下，設(shè)計(jì)了以下3 種情況。

第一種情況：合成互相關(guān)記錄DR,DS的速度模型為圖5a 所示的均勻模型，反演時(shí)計(jì)算MR,MS的速度模型也是圖5a 所示的均勻模型。在這種情況下，輸入模型與反演模型一致，且該模型為均勻模型，結(jié)構(gòu)比較簡單。

第二種情況：合成互相關(guān)記錄DR,DS的速度模型為圖5b 所示的不均勻速度模型，反演時(shí)計(jì)算MR,MS的速度模型也是圖5b 所示的不均勻速度模型。在這種情況下，輸入模型與反演模型一致，但速度模型為不均勻模型，結(jié)構(gòu)相對復(fù)雜一些。

第三種情況：合成互相關(guān)記錄DR,DS的速度模型為圖5b 所示的不均勻速度模型，反演時(shí)計(jì)算MR,MS的速度模型是圖5a 所示的均勻速度模型。在這種情況下，輸入模型與反演模型不一致，用于檢驗(yàn)在準(zhǔn)確的速度模型未知時(shí)反演噪聲源分布的穩(wěn)定性。

圖5 速度模型與臺站分布示意圖

根據(jù)圖6a～6b，在已知準(zhǔn)確速度模型的情況下，該方法可以準(zhǔn)確估計(jì)出噪聲源能量隨角度的分布情況，反演結(jié)果與輸入模型基本一致。

根據(jù)圖6c，在并不知道準(zhǔn)確的速度模型時(shí)，用一個(gè)近似的速度結(jié)構(gòu)去反演噪聲源分布時(shí)，大致的圖像還是相同的，但在細(xì)節(jié)上還是有一些差別。在實(shí)際情況中，一般并不知道準(zhǔn)確的地下速度結(jié)構(gòu)，則可以根據(jù)Tikhonov 反演方法對噪聲源能量分布進(jìn)行相對準(zhǔn)確的估計(jì)。這也是使用迭代方法去更新速度結(jié)構(gòu)和噪聲源能量分布的依據(jù)。

圖6 不同速度模型下噪聲源反演結(jié)果與輸入模型的對比圖

3.2 不同噪聲源分布模型

為了測試反演方法在不同的噪聲源能量隨方位角分布的模型下的穩(wěn)定性，在保證速度模型為圖2b所示的均勻模型、周期T=30 s等其他條件不變的情況下設(shè)計(jì)了3 種情況：①噪聲源均勻分布（圖7a）；②噪聲源一般不均勻分布，變化比較連續(xù)（圖7b）；③噪聲源特別不均勻分布，某些角度能量為1，某些角度能量為0（圖7c）。

圖7 不同噪聲源分布情況下噪聲源反演結(jié)果與輸入模型的對比圖

根據(jù)圖7a～7b，在噪聲源能量分布隨角度變化比較連續(xù)的情況下，該方法可以準(zhǔn)確估計(jì)出噪聲源能量隨方位角的分布情況，反演結(jié)果與輸入模型基本一致。

根據(jù)圖7c，在輸入的噪聲源模型特別不均勻的時(shí)候，反演結(jié)果大致的樣式與輸入模型還是相同的，但是反演結(jié)果在峰值和谷底會有一些起伏。

根據(jù)以上結(jié)果可以知道，Tikhonov 反演方法對于不同的噪聲源能量分布模型，能夠比較準(zhǔn)確地估計(jì)其分布；模型變化越平緩，反演估計(jì)出的解就越接近真解。

3.3 不同周期

為了測試反演方法在不同的周期下的穩(wěn)定性，在保證噪聲源分布（圖2a），且速度模型為均勻模型（圖2b）的情況下，考慮3 種周期T=10 s、T=30 s、T=50 s情況下噪聲源的反演結(jié)果（圖8）。

圖8 不同周期下噪聲源反演結(jié)果與輸入模型的對比以及相應(yīng)的臺站互相關(guān)函數(shù)圖

由圖8 可以看出，在3 種不同的周期下，反演估計(jì)出的噪聲源能量分布和輸入模型都基本一致。說明在多種周期下，利用Tikhonov 正則化方法均能對噪聲源方位角能量分布進(jìn)行準(zhǔn)確反演。

4 結(jié)論

噪聲源能量隨角度分布的反演問題本身是一個(gè)離散病態(tài)問題，直接利用最小二乘法，無法得到有效的解。針對該問題，應(yīng)用Tikhonov 正則化方法，在預(yù)設(shè)的各種情況下實(shí)現(xiàn)了對真解的近似，得到了與輸入模型相近的、有意義的解。

通過對不同條件的測試，Tikhonov 正則化方法在不同的速度模型、噪聲源分布模型以及周期下，對解的估計(jì)都是有效且準(zhǔn)確的。

噪聲源的不均勻分布，會影響相速度的測量，Tikhonov 正則化方法能夠在未知準(zhǔn)確的速度模型情況下，對噪聲源分布進(jìn)行反演。之后可以根據(jù)噪聲源反演結(jié)果，對測量的相速度進(jìn)行矯正，從而改善由于噪聲源的不均勻分布對背景噪聲成像的影響。

致謝研究工作得到北京大學(xué)高性能計(jì)算校級公共平臺支持；感謝參與討論并給出建議的蔣一然、溫景充、鮑鐵釗、石永祥、殷常陽等同學(xué)。