基于压缩感知和预条件随机梯度的航空电磁三维反演方法

本发明涉及地球物理电磁法反演，具体为基于压缩感知和预条件随机梯度的航空电磁三维反演方法。

背景技术：

1、我国国土面积近三分之二为地面人员难以接近的无人区，特别是广大西部地区，高山、沙漠和大面积森林覆盖导致这些地区地球物理勘探程度很低。航空电磁勘查技术是基于固定翼飞机、直升机或无人机等搭载平台的地球物理勘查技术。通过机载或吊挂不接地回线作为磁性发射源并采用磁传感器接收地下电磁信号以达到勘探的目的，具有采样速度快、勘探效率高、无需地面人员接近的特点。航空电磁法的这些优势完美地解决了无人区无法勘探的问题，因此发展航空电磁法是大势所趋。

2、航空电磁反演目的是通过在机载接收到的电磁响应来反推地下空间的真实地电模型，并以此辅助地质解释工作者得到更为准确的地质结构判断。对于航空电磁法，国内外学者做了大量的研究，由于航空电磁数据量庞大，数据反演解释方法通常采用成像或一维反演，但是成像精度低，一维反演具有极大的局限性，现如今已无法满足对地下精细结构的刻画。所各种优化算法被引用到三维航空电磁反演问题中，目前流行的包括：高斯牛顿法(mackie and madden，1993；newman and alumbaugh,2000；sasaki et al.,2013)、拟牛顿法(bae et al.,2012；刘云鹤等，2013)和非线性共轭梯度法(刘云鹤等，2013；kamm andpedersen，2014)等。

3、传统三维反演算法需要对测区或目标区所有测点数据进行计算，同时需要大量的伴随正演获得灵敏度矩阵，耗时长、效率低。随着机器学习理论的不断崛起，随机梯度法逐渐受到关注，该方法是在机器学习训练中，随机选取某一个或者一批样本进行训练，代替传统方法全批量计算，在达到训练目的的同时大大提高训练效率。但该方法由于数据量少使得梯度方向并不是最优的，导致反演收敛速度慢，同时会引入梯度噪声，导致其可能反复跳出最优解，而这些问题可以通过预条件技术得到有效改善。因此基于压缩感知和预条件随机梯度的时间域航空电磁三维反演方法可以保证与传统三维反演方法收敛速度相当，同时保证反演精度的同时极大的提升反演计算效率。

技术实现思路

1、本部分的目的在于概述本发明的实施方式的一些方面以及简要介绍一些较佳实施方式。在本部分以及本技术的说明书摘要和发明名称中可能会做些简化或省略以避免使本部分、说明书摘要和发明名称的目的模糊，而这种简化或省略不能用于限制本发明的范围。

2、鉴于上述和/或现有技术中存在的问题，提出了本发明。

3、因此，本发明的目的是提供基于压缩感知和预条件随机梯度的航空电磁三维反演方法，收敛速度可以媲美传统全批量数据三维反演，同时在保证反演精度的同时极大的提升反演计算效率。

4、为解决上述技术问题，根据本发明的一个方面，本发明提供了如下技术方案：

5、基于压缩感知和预条件随机梯度的航空电磁三维反演方法，包括以下步骤：

6、步骤一，筛选航空电磁实测数据中的可用数据用于反演；

7、步骤二，使用泊松圆盘采样法对测点进行随机采样，并对采样点利用有限体积法进行正演计算得到随机测点的预测数据，并使用压缩感知方法对所有测点的预测数据进行重构；

8、步骤三，将步骤一的实测数据和步骤二的预测数据做数据拟合，并结合模型粗糙度构建航空电磁三维正则化反演的目标函数；

9、步骤四，对正演方程两侧求导，利用伴随正演计算雅可比矩阵；

10、步骤五，利用灵敏度矩阵的转置与数据拟合差的乘积，加上模型电导率参数可计算梯度(目标函数一阶导数)，灵敏度矩阵结合模型协方差和数据协方差可得到海森矩阵(目标函数二阶导数)；

11、步骤六，建立预条件算子，形成预条件随机梯度-高斯-牛顿方法反演方程，计算模型更新量，得到新的迭代模型；

12、步骤七，重复执行步骤二到步骤六，直到达到最大迭代次数或者符合终止条件，得到最终的反演结果。

13、作为本发明所述的基于压缩感知和预条件随机梯度的航空电磁三维反演方法的一种优选方案，其中：所述步骤一中，将航空电磁数据进行天电噪声去除、背景场去除、运动噪声去除等预处理后，筛选可用的航空电磁数据用于反演。

