CN115186520B - 水平层状大地介质中矢量偶极源的时频电磁响应模拟方法 - Google Patents

水平层状大地介质中矢量偶极源的时频电磁响应模拟方法 Download PDF

Info

Publication number
CN115186520B
CN115186520B CN202211106686.3A CN202211106686A CN115186520B CN 115186520 B CN115186520 B CN 115186520B CN 202211106686 A CN202211106686 A CN 202211106686A CN 115186520 B CN115186520 B CN 115186520B
Authority
CN
China
Prior art keywords
formula
source
dipole source
layer
vector
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202211106686.3A
Other languages
English (en)
Other versions
CN115186520A (zh
Inventor
胡开颜
黄清华
胡文宝
唐新功
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Peking University
Original Assignee
Peking University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Peking University filed Critical Peking University
Priority to CN202211106686.3A priority Critical patent/CN115186520B/zh
Publication of CN115186520A publication Critical patent/CN115186520A/zh
Application granted granted Critical
Publication of CN115186520B publication Critical patent/CN115186520B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/11Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13Differential equations
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/15Correlation function computation including computation of convolution operations
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/16Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Computing Systems (AREA)
  • Operations Research (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明公布了一种水平层状大地介质中矢量偶极源的时频电磁响应模拟方法,通过建立在直角坐标系中通用的电型和磁型偶极子源引起的矢量势函数的传播系数矩阵,并递推求解得到地表的电磁场分量;分别计算水平和垂直偶极子源的电磁场,再合成得到任意矢量偶极子源的电磁场响应;可用于水平层状地中任意方向的电型偶极子源和磁型偶极子源;在得到偶极子源在均匀全空间势函数的振幅函数后进行计算,即可得到水平层状大地介质中偶极子源时频电磁响应。本发明可应用于模拟层状大地中地震前兆的电磁扰动、地下核污染源或煤矿巷道电型或磁型矢量偶极子源辐射的电磁波场,适用性好,应用价值高。

Description

水平层状大地介质中矢量偶极源的时频电磁响应模拟方法
技术领域
本发明提供一种模拟计算水平层状大地中存在任意矢量电偶极子源和磁偶极子源产生的电磁场的方法,具体涉及一种水平层状大地介质中矢量偶极子源的时频电磁响应的模拟方法,通过势函数传播系数矩阵进行模拟,求解时间域与频率域电磁场响应,属于地球电磁应用技术领域。
背景技术
目前已有地震电磁观测资料显示,地震在其孕育过程中可能产生电磁辐射而导致较大范围的地球电磁场扰动。对地震发生而导致破坏前产生的与之相关的电磁信号的解释对地震预测研究具有重要意义。地震破裂前辐射电磁波的源可视为电型或磁型偶极子源的叠加,基于此,模拟层状大地中矢量偶极子源产生的电磁响应是研究地震前兆电磁形象的基础。然而,目前用于计算层状介质中偶极子源电磁场响应的工作中未有适应任意方向的矢量偶极子源的全空间时间域和频率域的电磁场量计算的算法报道。Wait (1951)在文献(Wait J R. The magnetic dipole over the horizontally stratified earth. Can JPhys, 195l,29: 577-592.)中记载了利用水平电磁场量在界面的连续性求解了磁型源引起的矢量势函数的传播系数并实现了水平层状介质中磁偶极子源响应的计算。徐建华等(1994)在文献(徐建华,朱德怀,胡文宝,李鹏翔.多层介质中偶极子场的系数递推关系. 石油地球物理勘探,1994,29(1):69-74.)中采用与Wait类似的方式推导了磁型偶极子源引起的矢量势函数的层系数递推关系和传递矩阵,可用于求解多层介质中存在水平或垂直的偶极子产生的电磁场响应。已有的电磁场响应获取技术只能用于特定源或层结构,对于计算地球物理应用领域实际的测井中和地表的电磁场仍具有局限性,难以解决地震前兆电磁辐射和复杂偶极子源的电磁场响应的模拟计算。
发明内容
为了克服上述现有技术的不足,本发明提供一种水平层状大地介质中偶极子源时频电磁响应的模拟方法,分别计算了水平和垂直偶极子源的电磁场,然后合成得到任意矢量偶极子源的电磁场响应,可应用于解释层状大地中地震前兆中矢量偶极子源产生的电磁扰动现象,进一步可运用于大地震预测、识别地下核污染源或煤矿巷道中产生的电磁辐射。
本发明对现有的水平层状介质中磁偶极子源响应的计算方法和利用电型和磁型源引起的矢量势函数的层系数递推关系和传递矩阵求解多层介质中存在水平或垂直的偶极子产生的电磁场响应的方法加以改进,通过建立在直角坐标系中通用的电型和磁型偶极子源引起的矢量势函数的传播系数矩阵并递推求解得到地表的电磁场分量,可考虑用于水平层状地中任意方向的电型偶极子源和磁型偶极子源。对于不同类型的偶极子源,只需要得到偶极子源在均匀全空间势函数的振幅函数后进行计算即可得到水平层状大地介质中偶极子源时频电磁响应。
本发明的方法包括如下步骤:
A.建立水平层状大地介质中偶极子源电磁响应的一般表达式:
将由水平层状地中任意方向的电型偶极子源产生的电磁场的矢量势函数记为A=A z
Figure DEST_PATH_IMAGE001
Figure 252649DEST_PATH_IMAGE001
为直角坐标系中z轴方向的单位矢量;由磁型偶极子源产生的电磁场的矢量势函数记为F=F z
Figure 802580DEST_PATH_IMAGE001
即限定矢量势函数A和F只具有z分量A z F z
水平层状地中任意方向的电型偶极子源激发的电磁场的各分量分别表示为式1和式2:
Figure 385745DEST_PATH_IMAGE002
式1
Figure 679323DEST_PATH_IMAGE003
 式2
式中,复电导率
Figure 267431DEST_PATH_IMAGE004
=σ+iωε,i代表虚数单位,σ为介质电导率,ε为介电常数;E x 、E y 、E z 分别为直角坐标系下x、y、z轴方向的电场分量;H x H y 、H z 分别为x、y、z轴的磁场分量;其中H z 为TM极化场;
水平层状地中任意方向的磁型偶极子源引起的电磁场分量分别表示为式3和式4:
Figure 988262DEST_PATH_IMAGE005
     式3
Figure 419243DEST_PATH_IMAGE006
式4
其中,电场垂向分量E z 为零,E z 为TE极化场;
Figure 391879DEST_PATH_IMAGE007
为阻抗率,
Figure 693547DEST_PATH_IMAGE008
,i代表虚数单位,ω-信号角频率;
Figure 726225DEST_PATH_IMAGE009
为介质磁导率;式1和式4中k满足
Figure 644502DEST_PATH_IMAGE010
,是介质的波数;
B.计算电磁场的矢量势函数,分别求得无源区势函数一次场和含源层势函数的二次场的解;
分两种情况考虑,即无源区一次电磁场矢量势函数的解和含源层位二次电磁场矢量势函数的解;
B1. 对无源区势函数,可解齐次微分方程得:
Figure 14304DEST_PATH_IMAGE011
     式5
Figure 75539DEST_PATH_IMAGE012
     式6
式中,上波浪号代表对应势函数的傅里叶变换形式,A +F +代表对应的上行(+z方向)波解,A -F -则代表对应的下行(−z方向)波解,k x k y 分别是xy方向上的空间波数,u 2=k x 2+k y 2-k 2
B2. 求得含源层位势函数的通解:
对于含源层位,除互补解外,还需加上非齐次微分方程的特解。
对任意N层大地模型,z轴向上;设置源的中心为坐标原点,源所在的界面z s=0,该界面是虚拟界面,用于计算源的响应;含源界面上方即正z方向的层序号由下往上依次标记为i1,i2,…,iL;地层参数对应为电导率σ ii 、介电常数ε ii 和层厚度h ii ;含源界面下方即负z方向的层序号由上至下依次标记为j1,j2,…,jM,地层参数对应为电导率σ jj 、介电常数ε jj 和各层厚度h jj ;设A p F p 为源场的振幅,根据均匀全空间的平面波解,含源时非齐次微分方程的特解表示为式7和式8:
Figure 138173DEST_PATH_IMAGE013
        式7
Figure 153533DEST_PATH_IMAGE014
        式8
矢量势函数傅里叶变换形式的z分量
Figure 123763DEST_PATH_IMAGE016
Figure 540969DEST_PATH_IMAGE017
在源以上和源以下都产生衰减;N层大地中,n为层序号,含源的层位为n=i1n=j1
把式7和式8的源的特解与B1中的无源区势函数的解进行合并,即得到含源层的势函数的汉克尔积分形式,表示为式9和式10:
Figure 774504DEST_PATH_IMAGE018
式9
Figure 277161DEST_PATH_IMAGE019
   式10
式中,A n +A n - 、F n +F n -为待定系数,
Figure 254344DEST_PATH_IMAGE020
Figure 296030DEST_PATH_IMAGE021
Figure 231625DEST_PATH_IMAGE022
,用系数因子β n 区分有源和无源区,β n 取值为:
Figure 221578DEST_PATH_IMAGE023
C.根据地中垂直电偶极子源的势函数,计算得到势函数幅值:
将位于坐标原点处电极矩为Idl
Figure 736873DEST_PATH_IMAGE024
的垂直电偶极子源密度表示为:
Figure 987726DEST_PATH_IMAGE025
       式11
式中,
Figure 704009DEST_PATH_IMAGE026
是z轴方向的单位矢量;δ(x)、δ(y)、δ(z)为狄拉克函数,只有在坐标原点(0,0,0)为1,否则为0;垂直电偶极子源的势函数幅值表示为:
Figure 40313DEST_PATH_IMAGE027
          式12
由式9得任意层中含垂直电偶极子源的势函数在变换域的汉克尔积分形式:
Figure 359298DEST_PATH_IMAGE028
式13
传播系数A 0n +A 0n -由递推关系求得;递推关系包括:
C1. 由源和各层介质参数确定顶层和底层传递矩阵得到顶层和底层的传播系数,表示为式14和式15:
Figure 340024DEST_PATH_IMAGE029
 式14
Figure 351842DEST_PATH_IMAGE030
      式15
式中,q为2*2的传递矩阵Q的4个系数,系数值由各层介质参数决定;s i 为源向量S i 的2个系数
Figure 909863DEST_PATH_IMAGE031
s j 为源向量S j 的2个系数
Figure 406441DEST_PATH_IMAGE032
,系数值由源决定;
根据大地层状模型递推,得到源层的传递矩阵表达形式:
Figure 631886DEST_PATH_IMAGE033
       式16
Figure 549026DEST_PATH_IMAGE034
         式17
C2. 确定逐层传播系数;
在源的位置设置一个虚拟层界面,将大地模型分为源上方和源下方两种情况;
结合式16~17将源上方TM极化模式
Figure 867DEST_PATH_IMAGE035
传播系数写为:
Figure 927235DEST_PATH_IMAGE036
           式18
U i 的矩阵形式表示为:
Figure 616973DEST_PATH_IMAGE037
              式19
其中:
Figure 236174DEST_PATH_IMAGE038
       式20
Figure 768786DEST_PATH_IMAGE039
        式21
Figure 843052DEST_PATH_IMAGE040
      式22
Figure 777510DEST_PATH_IMAGE041
         式23
第i L 层为顶层底界面,向上为无限延伸空间,即A iL -=0无下行波,则有:
Figure 302033DEST_PATH_IMAGE042
     式24
源下方TM极化模式U j表示为矩阵形式:
Figure 695843DEST_PATH_IMAGE043
            式25
Figure 964013DEST_PATH_IMAGE044
               式26
Figure 752977DEST_PATH_IMAGE045
     式27
Figure 182822DEST_PATH_IMAGE046
        式28
Figure 565393DEST_PATH_IMAGE047
         式29
Figure 637254DEST_PATH_IMAGE048
      式30
类似地,TE极化模式的源层的传递矩阵形式相同,式中
Figure 280725DEST_PATH_IMAGE049
Figure 22416DEST_PATH_IMAGE050
A p F p
C3. 由顶层和底层的传播系数矩阵向模型内部逐层递推,由式16 ~式30得到源层的传播系数,可计算得到各层的传播系数A n +A n -n为层序号;
将传播系数代入式13求得各层的势函数;
将势函数代入式3~式6求得电磁场各分量核函数;
再应用数字滤波算法实现汉克尔变换,得到各场分量的频域响应;
D.建立地中水平电偶极子源的势函数:
水平偶极子源的电流源密度表示为:
Figure 751337DEST_PATH_IMAGE051
       式31
式中,
Figure 361310DEST_PATH_IMAGE052
x轴的方向单位矢量,δ(x)、δ(y)、δ(z)为狄拉克函数,只有在坐标原点(0,0,0)为1,否则为0;
任意层中含水平电偶极子源的势函数表示为:
Figure 265812DEST_PATH_IMAGE053
式32
Figure 37459DEST_PATH_IMAGE054
式33
式中,
Figure 519256DEST_PATH_IMAGE055
符号表示当z>0取–号时,z<0时取+号;
传播系数有:
Figure 306821DEST_PATH_IMAGE056
      式34
Figure 924884DEST_PATH_IMAGE057
      式35
Figure 867433DEST_PATH_IMAGE058
            式36
利用传播系数求得势函数,然后将式32~36代入式3~6求得电磁场各分量核函数,再应用数字滤波算法实现汉克尔变换,得到任意层n中水平电偶极子源激发的各场量频率域响应;
E.进行模拟计算得到地中任意方向的矢量电偶极子源场;包括:
E1. 场源分解;
由电偶极矩P e =Idl、方位角
Figure 711892DEST_PATH_IMAGE059
和倾角
Figure 663667DEST_PATH_IMAGE060
三个参数定义地中任意电偶极子,将电偶极矩分解为水平分量
Figure 136237DEST_PATH_IMAGE061
和垂直分量P z ,分别为:
Figure 656211DEST_PATH_IMAGE062
            式37
Figure 847021DEST_PATH_IMAGE063
             式38
对于层状介质而言任意方位的水平电偶极子源的场具有旋转对称性,可以通过坐标旋转得到任意方位水平电偶极子源的场;对于垂直电偶极子源,通过步骤C推导的垂直偶极子源响应的计算公式合成;
对于方位角为
Figure 336908DEST_PATH_IMAGE059
的水平电偶极子源,旋转后的坐标参数为:
Figure 929564DEST_PATH_IMAGE064
            式39
E2. 电磁场合成;
将电偶极矩为P e 时水平电偶极子源响应的场矢量(
Figure 620439DEST_PATH_IMAGE065
Figure 32966DEST_PATH_IMAGE066
)和垂直电偶极子源响应的场矢量(Ez和Hz)进行合成,得到地中任意电偶极子源的电磁场分量为:
Figure 983603DEST_PATH_IMAGE067
            式40
Figure 430765DEST_PATH_IMAGE068
            式41
E3. 对合成的电磁场做反傅里叶变换,得到时间域的电磁场分量;
由此实现对水平层状大地介质中偶极子源时频电磁响应的模拟计算。
与现有技术相比,本发明的有益技术效果:
本发明提供了水平层状大地中矢量偶极子源的时间域与频率域电磁场响应的计算方法,导出了适用于电型和磁型偶极子源和多层传播系数递推公式和求解算法;并利用多源组合实现线源和面源电磁响应的计算技术。本发明可应用于模拟层状大地中地震前兆的电磁扰动、地下核污染源或煤矿巷道电型或磁型矢量偶极子源辐射的电磁波场,可用于解释层状大地中地震前兆中矢量偶极子源产生的电磁扰动现象,进一步用于大地震预测、识别地下核污染源或煤矿巷道中产生的电磁辐射,具有较高的适用性和应用价值。
附图说明
图1是本发明提供的地中层状介质偶极子源的时频电磁响应模拟计算方法的流程框图。
图2是水平层状大地介质中的垂向偶极子源各层的示意图;
其中,AF表示由电型和磁型偶极子源引起的矢量势函数的传播系数,上角标+表示上行的电磁波,上角标-表示下行的电磁波。
图3是矢量电偶极源的分解示意图;
其中,
Figure 417175DEST_PATH_IMAGE059
Figure 723523DEST_PATH_IMAGE060
分别表示矢量偶极子源的方位角和倾角,P e 偶极源极距,可分解为水平分量
Figure 555212DEST_PATH_IMAGE061
和垂直分量P z
图4是本发明实施例中采用的4层大地模型的结构与参数示意图;
其中,
Figure 997826DEST_PATH_IMAGE069
代表介质磁导率,σ代表介质电导率,下角标代表层序号。
图5是地中矢量电偶极子源响应各场量10Hz幅值的x-y(水平)剖面分布图,
其中,(a)-(c)分别为电场频域响应的x、y和z分量;(d)-(f)分别为磁场频域响应的x、y和z分量。
图6是地中矢量电偶极子源响应各场量10Hz幅值的x-z(垂直)剖面分布图,
其中,(a)-(c)分别为电场频域响应的x、y和z分量;(d)-(f)分别为磁场频域响应的x、y和z分量。
图7是地中矢量电偶极子源0.1s各场量幅值的x-z(垂直)剖面分布图,
其中,(a)-(c)分别为电场响应的x、y和z分量;(d)-(f)分别为磁场响应的x、y和z分量。
具体实施方式
下面结合附图,通过实施例进一步描述本发明,但不以任何方式限制本发明的范围。
本发明提供一种水平层状大地介质中偶极子源时频电磁响应的模拟方法,通过导出适用于电型和磁型偶极子源和多层传播系数递推公式和求解算法,并利用多源组合实现线源和面源电磁响应的模拟计算,具有较高的适用性和应用价值。
具体实施时,设置空间大气层与地层模型,设置矢量电偶极子源的方位角和倾角以及源的坐标位置,此例设置源的位置(x sy sz s),并设置各层参数。本发明的水平层状大地介质中偶极子源时频电磁响应的模拟方法包括如下步骤:
A.建立水平层状大地介质中偶极子源电磁响应的一般表达形式:
由于不同水平层状大地介质的分界面分别与相应水平层状大地介质的等z平面重合,因此在解边值问题时可将地表电磁场的偏微分方程变换成关于等z平面的常微分方程。边值问题的通解是非齐次微分方程的特解和齐次方程的互补解的和。设由电型偶极子源产生的电磁场的矢量势函数记为A=A z
Figure 155138DEST_PATH_IMAGE070
Figure 683203DEST_PATH_IMAGE071
为直角坐标系中z轴方向的单位矢量,A z 为A的z分量
Figure 318583DEST_PATH_IMAGE072
;而由磁型偶极子源产生的电磁场的矢量势函数记为F=F z
Figure 5917DEST_PATH_IMAGE026
即限定矢量势函数A和F只具有z分量A z F z 。那么,电型偶极子源激发的电磁场的各分量分别写为:
Figure 708031DEST_PATH_IMAGE073
式1
Figure 848025DEST_PATH_IMAGE074
    式2
式中,复电导率
Figure 896884DEST_PATH_IMAGE075
=σ+iωε,σ-介质电导率(S/m);ε为介电常数,对一般应用,取ε=ε 0=8.854*10-12F/m。E x 、E y 、E z 分别为直角坐标系下x、y、z轴方向的电场分量,H x 、H y、 H z 分别为x、 y、z轴的磁场分量,其中由于H z (磁场垂向分量)为零,H z 被称为TM极化场。类似地,磁型偶极子源引起的电磁场分量分别表示为:
Figure 438724DEST_PATH_IMAGE076
     式3
Figure 813204DEST_PATH_IMAGE077
式4
其中,由于E z (电场垂向分量)为零,E z 被称为TE极化场。
Figure 174916DEST_PATH_IMAGE078
为阻抗率,
Figure 152099DEST_PATH_IMAGE079
i代表虚数单位;ω-信号角频率;
Figure 689391DEST_PATH_IMAGE080
为介质磁导率;对一般应用,取
Figure 93827DEST_PATH_IMAGE069
=
Figure 677255DEST_PATH_IMAGE069
0=4π*10-7 H/m。
B.计算电型偶极子源的电磁场和磁型偶极子源的电磁场的矢量势函数;
分两种情况,即无源区一次电磁场矢量势函数的解(齐次偏微分方程互补解)和含源层位二次电磁场矢量势函数的解(非齐偏微分方程特解);
B1. 对无源区势函数可解齐次微分方程得:
Figure 458129DEST_PATH_IMAGE011
     式5
Figure 348463DEST_PATH_IMAGE012
     式6
式中,上波浪号代表对应势函数的傅里叶变换形式,A +F +代表对应的上行(+z方向)波解,A -F -则代表对应的下行(−z方向)波解,k x k y 分别是xy方向上的空间波数,u 2=k x 2 +k y 2-k 2。对于层状大地介质的无源区,即除图2中i1或j1以外的区域,可在不同的层内区域求定解,与适当的边界条件结合可推得各层的幅值(常数)。
B2. 求解含源层位势函数的通解形式:
对于含源层位,除互补解外,还需加上非齐次微分方程的特解。考虑任意N层大地模型如图1所示,z轴向上。设置源的中心为坐标原点,源所在的界面z s=0,该界面是虚拟界面,主要用于计算源的响应。含源界面上方(正z方向)的层序号由下往上依次标记为i1,i2,…,iL,地层参数对应为电导率σ ii 、介电常数ε ii 和层厚度h ii ;含源界面下方(负z方向)的层序号由上至下依次标记为j1,j2,…,jM,地层参数对应为电导率σ jj 、介电常数ε jj 和各层厚度h jj 。设A p F p 为源场的振幅,则根据均匀全空间的平面波解,含源时非齐次微分方程的特解表示为:
Figure 923801DEST_PATH_IMAGE013
        式7
Figure 994525DEST_PATH_IMAGE014
        式8
矢量势函数傅里叶变换形式的z分量
Figure 188877DEST_PATH_IMAGE081
Figure 559815DEST_PATH_IMAGE082
在源以上和源以下都产生衰减,这里的A p F p 依源而定,n代表层位序号。如图2所示的N层大地模型中,含源的层位为n=i1n=j1。把式7和式8的源的特解与B1中的无源区势函数的解进行合并,即可得到含源层的势函数的汉克尔积分形式,表示为:
Figure 306054DEST_PATH_IMAGE018
式9
Figure 5020DEST_PATH_IMAGE019
   式10
式中,A n +A n - 、F n +F n -为待定系数,
Figure 127697DEST_PATH_IMAGE020
Figure 353142DEST_PATH_IMAGE021
Figure 145649DEST_PATH_IMAGE022
,这里用系数因子β n 区分有源和无源两种情况,系数β n 取值为:
Figure 190965DEST_PATH_IMAGE083
C.计算地中垂直电偶极子源的势函数:
一短载流导线,若其长度远小于与观测点的距离,可被等效为电偶极子源;若电偶极子源距离观测点较近,或是载流导线较长,可把载流导线发射源看作多个偶极子源的叠加。位于坐标原点处电极矩为Idl
Figure 851753DEST_PATH_IMAGE084
的垂直电偶极子源密度(A/m2)可表示为:
Figure 571186DEST_PATH_IMAGE085
     式11
式中,
Figure 659227DEST_PATH_IMAGE024
是z轴方向的单位矢量,δ(x)、δ(y)、δ(z)为狄拉克函数,只有在坐标原点(0,0,0)为1,否则为0。垂直电偶极子源的势函数幅值表示为:
Figure 191840DEST_PATH_IMAGE086
          式12
由式9得任意层中含垂直电偶极子源的势函数在变换域的汉克尔积分形式:
Figure 797265DEST_PATH_IMAGE087
式13
传播系数A n +A n -由递推关系求得,β n 参见步骤B2。递推关系包括:
C1. 由源和各层介质参数确定顶层和底层传递矩阵得到顶层和底层的传播系数:
Figure 997302DEST_PATH_IMAGE088
 式14
Figure 256245DEST_PATH_IMAGE089
      式15
式中,q为2*2的传递矩阵Q的4个系数,系数值由各层介质参数决定;s i 为源向量S i 的两个系数
Figure 151520DEST_PATH_IMAGE090
s j 为源向量S j 的两个系数
Figure 685269DEST_PATH_IMAGE091
,系数值由源决定。
根据图2的大地层状模型递推,得到源层的传递矩阵表达形式:
Figure 474234DEST_PATH_IMAGE033
       式16
Figure 779444DEST_PATH_IMAGE034
         式17
C2. 确定逐层传播系数递推关系式,在源的位置设置一个虚拟层界面,再将大地模型分为源上方和源下方两种情况考虑。
源上方TM极化模式表示为矩阵形式:
Figure 286649DEST_PATH_IMAGE036
           式18
Figure 358510DEST_PATH_IMAGE037
              式19
Figure 381742DEST_PATH_IMAGE092
       式20
Figure 513646DEST_PATH_IMAGE039
        式21
Figure 242567DEST_PATH_IMAGE093
      式22
Figure 993486DEST_PATH_IMAGE041
         式23
第i L 层为顶层底界面,向上为无限延伸空间,即A iL -=0无下行波,则有:
Figure 491463DEST_PATH_IMAGE042
     式24
源下方TM极化模式U j表示为矩阵形式:
Figure 528689DEST_PATH_IMAGE043
            式25
Figure 885852DEST_PATH_IMAGE044
               式26
Figure 33937DEST_PATH_IMAGE094
     式27
Figure 652000DEST_PATH_IMAGE046
        式28
Figure 735494DEST_PATH_IMAGE095
         式29
Figure 704587DEST_PATH_IMAGE096
      式30
类似地,TE极化模式得源层的传递矩阵形式相同,只需替换式中
Figure 656362DEST_PATH_IMAGE097
Figure 502833DEST_PATH_IMAGE098
A p F p
C3.由顶层(式14)和底层(式15)的传播系数矩阵向模型内部逐层递推,由式16 ~式30得到源层的传播系数,然后计算得到各层n的传播系数A n +A n -n为层序号;将计算得到的传播系数代入式13求得各层的势函数;最后将势函数代入式3~式6求得电磁场各分量核函数;再应用数字滤波算法实现汉克尔变换,得到各场分量的频域响应;
D.建立地中水平电偶极子源的势函数:
水平偶极子源的电流源密度可表示为:
Figure 881862DEST_PATH_IMAGE099
      式31
式中,
Figure 72672DEST_PATH_IMAGE052
x轴的方向单位矢量,δ(x)、δ(y)、δ(z)为狄拉克函数,只有在坐标原点(0,0,0)为1,否则为0。任意层中含水平电偶极子源的势函数可表示为:
Figure 703505DEST_PATH_IMAGE100
 式32
Figure 296160DEST_PATH_IMAGE101
 式33
式中,
Figure 846090DEST_PATH_IMAGE102
符号表示当z>0取–号时,z<0时取+号。传播系数(水平)有:
Figure 399562DEST_PATH_IMAGE103
      式34
Figure 693140DEST_PATH_IMAGE104
      式35
Figure 140302DEST_PATH_IMAGE105
            式36
与步骤C3类似处理,利用传播系数求得势函数然后代入式3~6求得电磁场各分量核函数,最后应用数字滤波算法实现汉克尔变换得到任意层n中水平电偶极子源激发的各场量频率域响应。
E.模拟计算得到地中任意方向的矢量电偶极子源场;包括:
E1. 场源分解
由电偶极矩P e =Idl、方位角
Figure 861133DEST_PATH_IMAGE059
和倾角
Figure 636323DEST_PATH_IMAGE060
三个参数定义地中任意电偶极子,如图3所示,可以将电偶极矩分解为水平分量
Figure 733592DEST_PATH_IMAGE061
和垂直分量P z ,分别为:
Figure 300839DEST_PATH_IMAGE106
            式37
Figure 566473DEST_PATH_IMAGE063
             式38
对于层状介质而言任意方位的水平电偶极子源的场具有旋转对称性,可以用已经导出的x方向场的计算公式基础上,直接通过坐标旋转得到任意方位水平电偶极子源的场。对于垂直电偶极子源,由于地面边界的作用,场的分布不具有旋转对称性,需要用步骤C推导的垂直偶极子源响应的计算公式合成。
对于方位角为
Figure 953592DEST_PATH_IMAGE059
的水平电偶极子源,旋转后的坐标参数为:
Figure 854552DEST_PATH_IMAGE064
            式39
E2. 电磁场合成
采用前面推导的公式用旋转后的坐标参数分别计算出电偶极矩为P e 时水平电偶极子源响应的场矢量(
Figure 151672DEST_PATH_IMAGE065
Figure 214306DEST_PATH_IMAGE066
)和垂直电偶极子源响应的场矢量(Ez和Hz),则可合成得到地中任意电偶极子源的场分量为:
Figure 88721DEST_PATH_IMAGE067
            式40
Figure 403159DEST_PATH_IMAGE068
            式41
E3. 对合成的电磁场做反傅里叶变换,得到水平层状大地介质中偶极子源时间域的电磁响应的x、y、z方向分量。
具体实施时,设置空间大气层与地层模型如图4所示,设置一方位角
Figure 679420DEST_PATH_IMAGE107
=60度,倾角
Figure 912955DEST_PATH_IMAGE108
=30度的矢量电偶极子源位于地中(0,0,z s)处,各层介质参数分别为:
Figure 415612DEST_PATH_IMAGE009
1=
Figure 392795DEST_PATH_IMAGE009
2=
Figure 523562DEST_PATH_IMAGE009
3=
Figure 567479DEST_PATH_IMAGE009
0σ 0 =10–12S/m;σ 1 =0.05 S/m、σ 2 =0.001 S/m、σ 3 =0.1 S/m。坐标设置如图所示,z轴向上,地面z=0;根据存在多层地下界面各层厚有h 1=2000 m,h 2=5000m,
Figure 416486DEST_PATH_IMAGE109
。垂直电偶极子源位于第2层中,z s=–4000 m。由步骤C和步骤D分别求得垂直电偶极子源和水平电偶极子源的场分量,再由步骤E1对电偶极矩进行分解和坐标旋转,再由步骤E2得到水平电偶极子源响应的场矢量(
Figure 197361DEST_PATH_IMAGE065
Figure 182634DEST_PATH_IMAGE066
)和垂直电偶极子源响应的场矢量(Ez和Hz)合成得到频率域的电磁场响应。根据本例设置的单位电偶极子源的倾角
Figure 898917DEST_PATH_IMAGE060
=30度和方位角
Figure 235221DEST_PATH_IMAGE059
=60度利用式37-41合成矢量源的各场量。选择频率为10 Hz为例显示如图4-5。其中图5取自近地表水平x-y剖面,图6取自x-z纵剖面(y=0m)。最后根据步骤E3对合成的电磁场频域响应进行反傅里叶变换得到时间域的电磁场分量的幅值,这里以时间为0.1s时的电磁场响应纵剖面为例图7中进行展示。
需要注意的是,公布实施例的目的在于帮助进一步理解本发明,但是本领域的技术人员可以理解:在不脱离本发明及所附权利要求的范围内,各种替换和修改都是可能的。因此,本发明不应局限于实施例所公开的内容,本发明要求保护的范围以权利要求书界定的范围为准。

Claims (3)

1.一种水平层状大地介质中矢量偶极源的时频电磁响应模拟方法,通过建立在直角坐标系中通用的电型和磁型偶极源引起的矢量势函数的传播系数矩阵,并递推求解得到地表的电磁场分量;分别计算水平和垂直偶极源的电磁场,再合成得到任意矢量偶极源的电磁场响应;在得到偶极源在均匀全空间势函数的振幅函数后进行计算,即可得到水平层状大地介质中偶极源时频电磁响应;包括如下步骤:
A.建立水平层状大地介质中偶极源电磁响应的一般表达式:
将由水平层状大地介质中任意方向的电型偶极源产生的电磁场的矢量势函数记为A=A z
Figure 633535DEST_PATH_IMAGE002
Figure 979197DEST_PATH_IMAGE003
为直角坐标系中z轴方向的单位矢量;由磁型偶极源产生的电磁场的矢量势函数记为F=F z
Figure 381360DEST_PATH_IMAGE004
,即限定矢量势函数A和F只具有z分量A z F z
水平层状大地介质中任意方向的电型偶极源激发的电磁场的各分量分别表示为式1和式2:
Figure DEST_PATH_IMAGE005
式1
Figure 112555DEST_PATH_IMAGE006
 式2
式中,复电导率
Figure DEST_PATH_IMAGE007
=σ+iωε,i代表虚数单位,σ为介质电导率,ε为介电常数;E x 、E y 、E z 分别为直角坐标系下x、y、z轴方向的电场分量;H x H y 、H z 分别为x、y、z轴的磁场分量;其中H z 为TM极化场;
水平层状大地介质中任意方向的磁型偶极源引起的电磁场分量分别表示为式3和式4:
Figure 119301DEST_PATH_IMAGE008
     式3
Figure DEST_PATH_IMAGE009
式4
其中,电场垂向分量E z 为零,E z 为TE极化场;
Figure 686680DEST_PATH_IMAGE010
为阻抗率,
Figure DEST_PATH_IMAGE011
,i代表虚数单位,ω为信号角频率;
Figure 486008DEST_PATH_IMAGE012
为介质磁导率;式1和式4中k满足
Figure 743814DEST_PATH_IMAGE013
,是介质的波数;
B.计算电磁场的矢量势函数,分别求得无源区势函数一次场和含源层势函数二次场的解;包括:
B1. 对无源区势函数,解齐次微分方程得:
Figure DEST_PATH_IMAGE014
     式5
Figure 986708DEST_PATH_IMAGE015
     式6
式中,上波浪号代表对应势函数的傅里叶变换形式,A +F +代表对应的上行即正z方向的波解,A -F -代表对应的下行即负z方向的波解,k x k y 分别是xy方向上的空间波数,u 2=k x 2 +k y 2-k 2
B2. 求得含源层位势函数的通解,包括:
对任意N层大地模型,z轴向上;设置源的中心为坐标原点,源所在的界面z s=0,该界面是虚拟界面,用于计算源的响应;含源界面上方即正z方向的层序号由下往上依次标记为i1,i2,…,iL;地层参数对应为电导率σ ii 、介电常数ε ii 和层厚度h ii ;含源界面下方即负z方向的层序号由上至下依次标记为j1,j2,…,jM,地层参数对应为电导率σ jj 、介电常数ε jj 和各层厚度h jj ;设A p F p 为源场的振幅,根据均匀全空间的平面波解,含源时非齐次微分方程的特解表示为式7和式8:
Figure DEST_PATH_IMAGE016
        式7
Figure 556230DEST_PATH_IMAGE017
        式8
nN层大地中的层序号;
将式7和式8的源的特解与步骤B1中的无源区势函数的解进行合并,即得到含源层的势函数的汉克尔积分形式,表示为式9和式10:
Figure DEST_PATH_IMAGE018
式9
Figure 907052DEST_PATH_IMAGE019
   式10
式中,A n +A n - 、F n +F n -为待定系数,
Figure DEST_PATH_IMAGE020
Figure 878419DEST_PATH_IMAGE021
Figure DEST_PATH_IMAGE022
,系数因子β n 用于区分有源区和无源区;
C.根据地中垂直电偶极源的势函数,计算得到势函数幅值:
由式9得任意层中含垂直电偶极源的势函数在变换域的汉克尔积分形式:
Figure 416848DEST_PATH_IMAGE023
式13
垂直电偶极源的势函数幅值表示为:
Figure DEST_PATH_IMAGE024
          式12
其中,位于坐标原点处电极矩为Idl
Figure DEST_PATH_IMAGE026
Idl
Figure 817873DEST_PATH_IMAGE026
的垂直电偶极源密度表示为:
Figure 365529DEST_PATH_IMAGE027
       式11
式中,
Figure DEST_PATH_IMAGE028
是z轴方向的单位矢量;δ(x)、δ(y)、δ(z)为狄拉克函数;
传播系数A 0n +A 0n -由递推关系求得;包括:
C1. 由源和各层介质参数确定顶层和底层传递矩阵得到顶层和底层的传播系数,表示为式14和式15:
Figure 942135DEST_PATH_IMAGE029
 式14
Figure DEST_PATH_IMAGE030
      式15
式中,q为传递矩阵Q的系数;s i 为源向量S i 的系数,s j 为源向量S j 的系数;
根据大地层状模型递推,得到源层的传递矩阵表达形式为:
Figure 776099DEST_PATH_IMAGE031
       式16
Figure DEST_PATH_IMAGE032
         式17
C2. 确定逐层传播系数;
在源的位置设置一个虚拟层界面,将大地模型分为源上方和源下方;
源上方TM极化模式
Figure 805366DEST_PATH_IMAGE033
传播系数表示为:
Figure DEST_PATH_IMAGE034
           式18
U i 的矩阵形式表示为:
Figure 891134DEST_PATH_IMAGE035
              式19
其中:
Figure 837093DEST_PATH_IMAGE036
       式20
Figure 514062DEST_PATH_IMAGE037
        式21
Figure DEST_PATH_IMAGE038
      式22
Figure 785554DEST_PATH_IMAGE039
         式23
第i L 层为顶层底界面,向上为无限延伸空间,即A iL -=0无下行波,则有:
Figure DEST_PATH_IMAGE040
     式24
源下方TM极化模式U j表示为矩阵形式:
Figure 534067DEST_PATH_IMAGE041
            式25
Figure 741057DEST_PATH_IMAGE042
               式26
Figure DEST_PATH_IMAGE043
     式27
Figure 933135DEST_PATH_IMAGE044
        式28
Figure 327208DEST_PATH_IMAGE045
         式29
Figure DEST_PATH_IMAGE046
      式30
TE极化模式的源层的传递矩阵形式为:将式25~30中的
Figure 879412DEST_PATH_IMAGE047
改为
Figure DEST_PATH_IMAGE048
A p 改为F p
C3. 由顶层和底层的传播系数矩阵向模型内部逐层递推,由式16 ~式30得到源层的传播系数,计算得到各层的传播系数A n +A n -
将传播系数代入式13求得各层的势函数;
将势函数代入式3~式6求得电磁场各分量核函数;
再应用数字滤波算法实现汉克尔变换,得到各场分量的频域响应;
D.建立地中水平电偶极源的势函数:
水平偶极源的电流源密度表示为:
Figure 19537DEST_PATH_IMAGE049
       式31
式中,
Figure DEST_PATH_IMAGE050
x轴的方向单位矢量;δ(x)、δ(y)、δ(z)为狄拉克函数;
任意层中含水平电偶极源的势函数表示为:
Figure 631784DEST_PATH_IMAGE051
 式32
Figure DEST_PATH_IMAGE052
 式33
式中,
Figure 247573DEST_PATH_IMAGE053
符号表示当z>0取–号时,z<0时取+号;
传播系数表示为:
Figure DEST_PATH_IMAGE054
      式34
Figure 354201DEST_PATH_IMAGE055
      式35
Figure DEST_PATH_IMAGE056
            式36
利用传播系数求得势函数,将式32~36代入式3~6求得电磁场各分量核函数,再应用数字滤波算法实现汉克尔变换,得到任意层n中水平电偶极源激发的各场量频率域响应;
E.进行模拟计算得到地中任意方向的矢量电偶极源场;包括:
E1. 场源分解;
由电偶极矩P e =Idl、方位角
Figure 863679DEST_PATH_IMAGE057
和倾角
Figure DEST_PATH_IMAGE058
定义地中任意电偶极子,将电偶极矩分解为水平分量
Figure 129051DEST_PATH_IMAGE059
和垂直分量P z ,分别为:
Figure DEST_PATH_IMAGE060
            式37
Figure 232136DEST_PATH_IMAGE061
             式38
对于方位角为
Figure 657301DEST_PATH_IMAGE057
的水平电偶极源,旋转后的坐标参数为:
Figure DEST_PATH_IMAGE062
            式39
E2. 电磁场合成;
将电偶极矩为P e 时水平电偶极源响应的场矢量
Figure 162232DEST_PATH_IMAGE063
Figure DEST_PATH_IMAGE064
,和垂直电偶极源响应的场矢量Ez和Hz进行合成,得到地中任意电偶极源的电磁场分量为:
Figure 335856DEST_PATH_IMAGE065
            式40
Figure DEST_PATH_IMAGE066
            式41
E3. 对合成的电磁场进行反傅里叶变换,得到时间域的电磁场分量;
由此实现对水平层状大地介质中偶极源时频电磁响应的模拟计算。
2.如权利要求1所述水平层状大地介质中矢量偶极源的时频电磁响应模拟方法,其特征是,步骤B2表示含源时非齐次微分方程的特解中,矢量势函数傅里叶变换形式的z分量
Figure 50871DEST_PATH_IMAGE067
Figure DEST_PATH_IMAGE068
在源以上和源以下均产生衰减;含源的层位为n=i1n=j1
3.如权利要求1所述水平层状大地介质中矢量偶极源的时频电磁响应模拟方法,其特征是,步骤B2中,β n 取值为:
Figure 499301DEST_PATH_IMAGE069
CN202211106686.3A 2022-09-09 2022-09-09 水平层状大地介质中矢量偶极源的时频电磁响应模拟方法 Active CN115186520B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211106686.3A CN115186520B (zh) 2022-09-09 2022-09-09 水平层状大地介质中矢量偶极源的时频电磁响应模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211106686.3A CN115186520B (zh) 2022-09-09 2022-09-09 水平层状大地介质中矢量偶极源的时频电磁响应模拟方法

Publications (2)

Publication Number Publication Date
CN115186520A CN115186520A (zh) 2022-10-14
CN115186520B true CN115186520B (zh) 2023-01-24

Family

ID=83524639

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211106686.3A Active CN115186520B (zh) 2022-09-09 2022-09-09 水平层状大地介质中矢量偶极源的时频电磁响应模拟方法

Country Status (1)

Country Link
CN (1) CN115186520B (zh)

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114239268A (zh) * 2021-12-16 2022-03-25 西北工业大学 一种基于Romberg获取水下双电偶极子阵列跨界面辐射场的方法

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8784704B2 (en) * 2012-08-28 2014-07-22 The United States Of America As Represented By The Secretary Of The Navy Broadband artificial dielectric with dynamically optical control
CN110376655A (zh) * 2019-07-25 2019-10-25 中南大学 层状地质条件下任意位置偶极源电磁场响应的计算方法
CN112285435B (zh) * 2020-10-29 2022-04-22 中国舰船研究设计中心 一种大功率磁场辐射源的等效模拟方法
CN113076508A (zh) * 2021-04-02 2021-07-06 北京环境特性研究所 基于垂直磁偶极子在半空间下的低频近场快速计算方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114239268A (zh) * 2021-12-16 2022-03-25 西北工业大学 一种基于Romberg获取水下双电偶极子阵列跨界面辐射场的方法

Also Published As

Publication number Publication date
CN115186520A (zh) 2022-10-14

Similar Documents

Publication Publication Date Title
Grayver et al. 3D inversion and resolution analysis of land-based CSEM data from the Ketzin CO 2 storage formation
CN110058317B (zh) 航空瞬变电磁数据和航空大地电磁数据联合反演方法
Huang et al. Bidimensional empirical mode decomposition (BEMD) for extraction of gravity anomalies associated with gold mineralization in the Tongshi gold field, Western Shandong Uplifted Block, Eastern China
CN110058315B (zh) 一种三维各向异性射频大地电磁自适应有限元正演方法
Teasdale Methods for understanding poorly exposed terranes: the interpretive geology and tectonothermal evolution of the western Gawler Craton
Saad Understanding gravity gradients—a tutorial
Wilson et al. Real-time 3D inversion of ultra-deep resistivity logging-while-drilling data
CN104992029B (zh) 一种多尺度非均匀月壤层内离散随机介质建模方法
CN111666534A (zh) 电性任意各向异性电磁场解耦方法
Papadopoulos et al. Electrical resistivity tomography for the modelling of cultural deposits and geomophological landscapes at Neolithic sites: a case study from Southeastern Hungary
Zaher et al. Integration of geophysical methods for groundwater exploration: A case study of El Sheikh Marzouq area, Farafra Oasis, Egypt
Zhang et al. Regional gravity survey and application in oil and gas exploration in China
Lee et al. Three-dimensional facies architecture and three-dimensional calcite concretion distributions in a tide-influenced delta front, Wall Creek Member, Frontier Formation, Wyoming
Zhang et al. New seismic evidence for continental collision during the assembly of the Central Asian Orogenic Belt
Ringrose et al. Permeability estimation functions based on forward modeling of sedimentary heterogeneity
CN115186520B (zh) 水平层状大地介质中矢量偶极源的时频电磁响应模拟方法
Ramdani et al. How in-place volumes of subsurface reservoir models are impacted by using 3D high-resolution outcrop analogue data. A case study using depositional architectural heterogeneity of stromatoporoid/coral buildups of the Hanifa Fm, Saudi Arabia
Chen et al. Induced polarization and magnetic responses of serpentinized ultramafic rocks from mid‐ocean ridges
Jol et al. Stratigraphic imaging of the Navajo Sandstone using ground-penetrating radar
Aghil et al. Delineation of electrical resistivity structure using Magnetotellurics: a case study from Dholera coastal region, Gujarat, India
EP1070970B1 (en) A method of three dimensional reconstructing a physical magnitude inside a borehole
Malik et al. Ground‐Penetrating Radar Investigations along Hajipur Fault: Himalayan Frontal Thrust—Attempt to Identify Near Subsurface Displacement, NW Himalaya, India
Thiel et al. Full 3D reservoir mapping using deep directional resistivity measurements
Abbaszadeh et al. Integrated geostatistical reservoir characterization of turbidite sandstone deposits in Chicontepec Basin, Gulf of Mexico
Hoversten et al. Subsalt imaging via seaborne electromagnetics

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant