跳到主要内容
Sandbox Physics
地震成像 · 边学边做实验

L2 · 第 5 章

地震定位的误差椭圆,以及它没算进去的误差

用来更新试探震源的那个矩阵,同样告诉我们随机的到时误差能把结果推多远。由它得到协方差矩阵,再由协方差画出一个椭圆,按设定的概率包住真实震中。这个承诺可以靠计数来检验:只有到时误差时,它成立;速度模型错了几个百分点,它就失效,而椭圆本身看不出任何迹象。

到时误差变成协方差矩阵

假设每个到时都带有独立的误差,标准差相同。在解附近,定位结果近似是到时的线性函数。到时误差服从高斯分布,定位误差也就服从高斯分布,协方差等于到时方差乘以上一章正规矩阵的逆。如果各个到时的质量不同,就按各自的方差给每一行加权;定位实验假定 S 到时的读数误差是 P 到时的两倍。

Klein 的 HYPOINVERSE-2000 手册在附录 4 用问题的特征值和特征向量写出了同一个矩阵。到时方差取读数误差 RDERR 的平方,加上系数 ERCOF 乘以该次地震残差均方根(RMS)的平方。ERCOF 取 1 时,拟合较差的地震会得到更大的椭圆。

对角元是发震时刻、两个水平坐标和深度的方差;非对角元表示它们之间的此消彼长,例如深度和发震时刻的关系。去掉发震时刻,剩下震源的 3 阶块;其中两个水平坐标的行和列,构成震中的 2 阶块。

C=σ2(GTG)−1=w2 VS−2VT,w2=RDERR2+ERCOF⋅RMS2\mathbf C=\sigma^2\left(\mathbf G^{\mathsf T}\mathbf G\right)^{-1}=w^2\,\mathbf V\mathbf S^{-2}\mathbf V^{\mathsf T},\qquad w^2=\mathrm{RDERR}^2+\mathrm{ERCOF}\cdot\mathrm{RMS}^2

从协方差到椭圆

取协方差中震中那个 2 阶块。它的特征向量定出椭圆两轴的方向,特征值开平方得到沿两轴的标准误差,由此画出 1 倍标准差椭圆。要达到 95%,每根轴都要乘一个系数,系数大小取决于这个区域包含几个坐标。

只有一个坐标时,95% 的系数是熟悉的 1.96;两个坐标是 2.45;三个坐标是 2.80。1 倍标准差区域内的概率也随维数下降:线段上 68.3%,椭圆内 39.3%,椭球内只有 19.9%。

Klein 的手册写道,震源有 32% 的概率落在 1 倍标准差椭球内,乘以 2.4 可得 95% 的椭球。这两个数都和三维高斯分布对不上,三维的正确值是 19.9% 和 2.80。2.4 倒是接近二维的系数。定位实验的震中椭圆用的是二维系数。

k95=χν, 0.952:ν=1→1.96,  ν=2→2.45,  ν=3→2.80;P1σ=68.3%, 39.3%, 19.9%k_{95}=\sqrt{\chi^2_{\nu,\,0.95}}:\quad \nu=1\to1.96,\ \ \nu=2\to2.45,\ \ \nu=3\to2.80;\qquad P_{1\sigma}=68.3\%,\ 39.3\%,\ 19.9\%

数一数椭圆包住真值的次数

95% 椭圆给出了一个可以检验的说法:对许多次到时噪声相互独立的地震,大约每 100 个真实震中里应有 95 个落在各自的椭圆内。定位实验从同一个震源反复生成固定种子的合成地震,逐个定位,然后计数。

环形台网、使用 S 到时、到时噪声 0.1 s 时,1000 个真实震中有 944 个落在椭圆内。1000 次计数本身大约有 0.7 个百分点的随机起伏,所以 94.4% 与承诺相符。

速度模型错了,答案会移动,椭圆却不变

现在告诉定位程序岩石比实际快 5%,其他一切不变,把同一组到时重新定位。平均震中误差从 0.31 km 增加到 0.42 km,平均深度从 9.99 km 移到 9.62 km。椭圆几乎没变,因为它只取决于台站几何和假定的到时噪声,与速度误差无关。现在 1000 个真实震中里只有 849 个落在椭圆内。

