地震层析成像LSQR算法的并行化_刘劲松
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
[ 3 , 4]
式中 x 为模型网格点上的速度扰动量 , A 为相应的 速度偏导数矩阵 , y 为研究区域内地震震源参数的 扰动量 , B 为相应的震源参数偏导数矩阵 , δ t 为走 时残差 . 下标表示地震序号 , L r 为区内地震总数 . 引 入正交投影算子 P B : PB = BB , 可将( 1) 式分解为两个方程分别求解 : ( I - P B) Ax =( I - PB ) δ t, By = P B ( δ t - Ax ) . ( 2) 式可写为下列形式 : A′ x =δ t′ ,
Abstract We discuss the LSQR algorithms used in earthquake travel time tomography . We keep the epicenter terms in the equation for regional events , and then use the orthogonal projection method to eliminate the epicenter terms . For tele-events , the classic smoothing process is used . The number of non-zero elements in the partial derivative matrix is increased by several times because of the orthogonal projection and smoothing processes . For a large scale inversion problem , the amount of non -zero elements can be dozens of Gigabytes or hundreds of Gigabytes . The huge amount of memory requirement becomes the bottle neck of LSQR algorithms . To solve this problem , we have studied the distribution property of non-zero elements in the partial derivative matrix , designed an efficient data structure for the sparse matrix , used a distributed memory and computation scheme for matrix computation , and implemented it on a multi-processor super-computer . We have derived an estimation formula of parallel efficiency and tested two real tomography models . Keywords Seismic tomography , Traveltime inversion , LSQR algorithms , Parallel computation
2
x0 = 0 , ρ 1 = 对角化 :
a 1 +λ , φ 1 =α 1 β1 ρ 1 . ( 8)
2
2
βi +1 ui +1 = Cv i -α i ui , α i+ 1 vi + 1 = C ui + 1 -βi +1 vi 正交变换 : ρ i = ρ si = βi +1 ρ i +βi +1 , fi = ρ i ρ i , i , ( 10) θ ρ i+ 1 = si α i+ 1 , i+ 1 =- fi α i+ 1 φ i = fi φ i , φ i+ 1 = si φ i . ( 11) 解矢量 : xi = x i -1 + ω iφ i ρ i , ω i+ 1 = vi+ 1 -ω iθ i+ 1 ρ i . ( 12)
地震层析成像 LSQR 算法的并行化
刘劲松 , 刘福田 , 刘 俊 , 郝天珧
1 中国科学院地质与地球物理研究所 , 北京 100029 2 中国科学院计算机网络信息中心超级计算中心 , 北京 100080 1 1 2 1
摘 要 讨论了地震层析成像的 LSQR 算法( 最小二乘 QR 分解). 在建立偏导数矩阵方程组时 , 对区内地震在方 程 中保留震源项 , 引入正交投影算子进行参数分离 , 对区外远震 采用传统 的平滑处 理方式 , 用 LSQR 法求解 联立的 方 程组 . 由于区内地震的正交分解处理和区外远震的平滑处理 , 使得偏导数矩阵中的非零元素成倍增加 , 对于大型反 演问题 , 这些非零元素常常达到几十 GB 到几百 GB 的数量级 , 巨量的内存占用成为 LSQR 算法的瓶颈 .针对这一问 题 , 本文研究了偏导数矩阵中非零元素的分布规律 , 设计 出合理 的存储 结构 , 采用分 布式存 储进行 矩阵计 算 , 提 出 了 LSQR 算法的并行化方案 , 并在联想深腾 6800 超级计算机上实现 .导出了 LSQR 算法的并行效率估算公式 .对两 个地区的实际地震层析成像数据进行了效率测试 . 关键词 地震层析成像 , 走时反演 , LSQR 算法 , 并行计算 文章编号 0001 -5733 ( 2006) 02 -0540 -06 中图分类号 P631 收稿日期 2005 -04 -06 , 2005 -12 -26 收修定稿
基金项目 中国科学院知识创新工程信息化建设重大专项( INF10 5 -SCE) 和国家自然科学基金( 40034010) 资助 . 作者简介 刘劲松 , 男 , 1965 年生 , 1988 年毕业于厦门大学物理系 , 1997 年在中 国地质大学 ( 北京) 获 博士学位 , 主要从事 地震层析 成像及深 部 地球物理研究 . E mail : jsliu @mai l . iggcas. ac . cn
第 49 卷 第 2 期 2006 年 3 月
地 球 物 理 学 报
CHINESE JOURNAL OF GEOPHYSICS
Vol . 49 , No . 2 Mar . , 2006
刘劲松 , 刘福田 , 刘 俊等 . 地震层析成像 LSQR 算法的 并行化 . 地球物理学报 , 2006 , 49( 2): 540 ～ 545 Liu J S , Liu F T , Liu J , et al . Parallel LSQR algorithms used in seismic tomography Chinese J .Geophys . ( in Chinese) , 2006 , 49( 2) : 540 ～ 545
[ 3] [ 8, 10]
( 2) ( 3) ( 4)
;
1984 年刘福田 引入正交投影算子 , 将震源参数与 速度参数分离 , 提出顺序正交三角化算法 , 减少了内 存占用量和 计算量 , 使得 最终的 方程 可用 SVD 法 求解 . 正交投影参数分离方法在原始的走时方程中保 留震源项 , 求解速度参数前引入正交投影算子 , 消去 震源项 , 得到只包含速度项的方程 . 由于在方程中 考虑了震源项 , 比重新定位后忽略震源项的方法更 接近实际情况 . 但是经过正交分解处理的矩阵方程 组 , 矩阵的稀疏性被降低 , 非零元素成倍增加 . 在用 LSQR 方法求解大型地震层析成像问题时 , 巨量的内 存占用成为瓶颈 . 本文针对这一问题 , 研究了偏导 数矩阵中非零元素的分布规律 , 采用分布式存储进 行矩阵计算 , 提出了 LSQR 算法的并行化方案 .
T T 1 T 2 T L T +
在模 型中 引入了 间断 面用 以描 述莫 霍面 和
[ 5 ～ 7] [ 8]
410 、 660km 等速度间断面 , 这一方法被广泛 应用于 国内的地震层析成像工作 ; 1992 年 Zhao et al 改进了射线追踪方法 , 使其适用于模型中的非水平 间断面 , 进一步接近了实际情况 . 在反演方法方面 , Aki 等最初采用阻尼最小二乘法 ; 随后各种迭代类 算法被用于地震层析成像反演 , 如 ART 方法 , SIRT [ 9] 方法等 ; 1982 年 Paige and Sanders 提出 LSQR 方法 后 , 该方 法被广泛 用于地震 层析成像 反演中
[ 11, 12]
2 正交投影地震层析成像反演法简介
用于速度图像重建的层析成像法
[ 3 , 4]
, 可归结 ( 1)
为求解一个矩阵方程组 . 对于区内地震有 Ax +By =δ t, 其中 A1 A2 A = B = AL r B1 B2 · · BL r ,
. 作为一种可选方案 ,
本文采 用 LSQR 方 法 求解 方 程 ( 5) . 与 其 他采 用 LSQR 法的地震层析成像方法不同的是 , 本文在区内 地震的方程中考虑震源项 , 引入正交投影算子进行 正交分解参数分离 , 消去震源项后求解 . 与 SVD 法相比 , LSQR 法具有计算量小的优点 . 对于中小规模的地震层析成像问题 , 与充分优化的 顺序正 交三角化 SVD 方法相比
[ 11 , 12]
, LBiblioteka BaiduQR 法所需
的内存量一般情况下均大于 SVD 法 . 设未知数和射
542
地 球 物 理 学 报( Chinese J . Geophys . )
49 卷
线数 分 别为 n , m , 用 SVD 法所 需 的 内存 量 约 为 4 n , LSQR 法所需的内存量约为 4 mns , s 为非 零元 素所占百分比 . 例如对于中国大陆地震成像数据 , 若采用 1° × 1° 的模型网格间距 , 方程未知数个数将 超过 10 万 ,( 5) 式中的矩阵 C 数据即使只存放非零 元素也将达到 60GB ; 如果采用 0. 5 ° × 0. 5° 的网格间 距 , 矩阵 C 数据量将超过 200GB , 采用常规 的办法 难以求 解这样 大型反 演问 题 , 必须采 用并 行计 算 方法 .
x =( x 1 , x 2 , …, xn ) , y =( y , y , … , y r). 对于远 震数 据 , 采 用平 滑处理 后可 消去 震源 项
[ 1]
, 直接得到形如( 4) 式的方程 . 区内地震和远震 Cx = b ( 5)
的方程联立 , 得到 : 矩阵 C 可写为 C1 C2 C = CL Ci 为经过正交分解参数分离或平滑处理后的每个 地震的偏导数矩阵 , L 为地震总数 . 文献[ 3] 采用先对方程 ( 5) 进行顺序正交三角 化 , 再进行 SVD 分解求解方程的办法 . 这一方法具 有节省内存 、 数值稳定 、 能够得到解的分辨和方差等 优点 . 但对于大型问题 , 即使采用并行计算 , 也会有 计算量大 、 耗时长的问题
2 期
刘劲松等 : 地震层析成像 LSQR 算法的并行化
541
δ t1
1 引 言
δ t = 自 1977 年 Aki 等的地震成像的开创性 工作以 来
[ 1]
δ t2 , δ tL r
, 体波地震成像方法得到了长 足的发展 . 在模
[ 2]
型的参数化方面 , 1982 年朱天飞和刘福田 采用网 格点插值模型取代了最初的分块常速模型 ; 刘福田 等
Parallel LSQR algorithms used in seismic tomography
LIU JinSong 1 , LIU FuTian1 , LIU Jun2 , HAO Tian -Yao 1
1 Institute of G eology and G eophysics , Chinese Academy of Sciences , Beijing 100029 , China 2 Supercomputing Center of Computer Network Information Center , Chinese Academy of Sciences , Beijing 100080 , China
式中 x 为模型网格点上的速度扰动量 , A 为相应的 速度偏导数矩阵 , y 为研究区域内地震震源参数的 扰动量 , B 为相应的震源参数偏导数矩阵 , δ t 为走 时残差 . 下标表示地震序号 , L r 为区内地震总数 . 引 入正交投影算子 P B : PB = BB , 可将( 1) 式分解为两个方程分别求解 : ( I - P B) Ax =( I - PB ) δ t, By = P B ( δ t - Ax ) . ( 2) 式可写为下列形式 : A′ x =δ t′ ,
Abstract We discuss the LSQR algorithms used in earthquake travel time tomography . We keep the epicenter terms in the equation for regional events , and then use the orthogonal projection method to eliminate the epicenter terms . For tele-events , the classic smoothing process is used . The number of non-zero elements in the partial derivative matrix is increased by several times because of the orthogonal projection and smoothing processes . For a large scale inversion problem , the amount of non -zero elements can be dozens of Gigabytes or hundreds of Gigabytes . The huge amount of memory requirement becomes the bottle neck of LSQR algorithms . To solve this problem , we have studied the distribution property of non-zero elements in the partial derivative matrix , designed an efficient data structure for the sparse matrix , used a distributed memory and computation scheme for matrix computation , and implemented it on a multi-processor super-computer . We have derived an estimation formula of parallel efficiency and tested two real tomography models . Keywords Seismic tomography , Traveltime inversion , LSQR algorithms , Parallel computation
2
x0 = 0 , ρ 1 = 对角化 :
a 1 +λ , φ 1 =α 1 β1 ρ 1 . ( 8)
2
2
βi +1 ui +1 = Cv i -α i ui , α i+ 1 vi + 1 = C ui + 1 -βi +1 vi 正交变换 : ρ i = ρ si = βi +1 ρ i +βi +1 , fi = ρ i ρ i , i , ( 10) θ ρ i+ 1 = si α i+ 1 , i+ 1 =- fi α i+ 1 φ i = fi φ i , φ i+ 1 = si φ i . ( 11) 解矢量 : xi = x i -1 + ω iφ i ρ i , ω i+ 1 = vi+ 1 -ω iθ i+ 1 ρ i . ( 12)
地震层析成像 LSQR 算法的并行化
刘劲松 , 刘福田 , 刘 俊 , 郝天珧
1 中国科学院地质与地球物理研究所 , 北京 100029 2 中国科学院计算机网络信息中心超级计算中心 , 北京 100080 1 1 2 1
摘 要 讨论了地震层析成像的 LSQR 算法( 最小二乘 QR 分解). 在建立偏导数矩阵方程组时 , 对区内地震在方 程 中保留震源项 , 引入正交投影算子进行参数分离 , 对区外远震 采用传统 的平滑处 理方式 , 用 LSQR 法求解 联立的 方 程组 . 由于区内地震的正交分解处理和区外远震的平滑处理 , 使得偏导数矩阵中的非零元素成倍增加 , 对于大型反 演问题 , 这些非零元素常常达到几十 GB 到几百 GB 的数量级 , 巨量的内存占用成为 LSQR 算法的瓶颈 .针对这一问 题 , 本文研究了偏导数矩阵中非零元素的分布规律 , 设计 出合理 的存储 结构 , 采用分 布式存 储进行 矩阵计 算 , 提 出 了 LSQR 算法的并行化方案 , 并在联想深腾 6800 超级计算机上实现 .导出了 LSQR 算法的并行效率估算公式 .对两 个地区的实际地震层析成像数据进行了效率测试 . 关键词 地震层析成像 , 走时反演 , LSQR 算法 , 并行计算 文章编号 0001 -5733 ( 2006) 02 -0540 -06 中图分类号 P631 收稿日期 2005 -04 -06 , 2005 -12 -26 收修定稿
基金项目 中国科学院知识创新工程信息化建设重大专项( INF10 5 -SCE) 和国家自然科学基金( 40034010) 资助 . 作者简介 刘劲松 , 男 , 1965 年生 , 1988 年毕业于厦门大学物理系 , 1997 年在中 国地质大学 ( 北京) 获 博士学位 , 主要从事 地震层析 成像及深 部 地球物理研究 . E mail : jsliu @mai l . iggcas. ac . cn
第 49 卷 第 2 期 2006 年 3 月
地 球 物 理 学 报
CHINESE JOURNAL OF GEOPHYSICS
Vol . 49 , No . 2 Mar . , 2006
刘劲松 , 刘福田 , 刘 俊等 . 地震层析成像 LSQR 算法的 并行化 . 地球物理学报 , 2006 , 49( 2): 540 ～ 545 Liu J S , Liu F T , Liu J , et al . Parallel LSQR algorithms used in seismic tomography Chinese J .Geophys . ( in Chinese) , 2006 , 49( 2) : 540 ～ 545
[ 3] [ 8, 10]
( 2) ( 3) ( 4)
;
1984 年刘福田 引入正交投影算子 , 将震源参数与 速度参数分离 , 提出顺序正交三角化算法 , 减少了内 存占用量和 计算量 , 使得 最终的 方程 可用 SVD 法 求解 . 正交投影参数分离方法在原始的走时方程中保 留震源项 , 求解速度参数前引入正交投影算子 , 消去 震源项 , 得到只包含速度项的方程 . 由于在方程中 考虑了震源项 , 比重新定位后忽略震源项的方法更 接近实际情况 . 但是经过正交分解处理的矩阵方程 组 , 矩阵的稀疏性被降低 , 非零元素成倍增加 . 在用 LSQR 方法求解大型地震层析成像问题时 , 巨量的内 存占用成为瓶颈 . 本文针对这一问题 , 研究了偏导 数矩阵中非零元素的分布规律 , 采用分布式存储进 行矩阵计算 , 提出了 LSQR 算法的并行化方案 .
T T 1 T 2 T L T +
在模 型中 引入了 间断 面用 以描 述莫 霍面 和
[ 5 ～ 7] [ 8]
410 、 660km 等速度间断面 , 这一方法被广泛 应用于 国内的地震层析成像工作 ; 1992 年 Zhao et al 改进了射线追踪方法 , 使其适用于模型中的非水平 间断面 , 进一步接近了实际情况 . 在反演方法方面 , Aki 等最初采用阻尼最小二乘法 ; 随后各种迭代类 算法被用于地震层析成像反演 , 如 ART 方法 , SIRT [ 9] 方法等 ; 1982 年 Paige and Sanders 提出 LSQR 方法 后 , 该方 法被广泛 用于地震 层析成像 反演中
[ 11, 12]
2 正交投影地震层析成像反演法简介
用于速度图像重建的层析成像法
[ 3 , 4]
, 可归结 ( 1)
为求解一个矩阵方程组 . 对于区内地震有 Ax +By =δ t, 其中 A1 A2 A = B = AL r B1 B2 · · BL r ,
. 作为一种可选方案 ,
本文采 用 LSQR 方 法 求解 方 程 ( 5) . 与 其 他采 用 LSQR 法的地震层析成像方法不同的是 , 本文在区内 地震的方程中考虑震源项 , 引入正交投影算子进行 正交分解参数分离 , 消去震源项后求解 . 与 SVD 法相比 , LSQR 法具有计算量小的优点 . 对于中小规模的地震层析成像问题 , 与充分优化的 顺序正 交三角化 SVD 方法相比
[ 11 , 12]
, LBiblioteka BaiduQR 法所需
的内存量一般情况下均大于 SVD 法 . 设未知数和射
542
地 球 物 理 学 报( Chinese J . Geophys . )
49 卷
线数 分 别为 n , m , 用 SVD 法所 需 的 内存 量 约 为 4 n , LSQR 法所需的内存量约为 4 mns , s 为非 零元 素所占百分比 . 例如对于中国大陆地震成像数据 , 若采用 1° × 1° 的模型网格间距 , 方程未知数个数将 超过 10 万 ,( 5) 式中的矩阵 C 数据即使只存放非零 元素也将达到 60GB ; 如果采用 0. 5 ° × 0. 5° 的网格间 距 , 矩阵 C 数据量将超过 200GB , 采用常规 的办法 难以求 解这样 大型反 演问 题 , 必须采 用并 行计 算 方法 .
x =( x 1 , x 2 , …, xn ) , y =( y , y , … , y r). 对于远 震数 据 , 采 用平 滑处理 后可 消去 震源 项
[ 1]
, 直接得到形如( 4) 式的方程 . 区内地震和远震 Cx = b ( 5)
的方程联立 , 得到 : 矩阵 C 可写为 C1 C2 C = CL Ci 为经过正交分解参数分离或平滑处理后的每个 地震的偏导数矩阵 , L 为地震总数 . 文献[ 3] 采用先对方程 ( 5) 进行顺序正交三角 化 , 再进行 SVD 分解求解方程的办法 . 这一方法具 有节省内存 、 数值稳定 、 能够得到解的分辨和方差等 优点 . 但对于大型问题 , 即使采用并行计算 , 也会有 计算量大 、 耗时长的问题
2 期
刘劲松等 : 地震层析成像 LSQR 算法的并行化
541
δ t1
1 引 言
δ t = 自 1977 年 Aki 等的地震成像的开创性 工作以 来
[ 1]
δ t2 , δ tL r
, 体波地震成像方法得到了长 足的发展 . 在模
[ 2]
型的参数化方面 , 1982 年朱天飞和刘福田 采用网 格点插值模型取代了最初的分块常速模型 ; 刘福田 等
Parallel LSQR algorithms used in seismic tomography
LIU JinSong 1 , LIU FuTian1 , LIU Jun2 , HAO Tian -Yao 1
1 Institute of G eology and G eophysics , Chinese Academy of Sciences , Beijing 100029 , China 2 Supercomputing Center of Computer Network Information Center , Chinese Academy of Sciences , Beijing 100080 , China