一维速度结构的地震波全波形反演理论测试

小二乘法的时间域全波形反演方法, 并建立全波形反演理论的基本框架。Tromp等[9]论证了全波形反演与有限频层析成像的关系, Chen等[10]反演了洛杉矶地区的三维地壳结构, Wang等[11]反演了比利牛斯山地区的地壳结构, Fichtner[12]反演了澳大得亚地区的壳幔速度结构, Zhu等[13]反演了欧洲地区的壳幔速度结构, Simutė等[14]反演了日本列岛地区的壳幔

[15]结构, Bozdağ等 开展了全球三维各向异性地壳和地幔结构的全波形反演。

尽管全波形反演已在很多区域得到应用, 并取得较好的结果, 然而该方法的局限性尚未得到很好的解决。由于初始模型、数据质量以及搜索范围等的限制, 模型的更新容易陷入错误的局部极小值,也很难对反演结果的质量做出较准确的评价。理论测试是反演研究的重要部分, 通过理论测试, 可以对反演过程中各个因素的影响有一个定性的认识,

[16]为实际数据的反演及评价提供基础。Kolb等 研究了一维速度结构反演中的多尺度以及噪声等问题。Burstedde等[17]基于一维理论测试, 研究了反演问题中的非线性优化算法、正则化和多尺度反演等

[18]问题。Aleardi等 研究了基于混合遗传算法的一维全波形反演的不确定性估计问题。本文研究一维介质的速度结构反演问题。以残差平方和作为目标函数, 利用伴随状态法推导计算目标函数梯度的公式, 并使用共轭梯度法, 在不同尺度上反演模型参数, 讨论不同尺度、震源、台站、初始模型以及噪声等因素对反演结果的影响。由于反演计算与空间维数无关, 本文的结果可以为更高维数以及实际数据的反演提供有效的参考。

1 研究方法及模型设计1.1 全波形反演的基本思想

一维声波方程满足us (0)  0,  us (0)  0, (1)  t   u 2 2 m  us s  f s,   t2 x2其中, m为模型参数, 与速度c的关系为m=1/c2 , us为震源fs引起的压力场。地震波全波形反演大致包含3 个步骤: 1) 选取一个目标函数J(m), 衡量理论波形(正演计算得到的波形) us与观测波形ds的接近程度; 2) 求得目标函数的导数∂J/∂m; 3) 根据目标函数及

其导数的值, 使用非线性优化算法求得目标函数的极小值minmj(m)及其对应的模型参数m。

1.2 目标函数

目标函数是反演过程中重要的一环, 合适的目标函数能有效地提高反演的收敛速度和结果质量。常用的目标函数有残差平方和、互相关到时差和时频波形差等。其中, 残差平方和具有计算速度快、适用范围较广的优点, 本文用其作为目标函数:

gi  g T  H ( m  m0) (13) i代入式(12), 即可得到Fletcher-reeves公式: g T g   i 1 i 1。 (14) i g T gi i

1.5 多尺度反演

反演问题的一个难点是初始模型的选取。由于求解反演问题时, 很容易陷入局部极小值, 因此初始模型的选取对反演结果的好坏有很大的影响。多尺度反演是减小局部极小值影响的有效手段, 此方法将模型离散为不同的尺度, 首先使用低频信息,在大尺度上得到一个较粗略的结果, 然后以此结果作为初始模型, 逐渐提高使用的频率, 并进行更小尺度的反演。由于在大尺度上未知参数少, 更容易找到全局极小值, 因此在大尺度上反演的结果可作为对初始模型的一个很好的猜测。与直接在小尺度上反演相比, 多尺度反演的结果更稳定, 对初始模型的依赖相对较少, 收敛速度更快。

本文在实现多尺度反演的过程中, 没有改变正演用到的网格数, 而是直接改变模型参数的数量,并通过插值的方法, 映射到正演所用的网格。这样做的好处是, 可以更自由地选择反演的尺度, 从而更好地避免局部极小值的影响。对于目标函数J(m)及模型参数m, 反演的步骤如下。第一步, 取低网格数模型m′。第二步, 选取一个合适的频率, 对观测及正演数据进行滤波。

第三步, 用共轭梯度法求得m′的最优解: 1) 在正演计算中, 将m′映射到与m相同的维数; 2) 计算

 J  J m梯度 ; 3) 用共轭梯度法解出m′。 m  m m第四步, 以m′的反演结果作为下一次反演初始模型, 提高m′的网格数及滤波频率, 重复计算直到得到m的最优解。

真实速度模型采用一维非均匀介质, 速度结构如图 1 所示。

正演计算中将模型划分为250个网格, 即Δx= 0.004 km。反演的模型参数为5~100个, 并通过线性插值映射到250个网格点。总记录时长为3 s, 划分为1500个时间步, 即dt=0.002 s。模型顶部采用自由表面条件, 底部为吸收边界条件[20]:

I i