拟合只变差了一点:平均残差均方根从 0.138 s 升到 0.173 s。定位程序可以通过挪动震源、调整发震时刻,把很大一部分速度误差吸收掉。手册的措辞因此很谨慎。Lee 和 Lahr 提醒,统计意义上的标准误差未必代表实际误差范围;Klein 写道,真正的误差椭圆需要知道地壳模型的不确定度,而这几乎不可能做到。

完整例题

把一个协方差换成椭圆

取震中协方差块:东向方差 0.040 平方千米,北向方差 0.010,协方差 0.010。

迹的一半是 0.025,对角元之差的一半是 0.015;把 0.015 和 0.010 各自平方、相加再开方,得 0.0180。两个特征值于是为 0.0430 和 0.00697 平方千米,沿两轴的标准误差是 0.2074 km 和 0.0835 km。

都乘以 2.448(即 5.991 的平方根),得到 95% 震中椭圆的两个半轴:0.508 km 和 0.204 km。长轴从正东向北偏转的角度,是正切为 0.020 比 0.030 的那个角的一半,即 16.8 度。

λ±=0.025±0.0152+0.0102=0.0430, 0.00697 km2;a95=2.448λ±=0.508, 0.204 km;θ=12arctan⁡0.0200.030=16.8∘\lambda_{\pm}=0.025\pm\sqrt{0.015^2+0.010^2}=0.0430,\ 0.00697\ \mathrm{km^2};\quad a_{95}=2.448\sqrt{\lambda_\pm}=0.508,\ 0.204\ \mathrm{km};\quad \theta=\tfrac12\arctan\frac{0.020}{0.030}=16.8^\circ

带着问题动手

先预测,再改变,最后解释

  1. 在定位实验中先选“环形台网:椭圆很小”,记下覆盖率读数。预测把假设波速的误差调到 5% 后,它会怎样变化。
  2. 把波速误差滑块调到 5%,或者直接选“速度模型错了”。
  3. 解释为什么椭圆大小几乎不变,覆盖率却下降了。

这些实验使用平坦、各向同性的弹性或声学岩石,波速是虚构但接近实际的数值,都不对应真实的观测。每个实验在模型说明里写明模型和适用范围,并可导出读数背后的数据。

试着自己解释

1. 1 倍标准差的震中椭圆要乘多大的系数才变成 99% 椭圆?

展开答案与解释

对两个坐标,落在 k 倍标准差内的概率是 1 减去 e 的负 k 平方除以 2 次方。令它等于 0.99,得 k 等于负 2 乘 0.01 的自然对数再开方,即 3.03。

k99=−2ln⁡0.01=3.03k_{99}=\sqrt{-2\ln 0.01}=3.03

2. 实验的覆盖率读数统计 40 次地震。如果椭圆是对的,计数会有多大的随机起伏?

展开答案与解释

计数服从二项分布,按百分比计的标准差是 0.95 乘 0.05 除以 40 再开方,约 3.4 个百分点。读数 90% 仍在两倍标准差以内;统计 1000 次时,起伏缩小到 0.7 个百分点。

3. 1000 次里有 849 次落在椭圆内,这与正确的 95% 椭圆相符吗?

展开答案与解释

不相符。期望值是 950 次,标准差约 6.9 次,849 比期望低了约 15 倍标准差。椭圆相对实际误差太小了。

可以放在旁边读的来源

  1. W. H. K. Lee and J. C. Lahr, HYPO71 (revised): a computer program for determining hypocenter, magnitude, and first motion pattern of local earthquakes, U.S. Geological Survey Open-File Report 75-311 (1975)

    第 4.1 节讲 Geiger 方法:把到时残差对发震时刻和三个坐标作一阶展开(式 2),解最小二乘正规方程(式 3 至 5),更新试探震源(式 6),然后重复。引言提醒:残差小、标准误差小,并不保证震源定得准。

  2. F. W. Klein, User’s guide to HYPOINVERSE-2000, a Fortran program to solve for earthquake locations and magnitudes, U.S. Geological Survey Open-File Report 02-171 (2002)

    附录 4 用定位问题的特征值和特征向量乘上到时方差,构造 4 行 4 列的协方差矩阵,再取其中只含空间坐标的 3 行 3 列得到误差椭球。附录指出,椭圆随假定的到时误差缩放;真正的误差还要算上地壳模型本身的不确定度,而这几乎无法知道。

反馈这个页面有错误、看不明白,或者缺了什么?告诉我们
反馈类型

提交时会附带本页地址、Lab ID、网站语言和设备类型。请不要填写个人信息。反馈最多保存 12 个月。隐私政策