伍 晗，宁杰远,2*

（1.北京大学地球与空间科学学院，北京 100871；2.河北红山地球物理国家野外科学观测研究站，河北 邢台 054000）

0 引言

背景噪声成像作为一种不依赖于天然地震分布、对浅部结构更敏感的成像技术，近些年来得到了广泛关注[1-18]。然而背景噪声成像的准确度，受到各种因素的影响。在许多研究中[3-5]，假定背景噪声源是均匀分布的。然而，噪声源的实际分布既不是各向同性的，也不是始终稳定分布的，而是在方位和时间上都有变化的[15]；而且周期不同，噪声源的分布也会有所差异[19]。

研究噪声源的时空变化，可以用于监测地壳运动、海洋运动和人为活动的变化。此外，对结构的精细成像也需要考虑噪声源的不均匀分布带来的测量误差。有些学者定量地从理论和实际数据出发，对源的不均匀分布对背景噪声成像结果的影响进行了研究[11,15,20-21]，多数人认为噪声源的不均匀分布所导致相速度测量偏差不大，但需进一步确认[15,21]。噪声源不均匀也会影响从背景噪声互相关中提取用于研究介质衰减的振幅信息。在定量描述噪声源的不均匀分布造成的上述偏差时，首先要对噪声源分布有定量的估计。

本文将在速度模型和台阵内部台站对地震波形互相关函数均给定的条件下，尝试基于Tikhonov正则化[22]方法，利用互相关函数和噪声源分布的线性关系，反演背景噪声源的分布。该方法相较于前人的方法[11,15]，在噪声源分布反演时可以得到更为稳定的解，有利于通过迭代得到更加正确的噪声源分布及波速反演结果。

1 台站对互相关函数的推导及合成

1.1 平面波假设

当主要噪声源和台阵之间的距离远大于台阵的规模时，如果不考虑台阵内或附近的不均匀结构引起的局部散射，也不考虑台阵附近的噪声源对背景噪声干涉的贡献，可以假设弹性波能量是以平面波的形式传播。在区域台阵面波分析和背景噪声干涉理论推导中，为了简化，平面波近似被广泛使用[8,11,15,21]。在平面波假设下，噪声源的能量分布仅仅与周期和入射角度有关，用E(ω,θ)表示。其中：ω表示频率，θ表示平面波入射角度和正东方向的夹角（图1）。

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

1.2 互相关函数的推导

对于某一个固定频率 ω，入射角固定的二维平面波可以表示为：

式中：A表 示振幅；ω表示频率；表示平面波的传播方向；表示观测点的位置。

对于台站对AB，从某一方向入射平面波的互相关函数为

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

式中：δt是入射角为 θ的平面波到台站对AB的走时差。由图1 可知，在射线近似下两者的时间差 δt满足

式（3）是对于某一个角度入射平面波对于互相关函数的贡献，那么对 θ进行积分就能得到AB台站对的互相关函数：

根据Yao et al.[15]，对台站对的互相关函数离散化，可以表示为：

式中：Cn(ω,t)表示互相关函数；Δn表示第n个台站对的台间距；Wn(t,Δn)表示窗函数；Em(θm)是角度为θm频率为ω的噪声源能量；δtnm是第m个入射角对应的平面波在第n个台站对的走时差；Hnm(t,δtnm)是taper函数。

相应的taper 系数和窗函数定义为：

式中：T为平面波的周期；T*是Hnm(t,δtmn)的宽度，这里取T*=5T；vminvmax表示最大和最小的参考速度，用于确定主要的面波时窗，这里vmin=2 km/s,vmax=5 km/s。

基于公式（6），在已知速度模型和噪声源分布的情况下，可以理论地合成固定频率T=30 s时台站对的互相关函数。当噪声源随方位角分布的输入模型如图2a 所示，速度模型和台站分布如图2b 所示，可以计算出台站对所对应的互相关函数结果（图2c）。从图2 中可以知道，即使速度结构是均匀的，噪声源的不均匀分布也会导致互相关函数有明显的不同，会对到时和振幅的测量产生不可忽略的影响。

