• / 7
  • 下载费用:1 下载币  

地震学百科知识(三)——地震学反演问题

关 键 词:
地震学 百科 知识 反演 问题
资源描述:
第4期(总第412期) 2013年4月 国际地震动态 o.4(o.412) 013 地震学百科知识(三) 地震学反演问题* 许忠淮 (中国地震局地球物理研究所,北京100081) 中图分类号:1; 文献标识码: A;0.3969/J.235—4975.2013.04.009 引言 求解地震学反演问题是一般数学物理反 演问题的分析方法在地震学研究中的应用。 一般反演问题的提出可表述如下。 当人们研究一个物理对象时,首先要选 定一些描述对象特征的参量X,根据已知的 物理规律列出 变化所遵从的数量联系,即 遵从的方程 F( )===0 (1) 函数地震波的波动方程、 地震波的走时方程等)由研究对象运动变化 的具体规律决定。矢量 一( ,32 ,…; ) 表示共有标函 数矢量 ,F。,…,F ) 表示式(1)代 表含(z1,,32L)一0 j …… (2) 1 J FⅣ(2,…,0 对许多问题,人们还没有完全认识有关 的规律,只是在观测和实验的基础上提出表 示各参量之间可能物理联系的假设关系式来 试验性地分析问题,这时象的模型。参量部分 是未知的构建模型用的模型参数m;另一部 收稿日期:2012—11是可观测的参数d。观测量是会有误差 的,但在一定的认识条件和认识水平下,有 些量的误差可认为很小,它们对问题分析结 果的影响可忽略,这些量就被当作模型中的 常量。考虑模型参数和观测参数的划分后, 式(1)可改写为 F(m,d)=0 (3) 对有些实际问题,式(3)可以写成以下观测 方程的形式口]: G(m)=d (4) 根据设定的模型关系式,对给定的模型 参量,通过分析演绎,预测会有怎样的观测 量的问题,称为正演问题;由已知的观测 值,根据模型推断可能有怎样的模型参数值 的问题称为反演问题。 在求解反演问题之前,人们必须先能求 解正演问题,即需要能理解产生观测结果的 物理过程,以便建立该过程的可靠的数学 模型。 在确定模型F(或G)和已知数据,反演问题就归结为根据式(3)或式(4)求 解模型参数优的数学问题。对不少实际问 题,函数F(m,d)或G(非线性函数, 这时要求解的数学问题是个非线性反演问 题;当它们是线性函数时,则要求解的是线 性反演问题。对于许多实际问题,这一求解 过程常常是比较复杂的,目前已发展了许多 36 国际地震动态 求解反演问题的数学方法。 在地震学中,人们已经研究了许多“正 演问题”的解,即对给定的地球结构和地震 震源的模型,人们已计算出了可观测的地震 运动的特性,例如,各种地震波的走时、地 震波的频散曲线、地震引起的强地震动与近 场地震动的频谱、地球自由震荡周期、远场 地震波的波形和完整的地震图,等等。“反 演问题”就是将这些计算结果与实际观测数 据进行比较,使计算结果能最好地拟合观测 数据,以便确定描述地球结构和地震震源的 各种模型参数。 1近震定位问题 作为地震学反演问题的例子,现分析一 个最简单的近地震的定位问题。 设均匀水平地壳内 已知量。选取一原点 轴向东、地壳内 在时刻r(发震时刻)、深度为 处发生一个 小地震,其震中位置的坐标为( , )。设地 震周围有方位不同、离地震远近不同的N 个地震台接收到地震发出的直达时间t ( 一1,2,…,N),各台的空间位 置( 2= )可认为足已知量。根据所用 的均匀地壳模型,各台站的理论到时T,应 等于发震时刻Tj(r )一r+ 一r+ (z—z ) +(3,一Y ) +(z— ) P (,…,N) (5) 式中R 表示地震至台站们求 解的未知地震参数(r, ,Y,z)应能使各台 的计算到时与各台的观测到时相符: 丁J(r,z,Y,z)一t , ( 一1,2,…,N) (6) 可将上式的(m):d (7) 式中 一(r,z,Y,z) 是4维的模型参数 矢量,T(m)===(T1(m),…,2)) 是N 维的计算到时矢量,, ,…, )是 (7)所示的地震定 位问题正是一般反演问题公式(4)的一个具 体表达形式。式(6)或(7)可称为地震定位问 题的观测方程组。由式(5)可见,各台的计 算走时T,是模型参数的非线性函数,即 T(m)代表一组非线性函数,因此,这里的地 震定位问题是个非线性反演问题。 由于定位模型的近似性和观测存在误差 等因素,一般是不能找到使式(6)或(7)严格 成立的模型参数的。方程组(7)中未知数是4 个,一般台站数式(7) 是个超定方程组,无法求数学上严格成立的 解答。求解反演问题通常是求能使残差矢量 ,一T(m)一d (8) 尽量小的模型参数。Tj(r, z,Y,z)一t 是第的差异,即第们所要 求的解答是要使各台的计算到时与观测到时 总体上能尽量符合。 可以有不同的标准来衡量计算到时与观 测到时的符合程度,一种常用方法是使以下 的失配函数Q(亦称目标函数)取极小值: Q(r, ,Y, )一ll r N ∑ 一9) ll r r}+r;+…+r 是矢量几 里德)长度。由于上式是使各台的计算到时 与观测到时之差的平方和取极小,所以用式 (9)求解答的方法也称最小二乘法。也可取 其他形式的失配函数,例如使各台的计算到 时与观测到时之差的绝对值之和取极小: N Q 一>:1r l: N l (r,z,Y,z)一t f—=1 选取不同的失配函数常会使最后的解答有些 第4期 许忠淮等:地震学百科知识(三)——地震学反演问题 37 差异。 2反演算法 有多种从式(9)求解答m一(r,z,Y, ) 的方法,大致可分为两类。一类是只计算并 比较失配函数值的搜索法,这类方法多用来 求解非线性反演问题;另一类是需要计算失 配函数对未知参数的导数的方法,或称求导 法,使用该法的前提是存在失配函数对未知 量的一定阶数的导数。 (1)搜索法也有多种。例如随机尝试 法,亦称蒙特卡罗法,即在(r, ,Y,z)的可 能取值范围内用很多组(数量可很大)随机数 构成尝试解答,由式(9)计算每组尝试解的 后将与最小接受的解答。又如应用较广的遗传算法, 这是模仿生物进化过程而设计的一种不断改 变和更新一组尝试模型参数、逐步找到逼近 观测结果的优化模型参数的算法。 (2)求导法中简单的一种是只求一阶导 数的方法,现以上述近震定位的例子说明 如下。 首先要选定一组近似的初始模型参数 m 7r。,o,2o) 。例如,在台站分布 较密时,可将初至波到时最早的台站位置设 为震中的初始位置,根据以往对该地区地壳 地震的了解设定一个可能的震源深度值,进 而由各台波的到时减去由初始震源位置计算 的到时,可估算出一个平均的初始发震时 刻。也可用其他方法(如随机尝试法等)选定 初始模型参数。 由式(5)知,各台的计算到时T,是待求 模型参数的非线性函数,这样导致式(8)是 一个难以直接求解的非线性方程组。实际求 解的方法之一是将T,(m)在初值m∞ 处作泰 勒展开,并只保留一阶导数项 ]: 3个‘ T,(m)≈T,(m∞ )+ 上l(r— 等l。(70)+等l。( l( 一 一J & (10) 式中偏导数的下标0表示取初值m∽ 处的导 数值;后式中用A (是一1,2,3,4)表示了 前式中的4个偏导数值;后式中的3r—等等。对各个T,(m)如此作泰勒展开近 似后,求解非线性方程组(7)的问题遂可转 化为求解近似的线性观测方程组 Ax—b (11) 的问题,式中X===(3r,3x,3y,3z)一 一 m—m∞ 是待求的对初始模型参数m‘。 的修 正量;矩阵 A== aT ar v A“…21…~1…A~4 x 3x 3x y 3y 3y z 3z 3z (12) 后项矩阵的下标0表示矩阵内各偏导数皆取 初始模型参数 ∞ 处的数值。式(11)中的数 据矢量b—T(m∞ )一d。与观测方程(11)相 应的残差矢量r—A —b。 用最小二乘法求式(11)问题的解答,就 是要求使 Q—ll,( 丁A b+b b (13) 取极小值的解答 n’。分别求4个 分量的偏导数,并令这4个偏导数的表达式 等于0后,可得到 2A 28 国际地震动态 即 b (14) 该方程组称为线性方程组Ax—组,由于未知量是4个,矩阵A 是个 4×4的对称方阵,其本征值都是非负实数。 求矩阵正值 X‘ =(A ) A b (15) 将初始模型参数加上这一修正量,即得待求 的地震定位参数m 一 ‘。 +X“ 。通过迭代 法,每次将新求得的参数再当作初始模型参 数,重复上述求解过程,可求得第结果m ’一m +X ,直至m ’与 m 的差异已在容许误差之内,即得最终 的地震定位结果。 上述所用的求导法只在式(10)的计算到 时泰勒展开中保留了一阶导数项,对于非线 性较强的复杂问题,一阶导数近似法常收敛 很慢,有时甚至得不到合适解答;于是需要 在泰勒展开中保留2阶导数项,相应地已发 展了专门算法。 对于实际地震定位问题,需要考虑地壳 地震波速的分层结构,或波速随深度连续变 化的结构,甚至还要考虑波速的横向变化, 这时理论到时的计算方法要远比式(5)表达 的公式复杂得多,因而通常的地震定位问题 在数学上也是个高度非线性的问题。 3解的不唯一性 反演问题的解答一般都是不唯一的。这 主要是由于观测数据中缺少全部未知模型参 数的足够信息,有限的观测结果造成了数学 问题的欠定性。以求解式(11)所示的线性反 演问题 b 为例,对上述简化的地震定位问题,未知量 常台站数会多于4个, 或有些台可有不同震相的到时数据,即观测 方程数多于4个。该方程组是超定的,一般 应能求出使A ≈6的解答来。但实际问题常 使形式上的超定方程组实际是欠定的。例 如,当缺少靠近地震震中的近距离台站时, 单靠震中距较大的台站据不能将震源深度约束住,即观测矢量b 中缺少未知量深度的信息,在数学上表现为 式(11)中有些线性方程是近似相关的,真正 独立的方程少于4个。又如,当用分布在一 个小区域中的台网去测定在台网外较远处的 地震的位置时,若只用各台的,考虑观测误差后,各台的到时几乎相当 于一个地震台的观测结果,这时待解方程组 将是高度欠定的。 对其他未知模型参数很多的地震学反演 问题,一组观测值常常对一些参数是超定 的,而对另一些参数却是欠定的。例如,反 演多台记录的多个地震的震相到时数据,以 求某地区地震波速度局部变化的地震层析成 像问题,需要将研究区分成许多小的区块, 未知数是各区块的波速相对于全区平均速度 模型(初始模型)的扰动量。由于台站和地震 位置的局限性,有些区块可能没有地震射线 穿过,或只有稀少的、方向比较单一的射线 穿过,于是到时观测数据中就不含有这些区 块的波速信息,造成反演的数学问题对这些 区块的波速扰动量是欠定的。 此外,采用不同的计算方法,有时也会 得到不大一样的解答。例如,选择不同的失 配函数,有可能使我们所求的最优化解答结 果不同。对复杂的非线性反演问题,失配函 数常常是模型参数的多极值函数,当用搜索 法求解时,有时会将与失配函数局部极小值 (而不是全局极小值)对应的模型参数确定为 解答,这也会造成解答的不唯一。 以上所述是指对给定的物理模型,在求 解反演的数学问题中出现的解答的不唯一 性。物理模型的不确定性也会使问题的解答 含有不定性。例如,基于初步认识,人们对 地震定位问题选择了速度均匀的地壳模型, 第4期 许忠淮等:地震学百科知识(三)——地震学反演问题 39 如果另一些研究人员采用了地壳分层的速度 模型,则地震的定位结果将会有所不同。这 种模型的近似性可能会影响与之相关的反演 的数学问题是否适定,特别是不符合实际的 模型是难以从实际观测数据中找到合理的支 持证据的;但反演问题的不适定性主要是从 分析数学问题中来讨论的。 4算法的选择 在求解反演问题时,选择合适的计算方 法是非常重要的。 例如,对于经常遇到的待解的线性方程 组Ax:b,当论上可 通过求_。b。 但对于许多实际反演问题,矩阵是非常大的,这种情况下,即使,也不能用经典线性代数中的克莱姆法则 求其逆矩阵因为这需要耗费大量机 时,甚至不可能算出。但使用高斯消元法则 可很快求出解答。这说明解反演问题需要使 用合适的算法。 即使能顺利求出Ax=还有 解答是否稳定的问题。例如,对于问题_2] ( ( 有解答 —f 1,矩阵的行列式值为1.0。但 \0/ 上式右端的“观测数据”若稍有变动为 ( : ),则解就变为 一(一 :; ),解答就 完全不同了。实际观测数据总是带有一定误 差的,为了解问题Ax=,需要分析模型矩阵在实际反演问题中,数的矩阵,无法使用消元法之类的解法。 在上述地震定位问题中采用了最小二乘法来 求解X,导致求解正规方程组(14)而得到解 答式(15)。这一算法的优点是求解的式(14) 是一个低阶的适定方程组,矩阵 称方阵,其本征值全是非负实数,完全可 以求逆矩阵。但这一算法也有明显缺点:① 用计算机解正规方程时,计算精度需为求解 原始方程组的两倍;②组成矩阵 A ’坏了原始方程中的某些信息。 5奇异值分解法 对于一般的线性方程组 A × X ×l—b ×l (16) 式中各量下标表示行数乘列数,由于 ≠ , 而无法求逆矩阵A~。 数学家兰克卓斯(C.出了一种A 的奇异值分解方法口],由此可计算出逆矩阵。奇异值分解法将矩阵3个矩阵的乘积: A 一U 0 ···0 O i 0 …8p j } l J 2 J r I I 一一 1一一 一一 一一 一一 一一 式中, 阵,…,S )构成的对 角矩阵;矩阵u 由 , …,U )构成,每个矢量含 个分量;矩阵 y 由, ,…, )构成, 每个矢量含 个分量。这一分解方法的数学 推导可参阅专业书籍,例如《定量地震学》 12.3节。早有人已编制好奇异值分解法的 计算程序,下面给出一个具体分解实例,可见由右侧3个矩阵的乘积确 实可得到左侧矩阵的具体数值,并可见到矩 40 国际地震动态 阵U 和V 列向量的相互正交性(注意下式 最右端矩阵的行向量即是y 的列向量)。 1 2 —1 2 3 —1 1 3 0 —2 2 —2 —1 —1 0 一O。337 5 一O.098 7 O.567 9 —0.260 0 —0.697 3 数学上有人提出了一种关于广义逆矩阵 的定义如下。 设A为 x 的实矩阵,若某个 × 的矩阵,=H,(一 月一18) 则将这样的矩阵义逆矩阵。根据奇异值分解表示式(17),并注意 u;p,V; (19) I 是P×容易验证矩阵 A ===V S U; (20) 是满足式(18)中的诸条件的,因而A 就称 为矩阵式中的矩阵由矩阵 ,…,S )的倒 1 1 1 数f÷, ,…, 1构成的对角方阵。 1 52 P, 先后用矩阵A 和V 左乘Ax—,利用式(17)和式(2O)的分解表示,并注 意式(19)结果,可得到 V )一0 于是可求得方程组(16)中模型矢量 的广义 逆解为 X 一A b (21) 在求解线性方程组Ax=不存在 ,就可取式(21) 表达的广义逆解作为问题的解答。下面看看 该解答的一些特性。 (1)分辨率矩阵口 根据式(16)有b=入式(21)并考 O。326 3 —0.727 3 —0.364 6 一O.443 7 —0.186 4 一O.834 0 ——0.045 8 ——0.541 6 —0.09 5 1 0.004 4 虑式(17)对得 g S u;U ; === V V;x (22) 矩阵R=V V 是个m×是模型 参数的数目),称为广义逆解答的模型空间 的分辨率矩阵,该矩阵决定了广义逆解X 反映“真解”好的情况是近于单位矩阵J,这时广义逆解可很好地 表示真实解答。但对一般反演问题,它与单 位矩阵是会有较大偏离的。在上述地震定位 问题中,如果有包含近台的包围地震的台站 分布,甚至可利用清晰的深度震相的到时记 录,并且已有合理的震源参数的初始值,在 这样的条件下是可以将4个震源参数的修订 值定得比较准确的,在数学上就表现为由A 矩阵经奇异值分解所求得的分辨矩阵近J。反之,如果定位条件不好,则得出的 矩阵非对角元 素也有非零数值了。这时广义逆解中某个具 体模型参数将是包含其他真实参数的加权平 均结果,即解答被模糊了,不同参数不能被 很好地分辨清楚,甚至模型中的有些参数根 本无法“看见”。例如,在用地震层析成像方 法反演地震波走时求解地下地震波速度分布 时,如果某些小区域没有地震射线通过,则 我们的模型对这些区域的波速高低是没有任 何分辨能力的。 (2)解的误差估计测数据的误差会带来解答的误差。待 求解的方程组实际是6±这里假 4 O 6 5 7 9 5 1 7 1 7 6 O 0 O 一 一 2 4 3 O 6 1 6 3 1 6 4 6 O 0 O 一 一 9 6 5 4 3 5 3 4 0 7 5 4 O 0 O 一 一 一 ,,,...........,..................1 1●●●●●●●●●●●●J 9 2 0 0鼹 L 3 6 0 0 4 O 5 O O 5 ,, ........... ....................1 第4期 许忠淮等:地震学百科知识(三)——地震学反演问题 41 定于广义逆解,已有数学家 证明 ,由可由下式给出 ≤ 式中 —s 。 是矩阵值的比值,并称为矩阵1.1I 表示矢量的范数,等于各分量平方和的平方 根。式(23)表明,观测量相 对误差放大量的上界。式(23)说明,解的误 差不但与观测误差有关,也与模型矩阵构有关。 [1] [2] [3] 若矩阵>10 ),显著的奇异性,即不能求得稳定的解答。 对某个实际反演问题,如果矩阵值(这里的“零”是数值计算意义上的零, 即虽有很小的数值,但已与计算误差无法区 别),则不能用上述方法求得广义逆解,需 用另法(例如阻尼最小二乘法)求取反演问题 的可能解答。 参考文献 (作者电子信箱,许忠淮:美]安艺敬一,P 量地震学:理论和方法.北京:地震出版社,1987 [美]李泫鉴,S 震台网的原理及应用.唐美华,译;叶世元,赵仲和,校.北京:地震出 版社,1984 E, A, B.977 (上接第29页) of PS i,00036,he 7 iS PS of 999 to 0 1 1 e by n we of 7
展开阅读全文
  石油文库所有资源均是用户自行上传分享,仅供网友学习交流,未经上传用户书面授权,请勿作他用。
0条评论

还可以输入200字符

暂无评论,赶快抢占沙发吧。

关于本文
本文标题:地震学百科知识(三)——地震学反演问题
链接地址:http://www.oilwenku.com/p-63559.html
关于我们 - 网站声明 - 网站地图 - 资源地图 - 友情链接 - 网站客服客服 - 联系我们
copyright@ 2016-2020 石油文库网站版权所有
经营许可证编号:川B2-20120048,ICP备案号:蜀ICP备11026253号-10号
收起
展开