一种基于移动脚印驱动策略的航空电磁随机近似反演方法

本发明涉及地球物理航空电磁的反演，具体为一种基于移动脚印驱动策略的航空电磁随机近似反演方法。

背景技术：

1、航空电磁法(aem)是一种高效的地球物理勘探方法，一次勘查可覆盖数百线公里。海量的数据集提供了丰富的地下电性信息，同时也给用于获取高精度电导率图像和精细结构解释的三维(3d)反演带来了巨大挑战。基于梯度下降的确定性反演如高斯牛顿(gn)、非线性共轭梯度(nlcg)和拟牛顿法等可以有效实现三维反演，但其只能提供一个地下电性结构结果，需要更多额外信息进行结果可靠性评价。基于统计的概率类反演如贝叶斯方法可以提供结果的不确定信息，但其大量的模型统计计算使得这类方法目前无法实现三维运算。如何快速地获取具有可靠性信息的三维反演结果对航空电磁高精度精细结构解释具有重要意义。

2、在3d确定性反演中，正演模拟和雅可比矩阵或梯度被认为是最耗时的两个部分。首先，在时域aem正演模拟中，需要对所有测量点进行数百个时间步长的运算。然后，通过采用伴随正演方法来获得灵敏度矩阵或梯度，但这一过程仍然需要求解大量线性方程组。为提高3d反演效率，基于预条件的随机梯度方法结合压缩感知理论的航空电磁三维反演已经取得一定突破，但这些方法仍然只能提供一个地电结构结果，难以进行结果可靠性评价。

3、在随机优化反演中，随机近似(sa)方法对数据进行k组随机欠采样，获得k个模型反演结果并通过平均计算得到最终的地下电性结果，很好地结合了随机优化和确定性反演，避免了仅恢复一个模型带来的偏差。sa方法与随机梯度相比，由于采样的数据较多，仍然面临收敛速度慢、计算量大的问题。wang等(2021)表明，使用压缩感知(cs)方法，30％～40％的随机采样可以加快3d正演模拟约70％，ren等(2023)进一步提出预处理的随机梯度方法实现3d反演，大幅度提高计算效率。然而，当采用传统sa法进行k组数据集随机采样时，可能会导致大量数据参与反演并降低反演速度。同时，由于航空电磁是发射-接收紧凑系统，每个发射源下方的影响范围(移动脚印)远小于测区大小，导致传统的sa中简单k个模型平均获取最终电导率结果已经不再适用。

4、我们提出一种基于移动脚印驱动策略的航空电磁三维反演算法。在正演中基于cs理论和有限体积法快速实现三维正演，这里需要随机采样k组数据，而这k组数据总体分布必须满足预设的采样率条件，从而保证所有测点3d正演响应的重建精度。k组数据集将在每次迭代中进行重新采样。此外，考虑到具有紧凑型的aem系统，sa中k个反演模型的常规平均操作不再适用，为此提出一种基于移动脚印驱动的方法估计最终恢复的电导率模型。主要思想是，我们首先分析每个剖分的网格单元位置，根据数据集的脚印分布判断是否包括了这一网格，进而判断其对应的模型是否被添加到统计数据中，然后从k个模型中得到可接受的m2个模型(m2<＝k)，然后我们使用m2个恢复模型的平均值和标准差计算每次迭代中的最终电导率更新量。

技术实现思路

1、(一)解决的技术问题

2、针对现有技术的不足，本发明提供了一种基于移动脚印驱动策略的航空电磁随机近似反演方法，解决了航空电磁三维反演计算效率低，很难提供不确定性分析的问题。

3、(二)技术方案。

4、本发明为了实现上述目的具体采用以下技术方案：

5、一种基于移动脚印驱动策略的航空电磁随机近似反演方法，其特征在于：包括以下步骤：

6、1)、将航空电磁数据进行天电噪声去除、背景场去除、运动噪声去除等预处理后，筛选可用的航空电磁数据用于反演；

7、2)、对预处理后数据的所有测点进行随机采样和重构，通过重构精度评价获取最佳采样率，并作为预设采样率；

8、3)、使用泊松采样方法以预设的采样率获取随机测点；

9、4)、对3)中获取的采样数据利用随机数生成方式进行重采样，得到每组中包含l个测点的随机k组数据集；

10、5)、对采样点在初始模型下进行有限体积法正演得到随机测点的正演数据，并使用高精度的压缩感知重构得到所有时间道上所有测点的响应，利用重构的响应和实测数据进行目标函数和rms的计算；

11、6)、通过统计方法，可以在每次迭代中获取不重叠的采样测点，并生成列缺失的灵敏度矩阵，这里的列缺失表示某些数据未被选择，因此相应的列为零。这样可以一次计算出针对所有k组数据集的灵敏度矩阵；

12、7)、对灵敏度矩阵进行k次重组，快速计算k组数据集的预条件梯度和海森矩阵，并组装成基于预条件算子的高斯-牛顿k组数据集的反演方程；

13、8)、通过解反演方程得到k组模型，并使用基于移动脚印驱动平均策略形成一个最终的模型更新项；