14、作为本发明所述的基于压缩感知和预条件随机梯度的航空电磁三维反演方法的一种优选方案，其中：所述步骤二中，泊松圆盘采样将所有测点的区域划分为若干相同半径且不重叠的圆圈区域，在每个圆圈里只取一个点，完成对所有测点的随机采样，避免传统高斯采样中局部样点过于密集或稀疏的问题。针对随机欠采样点，利用有限体积法进行三维正演计算，电场的双旋度方程可表示为：

15、

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

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

18、

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

20、对方程(2)采用无条件稳定的后推欧拉法进行时间离散，考虑随机采样获得的所有发射源、所有分量和所有时间道，将其整理为统一形式：

21、

22、其中

23、

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

25、对随机测点的正演响应计算结束后，采用压缩感知对所有测点的数据进行重构，重构是一个求解最优化的过程，即求解：

26、

27、其中ε代表数据噪声，s为采样矩阵获得的测量结果，ψ为稀疏变换矩阵。即通过l1范数最小化，不断寻找一个更为稀疏的再利用稀疏逆变换得到测区所有测点正演结果。

28、作为本发明所述的基于压缩感知和预条件随机梯度的航空电磁三维反演方法的一种优选方案，其中：所述步骤三中，根据航空电磁观测数据、正演得到的预测电磁数据构建正则化反演的目标函数为：

29、

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

31、为数据拟合项，具体形式为：

32、

33、其中，f(m)表示模型参数m对应正演响应，cd表示航空电磁数据的协方差矩阵，t为转置算符。

34、为模型约束项，具体形式为：

35、

36、其中，m0表示模型的初始参数，cm是用来约束反演迭代过程中模型参数变化的模型协方差。

37、作为本发明所述的基于压缩感知和预条件随机梯度的航空电磁三维反演方法的一种优选方案，其中：所述步骤四中，计算出的灵敏度矩阵是缺失的，计算灵敏度矩阵的过程为：

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

39、

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

41、

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

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

44、

45、

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

47、记中间变量故：

48、

49、

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

51、作为本发明所述的基于压缩感知和预条件随机梯度的航空电磁三维反演方法的一种优选方案，其中：所述步骤五中，航空电磁正则化反演目标函数对模型参数求一次导数得到梯度：

52、

53、其中，为中间变量，且

54、航空电磁正则化反演目标函数对模型参数的二阶导数可得到海森矩阵：

55、

56、其中海森矩阵的第一项对应矩阵是满阵且非正定，求解比较复杂，故在高斯-牛顿反演方法中，该项可忽略，即海森矩阵表达式为：

57、

58、作为本发明所述的基于压缩感知和预条件随机梯度的航空电磁三维反演方法的一种优选方案，其中：所述步骤六中，预条件随机梯度的高斯-牛顿方法的反演方程式为：

59、

60、其中表示随机采样的测点计算出来的随机梯度，p为预条件算子，是一个对角且正定的矩阵，其形式可以表示为：

61、pij＝α/β+uij                     (18)

62、其中β是采样率，α为权重系数，uij为噪声。

63、最后通过求解反演方程(17)得到模型更新量。

64、作为本发明所述的基于压缩感知和预条件随机梯度的航空电磁三维反演方法的一种优选方案，其中：所述步骤七中，在迭代的过程中，若没有满足迭代终止条件，则判断相邻两次迭代数据拟合的均方根误差的减小量是否小于减小量阈值，若小于减小量阈值则更新正则化因子，继续迭代计算。

65、作为本发明所述的基于压缩感知和预条件随机梯度的航空电磁三维反演方法的一种优选方案，其中：所述步骤七中，迭代终止条件包括：(1)数据拟合的均方根误差小于误差阈值；(2)航空电磁正则化反演目标函数收敛；(3)达到预设最大迭代次数，只要满足任何一个迭代终止条件就迭代终止。

66、与现有技术相比：本发明提供的基于压缩感知和预条件随机梯度的时间域航空电磁快速三维反演方法中，反演采用的数值模拟算法为基于压缩感知的三维时间域有限体积法快速正演算法，能够模拟任意电导率结构下的电磁场信号，其具有较快速度且较高精度模拟复杂地质地形结构并节省内存的优点；同时反演方法基于预条件随机梯度的高斯-牛顿法进行求解，并提出了有效的预条件算子估计方法、正则化约束方法及灵敏度求解方法，能够在保证精度和收敛速度的情况下，极大的提升三维反演计算效率。