2007_Mogi模型在长白山天池火山区的应用_胡亚轩
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数才可达到较优的结果 (图 2 a)。而公式 (1 )、(2 )和 (3 )中 R 与位移有很好的比例关系 , 如果用 成倍的等效源半径 (如 RS 取 600m , 倍数为 10 )作为可能的理论值 , 相当于增大偏导数矩阵 (9 ) 的值 , 在做残差比较时再重新考虑 R 与各方向位移的比例关系 , 这样通过增大每次叠代时各模 型参数的改正向量的估计值 Δ b^, 减少叠代次数 , 在其它参数相同的情况下可以较快地得到最 优化的反演结果 (图 2 b)。另外 , 从图 2 还可以看出 , 反演时 2 种方法都存在 2 个区 , 其中 Ⅰ 区 说明当反演达到一定程度时 , 增加叠代次数 , 反演结果改进不很明显 ;Ⅱ区说明若给定此范围内 的任意值作为初始值 , 此范围内的极值解对给定的初始值依赖性很强 。 故应用此方法进行计算 时 , 注意形变量的大小和各参数的关系 , 合理选择参数值对反演结果极为重要 。
M og i模型对不同方向的位移变化响应有一定的局限性 。 D ie te rich等 (1975 )应用有限元模
型研究得出地表垂直位移对压力源的中心位置比较敏感 , 而水平位移对压力源的形状和深度有
所反映 。故分析火山区压力源参数时 , 水平变形和垂直变形两者都是必须的 。下面我们分析如
何应用垂直和水平形变数据联合反演火山区岩浆压力源的参数 。
为岩浆源中心深度 。
图 1 M ogi模型 坐标示意图 F ig. 1 Coord inate ske tch m ap o fM ogi mode.l
若岩浆压力源中心坐标为 (0, 0, - d), 依据式 (1) ～ (3), 当 r =0时 , 有最大垂直位移量 ;当
r =d / 2时 , 有最大水平位移量 , 且水平位移呈对称放射状分布 。
1期
胡亚轩等 :M og i模型在长白山天池火山区的应用
145
达式 。 假定地面直角坐标系为 (x, y, z), 岩浆压力源中心坐标为 (x0, y0 , - d )(图 1), 如果用等 效半径 R 表示岩浆流入或流出的体积 , 则根据 M og i模型 , 岩浆房膨胀所引起的地表位移表达式
分别为 (王勇 , 1996;胡亚轩等 , 2003;施
从公式 (1 )、(2 )和 (3 )可以看出 , 用等效球状源的参数 R (等效半径 )及源的位置参数 (x0 , y0 , d0 )作为反演参数应用阻尼最小二乘反演方法联合反演时 , 符合的条件为几何形变实测值与 模型理论值的残差平方和最小 , 如果考虑点位高程的影响 , 此条件可以表示为
n
m
∑ ∑ S =Qm in = [ ΔU zi - ΔUz (x, y, z;x0 , y0 , d0 , R )] 2 + [ ΔU xi -
1期
胡亚轩等 :M og i模型在长白山天池火山区的应用
147
M og i模型 , 当等效压力源体积变化较小时 , 反演时会出 现困难 。 从 3 个位移公式 (1 )、(2 )和 (3 )可以看出 , 压力源变化的大小 (等效半径 R )对 3 个方向的位移影响较大 。 不同大小的 R 值 对反演过程也有一定的影响 。 这里假定压力源等效半径 RS 为 60m , 位置参数 xS , yS , dS 分别为 3 000m , 4 500m , 5 000m , 如果给出地面形变监测点的坐标 x, y, z(这里取 z =0 ), 依据 M og i模型 公式可以计算出位移值 , 相当于实际工作中的观测值 。反演时假定 1 组初始的压力源参数 x0 , y0 , d0 , R0 分别为 300m、800m 、2 000m 和 300m , 计算得出的位移值为理论初始值 。应用以上方 法进行反演 。
2002年 6月汪清地震后 , 地震 、地化和形变等资料均表明长白山天池火山区经历了一次比 较强烈的岩浆上涌活动过程 , 而且火山变形明显地出现在天池旧火山口附近 , 水平位移接近放 射状 , 与 M og i模型的理论位移形状相符 (M og i, 1958;Yang et a l. , 1992)。由于岩浆的相对变化 较小 , ATA 是奇异矩阵 。 本文在理论证明阻尼最小二乘方法可行性的前提下 , 应用此方法反演 了长白山天池火山区压力源的位置和等效半径等参数 , 分析了近年来岩浆可能的活动特征 。
z为观测点的高程 , d 为压力源深度 。ΔU zi 、ΔU xi 和 ΔU yi 分别为实测垂直形变量 , 实测 X 方向的 水平分量及 Y 方向的水平分量 ;ΔUz 、ΔUx 和 ΔUy 分别为 M og i模型计算的各方向理论值 。
3 阻尼最小二乘反演初始值的确定
进行阻尼最小二乘反演时 , 一般依据其它资料给出合适的初始值 (胡亚轩等 , 2003 )。 应用
中图分类号 :P 317
文献标识码 :A
文章编号 :0253 - 4967(2007)01 - 0144 - 08
0 引言
M og i模型是模拟研究火山地表形变的最简单适用的模型 。该模型已用来模拟和解释大量 的火山区 地 面变 形 (Yokoyam a, 1971;Decke r et a l. , 1983;Berrino et a l. , 1984;Yang et a l. , 19 92)。 应 用火 山区 的形 变测 量数 据来 反 演 岩 浆 压 力源 参 数 对 火 山 的 危 险 性 评 价 有 很 大 的 意 义 。反演方法很多 , 经典的最小二乘方法在很多模型中有很好的应用 , 但是该模型有一定的局 限性 。 如果所求函数在某一初始值附近作泰勒展开形成的偏导矩阵为 A , 则当 ATA 是奇异矩 阵或可近似看作奇异矩阵时出现发散困难 。阻尼最小二乘方法可用来解决此难题 , 该方法已在 其它模型中得到了成功的应用 (王庆良等 , 1998)。
行觉等 , 2004)
Uz
= (r2
R3d +d2
)3 2
Ux
=R3 (x (r2 +d
x0 )
) 2
3 2
(1) (2)
Uy
=R3 (y (r2 +d
y0 )
) 2
3 2
(3)
式 (1) ～ (3)中 , Uz 为 z分 量位 移 , Ux
为 x分量位移 , Uy 为 y 分量位移 。 r = (x - x0 )2 +(y - y0 )2 为 径 向 半 径 , d
2 M og i模型的阻尼最小二乘原理
假定某一模型函数可表示为
U =U (X , B )
(4)
式 (4)中 , U 是参数 B 一般形式的函数 , B 可以是单个变量 , 也可以是 m 个变量 , 即 B =(b1 , b2, … , bm );X 可以是单个变量 , 也可以是 p个变量 , 即 X =(x1, x2, …, xp )。 对于 M og i模型 , 当用形 变资料反演压力源参数时 , 此处的 U 函数即为 (1)、(2)和 (3)式 。 X 表示观测点的坐标 , B 表示
S =εT ε=(ΔU - AΔb)T (ΔU - AΔb) =m in
(12)
取极小值 , 以获得模型参数的改正向量 Δb 的估计值 Δb^。式 (12)中 , T 表示转置 。按最小二乘
准则得到下列的正则方程 :
ATA Δb^=AT ΔU
(13)
源自文库此时 , 如果 ATA 是非奇异矩阵 , 则可得 Δb^:
压力源参数 , U 表示对应的形变观测量 。在实际问题中 , 绝大多数情况下模型参数与观测数据
都构成非线性关系
。
用
U
1 i
表示第
i个观测值 , 不失一般性 , 设以
Ui
=U i (X , B )表示相应的理论
计算值 , 则残差为
ε1i
=U
1 i
- Ui ,
i
=1,
2,
…,
n
(5)
式 (5)中 , U i 是参数 bj 的函数 , 将其在某一初始值 b0j 处作泰勒展开 , 然后略去 2次和 2次以上
14 6
地 震 地 质
29 卷
则 εi 可以解释为观测的残差减去理论的残差 : εi =ΔU i - A ijΔbj
用矩阵符号 , 可以把上式表示为
(10 )
ε=ΔU - A Δb
(11 )
若观测方程的个数 n 大于模型参数 bj 的个数 m , 则方程组 (8)便是不相容方程组 。 在这种 情况下 , 可以按最小二乘准则求解 , 这就是使残差平方和
Δb^=(A TA )-1 A T ΔU
(14)
如果 ATA 是奇异矩阵或可近似看作奇异矩阵则出现发散困难 。
应用 M ogi模型 , 对于岩浆等效体积变化较小的火山区来说 , 如果采用经典的最小二乘反
演 , 由于对压力源各参数偏导形成的矩阵 A 各量较小 , ATA 是奇异矩阵或可近似看作奇异矩
阵 。此处采用阻尼最小二乘法解决 , 则 Δb^为
和水平位移的联 合反演 , 数据为 2002— 2005年共 4 期的水准 和 G PS观测 资料 , 最后结 合其它 观测资
料 , 分析近年来岩浆可能的活动特征 。结果显 示岩浆 的位置在 发生变 化 , 体积增 量逐年 减小 , 表明火
山岩浆在 2002— 2005年的活动逐年减弱 。
关键词 阻尼 最小二乘法 联合反演 长白山天池火山 岩浆活动特征
摘 要 文中首先介绍 了 M og i模型的阻尼最小二乘反 演 , 得出应用形变资料 联合反演压 力源参数
的公式 , 再指出用等效源半径反演时不同初始值对反演结果的影响 , 说明选择合适的初始值 可以方便
地得到可靠的反 演结果 。 然后用此方法反演长白山天池火山区压力源的大小和位置参 数 。 采用垂直
Δb^=(A TA +θ2 I)-1 A T ΔU
(15)
式 (15)中 , θ2 为阻尼系数 , I为单位矩阵 。
进行反演计算时 , 给定各参数的初值 b0, 求得新值 b =b0 +Δb^。 然后以新值作为下次的初
值 , 逐次叠代 , 直至参数满足预先给定的收敛准则 (王庆良等 , 1998 ;陈运泰等 , 1979 )。
图 2 等效半径的初始值与反演叠代次数的关系 F ig. 2 Re lation of in itia lR va lue and itera tive num be r.
第 29卷 第 1期 2007年 3月
地 震 地 质
SE ISMOLOGY AND GEOLOGY
V o .l 29, N o. 1 M ar., 2007
Mog i模型在长白山天池火山区的应用
胡亚轩 1) 王庆良1) 崔笃信 1) 王文萍1) 李 克 2) 陈红卫 1)
1)中国地震局第二监测中心 , 西安 710054 2)吉林省地震局 , 长春 130022
高次项
式 (6)中 ,
∑m
ε1i ≈ εi
=U
1 i
- [ U0i (X , B ) +
j=1
Ui bj
Δbj ]
b0j
(6)
Δbj =bj - b0j
(7)
若以 ΔU i 表示观测的残差 :
ΔUi
=U
1 i
-
U
0 i
(X
,
B
)
(8)
以 A ij表示偏导数矩阵 :
Aij =
Ui bj b0j
(9 )
1 M og i模型
M og i模型是将岩浆压力源置于均匀弹性半空间中 , 在源的半径远小于源的深度时 , 可得出 由爆炸性或收缩性压力源引起的弹性半空间地表的水平变形和垂直变形与源的参数之间的表
〔收稿日期 〕 2006 - 07 - 05收稿 , 2006 - 11 - 10改回 。 〔基金项目 〕 国家自然科学基金 (40574041)和地震科学联合基金 (605004)共同资助 。
最小二乘反演叠代终止的条件为相邻 2 次叠代得到的残差平方和 S(k)和 S(k +1 )满足不等式
r1
=
S (k) - S (k +1) S (k)
≤γ
(17)
γ为 1 个事先给定的 <1 的正数 。 反演时 , 如果直接采用 RS =60m 作为反演需达到的理论值 , 可以看到需要叠代相当多的次
i
i
m
∑ ΔU x(x, y , z;x0 , y0 , d0 , R )] 2 + [ ΔUyi - ΔUy (x, y , z;x0 , y0 , d0 , R )] 2 i
(16 )
式 (16 )中 , n为地面垂直形变观测点数 , m 为地面水平形变观测点数 , x, y 为地面点的水平坐标 ,
M og i模型对不同方向的位移变化响应有一定的局限性 。 D ie te rich等 (1975 )应用有限元模
型研究得出地表垂直位移对压力源的中心位置比较敏感 , 而水平位移对压力源的形状和深度有
所反映 。故分析火山区压力源参数时 , 水平变形和垂直变形两者都是必须的 。下面我们分析如
何应用垂直和水平形变数据联合反演火山区岩浆压力源的参数 。
为岩浆源中心深度 。
图 1 M ogi模型 坐标示意图 F ig. 1 Coord inate ske tch m ap o fM ogi mode.l
若岩浆压力源中心坐标为 (0, 0, - d), 依据式 (1) ～ (3), 当 r =0时 , 有最大垂直位移量 ;当
r =d / 2时 , 有最大水平位移量 , 且水平位移呈对称放射状分布 。
1期
胡亚轩等 :M og i模型在长白山天池火山区的应用
145
达式 。 假定地面直角坐标系为 (x, y, z), 岩浆压力源中心坐标为 (x0, y0 , - d )(图 1), 如果用等 效半径 R 表示岩浆流入或流出的体积 , 则根据 M og i模型 , 岩浆房膨胀所引起的地表位移表达式
分别为 (王勇 , 1996;胡亚轩等 , 2003;施
从公式 (1 )、(2 )和 (3 )可以看出 , 用等效球状源的参数 R (等效半径 )及源的位置参数 (x0 , y0 , d0 )作为反演参数应用阻尼最小二乘反演方法联合反演时 , 符合的条件为几何形变实测值与 模型理论值的残差平方和最小 , 如果考虑点位高程的影响 , 此条件可以表示为
n
m
∑ ∑ S =Qm in = [ ΔU zi - ΔUz (x, y, z;x0 , y0 , d0 , R )] 2 + [ ΔU xi -
1期
胡亚轩等 :M og i模型在长白山天池火山区的应用
147
M og i模型 , 当等效压力源体积变化较小时 , 反演时会出 现困难 。 从 3 个位移公式 (1 )、(2 )和 (3 )可以看出 , 压力源变化的大小 (等效半径 R )对 3 个方向的位移影响较大 。 不同大小的 R 值 对反演过程也有一定的影响 。 这里假定压力源等效半径 RS 为 60m , 位置参数 xS , yS , dS 分别为 3 000m , 4 500m , 5 000m , 如果给出地面形变监测点的坐标 x, y, z(这里取 z =0 ), 依据 M og i模型 公式可以计算出位移值 , 相当于实际工作中的观测值 。反演时假定 1 组初始的压力源参数 x0 , y0 , d0 , R0 分别为 300m、800m 、2 000m 和 300m , 计算得出的位移值为理论初始值 。应用以上方 法进行反演 。
2002年 6月汪清地震后 , 地震 、地化和形变等资料均表明长白山天池火山区经历了一次比 较强烈的岩浆上涌活动过程 , 而且火山变形明显地出现在天池旧火山口附近 , 水平位移接近放 射状 , 与 M og i模型的理论位移形状相符 (M og i, 1958;Yang et a l. , 1992)。由于岩浆的相对变化 较小 , ATA 是奇异矩阵 。 本文在理论证明阻尼最小二乘方法可行性的前提下 , 应用此方法反演 了长白山天池火山区压力源的位置和等效半径等参数 , 分析了近年来岩浆可能的活动特征 。
z为观测点的高程 , d 为压力源深度 。ΔU zi 、ΔU xi 和 ΔU yi 分别为实测垂直形变量 , 实测 X 方向的 水平分量及 Y 方向的水平分量 ;ΔUz 、ΔUx 和 ΔUy 分别为 M og i模型计算的各方向理论值 。
3 阻尼最小二乘反演初始值的确定
进行阻尼最小二乘反演时 , 一般依据其它资料给出合适的初始值 (胡亚轩等 , 2003 )。 应用
中图分类号 :P 317
文献标识码 :A
文章编号 :0253 - 4967(2007)01 - 0144 - 08
0 引言
M og i模型是模拟研究火山地表形变的最简单适用的模型 。该模型已用来模拟和解释大量 的火山区 地 面变 形 (Yokoyam a, 1971;Decke r et a l. , 1983;Berrino et a l. , 1984;Yang et a l. , 19 92)。 应 用火 山区 的形 变测 量数 据来 反 演 岩 浆 压 力源 参 数 对 火 山 的 危 险 性 评 价 有 很 大 的 意 义 。反演方法很多 , 经典的最小二乘方法在很多模型中有很好的应用 , 但是该模型有一定的局 限性 。 如果所求函数在某一初始值附近作泰勒展开形成的偏导矩阵为 A , 则当 ATA 是奇异矩 阵或可近似看作奇异矩阵时出现发散困难 。阻尼最小二乘方法可用来解决此难题 , 该方法已在 其它模型中得到了成功的应用 (王庆良等 , 1998)。
行觉等 , 2004)
Uz
= (r2
R3d +d2
)3 2
Ux
=R3 (x (r2 +d
x0 )
) 2
3 2
(1) (2)
Uy
=R3 (y (r2 +d
y0 )
) 2
3 2
(3)
式 (1) ～ (3)中 , Uz 为 z分 量位 移 , Ux
为 x分量位移 , Uy 为 y 分量位移 。 r = (x - x0 )2 +(y - y0 )2 为 径 向 半 径 , d
2 M og i模型的阻尼最小二乘原理
假定某一模型函数可表示为
U =U (X , B )
(4)
式 (4)中 , U 是参数 B 一般形式的函数 , B 可以是单个变量 , 也可以是 m 个变量 , 即 B =(b1 , b2, … , bm );X 可以是单个变量 , 也可以是 p个变量 , 即 X =(x1, x2, …, xp )。 对于 M og i模型 , 当用形 变资料反演压力源参数时 , 此处的 U 函数即为 (1)、(2)和 (3)式 。 X 表示观测点的坐标 , B 表示
S =εT ε=(ΔU - AΔb)T (ΔU - AΔb) =m in
(12)
取极小值 , 以获得模型参数的改正向量 Δb 的估计值 Δb^。式 (12)中 , T 表示转置 。按最小二乘
准则得到下列的正则方程 :
ATA Δb^=AT ΔU
(13)
源自文库此时 , 如果 ATA 是非奇异矩阵 , 则可得 Δb^:
压力源参数 , U 表示对应的形变观测量 。在实际问题中 , 绝大多数情况下模型参数与观测数据
都构成非线性关系
。
用
U
1 i
表示第
i个观测值 , 不失一般性 , 设以
Ui
=U i (X , B )表示相应的理论
计算值 , 则残差为
ε1i
=U
1 i
- Ui ,
i
=1,
2,
…,
n
(5)
式 (5)中 , U i 是参数 bj 的函数 , 将其在某一初始值 b0j 处作泰勒展开 , 然后略去 2次和 2次以上
14 6
地 震 地 质
29 卷
则 εi 可以解释为观测的残差减去理论的残差 : εi =ΔU i - A ijΔbj
用矩阵符号 , 可以把上式表示为
(10 )
ε=ΔU - A Δb
(11 )
若观测方程的个数 n 大于模型参数 bj 的个数 m , 则方程组 (8)便是不相容方程组 。 在这种 情况下 , 可以按最小二乘准则求解 , 这就是使残差平方和
Δb^=(A TA )-1 A T ΔU
(14)
如果 ATA 是奇异矩阵或可近似看作奇异矩阵则出现发散困难 。
应用 M ogi模型 , 对于岩浆等效体积变化较小的火山区来说 , 如果采用经典的最小二乘反
演 , 由于对压力源各参数偏导形成的矩阵 A 各量较小 , ATA 是奇异矩阵或可近似看作奇异矩
阵 。此处采用阻尼最小二乘法解决 , 则 Δb^为
和水平位移的联 合反演 , 数据为 2002— 2005年共 4 期的水准 和 G PS观测 资料 , 最后结 合其它 观测资
料 , 分析近年来岩浆可能的活动特征 。结果显 示岩浆 的位置在 发生变 化 , 体积增 量逐年 减小 , 表明火
山岩浆在 2002— 2005年的活动逐年减弱 。
关键词 阻尼 最小二乘法 联合反演 长白山天池火山 岩浆活动特征
摘 要 文中首先介绍 了 M og i模型的阻尼最小二乘反 演 , 得出应用形变资料 联合反演压 力源参数
的公式 , 再指出用等效源半径反演时不同初始值对反演结果的影响 , 说明选择合适的初始值 可以方便
地得到可靠的反 演结果 。 然后用此方法反演长白山天池火山区压力源的大小和位置参 数 。 采用垂直
Δb^=(A TA +θ2 I)-1 A T ΔU
(15)
式 (15)中 , θ2 为阻尼系数 , I为单位矩阵 。
进行反演计算时 , 给定各参数的初值 b0, 求得新值 b =b0 +Δb^。 然后以新值作为下次的初
值 , 逐次叠代 , 直至参数满足预先给定的收敛准则 (王庆良等 , 1998 ;陈运泰等 , 1979 )。
图 2 等效半径的初始值与反演叠代次数的关系 F ig. 2 Re lation of in itia lR va lue and itera tive num be r.
第 29卷 第 1期 2007年 3月
地 震 地 质
SE ISMOLOGY AND GEOLOGY
V o .l 29, N o. 1 M ar., 2007
Mog i模型在长白山天池火山区的应用
胡亚轩 1) 王庆良1) 崔笃信 1) 王文萍1) 李 克 2) 陈红卫 1)
1)中国地震局第二监测中心 , 西安 710054 2)吉林省地震局 , 长春 130022
高次项
式 (6)中 ,
∑m
ε1i ≈ εi
=U
1 i
- [ U0i (X , B ) +
j=1
Ui bj
Δbj ]
b0j
(6)
Δbj =bj - b0j
(7)
若以 ΔU i 表示观测的残差 :
ΔUi
=U
1 i
-
U
0 i
(X
,
B
)
(8)
以 A ij表示偏导数矩阵 :
Aij =
Ui bj b0j
(9 )
1 M og i模型
M og i模型是将岩浆压力源置于均匀弹性半空间中 , 在源的半径远小于源的深度时 , 可得出 由爆炸性或收缩性压力源引起的弹性半空间地表的水平变形和垂直变形与源的参数之间的表
〔收稿日期 〕 2006 - 07 - 05收稿 , 2006 - 11 - 10改回 。 〔基金项目 〕 国家自然科学基金 (40574041)和地震科学联合基金 (605004)共同资助 。
最小二乘反演叠代终止的条件为相邻 2 次叠代得到的残差平方和 S(k)和 S(k +1 )满足不等式
r1
=
S (k) - S (k +1) S (k)
≤γ
(17)
γ为 1 个事先给定的 <1 的正数 。 反演时 , 如果直接采用 RS =60m 作为反演需达到的理论值 , 可以看到需要叠代相当多的次
i
i
m
∑ ΔU x(x, y , z;x0 , y0 , d0 , R )] 2 + [ ΔUyi - ΔUy (x, y , z;x0 , y0 , d0 , R )] 2 i
(16 )
式 (16 )中 , n为地面垂直形变观测点数 , m 为地面水平形变观测点数 , x, y 为地面点的水平坐标 ,