14、9)、重复执行步骤3～步骤8，直到达到最大迭代次数或者符合终止条件，

15、得到最终的反演参考模型。

16、在步骤2和步骤5中，有限体积法正演计算时，电场的双旋度方程可表示为：

17、

18、其中，eb和es分别代表背景场和二次场的电场强度，t表示时间，σ、ε和μ分别表示电导率、介电常数和磁导率，表示哈密顿算子，时间二阶导数项为位移电流项，此项远小于传导电流，可忽略。

19、对方程(1)进行规则网格有限体积法离散，随后对每个控制体积进行积分，可得到方程(1)的积分形式：

20、

21、其中，表示任意控制体积v的表面积，n代表法向。

22、考虑所有发射源、所有分量和所有时间道，将其整理为统一形式：

23、

24、其中

25、

26、其中nc为航空电磁空间随机采样的测点数量，tn为程序内部计算使用的时间道数，c和d为与网格尺寸相关的参数，和表示由随机采样测点得到的系数、电场值和右端项。

27、三维反演目标函数计算方式可表示为

28、

29、其中第一项和第二项分别是数据拟合项和模型约束项，m＝(m1,m2,…mm)t是m维模型电导率参数向量，d＝(d1,d2,…,dn)t是n维数据向量，且n＝nst×n,nst和n分别表示航空电磁测点数量和时间道个数，γ为正则化因子。

30、rms的计算公式为：

31、

32、其中，εi为数据误差，n为数据点个数。

33、进一步地，所述步骤5中，重构是一个求解最优化的过程，即求解：

34、

35、其中ε代表数据噪声，s为采样矩阵获得的测量结果，ψ为稀疏变换矩阵。即通过l1范数最小化，不断寻找一个更为稀疏的

36、进一步地，所述步骤6中，使用伴随正演求解灵敏度矩阵的过程为：

37、首先对正演方程(3)两侧对模型参数m取微分，经过简单变换，得到：

38、

39、模型响应对模型参数的偏导数作为灵敏度矩阵j：

40、

41、其中，fb和fs分别代表背景场响应和异常体产生的异常场响应。

42、在背景场为一维的半空间或层状介质响应时，不随大地任意一点电导率变化而变化，故：

43、

44、

45、其中，表示列缺失数据对应的正演系数矩阵，是由系数矩阵和二次电场构成的中间参数，l包含了由正演求解的各位置所有时刻电场直接获得接收点处各实测时间道响应的插值过程，相对于传统方法计算的灵敏度矩阵j是按列缺失的。

46、记中间变量故：

47、

48、

49、所以通过解伴随正演方程组(12)便可得到v，再将v套入方程(10)，便可求得灵敏度矩阵的转置。

50、进一步地，所述步骤7中，预条件梯度和海森矩阵可表示为：

51、

52、

53、预条件的高斯-牛顿k组数据集的反演方程表示为：

54、

55、其中，p＝diag[α·γ-1]π(θ)，γ是采样率矩阵，α是权重系数，π是考虑梯度分布和梯度截断的算子矩阵。为第k个数据集的预条件算子，与预条件项关系为：和是根据第k个采样数据集计算的第k个梯度和海森矩阵。

56、进一步地，所述步骤8中，基于脚印驱动平均策略可表示为：

57、

58、

59、m1＝count|mk(i',j',k')∈ft (18)

60、

61、式中，nx、ny、nz分别为地下电导率x、y、z方向的单元数；ft代表反演数据的脚印。

62、进一步地，所述公式(16)～(18)表示：对于每一个离散化的单元，我们需要判断它是否位于第k个数据集的所有测量点的脚印中，如果是，我们将其对应的反演结果添加到有效模型的统计中，否则意味着第k个数据集在这一单元没有灵敏度，则其结果不计入统计，然后我们可以获得对地下电导率有贡献的m1个模型(m1<＝k)。这些m1个电导率可能会有很大差异，因为不同位置的数据具有不同的灵敏度。因此，我们进一步根据电导率的平均值和标准差(公式(19))将电导率值划分为一系列范围，并统计位于每个电导率范围内的模型数量。假设本次计数过程中最大的数为m2，则通过对m2电导率进行平均，即可得到每次最终更新的电导率分布结果，如式(16)和(17)所示。

63、(三)有益效果

64、与现有技术相比，本发明提供了一种基于移动脚印驱动策略的航空电磁随机近似反演方法，具备以下有益效果：

65、本发明基于移动脚印驱动策略的航空电磁随机近似反演方法中，结合了随机优化和确定性反演来进行时域aem数据反演，我们在正演模拟中采用压缩感知(cs)和三维时间域有限体积法，并使用随机采样来获得完全满足预设采样率的k组(k≥1)数据，k集将在每次迭代中重新采样，从而形成k组预条件反演方程并利用共轭梯度法进行求解，得到的k组模型可用于评估结果可靠性。模型更新是通过脚印驱动的平均技术获得的，该技术考虑了k组数据在每次更新中对同一位置的贡献，并使用统计平均算法来实现最终的模型更新。本发明既保留了确定性反演的收敛速度快，也保存了随机优化能进行不确定分析的优势。