CN114547938A - 一种基于有理Krylov子空间的三维多频可控源电磁反演方法及其系统 - Google Patents
一种基于有理Krylov子空间的三维多频可控源电磁反演方法及其系统 Download PDFInfo
- Publication number
- CN114547938A CN114547938A CN202210162406.4A CN202210162406A CN114547938A CN 114547938 A CN114547938 A CN 114547938A CN 202210162406 A CN202210162406 A CN 202210162406A CN 114547938 A CN114547938 A CN 114547938A
- Authority
- CN
- China
- Prior art keywords
- frequency
- krylov subspace
- controllable source
- inversion
- subspace
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 81
- 238000004364 calculation method Methods 0.000 claims abstract description 40
- 230000006870 function Effects 0.000 claims abstract description 34
- 230000008569 process Effects 0.000 claims abstract description 30
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 23
- 238000005259 measurement Methods 0.000 claims abstract description 8
- 239000011159 matrix material Substances 0.000 claims description 62
- 238000010276 construction Methods 0.000 claims description 19
- 230000005684 electric field Effects 0.000 claims description 14
- 239000013598 vector Substances 0.000 claims description 14
- 238000004590 computer program Methods 0.000 claims description 8
- 230000011218 segmentation Effects 0.000 claims description 8
- 230000035945 sensitivity Effects 0.000 claims description 8
- 238000012937 correction Methods 0.000 claims description 4
- 238000012546 transfer Methods 0.000 claims description 4
- 230000017105 transposition Effects 0.000 claims description 4
- 230000009467 reduction Effects 0.000 abstract description 7
- 238000000354 decomposition reaction Methods 0.000 description 10
- 238000004891 communication Methods 0.000 description 7
- 230000000694 effects Effects 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 238000013459 approach Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 239000006185 dispersion Substances 0.000 description 2
- 230000002349 favourable effect Effects 0.000 description 2
- 229910052500 inorganic mineral Inorganic materials 0.000 description 2
- 239000011707 mineral Substances 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 230000002829 reductive effect Effects 0.000 description 2
- 230000003252 repetitive effect Effects 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 2
- 208000005647 Mumps Diseases 0.000 description 1
- 230000002159 abnormal effect Effects 0.000 description 1
- 230000001133 acceleration Effects 0.000 description 1
- 230000004075 alteration Effects 0.000 description 1
- 230000003321 amplification Effects 0.000 description 1
- 230000002547 anomalous effect Effects 0.000 description 1
- 238000003491 array Methods 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 230000002452 interceptive effect Effects 0.000 description 1
- 230000000670 limiting effect Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 208000010805 mumps infectious disease Diseases 0.000 description 1
- 238000003199 nucleic acid amplification method Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 230000036961 partial effect Effects 0.000 description 1
- 230000002093 peripheral effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/12—Simultaneous equations, e.g. systems of linear equations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/17—Function evaluation by approximation methods, e.g. inter- or extrapolation, smoothing, least mean square method
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/30—Assessment of water resources
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)
- Complex Calculations (AREA)
Abstract
本发明公开了一种基于有理Krylov子空间的三维多频可控源电磁反演方法及其系统，该方法包括：构建三维多频可控源电磁正演问题的多频共用Krylov子空间；其中，基于Krylov子空间用以得到正演结果；构建三维多频可控源电磁伴随正演问题中每个频率的Krylov子空间；其中，基于每个频率的Krylov子空间得到每个频率的伴随正演结果；基于反演区域对应的正演结果、伴随正演结果以及观测数据构建反演目标函数及其梯度，应用梯度类算法得到测量区域的反演结果。本发明引入基于有理Krylov子空间的模型降阶算法引入到了可控源正演方程的近似求解以及多频伴随正演方程的求解，提高了反演解释过程的计算效率，节省运算内存。
Description
技术领域
本发明属于勘探地球物理领域，具体涉及一种基于有理Krylov子空间的三维多频可控源电磁反演方法及其系统。
背景技术
近年来，频率域可控源电磁法被广泛应用于矿产勘查，地下水和油气勘探等领域。其中，野外数据的解释大多依赖于一维反演，这不可避免的造成严重偏差。而三维反演解释主要取决于三维正演和伴随正演的计算精度和速度。
其中，三维反演的基础是三维正演。由于三维地质体的复杂性，三维正演的偏微分方程是无法得到解析解的，因此只能使用有限单元法等数值方法求取三维正演问题的数值解。因此，有限单元法的使用，使得空气和大地被剖分成尺度较小的网格单元，这些网格单元的电阻率或电导率是反演的模型参数。对于频率域可控源电磁问题，在进行有限元离散之后会形成一个大型的线性方程组。频率的个数决定了方程的求解次数，如何高效的求解这些大规模方程是限制可控源三维反演的主要因素。由于各个频点之间的大规模方程是独立的，因此大部分三维反演采用基于频点的MPI并行，这种加速策略严重依赖于计算机的硬件水平，随着频点数量的增加，CPU的数量和内存的大小限制了三维反演的计算速度。
因此，本发明致力于从快速计算多频CSEM的正演问题和伴随正演问题来实现多频CSEM的快速反演。
发明内容
本发明的目的是针对三维反演需要大量运算内存，三维反演解释的计算速度问题，提供一种基于有理Krylov子空间的三维多频可控源电磁反演方法及其系统。本发明所述方法从维多频可控源电磁正演和维多频可控源电磁伴随正演两个角度优化了计算效率，引入基于有理Krylov子空间的模型降阶算法引入到了可控源正演方程的近似求解以及多频伴随正演方程的求解，整体上提高了反演解释过程的计算效率，大大节省运算内存。
一方面，本发明提供的一种基于有理Krylov子空间的三维多频可控源电磁反演方法，包括以下步骤：
构建三维多频可控源电磁正演问题的多频共用Krylov子空间；其中，基于所述多频共用Krylov子空间用以得到三维多频可控源电磁正演问题的正演结果；
构建三维多频可控源电磁伴随正演问题中每个频率的Krylov子空间；其中，基于每个频率的Krylov子空间用以得到每个频率的伴随正演结果；
基于反演区域对应的正演结果、伴随正演结果以及观测数据构建反演目标函数及其梯度，应用梯度类算法得到测量区域的反演结果。
本发明所述方法考虑到了现有技术的计算效率问题，将基于有理Krylov子空间的模型降阶算法引入到了可控源正演方程的近似求解，以构建多频率共用的Krylov子空间，进而将大型稀疏的系数矩阵投影到子空间上，得到较小维度的投影矩阵。接着使用一些通用的直接求解器可快速求取投影矩阵的逆问题，这一过程不依赖于频率的个数，一旦构建了Krylov子空间，那么便可以快速得到多个频率的正演方程的近似解，可大大节省运算内存，提高计算效率。
另一方面，考虑到了在反演过程中，常使用伴随正演的方法求取梯度，该过程和正演过程类似，但其右端项是与频率无关的。在进行多频伴随正演问题的求解时，由于不同频率的右端项呈线性无关，不能再构造同一个子空间来进行多频计算，进而本发明针对每个频率分别构建了一个Krylov子空间。
应当理解，由于本发明中每个频率的伴随正演子空间是相互独立的，因此，优选使用OpenMP并行构建多个频率的伴随正演子空间。
进一步可选地，构建所述多频共用Krylov子空间时，若宽频超过了预设阈值，采用频率分段策略将计算频率分段，每个频段构建单个最优化重复极点的所述多频共用Krylov子空间。
本发明考虑到宽频计算问题，即当计算频率范围较大时，需要较大的子空间才能保证所有频率的方程的解都满足误差要求，这不利于子空间的构建。为此，本发明使用频率分段策略将计算频率分段，并优先使用MPI技术结合有理Lanczos算法并行构建每个频段的子空间。
应当理解，预设阈值的取值范围是经验值，依据需求以及运算效率而定。针对宽频分为多少个频段也是依据需求以及运算效率而定，本发明对此不进行具体的限定。
进一步可选地，每个频段对应的所述多频共用Krylov子空间中单个最优化重复极点为最优实数极点，表示为：
式中，ξj,ξ0为第j个极点、第0个极点；ωmin和ωmax分别为一个频段内的频率最小值和频率最大值。
本发明研究发现，按照上述取值，可以使得整个频段范围内的收敛函数的最小值最大，即具有最优收敛率。
进一步可选地，若正演结果为电场强度，基于所述多频共用Krylov子空间得到三维多频可控源电磁正演问题的正演结果的过程为：利用Krylov子空间构造了电场强度的有理Lanczos逼近式，最后基于所述有理Lanczos逼近式求解电场强度，所述有理Lanczos逼近式表示为：
E≈||X||M Vm+1gω(Am+1)e1
式中，
T为矩阵转置符号，E为电场强度，X＝M-1b，M为质量矩阵，b为可控源电磁的有限元线性方程组的右端源项，Vm+1＝[v1,..,vm+1]表示一组Krylov子空间的正交基向量，v1，vm+1为正交基向量Vm+1中的第1个、第m+1个正交基，m为子空间的维数，gω(Am+1)为传递函数，Am+1是A使用正交基Vm+1在Krylov子空间上的正交投影矩阵，投影矩阵为：
A＝M-1C，C为旋度矩阵。
进一步可选地，所述三维多频可控源电磁伴随正演问题表示为：
(C+iωM)x＝LTr
式中，C为旋度矩阵，M为质量矩阵，L为插值算子，T为矩阵转置符号，ω为频率，i为虚数单位，x为伴随正演问题求解的目标量，r表示预测数据与观测数据之间的残差向量；
其中，基于每个频率的Krylov子空间用以得到每个频率的伴随正演结为：将每个频率的伴随正演线性方程组投影到对应的Krylov子空间上得到目标量的最终解。
进一步可选地，所述目标量的最终解的获取过程是采用多次投影校正方法得到最终解，具体为：
将每个频率的伴随正演线性方程组投影到对应的Krylov子空间上得到目标量的初解x0以及相对残差和残差向量r；
再判断相对残差是否满足误差要求，若不满足误差要求，将得到的残差向量替代右端项LTr，继续构造Krylov子空间进行求解得到新解x1，并修正目标量的解x＝x0+x1，进而更新残差r1和相对残差，重复上述过程，直至满足误差要求。
进一步可选地，所述目标函数和目标函数的梯度表示为：
其中，
为目标函数，
为目标函数的梯度，
为数据拟合项，
为模型拟合项，数据拟合项的梯度
表示为：
其中，m为反演的模型参数，r＝DD(dobs-dpre)*，表示预测数据与观测数据之间的残差；
为灵敏度矩阵，dobs为观测数据，dpre为预测数据，灵敏度矩阵与伴随正演结果的关系为：
JTr＝GTx
其中，x为伴随正演问题求解的目标量，参量G满足：
M为质量矩阵，C为旋度矩阵，E为电场强度，ω为频率，i为虚数单位。
第二方面，本发明提供的一种基于所述三维多频可控源电磁反演方法的系统，其包括：
正演问题的Krylov子空间构建模块，用于构建三维多频可控源电磁正演问题的多频共用Krylov子空间；
正演模块，用于基于所述多频共用Krylov子空间得到三维多频可控源电磁正演问题的正演结果；
伴随正演问题的Krylov子空间构建模块，用于构建三维多频可控源电磁伴随正演问题中每个频率的Krylov子空间；
伴随正演模块，用于基于每个频率的Krylov子空间得到每个频率的伴随正演结果；
反演模块，用于基于反演区域对应的正演结果、伴随正演结果以及观测数据构建反演目标函数及其梯度，应用梯度类算法得到测量区域的反演结果。
第三方面，本发明提供一种电子终端，其至少包括：
一个或多个处理器；
存储了一个或多个计算机程序的存储器，所述处理器调用所述计算机程序以实现：
一种基于有理Krylov子空间的三维多频可控源电磁反演方法的步骤。
第四方面，本发明提供一种可读存储介质，其存储了计算机程序，所述计算机程序被处理器调用以实现：
一种基于有理Krylov子空间的三维多频可控源电磁反演方法的步骤。
有益效果
1.本发明提供的一种基于有理Krylov子空间的三维多频可控源电磁反演方法，其从三维多频可控源电磁正演和三维多频可控源电磁伴随正演两个角度优化了计算效率，引入基于有理Krylov子空间的模型降阶算法引入到了可控源正演方程的近似求解，以构建多频率共用的Krylov子空间，进而将大型稀疏的系数矩阵投影到子空间上，得到较小维度的投影矩阵。接着使用一些通用的直接求解器可快速求取投影矩阵的逆问题，这一过程不依赖于频率的个数，一旦构建了Krylov子空间，那么便可以快速得到多个频率的正演方程的近似解，可大大节省运算内存，提高计算效率。另一方面，在进行多频伴随正演问题的求解时，由于不同频率的右端项呈线性无关，不能再构造同一个子空间来进行多频计算，进而本发明针对每个频率分别构建了一个Krylov子空间，使得上述基于有理Krylov子空间的模型降阶算法可以成功引入伴随正演解释过程，进而整体上提高了反演解释过程的计算效率，大大节省运算内存。
2.本发明进一步的优选方案中，考虑到宽频计算问题，即当计算频率范围较大时，需要较大的子空间才能保证所有频率的方程的解都满足误差要求，这不利于子空间的构建。为此，本发明使用频率分段策略将计算频率分段，并优先使用MPI技术结合有理Lanczos算法并行构建每个频段的子空间。
3.本发明进一步的优选方案中，针对每个频段，单极点的选择使得原始复数矩阵变为实数矩阵,只需一次矩阵分解便可快速构建一个频段的Krylov子空间的正交基，大大节省计算内存。
4.本发明进一步的优选方案中，针对多频伴随正演问题的计算，提出多次投影校正算法，在快速求解多个频率的伴随正演结果的基础上，进一步保证计算精度。
综上，本发明所述方法在内存消耗较小的前提下，可以快速实现多频可控源电磁的三维正演和伴随正演计算，进而快速实现多频可控源电磁的快速反演。
附图说明
图1是收敛率函数随频率范围的变化示意图。其中，随着频率范围的增大，远离极点的频率的收敛率接近线性收敛，不利于子空间的构建，因此使用频率分段策略能够提高每一段频率的收敛率。
图2为反演的理论模型。其中，异常体的大小为1km×1km×0.5km，电阻率为10Ωm。将其嵌入在100Ωm的均匀半空间中，中心埋深坐标为(0km，5km，0.55km)。场源的长度为100m，沿着x轴放置，场源中心放置在坐标原点。频率范围为[1,8192]Hz，共30个。
图3为反演的rms下降曲线示意图。
其中，整个正演区域的外边界沿着x,y,z三个方向延伸至±50km，网格单元数为178,176，棱边数为553,834。共441个测点，给正演数据添加3％的高斯噪声作为观测数据进行反演。反演的初始模型设为100Ωm的均匀半空间。图3给出了反演的rms下降情况。共进行了29次NLCG迭代，终止拟合误差为1.04，总耗时仅1.9h，最大内存消耗为25.7G。整个过程共进行了2610次子正演，平均一次子正演仅仅耗时3.3s，而传统方法一次子正演耗时25s，内存占用8.5G，若采用常规的多频并行求解，内存占用将大大增加。可见本发明方法的提速效果非常明显。
图4为反演的结果，即反演结果的简单切片和三维展示，图4(a)为y＝5000m的垂直切片，图4(b)为x＝0m的垂直切片，图4(c)为z＝500m的水平切片。图4(d)为小于50Ωm的三维结果分布。其中，反演结果表明，本发明提供的方法能够获得较好的反演效果，反演结果与真实异常体吻合较好。
图5是本发明提供的一种基于有理Krylov子空间的三维多频可控源电磁反演方法的流程示意图。
具体实施方式
本发明提供的一种基于有理Krylov子空间的三维多频可控源电磁反演方法，是应用于矿产勘查，地下水和油气勘探等领域，反演的过程是建立一个关于地下导电性参数、观测数据的目标函数，进而推演出地下介质的导电性参数。本发明为了提高反演计算的效率以及降低内存运行要求，引入基于有理Krylov子空间的模型降阶算法引入到了可控源正演方程的近似求解以及多频伴随正演方程的求解中。下面将结合实施例对本发明做进一步的说明。
实施例1：
本实施例提供一种基于有理Krylov子空间的三维多频可控源电磁反演方法，包括以下步骤：
构建三维多频可控源电磁正演问题的多频共用Krylov子空间；其中，基于所述多频共用Krylov子空间用以得到三维多频可控源电磁正演问题的正演结果；
构建三维多频可控源电磁伴随正演问题中每个频率的Krylov子空间；其中，基于每个频率的Krylov子空间用以得到每个频率的伴随正演结果；
基于反演区域对应的正演结果、伴随正演结果以及观测数据构建反演目标函数及其梯度，应用梯度类算法得到测量区域的反演结果。
关于三维多频可控源电磁正演问题：
本实施例使用六面体网格进行空间离散，进而得到三维可控源电磁正演问题的有限元线性方程组为：
(C+iωM)E＝-iωb (1)
其中，b为右端源项，M为质量矩阵C为旋度矩阵，E为正演问题中待求的电场强度，i为虚数单位，ω为角频率。
本实施例针对宽频计算问题，使用频率分段策略将计算频率分段，使用MPI技术结合有理Lanczos算法并行构建每个频段的子空间。本实施例中，涉及的计算频率超过3～4个数量级时，在频段正中间分段，分为两个频段。若频率范围本身就很小(1～2个数量级)，则不需要分段。
关于每个频段的子空间，本发明使用Lanczos算法结合单个最优化重复极点构建有理Krylov子空间，表示为：
式中，
表示有理Krylov子空间，k为子空间的维度。J表示Krylov子空间的第j个向量。span{}表示扩张空间。ξ0表示第0个极点，ξj表示第j个极点。从式(5)可知，ξj的个数决定了矩阵方程的分解次数。为了减少矩阵方程的分解，本发明使用单个最优化重复实数极点
(ωmin和ωmax分别为计算频段范围内的最小值和最大值)，则针对该频段只需一次矩阵分解便可构造子空间。继而将原始的系数矩阵投影到子空间上，大大减小系数矩阵的维数。
应当说明，本实施例选择单个最优化重复实数极点是考虑到：在子空间构造之前，极点的选择影响着整个频段的收敛性能。当给定计算频率范围，单极点算法的收敛率函数R为：
其中，ζ为中间变量，ξ为极点参数。当ξ0为最优实数极点时：
时，可以使得整个频率范围内的收敛函数的最小值最大。收敛函数的最小值表示为：
从上式可知，当计算频率范围越宽，c值越大，而收敛函数的最小值Rmin慢慢减小到1。当c＝1时，即当ωmin＝ωmax或ω＝-ξ0时，具有最优收敛率
即频率越靠近极点，收敛的越快。因此，本实施例针对每个频段，分别利用上述公式确定该频段对应的单重复极点。
图1给出了收敛率函数随频率范围的变化，可知随着频率范围的增大，接近线性收敛，需要构建更大的子空间才能保证误差收敛，这不利于子空间的构建。为了提高收敛率，减小子空间的维度，本发明对频率范围进行分段。
为了简化子空间的构造，令式(5)中的A＝M-1C，X＝M-1b，构造传递函数gω(A)＝-iω(A+iωI)-1，I为单位矩阵，由于传递函数gω(A)是关于大型稀疏矩阵的逆矩阵，直接求解相当困难，因此Krylov子空间方法可用于构造此类矩阵的有效逼近。
基于上述Krylov子空间，使用Lanczos算法构建有理Krylov子空间的正交基Vm+1＝[v1,..,vm+1]，其递归表达式如下：
vj+1＝((A-ξjI)-1vj-αjvj-βj-1vj-1)/βj
＝((C-ξjM)-1Mvj-αjvj-βj-1vj-1)/βj (5)
其中，αj＝(Mvj)T(C-ξjM)-1Mvj，β0v0＝0，βj＝||(C-ξjM)-1Mvj-αjvj-βj-1vj-1||M。本发明通过这种方法构建的标准正交矩阵满足：
最终构造了电场强度的有理Lanczos逼近式：
E≈||X||MVm+1gω(Am+1)e1 (6)
式中，
投影矩阵表示为：
Am+1是A使用正交基Vm+1在Krylov子空间上的正交投影矩阵。Am+1的维度远远小于A的维度，可以利用PARDISO或MUMPS等求解器直接求逆，进而大大降低了解方程的运算量。X的计算实质上是一个方程的求解，即MX＝b。虽然涉及M矩阵的求逆，但由于M矩阵是正定对称的稀疏矩阵，使用Pardiso求解器计算常常只需要几秒钟。Vm+1是一个长方阵，在子空间的构造过程中逐步扩大，直接存储在内存当中，方便调用。而得到Am+1后，进而使用Pardiso求解器求解gω(Am+1)。其余的运算均是向量和矩阵的乘积运算，可以使用自编程序、编译器自带的内部函数、MKL库函数等接口进行快速求解。
应当理解，按照上述思路，针对每个频段构建了子空间，进而按照公式(6)进行正演推算得到电场强度。
关于三维多频可控源电磁伴随正演问题：
为了清楚说明伴随正演的必要性，本实施例构建的反演的目标函数
为：
其中，m为反演的模型参数，
为数据拟合项，
为模型拟合项，λ为正则化因子。D为数据协方差矩阵，W为模型协方差矩阵，T*为共轭转置操作符，m0为模型先验信息。λ为正则化因子，用来权衡数据拟合项和模型拟合项的比重。dobs为观测数据，即待反演的野外数据或者是合成数据。dpre为预测数据，即正演问题的解，本实施例使用有理Krylov子空间模型降阶算法来进行快速求解，即式(6)。可以求得目标函数的梯度为：
其中，数据拟合项的梯度表示为：
其中，r＝DD(dobs-dpre)*，表示预测数据与观测数据之间的残差。
为灵敏度矩阵，由于预测数据与模型参数之间是非线性的隐函数关系，很难显式求取灵敏度矩阵。一般在一维反演中通过扰动法求灵敏度矩阵，但这对于三维反演来说灵敏度矩阵J的计算量是巨大的。在NLCG反演中，常使用伴随正演来求取数据拟合项的梯度。伴随正演可表示为如下过程：
(C+iωM)x＝LTr (10)
JTr＝GTx (11)
其中，L为插值算子。一旦给定m，(C+iωM)是有限元分析所得，是已知项，而电场强度E也已经在正演过程中求得，那么
是可以解析求得的。数据拟合项梯度的求取可以转化为先求解式(10)形式的伴随正演。在进行伴随正演计算时，由于不同频率的伴随正演所涉及的线性方程组的右端项相互独立，无法再构建一个所有频率共用的Krylov子空间。因此，本发明提出分别构建每个频率的小维度Krylov子空间，由于每个频率之间完全独立，可以使用OpenMp并行构建。
由于伴随正演过程在可控源电磁正演过程之后，由于子空间的构建过程类似(只是右端项不同)，从而不需要再进行矩阵分解，可以直接使用正演过程中矩阵分解的结果。在得到每个频率对应的子空间后，将每个频率的伴随正演线性方程组投影到对应的小Krylov子空间上，得到初解x0和相对残差和残差向量r0。相对残差errori表示为：
其中，nf为频率的个数，ωi为第i个频率，ri为第i个频率的数据残差，
为第i个频率方程(10)的模型降阶结果，即对应公式(10)的x。求解过程为：
其中，m为子空间的维数，X＝M-1LTri，A＝M-1C，gω(A)＝(A+iωI)-1，
T表矩阵转置。投影矩阵为：
若相对残差满足误差要求，则得到该频率的伴随正演结果
否则将残差向量代替右端项LTr，继续构造子空间，进而得到新解x1，目标解被修正为x＝x0+x1，得到残差r1和相对残差。如此往复，直到相对残差满足要求，即可得到所有频率的伴随正演结果。
关于反演过程
按照前述方法得到正演和伴随正演的结果，即可得到目标函数和目标函数的梯度，那么梯度类算法便可以用来求解非线性反问题的极值问题，本实施例选择经典的非线性共轭梯度算法来极小化目标函数。
关于本实施例的模拟过程，其参数设定以及系统配置为：
(1)对计算频率进行分段。如果涉及的计算频率超过3～4个数量级时，在频段正中间分段，分为两个频段。若频率范围本身就很小(1～2个数量级)，则不需要分段。
(2)MPI并行构建子空间。首先将主程序命名为CPU0，在CPU0上进行数据整合及构建第一个频段Krylov子空间。在CPU1上进行第二个频段的子空间构建。整个过程只进行两次矩阵分解，并且将矩阵分解的结果保存，以便伴随正演直接使用。
(3)OpenMP并行构建伴随正演子空间。同理，按照正演的频率分段，将伴随正演的频率分别分配到CPU0和CPU1上，伴随正演虽然需要构造多个频率的子空间，但仍然只需两次矩阵分解，且矩阵分解的结果已经在正演获得。
(4)多次投影校正算法。当频率较多时，伴随正演的过大的子空间需要占用大量的内存，因此只构造5-20维的子空间以节省内存。模拟分析表明，本方法能够在内存消耗较小的前提下，快速获得三维CSEM多频反演的结果。
实施例2：
本实施例提供一种基于所述三维多频可控源电磁反演方法的系统，其包括：
正演问题的Krylov子空间构建模块，用于构建三维多频可控源电磁正演问题的多频共用Krylov子空间；
正演模块，用于基于所述多频共用Krylov子空间得到三维多频可控源电磁正演问题的正演结果；
伴随正演问题的Krylov子空间构建模块，用于构建三维多频可控源电磁伴随正演问题中每个频率的Krylov子空间；
伴随正演模块，用于基于每个频率的Krylov子空间得到每个频率的伴随正演结果；
反演模块，用于基于反演区域对应的正演结果、伴随正演结果以及观测数据构建反演目标函数及其梯度，应用梯度类算法得到测量区域的反演结果。
其中，各个单元模块的具体实现过程请参照前述方法的对应过程。应当理解，上述单元模块的具体实现过程参照方法内容，本发明在此不进行具体的赘述，且上述功能模块单元的划分仅仅是一种逻辑功能的划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。同时，上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实现。
实施例3：
本实施例提供一种电子终端，其包括：一个或多个处理器以及存储了一个或多个计算机程序的存储器，存储器中的计算机程序被处理器调用以实现：一种基于有理Krylov子空间的三维多频可控源电磁反演方法的步骤。
该终端还包括：通信接口，用于与外界设备进行通信，进行数据交互传输。
其中，存储器可能包含高速RAM存储器，也可能还包括非易失性除颤器，例如至少一个磁盘存储器。
如果存储器、处理器和通信接口独立实现，则存储器、处理器和通信接口可以通过总线相互连接并完成相互间的通信。所述总线可以是工业标准体系结构总线，外部设备互联总线或扩展工业标准体系结构总线等。所述总线可以分为地址总线、数据总线、控制总线等。
可选的，在具体实现上，如果存储器、处理器和通信接口集成在一块芯片上，则存储器、处理器即通信接口可以通过内部接口完成相互之间的通信。
各个步骤的具体实现过程请参照前述方法的阐述。
应当理解，在本发明实施例中，所称处理器可以是中央处理单元(CentralProcessing Unit，CPU)，该处理器还可以是其他通用处理器、数字信号处理器(DigitalSignal Processor，DSP)、专用集成电路(Application Specific Integrated Circuit，ASIC)、现成可编程门阵列(Field-Programmable Gate Array，FPGA)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。存储器可以包括只读存储器和随机存取存储器，并向处理器提供指令和数据。存储器的一部分还可以包括非易失性随机存取存储器。例如，存储器还可以存储设备类型的信息。
实施例4：
本实施例提供一种可读存储介质，其存储了计算机程序，所述计算机程序被处理器调用以实现：一种基于有理Krylov子空间的三维多频可控源电磁反演方法的步骤。
各个步骤的具体实现过程请参照前述方法的阐述。
所述可读存储介质为计算机可读存储介质，其可以是前述任一实施例所述的控制器的内部存储单元，例如控制器的硬盘或内存。所述可读存储介质也可以是所述控制器的外部存储设备，例如所述控制器上配备的插接式硬盘，智能存储卡(Smart Media Card，SMC)，安全数字(Secure Digital，SD)卡，闪存卡(Flash Card)等。进一步地，所述可读存储介质还可以既包括所述控制器的内部存储单元也包括外部存储设备。所述可读存储介质用于存储所述计算机程序以及所述控制器所需的其他程序和数据。所述可读存储介质还可以用于暂时地存储已经输出或者将要输出的数据。
基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分，或者该技术方案的全部或部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的可读存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。
需要强调的是，本发明所述的实例是说明性的，而不是限定性的，因此本发明不限于具体实施方式中所述的实例，凡是由本领域技术人员根据本发明的技术方案得出的其他实施方式，不脱离本发明宗旨和范围的，不论是修改还是替换，同样属于本发明的保护范围。
Claims (10)
1.一种基于有理Krylov子空间的三维多频可控源电磁反演方法，其特征在于：包括以下步骤：
构建三维多频可控源电磁正演问题的多频共用Krylov子空间；其中，基于所述多频共用Krylov子空间用以得到三维多频可控源电磁正演问题的正演结果；
构建三维多频可控源电磁伴随正演问题中每个频率的Krylov子空间；其中，基于每个频率的Krylov子空间用以得到每个频率的伴随正演结果；
基于反演区域对应的正演结果、伴随正演结果以及观测数据构建反演目标函数及其梯度，应用梯度类算法得到测量区域的反演结果。
2.根据权利要求1所述的三维多频可控源电磁反演方法，其特征在于：构建所述多频共用Krylov子空间时，若宽频超过了预设阈值，采用频率分段策略将计算频率分段，针对每个频段构建单个最优化重复极点的所述多频共用Krylov子空间。
3.根据权利要求2所述的三维多频可控源电磁反演方法，其特征在于：每个频段对应的所述多频共用Krylov子空间中单个最优化重复极点为最优实数极点，表示为：
式中，ξj,ξ0为第j个极点、第0个极点；ωmin和ωmax分别为一个频段内的频率最小值和频率最大值。
4.根据权利要求1所述的三维多频可控源电磁反演方法，其特征在于：若正演结果为电场强度，基于所述多频共用Krylov子空间得到三维多频可控源电磁正演问题的正演结果的过程为：利用Krylov子空间构造了电场强度的有理Lanczos逼近式，最后基于所述有理Lanczos逼近式求解电场强度，所述有理Lanczos逼近式表示为：
E≈||X||MVm+1gω(Am+1)e1
式中，
T为矩阵转置符号，E为电场强度，X＝M-1b，M为质量矩阵，b为可控源电磁的有限元线性方程组的右端源项，Vm+1＝[v1,..,vm+1]表示一组Krylov子空间的正交基向量，v1，vm+1为正交基向量Vm+1中的第1个、第m+1个正交基，m为子空间的维数，gω(Am+1)为传递函数，Am+1是A使用正交基Vm+1在Krylov子空间上的正交投影矩阵，投影矩阵为：
A＝M-1C，C为旋度矩阵。
5.根据权利要求1所述的三维多频可控源电磁反演方法，其特征在于：所述三维多频可控源电磁伴随正演问题表示为：
(C+iωM)x＝LTr
式中，C为旋度矩阵，M为质量矩阵，L为插值算子，T为矩阵转置符号，ω为频率，i为虚数单位，x为伴随正演问题求解的目标量，r表示预测数据与观测数据之间的残差向量；
其中，基于每个频率的Krylov子空间用以得到每个频率的伴随正演结为：将每个频率的伴随正演线性方程组投影到对应的Krylov子空间上得到目标量的最终解。
6.根据权利要求5所述的三维多频可控源电磁反演方法，其特征在于：所述目标量的最终解的获取过程是采用多次投影校正方法得到最终解，具体为：
将每个频率的伴随正演线性方程组投影到对应的Krylov子空间上得到目标量的初解x0以及相对残差和残差向量r；
再判断相对残差是否满足误差要求，若满足误差要求，得到的目标量为对应频率下的伴随正演结果，若不满足误差要求，将得到的残差向量替代右端项LTr，继续构造Krylov子空间进行求解得到新解x1，并修正目标量的解x＝x0+x1，进而更新残差r1和相对残差，重复上述过程，直至满足误差要求。
7.根据权利要求1所述的三维多频可控源电磁反演方法，其特征在于：所述目标函数和目标函数的梯度表示为：
其中，
为目标函数，
为目标函数的梯度，
为数据拟合项，
为模型拟合项，数据拟合项的梯度
表示为：
其中，m为反演的模型参数，r＝DD(dobs-dpre)*，表示预测数据与观测数据之间的残差；
为灵敏度矩阵，dobs为观测数据，dpre为预测数据，灵敏度矩阵与伴随正演结果的关系为：
JTr＝GTx
其中，x为伴随正演问题求解的目标量，参量G满足：
M为质量矩阵，C为旋度矩阵，E为电场强度，ω为频率，i为虚数单位。
8.一种基于权利要求1-7任一项所述三维多频可控源电磁反演方法的系统，其特征在于：包括：
正演问题的Krylov子空间构建模块，用于构建三维多频可控源电磁正演问题的多频共用Krylov子空间；
正演模块，用于基于所述多频共用Krylov子空间得到三维多频可控源电磁正演问题的正演结果；
伴随正演问题的Krylov子空间构建模块，用于构建三维多频可控源电磁伴随正演问题中每个频率的Krylov子空间；
伴随正演模块，用于基于每个频率的Krylov子空间得到每个频率的伴随正演结果；
反演模块，用于基于反演区域对应的正演结果、伴随正演结果以及观测数据构建反演目标函数及其梯度，应用梯度类算法得到测量区域的反演结果。
9.一种电子终端，其特征在于：其至少包括：
一个或多个处理器；
存储了一个或多个计算机程序的存储器，所述处理器调用所述计算机程序以实现：
权利要求1-7任一项所述基于有理Krylov子空间的三维多频可控源电磁反演方法的步骤。
10.一种可读存储介质，其特征在于：其存储了计算机程序，所述计算机程序被处理器调用以实现：
权利要求1-7任一项所述基于有理Krylov子空间的三维多频可控源电磁反演方法的步骤。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210162406.4A CN114547938B (zh) | 2022-02-22 | 一种基于有理Krylov子空间的三维多频可控源电磁反演方法及其系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210162406.4A CN114547938B (zh) | 2022-02-22 | 一种基于有理Krylov子空间的三维多频可控源电磁反演方法及其系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114547938A true CN114547938A (zh) | 2022-05-27 |
| CN114547938B CN114547938B (zh) | 2024-10-29 |
Family
ID=
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115099089A (zh) * | 2022-06-23 | 2022-09-23 | 中国人民解放军国防科技大学 | 均匀背景下的te极化快速互相关对比源电磁反演方法 |
| CN115114821A (zh) * | 2022-06-23 | 2022-09-27 | 中国人民解放军国防科技大学 | 均匀背景下的三维快速互相关对比源电磁反演方法 |
| CN115146573A (zh) * | 2022-07-05 | 2022-10-04 | 北京华大九天科技股份有限公司 | 模拟电路中小信号频率响应特性的仿真方法及相关产品 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20090083006A1 (en) * | 2007-09-20 | 2009-03-26 | Randall Mackie | Methods and apparatus for three-dimensional inversion of electromagnetic data |
| CN107305600A (zh) * | 2016-04-21 | 2017-10-31 | 新疆维吾尔自治区煤炭科学研究所 | 最小二乘法电阻率三维近似反演技术 |
| CN110058315A (zh) * | 2019-05-29 | 2019-07-26 | 中南大学 | 一种三维各向异性射频大地电磁自适应有限元正演方法 |
| CN110826283A (zh) * | 2019-11-15 | 2020-02-21 | 中南大学 | 一种预处理器及基于此预处理器的三维有限差分电磁正演计算方法 |
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20090083006A1 (en) * | 2007-09-20 | 2009-03-26 | Randall Mackie | Methods and apparatus for three-dimensional inversion of electromagnetic data |
| CN107305600A (zh) * | 2016-04-21 | 2017-10-31 | 新疆维吾尔自治区煤炭科学研究所 | 最小二乘法电阻率三维近似反演技术 |
| CN110058315A (zh) * | 2019-05-29 | 2019-07-26 | 中南大学 | 一种三维各向异性射频大地电磁自适应有限元正演方法 |
| CN110826283A (zh) * | 2019-11-15 | 2020-02-21 | 中南大学 | 一种预处理器及基于此预处理器的三维有限差分电磁正演计算方法 |
Non-Patent Citations (1)
| Title |
|---|
| 张继锋等: "三维陆地可控源电磁法有限元模型降阶快速正演", 地球物理学报, vol. 63, no. 9, 30 September 2020 (2020-09-30), pages 3520 - 3533 * |
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115099089A (zh) * | 2022-06-23 | 2022-09-23 | 中国人民解放军国防科技大学 | 均匀背景下的te极化快速互相关对比源电磁反演方法 |
| CN115114821A (zh) * | 2022-06-23 | 2022-09-27 | 中国人民解放军国防科技大学 | 均匀背景下的三维快速互相关对比源电磁反演方法 |
| CN115114821B (zh) * | 2022-06-23 | 2024-04-05 | 中国人民解放军国防科技大学 | 均匀背景下的三维快速互相关对比源电磁反演方法 |
| CN115099089B (zh) * | 2022-06-23 | 2024-04-12 | 中国人民解放军国防科技大学 | 均匀背景下的te极化快速互相关对比源电磁反演方法 |
| CN115146573A (zh) * | 2022-07-05 | 2022-10-04 | 北京华大九天科技股份有限公司 | 模拟电路中小信号频率响应特性的仿真方法及相关产品 |
| CN115146573B (zh) * | 2022-07-05 | 2024-08-20 | 北京华大九天科技股份有限公司 | 模拟电路中小信号频率响应特性的仿真方法及相关产品 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Kampolis et al. | CFD-based analysis and two-level aerodynamic optimization on graphics processing units | |
| CN103617150B (zh) | 一种基于gpu的大规模电力系统潮流并行计算的系统及其方法 | |
| CN108170639B (zh) | 基于分布式环境的张量cp分解实现方法 | |
| Dang et al. | An efficient graphics processing unit‐based parallel algorithm for pricing multi‐asset American options | |
| Agathos et al. | An adapted deflated conjugate gradient solver for robust extended/generalised finite element solutions of large scale, 3D crack propagation problems | |
| CN114547938A (zh) | 一种基于有理Krylov子空间的三维多频可控源电磁反演方法及其系统 | |
| CN114547938B (zh) | 一种基于有理Krylov子空间的三维多频可控源电磁反演方法及其系统 | |
| Datta et al. | Parallel and large scale matrix computations in control: some ideas | |
| CN109408927A (zh) | 一种基于黑盒传输线模型的二维静磁场并行有限元加速方法 | |
| CN111614078B (zh) | 电力系统小干扰稳定性分析方法、装置、设备及存储介质 | |
| CN113779498A (zh) | 离散傅里叶矩阵重构方法、装置、设备和存储介质 | |
| Cupertino et al. | LU decomposition on GPUs: The impact of memory access | |
| Rodohan et al. | A distributed implementation of the finite difference time‐domain (FDTD) method | |
| Tan et al. | Parallel particle swarm optimization algorithm based on graphic processing units | |
| CN113656976B (zh) | 一种二维磁梯度张量快速数值模拟方法、装置和设备 | |
| Garcia et al. | GPU-accelerated Poincaré map method for harmonic-oriented analyses of power systems | |
| Magaña-Lemus et al. | Periodic steady state solution of power systems by selective transition matrix identification, LU decomposition and graphic processing units | |
| Naik et al. | Implementing the gauss seidel algorithm for solving eigenvalues of symmetric matrices with CUDA | |
| Almeer et al. | The performance of the CUDA implementations of HABC and CFS-PML ABCs on GPU hardware | |
| Carpentieri et al. | Fast block Krylov subspace methods for solving sequences of dense MoM linear systems with multiple right-hand sides | |
| CN114547542A (zh) | 一种三维可控源电磁正演方法及系统 | |
| CN107562690A (zh) | 一种基于gpu的aim‑po混合方法 | |
| Ritter et al. | A stable subgridding finite difference time domain method on multi-GPU cluster | |
| Guo et al. | High performance lattice Boltzmann algorithms for fluid flows | |
| Yamada et al. | Communication avoiding Neumann expansion preconditioner for LOBPCG method: convergence property of exact diagonalization method for Hubbard model |
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 |