本发明涉及地球物理探测，尤其涉及一种基于模型降阶算法的瞬变电磁全空间快速正演方法。
背景技术：
1、瞬变电磁法(transient electromagnetic method，简称tem)是一种重要的地球物理勘探方法，以其显著的优点(经济、无损、快速、分辨率高等)广泛应用于金属矿产勘探、煤田灾害探测、地下断层划分、地下水探测和隧道超前地质预报等领域。
2、瞬变电磁的三维正演方法能为勘探的最优化设计提供强有力的帮助，是三维反演和成像技术发函的前提。从发射源的物理性质来看，发射源可概括为回线源和接地导线源两类。由于目前时间域的瞬变电磁三维正演方法提供需要求解大型线性方程组和隐式时间步迭代，即使结合并行求解器，仍然需要若干次的矩阵分解和回代求解，计算成本较大，无法为大规模的三维反演和成像技术提供基础，因此如何提高瞬变电磁的正演效率一直是三维正演与反演研究的一个具有挑战性的课题。按照空间离散的不同处理方法，模拟复杂地层模型的瞬变电磁三维正演算法主要包括有限差分法，有限体积法，有限元法。上述多数计算方法都是先计算频率域的瞬变电磁场，然后通过各种时频转换方法得到时间域的瞬变电磁响应，计算效率较低。直接时间域的三维矢量有限元正演，相比于以往方法，虽然直接在时间域求解，但是一次正演需要几十次的矩阵分解和上千次的回代，正演时间也较长。
3、一般来说，提高正演效率的方法有两种。第一种是通过降低系数矩阵的阶数提高正演效率，例如八叉树网格和非结构化网格；另一种方法是优化时间域积分方程提高正演效率。因此，本发明根据最优偏移思想：该方法对于给定的计算时间范围，计算出最优的偏移量与最小子空间阶数。使用这种方法，在保证计算精度的前提下，可以极大地提高瞬变电磁三维正演的效率。
技术实现思路
1、本发明的目的是提供一种基于模型降阶算法的瞬变电磁全空间快速正演方法，基于sai-krylov子空间算法的矢量有限元方法，在保证精度的前提下，极大地加快了正演的计算时间，这为瞬变电磁的反演及成像工作奠定了基础，实现对形状复杂模型的三维快速正演。
2、为实现上述目的，本发明提供了一种基于模型降阶算法的瞬变电磁全空间快速正演方法，包括以下步骤：
3、s1、基于sai-krylov子空间方法，求解瞬变电磁法的控制方程；
4、s2、对回线源和电性源初始电场进行求取；
5、s3、进行子空间阶数与最优偏移的求取；
6、s4、在使用sai-krylov子空间算法时，进行m次方程求解。
7、优选的，在步骤s1中，基于sai-krylov子空间方法，构建忽略位移电流的时间域的maxwell控制方程，如下所示：
8、
9、▽×μ-1b-σe＝s (2)
10、其中，e为电场强度，b是磁感应强度，t是时间，σ是电导率张量，μ是磁导率，s是外加电流密度。
11、优选的，在步骤s1中，瞬变电磁法的控制方程满足第一类边界条件，如下所示：
12、▽×e|γ＝0 (3)
13、▽×b|γ＝0 (4)。
14、优选的，在步骤s1中，基于galerkin法，二次场的控制方程，如下所示：
15、
16、其中，f为矢量基函数；
17、对上式离散后，得：
18、
19、其中，ni和nj为x,y,z三个方向的插值基函数；
20、令m1＝∫v(▽×ni)·(▽×nj)dv；m2＝∫vμ0σni·njdv，则上述式(6)化简为：
21、
22、令则最终求解问题，如下所示：
23、
24、e(t＝0)＝e0 (9)
25、其中，e为待求电场，e0为棱边上的初始电场。
26、优选的，e(t)的近似求解公式，如下所示：
27、
28、其中，γ为sai-krylov子空间的最优偏移量，m为子空间的阶数，vm是有m列的标准正交基，hm,m是hessenberg矩阵，e1是长度为m的单位向量。
29、优选的，在步骤s2中，求取回线源初始电场，具体过程如下：
30、对于回线源的瞬变电磁接收装置，在电流关断前，空间中只存在稳定电流产生的静态磁场，即只存在静态的磁矢势，该磁矢势分布与模型的电阻率无关，该磁矢势即为t＝0时刻整个空间中的分布a0＝a(0)；设地下模型的磁导率为空气磁导率，均和真空磁导率μ0相同，则柱坐标系下的磁矢势的解析表达式，如下所示：
31、
32、其中，a为回线源的半径；回线源的中心点为(x0,y0,z0)；接受点位置为(xr,yr,zr)，且x＝xr-x0，y＝yr-y0，z＝zr-z0，r＝(x2+y2)12，
33、分别为第一类与第二类完全椭圆积分；
34、将柱坐标中的磁矢势转化到直角坐标系中，如下所示：
35、
36、在计算关断后的初始电场时，电场与磁矢势的控制方程，如下所示：
37、▽×▽×a-μ0σe0＝0 (16)
38、根据加权余量法，采用第一类边界条件▽×a|γ＝0，上式最终可以写为：
39、∫v(▽×ni)·(▽×nj)dva-μ0σ∫vni·njdve0＝0 (17)
40、其中，ni,nj为witney型插值基函数，a是磁矢势，e0是棱边中点处的初始电场；将上述式(17)写为ax＝b的形式，如下所示：
41、μ0σ∫vni·njdve0＝∫v(▽×ni)·(▽×nj)dva (18)
42、其中，μ0σ∫vni·njdv为系数矩阵，∫v(▽×ni)·(▽×nj)dva为上述线性方程组的右端项，通过求解上式，得到棱边中点处的初始电场e0。
43、优选的，在步骤s2中，求取电性源初始电场，具体过程如下：
44、当发射波形为阶跃波且发射源为电性源时，在t＝0时刻之前，由于电极持续供电，会在介质中产生稳定的直流电场edc；此时，初始电场分为两部分，如下所示：
45、e0＝esource+edc (19)
46、其中，通过欧姆定律计算电场esource，通过电势计算edc。
47、计算点电极产生的电势是一个三维直流电阻率问题，根据欧姆定律的微分形式，电流密度矢量j，如下所示：
48、j(r)＝σ(r)e(r) (20)
49、电场强度e为电势φ的负梯度，如下所示：
50、e(r)＝-▽φ(r) (21)
51、使用总势方法求解上式，在外边界上，对电势φ施加第一类dirichlet边界条件，即
52、
53、总势方法能够保证空气中电场切向分量不为0，dirichlet边界条件能够保证在外边界上电场切向分量为0，因此直流电场edc和初始电场e0在边界上保持一致；
54、使用加权余量法推导标量有限元的电磁场问题，对于每个六面体单元，电势可近似表示为：
55、
56、其中，li是标量基函数，φi是电势，n是单元自由度的个数；
57、由galerkin方法得：
58、
59、将上式写成矩阵表达式，如下所示：
60、aφ＝s (25)
61、其中，s＝∫ωilδ(r-r0)dv。
62、优选的，在步骤s3中，对于一个确定的γi＝0.026ti，首先计算它在所有时间点上对应的子空间阶数，公式如下：
63、
64、其中，为γi在第j个时刻的子空间阶数；
65、然后，计算sai-krylov子空间阶数的最优值，如下所示：
66、
67、其中，对应的γi即为γopt。
68、优选的，在步骤s4中，使用sai-krylov子空间算法时，对如下方程进行m次求解：
69、(m2+γm1)wj＝m2vj (28)
70、将上式写成ax＝b；其中，a＝m2+γm1，x＝wj，b＝m2vj。
71、因此，本发明采用上述一种基于模型降阶算法的瞬变电磁全空间快速正演方法，基于sai-krylov子空间算法的矢量有限元方法，在保证精度的前提下，极大地加快了正演的计算时间，这为瞬变电磁的反演及成像工作奠定了基础，实现对形状复杂模型的三维快速正演。
72、下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。
1.一种基于模型降阶算法的瞬变电磁全空间快速正演方法，其特征在于，包括以下步骤：
2.根据权利要求1所述的一种基于模型降阶算法的瞬变电磁全空间快速正演方法，其特征在于，在步骤s1中，基于sai-krylov子空间方法，构建忽略位移电流的时间域的maxwell控制方程，如下所示：
3.根据权利要求1所述的一种基于模型降阶算法的瞬变电磁全空间快速正演方法，其特征在于，在步骤s1中，瞬变电磁法的控制方程满足第一类边界条件，如下所示：
4.根据权利要求1所述的一种基于模型降阶算法的瞬变电磁全空间快速正演方法，其特征在于，在步骤s1中，基于galerkin法，二次场的控制方程，如下所示：
5.根据权利要求4所述的一种基于模型降阶算法的瞬变电磁全空间快速正演方法，其特征在于，e(t)的近似求解公式，如下所示：
6.根据权利要求1所述的一种基于模型降阶算法的瞬变电磁全空间快速正演方法，其特征在于，在步骤s2中，求取回线源初始电场，具体过程如下：
7.根据权利要求1所述的一种基于模型降阶算法的瞬变电磁全空间快速正演方法，其特征在于，在步骤s2中，求取电性源初始电场，具体过程如下：
8.根据权利要求1所述的一种基于模型降阶算法的瞬变电磁全空间快速正演方法，其特征在于，在步骤s3中，对于一个确定的γi＝0.026ti，首先计算它在所有时间点上对应的子空间阶数，公式如下：
9.根据权利要求1所述的一种基于模型降阶算法的瞬变电磁全空间快速正演方法，其特征在于，在步骤s4中，使用sai-krylov子空间算法时，对如下方程进行m次求解：