图2 互相关函数合成示意图

2 利用Tikhonov 正则化对噪声源进行反演

2.1 利用Tikhonov 正则化对噪声源分布进行估计

在已知台站对的互相关函数和参考速度模型时，可以线性反演出噪声源随方位角的分布。为了方便计算，将时间域的记录转换到频率域。

利用傅里叶变换，将（6）式变换到频率域：

取傅里叶变换后，频率 ω对应的实部和虚部，式（9）即可用矩阵的方式表达为

式中：

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

根据以上描述，可以将互相关函数、噪声源分布、速度模型参数用线性方程组联系起来。在已知互相关函数和速度模型的情况下，对噪声源分布进行线性反演。该问题由于矩阵非零的最大和最小奇异值之比很大，是离散病态问题，不能通过标准的方法直接求出有意义的解，可以利用更复杂的正则化方法，得到有意义的解。本文将采用Tikhonov正则化方法[22]，对噪声源的能量分布进行反演。

Tikhonov 正则化针对离散病态问题，能够估计出稳定的解。定义目标函数χ(E)：

式中：λ是约束系数，由L-曲线准则[23]确定。

对目标函数求极小值，可以得到反演估计的解E*的表达式：

式中：I是单位矩阵。

至此，本文利用Tikhonov 正则化，通过（12）式可以对噪声源方位角分布进行估计。

下面展示一个例子，来看噪声源分布的反演结果与输入模型的比较。噪声源的输入模型和速度模型都如图2 所示，均匀模型、噪声源能量随角度不均匀分布，反演时用到的速度模型与输入的速度模型一致。利用（12）式进行反演，得到的反演结果如图3 所示。噪声源分布的反演结果与输入模型基本完全吻合，说明在利用Tikhonov 正则化能够得到一个稳定的，与真解相近的解。

图3 噪声源分布反演结果与输入模型对比图

2.2 对背景噪声源反演结果的评价

由于噪声源分布反演问题的解具有非唯一性，在完成解的估计之后，还需要对噪声源反演的结果进行评价。

定义正演算子A：

定义反演算子A*：

那么相应的模型分辨率矩阵Rmodel，数据分辨率矩阵Rdata，单位协方差矩阵Cov(E*)为：

在图3 的例子中，可以计算出相应的模型分辨率矩阵、数据分辨率矩阵、单位协方差矩阵（图4）。

模型分辨率矩阵展现了模型的估计值对模型真实值的逼近程度，如果Rmodel等于单位矩阵I，说明模型可以完全分辨，否则模型参数的估计值就是真实模型参数的加权平均值；如图4a 所示，模型分辨率矩阵的每行有中心位于对角线上的窄峰，接近于I，说明在反演中可以很好地分辨模型。数据分辨率矩阵描述了预测数据与观测数据的拟合程度，如果Rdata等于单位矩阵I，说明数据预测的误差为0；如图4b 所示，数据分辨率矩阵的每行有一个尖锐的极大值，其中心位于对角线上，说明数据在反演中也得到了很好的拟合。数据的单位协方差矩阵反映了数据中误差的放大程度，如果数据中各分量是独立互不相关的，那么协方差矩阵就是一个单位矩阵；如图4c 所示，单位协方差矩阵不是单位矩阵，说明数据中的分量是相关的，单位协方差矩阵的值主要集中在对角线附近，说明数据中距离相近的分量相关性比较强，单位协方差矩阵的幅值表明了误差的放大程度；图4c 所示的幅值最大为8 左右，说明反演方法对数据误差的响应并不大，反演结果比较稳定。

3 对反演方法的稳定性测试

为了测试上述反演方法的稳定性，设计了一系列不同的模型来展示在不同情况下解的估计。主要是从速度模型、不同噪声源分布模型、周期和增加噪声的角度去测试。

3.1 不同速度模型

为了测试反演方法在不同的速度模型下的稳定性，在保证噪声源分布如图2a 所示，周期T=30 s等其他条件不变的情况下，设计了以下3 种情况。

第一种情况：合成互相关记录DR,DS的速度模型为图5a 所示的均匀模型，反演时计算MR,MS的速度模型也是图5a 所示的均匀模型。在这种情况下，输入模型与反演模型一致，且该模型为均匀模型，结构比较简单。

第二种情况：合成互相关记录DR,DS的速度模型为图5b 所示的不均匀速度模型，反演时计算MR,MS的速度模型也是图5b 所示的不均匀速度模型。在这种情况下，输入模型与反演模型一致，但速度模型为不均匀模型，结构相对复杂一些。

第三种情况：合成互相关记录DR,DS的速度模型为图5b 所示的不均匀速度模型，反演时计算MR,MS的速度模型是图5a 所示的均匀速度模型。在这种情况下，输入模型与反演模型不一致，用于检验在准确的速度模型未知时反演噪声源分布的稳定性。

图5 速度模型与台站分布示意图

根据图6a～6b，在已知准确速度模型的情况下，该方法可以准确估计出噪声源能量随角度的分布情况，反演结果与输入模型基本一致。

根据图6c，在并不知道准确的速度模型时，用一个近似的速度结构去反演噪声源分布时，大致的图像还是相同的，但在细节上还是有一些差别。在实际情况中，一般并不知道准确的地下速度结构，则可以根据Tikhonov 反演方法对噪声源能量分布进行相对准确的估计。这也是使用迭代方法去更新速度结构和噪声源能量分布的依据。

图6 不同速度模型下噪声源反演结果与输入模型的对比图

3.2 不同噪声源分布模型

为了测试反演方法在不同的噪声源能量随方位角分布的模型下的稳定性，在保证速度模型为图2b所示的均匀模型、周期T=30 s等其他条件不变的情况下设计了3 种情况：①噪声源均匀分布（图7a）；②噪声源一般不均匀分布，变化比较连续（图7b）；③噪声源特别不均匀分布，某些角度能量为1，某些角度能量为0（图7c）。

图7 不同噪声源分布情况下噪声源反演结果与输入模型的对比图

根据图7a～7b，在噪声源能量分布随角度变化比较连续的情况下，该方法可以准确估计出噪声源能量随方位角的分布情况，反演结果与输入模型基本一致。

根据图7c，在输入的噪声源模型特别不均匀的时候，反演结果大致的样式与输入模型还是相同的，但是反演结果在峰值和谷底会有一些起伏。

根据以上结果可以知道，Tikhonov 反演方法对于不同的噪声源能量分布模型，能够比较准确地估计其分布；模型变化越平缓，反演估计出的解就越接近真解。

3.3 不同周期

为了测试反演方法在不同的周期下的稳定性，在保证噪声源分布（图2a），且速度模型为均匀模型（图2b）的情况下，考虑3 种周期T=10 s、T=30 s、T=50 s情况下噪声源的反演结果（图8）。

图8 不同周期下噪声源反演结果与输入模型的对比以及相应的台站互相关函数图

由图8 可以看出，在3 种不同的周期下，反演估计出的噪声源能量分布和输入模型都基本一致。说明在多种周期下，利用Tikhonov 正则化方法均能对噪声源方位角能量分布进行准确反演。

4 结论

噪声源能量随角度分布的反演问题本身是一个离散病态问题，直接利用最小二乘法，无法得到有效的解。针对该问题，应用Tikhonov 正则化方法，在预设的各种情况下实现了对真解的近似，得到了与输入模型相近的、有意义的解。

通过对不同条件的测试，Tikhonov 正则化方法在不同的速度模型、噪声源分布模型以及周期下，对解的估计都是有效且准确的。

噪声源的不均匀分布，会影响相速度的测量，Tikhonov 正则化方法能够在未知准确的速度模型情况下，对噪声源分布进行反演。之后可以根据噪声源反演结果，对测量的相速度进行矫正，从而改善由于噪声源的不均匀分布对背景噪声成像的影响。

致谢研究工作得到北京大学高性能计算校级公共平台支持；感谢参与讨论并给出建议的蒋一然、温景充、鲍铁钊、石永祥、殷常阳等同学。