电性参数分块连续的大地电磁二维有限元数值模拟
|
摘要
大地电磁测深法(MT)属于电磁法勘探中的频率域方法,利用天然交变场源对地球深部岩石的电性参数进行研究。目前,MT法的二维正演问题已基本解决,在使用有限单元法,有限差分法和积分方程法等数值模拟方法解决二维MT正演问题方面有广泛的研究结果,但在网格剖分方式,电性参数设定,辅助场定义和起伏地形模拟等方面,仍有改进之处。本文研究电性参数分块连续变化情况下的MT二维正演问题,推导出二维水平和起伏地形条件下大地电磁法有限元数值模拟算法,以提高正演模拟精度,为反演解释工作的进行提供良好的基础。
本文以频率域Maxwell方程组为基础,并结合大地电磁测深法的边界条件探讨,推导出二维MT法的变分方程。在有限单元法网格剖分方式上,采用矩形网格内剖分三角形网格的方案,这种剖分方式便于对倾斜异常体和倾斜地形的模拟,以适应各种水平或起伏地形情况。考虑到实际地层中的岩石,矿物体等在水平方向和垂直方向上电性参数是连续变化的,而在一些反演方法中,反演结果的电性参数也是连续变化,故将网格单元内的电性参数设定为线性变化。根据单元节点主场值和线性插值形函数间的关系,计算出单元节点辅助场值。在方程组的求解方面,采用变带宽存储解决含有大量零元素的大型稀疏矩阵的存储和方程组的求解问题,以节约内存使用量和提高计算速度。根据起伏地形情况下实测电磁场分量的特征,定义TE和TM两种模式下的视电阻率和阻抗相位计算公式。编制一套二维大地电磁法的正演程序,实现了使用高程信息自动生成网格和使用不同的模型参数计算指定频率的视电阻率和相位值。
通过对多种模型的验证,不同地电断面的水平地形和起伏地形的正演模拟结果与前人模拟结果一致,模型参数基本吻合。结合两种极化模式的视电阻率和相位信息的横向和纵向分辨率特点,对地下异常体的深度定位,规模大小和倾斜方位判定表现出良好的效果。山峰、山谷和斜坡等非水平地形情况下的正演响应结果也与正演模型基本符合,验证了本文方法的正确性和有效性。
Magnetotelluric sounding (MT) is one of the electromagnetic exploration methods which belong to frequency domain. The natural alternating field source is used to study electrical parameters of rocks in deep earth with this method. Currently, the two-dimensional forward problem in MT is basically solved. The use of finite element method, finite difference method and integral equation method and other numerical methods in solving the problem of MT two-dimensional forward has a widespread research. However there are still some improvements to study in the area of mesh generation method, electrical parameters set, the auxiliary field definition and topographic simulation, etc. In this dissertation, continuous electrical parameters with each block is involved in two-dimensional forward modeling, numerical simulation algorithms of magnetotelluric sounding in the condition of two-dimensional flat and topographic terrain are derived. All the studies above can help to improve the accuracy of forward simulation and make a good basis for work of the inversion interpretation.
In the beginning of this dissertation, the two-dimensional variational equation of MT method is derived from frequency domain Maxwell equations and the combination with discussion on boundary conditions. The rectangular grid within the subdivision triangle mesh is introduced in the FEM mesh generation method, which facilitates to simulate the sloping abnormal body and terrain to adapt the flat or topographic conditions. Accounting to the actual situation, that the rocks, minerals and other bodies in the horizontal direction and vertical direction usually have continuous variation of electrical parameters, and in some inversion methods, the inversion results of the electrical parameters have continuous variation, so the electrical parameters in each element is set to linear. The auxiliary field is calculated based on the relationship between the main field of the nodes and linear shape function in each element. In the aspect of equations solving, the variable bandwidth storage is used to solve storage problems of the large sparse matrix with large number of zero elements, which can help to save memory usage and improve the computing speed. The definition of apparent resistivity and impedance phase formula in TE and TM mode is based on the characteristics of the measured electromagnetic field components of flat and topographic condition. A set of two-dimensional MT forward procedure is created, which has the function of auto grid generating and using elevation information of the model parameters. It’s able to calculate apparent resistivity and impedance phase values with different specified frequency.
Through the verification on different geoelectric section models in the condition of flat and topographic terrains, it is proved that forward results are consistent with model parameters. Combining with the different resolution on horizontal and vertical direction of the two polarization modes of the apparent resistivity and phase characteristics, it shows a good result in the determination of underground abnormal body depth, size and sloping direction. The forward response of peaks, valleys and slopes and other topographic terrain are also consistent with their models, which verifies the correctness and validity of the method in this dissertation.
本文以频率域Maxwell方程组为基础,并结合大地电磁测深法的边界条件探讨,推导出二维MT法的变分方程。在有限单元法网格剖分方式上,采用矩形网格内剖分三角形网格的方案,这种剖分方式便于对倾斜异常体和倾斜地形的模拟,以适应各种水平或起伏地形情况。考虑到实际地层中的岩石,矿物体等在水平方向和垂直方向上电性参数是连续变化的,而在一些反演方法中,反演结果的电性参数也是连续变化,故将网格单元内的电性参数设定为线性变化。根据单元节点主场值和线性插值形函数间的关系,计算出单元节点辅助场值。在方程组的求解方面,采用变带宽存储解决含有大量零元素的大型稀疏矩阵的存储和方程组的求解问题,以节约内存使用量和提高计算速度。根据起伏地形情况下实测电磁场分量的特征,定义TE和TM两种模式下的视电阻率和阻抗相位计算公式。编制一套二维大地电磁法的正演程序,实现了使用高程信息自动生成网格和使用不同的模型参数计算指定频率的视电阻率和相位值。
通过对多种模型的验证,不同地电断面的水平地形和起伏地形的正演模拟结果与前人模拟结果一致,模型参数基本吻合。结合两种极化模式的视电阻率和相位信息的横向和纵向分辨率特点,对地下异常体的深度定位,规模大小和倾斜方位判定表现出良好的效果。山峰、山谷和斜坡等非水平地形情况下的正演响应结果也与正演模型基本符合,验证了本文方法的正确性和有效性。
Magnetotelluric sounding (MT) is one of the electromagnetic exploration methods which belong to frequency domain. The natural alternating field source is used to study electrical parameters of rocks in deep earth with this method. Currently, the two-dimensional forward problem in MT is basically solved. The use of finite element method, finite difference method and integral equation method and other numerical methods in solving the problem of MT two-dimensional forward has a widespread research. However there are still some improvements to study in the area of mesh generation method, electrical parameters set, the auxiliary field definition and topographic simulation, etc. In this dissertation, continuous electrical parameters with each block is involved in two-dimensional forward modeling, numerical simulation algorithms of magnetotelluric sounding in the condition of two-dimensional flat and topographic terrain are derived. All the studies above can help to improve the accuracy of forward simulation and make a good basis for work of the inversion interpretation.
In the beginning of this dissertation, the two-dimensional variational equation of MT method is derived from frequency domain Maxwell equations and the combination with discussion on boundary conditions. The rectangular grid within the subdivision triangle mesh is introduced in the FEM mesh generation method, which facilitates to simulate the sloping abnormal body and terrain to adapt the flat or topographic conditions. Accounting to the actual situation, that the rocks, minerals and other bodies in the horizontal direction and vertical direction usually have continuous variation of electrical parameters, and in some inversion methods, the inversion results of the electrical parameters have continuous variation, so the electrical parameters in each element is set to linear. The auxiliary field is calculated based on the relationship between the main field of the nodes and linear shape function in each element. In the aspect of equations solving, the variable bandwidth storage is used to solve storage problems of the large sparse matrix with large number of zero elements, which can help to save memory usage and improve the computing speed. The definition of apparent resistivity and impedance phase formula in TE and TM mode is based on the characteristics of the measured electromagnetic field components of flat and topographic condition. A set of two-dimensional MT forward procedure is created, which has the function of auto grid generating and using elevation information of the model parameters. It’s able to calculate apparent resistivity and impedance phase values with different specified frequency.
Through the verification on different geoelectric section models in the condition of flat and topographic terrains, it is proved that forward results are consistent with model parameters. Combining with the different resolution on horizontal and vertical direction of the two polarization modes of the apparent resistivity and phase characteristics, it shows a good result in the determination of underground abnormal body depth, size and sloping direction. The forward response of peaks, valleys and slopes and other topographic terrain are also consistent with their models, which verifies the correctness and validity of the method in this dissertation.
引文
1) Michael S. Zhdanov.Electromagnetic geophysics:Note from the past and the road ahead [J]. Geophysics,2010,75/5:75A49-75A66.
2) Philip E. Wannamaker,John A. Stodt,and Luis Rijo.Two-dimensional topographic responses in magnetotellurics modeled using finite elements[J]. Geophysics, 1986,51/11:2131-2144.
3) Micheal S. Zhdanov,George V. Keller.The Geoelectrical Methods in Geophysical Exploration (Methods in Geochemistry and Geophysics)[M].New York:Elsevier Science,1994.
4) Michael S. Zhdanov.Geophysical Electromagnetic Theory and Methods, Volume 43 (Methods in Geochemistry and Geophysics)[M](1).New York:Elsevier Science,2009.
5) A.A.考夫曼,G.V.凯勒,王建谋译.频率域与时间域电磁测深[M].北京:地质出版社,1987.
6)陈乐寿,王光锷.大地电磁测深法[M].北京:地质出版社,1990.
7)陈小斌.MT二维正演计算中地形影响的研究[J].石油物探,2000,39/3:112-120.
8)陈小斌.有限元直接迭代算法中一种基本结构的研究[J].物探化探计算技术,2000,22/2:132-136.
9)陈小斌,胡文宝.有限元直接迭代算法及其在线源频率域电磁响应计算中的应用[J].地球物理学报,2002,45/1:120-129.
10)陈小斌,张翔,胡文宝.有限元直接迭代算法在MT二维正演计算中的应用[J].石油地球物理勘探,2000,35/4:487-496.
11)程志平.电法勘探教程[M](1).北京:冶金工业出版社,2007.
12)刘德贵,费景高,于江永等.Fortran算法汇编(第一分册)[M].北京:国防工业出版社,1980.
13)刘云,王绪本.大地电磁二维自适应地形有限元正演模拟[J].地震地质,2010,32/3:382-391.
14)李予国,徐世浙,刘斌等.电导率分块连续变化的二维MT有限元模拟(II)[J].高校地质学报,1996,2/4:448-452.
15)马为.MT二维正演中辅助场计算新方法研究[D].北京:中国地震局,2007.
16)彭国伦.Fortran 95程序设计[M](第一版).北京:中国电力出版社,2002.
17)阮百尧.三角单元部分电导率分块连续变化点源二维电场有限元数值模拟[J].广西科学,2001,8/1:1-3.
18)阮百尧,熊彬.大型对称变带宽方程组的Cholesky分解法[J].物探化探计算技术,2000,22/4:361-363.
19)阮百尧,徐世浙.线元二维时间域电磁响应的有限差分解[J].高校地质学报,1996,2(4):437-447.
20)史明娟,徐世浙,刘斌.大地电磁二次函数插值的有限元法[J].地球物理学报,1997,40/3:421-429.
21)蔚宝强,胡文宝.MT资料的遗传法反演[J].江汉石油学院学报,1999,21(3):25-29.
22)王建谋等.大地电磁测深译文集.第一集[M](1).北京:地质出版社,1987.
23)熊彬,罗延钟.电导率分块均匀的瞬变电磁2.5维有限元数值模拟[J].地球物理学报,2006,49/2:590-597.
24)谢飞.带地形大地电磁场二维有限元数值模拟研究[D].北京:中国地质大学,2006.
25)徐世浙,于涛,李予国等.电导率分块连续变化的二维MT有限元模拟(I)[J].高校地质学报,1995,1/2:65-73.
26)徐士良.计算机常用算法[M](第二版).北京:清华大学出版社,1995.
27)徐世浙.地球物理中的有限单元法[M](第一版).北京:科学出版社,1994.
28)曾国.大地电磁二维有限元正演数值模拟[D].长沙:中南大学,2008.
29)周熙襄,钟本善.电法勘探数值模拟技术[M](第一版).成都:四川科学技术出版社,1986.
30)王建谋,陈乐寿等译.大地电磁测深译文集(第一集) [M](第一版).北京:地质出版社,1987.
31) Fortran. http://zh.wikipedia.org/wiki/fortran.
32)石应骏、刘国栋、王家映等.大地电磁测深法教程[M].成都:成都理工大学信息工程学院,1985.
33) Coggon J H.Electromagnetic and electrical modeling by the finite element method[J]. Geophysics, 1971, 36(1):132-155.
34) William L, Rodi. A technique for improving theaccuracy of finite element solution for magnetotelluric data[J].Geophysics, 1976,50 (7) :483-506.
35) Wannamaker P E, Stodt J A,Rijo L. A stable finite element solution for two-dimensional magnetotellnric modeling[J].Geophysics, 1987, 88 (5):277 -296.
36) Lugao P,Wannalnaker P E. Calculating the two-dimensional magnetotelluric Jacohian in finite elements using reciprocity[J].Geophysical Journal International, 1996, 127(2):806-810
37) Frank A, Ralph B, Klaus S. 2D finite element modeling of plane-wave diffusive time-harmonic electromagnetic fields using adaptive unstructured grids[A].Proceeding of the 17th Workshop,2004, S2:1-6.
38)简兴祥,大地电磁二三维建模及正演模拟[D].成都理工大学.2009
2) Philip E. Wannamaker,John A. Stodt,and Luis Rijo.Two-dimensional topographic responses in magnetotellurics modeled using finite elements[J]. Geophysics, 1986,51/11:2131-2144.
3) Micheal S. Zhdanov,George V. Keller.The Geoelectrical Methods in Geophysical Exploration (Methods in Geochemistry and Geophysics)[M].New York:Elsevier Science,1994.
4) Michael S. Zhdanov.Geophysical Electromagnetic Theory and Methods, Volume 43 (Methods in Geochemistry and Geophysics)[M](1).New York:Elsevier Science,2009.
5) A.A.考夫曼,G.V.凯勒,王建谋译.频率域与时间域电磁测深[M].北京:地质出版社,1987.
6)陈乐寿,王光锷.大地电磁测深法[M].北京:地质出版社,1990.
7)陈小斌.MT二维正演计算中地形影响的研究[J].石油物探,2000,39/3:112-120.
8)陈小斌.有限元直接迭代算法中一种基本结构的研究[J].物探化探计算技术,2000,22/2:132-136.
9)陈小斌,胡文宝.有限元直接迭代算法及其在线源频率域电磁响应计算中的应用[J].地球物理学报,2002,45/1:120-129.
10)陈小斌,张翔,胡文宝.有限元直接迭代算法在MT二维正演计算中的应用[J].石油地球物理勘探,2000,35/4:487-496.
11)程志平.电法勘探教程[M](1).北京:冶金工业出版社,2007.
12)刘德贵,费景高,于江永等.Fortran算法汇编(第一分册)[M].北京:国防工业出版社,1980.
13)刘云,王绪本.大地电磁二维自适应地形有限元正演模拟[J].地震地质,2010,32/3:382-391.
14)李予国,徐世浙,刘斌等.电导率分块连续变化的二维MT有限元模拟(II)[J].高校地质学报,1996,2/4:448-452.
15)马为.MT二维正演中辅助场计算新方法研究[D].北京:中国地震局,2007.
16)彭国伦.Fortran 95程序设计[M](第一版).北京:中国电力出版社,2002.
17)阮百尧.三角单元部分电导率分块连续变化点源二维电场有限元数值模拟[J].广西科学,2001,8/1:1-3.
18)阮百尧,熊彬.大型对称变带宽方程组的Cholesky分解法[J].物探化探计算技术,2000,22/4:361-363.
19)阮百尧,徐世浙.线元二维时间域电磁响应的有限差分解[J].高校地质学报,1996,2(4):437-447.
20)史明娟,徐世浙,刘斌.大地电磁二次函数插值的有限元法[J].地球物理学报,1997,40/3:421-429.
21)蔚宝强,胡文宝.MT资料的遗传法反演[J].江汉石油学院学报,1999,21(3):25-29.
22)王建谋等.大地电磁测深译文集.第一集[M](1).北京:地质出版社,1987.
23)熊彬,罗延钟.电导率分块均匀的瞬变电磁2.5维有限元数值模拟[J].地球物理学报,2006,49/2:590-597.
24)谢飞.带地形大地电磁场二维有限元数值模拟研究[D].北京:中国地质大学,2006.
25)徐世浙,于涛,李予国等.电导率分块连续变化的二维MT有限元模拟(I)[J].高校地质学报,1995,1/2:65-73.
26)徐士良.计算机常用算法[M](第二版).北京:清华大学出版社,1995.
27)徐世浙.地球物理中的有限单元法[M](第一版).北京:科学出版社,1994.
28)曾国.大地电磁二维有限元正演数值模拟[D].长沙:中南大学,2008.
29)周熙襄,钟本善.电法勘探数值模拟技术[M](第一版).成都:四川科学技术出版社,1986.
30)王建谋,陈乐寿等译.大地电磁测深译文集(第一集) [M](第一版).北京:地质出版社,1987.
31) Fortran. http://zh.wikipedia.org/wiki/fortran.
32)石应骏、刘国栋、王家映等.大地电磁测深法教程[M].成都:成都理工大学信息工程学院,1985.
33) Coggon J H.Electromagnetic and electrical modeling by the finite element method[J]. Geophysics, 1971, 36(1):132-155.
34) William L, Rodi. A technique for improving theaccuracy of finite element solution for magnetotelluric data[J].Geophysics, 1976,50 (7) :483-506.
35) Wannamaker P E, Stodt J A,Rijo L. A stable finite element solution for two-dimensional magnetotellnric modeling[J].Geophysics, 1987, 88 (5):277 -296.
36) Lugao P,Wannalnaker P E. Calculating the two-dimensional magnetotelluric Jacohian in finite elements using reciprocity[J].Geophysical Journal International, 1996, 127(2):806-810
37) Frank A, Ralph B, Klaus S. 2D finite element modeling of plane-wave diffusive time-harmonic electromagnetic fields using adaptive unstructured grids[A].Proceeding of the 17th Workshop,2004, S2:1-6.
38)简兴祥,大地电磁二三维建模及正演模拟[D].成都理工大学.2009