本发明属于地质探测领域，具体涉及一种基于混合正则化的航空电磁数据反演方法。
背景技术：
1、随着经济快速发展，对于自然资源的依存度与日俱增，给地球物理勘探技术带来了严峻的挑战。实现对近地表精细化探测是矿产资源开发和地质灾害预测等领域的基础，关系到人类安全和国民经济的可持续发展。航空电磁勘探法适合于在高山，沙漠和湖泊等地形复杂、环境恶劣的地区进行矿产、油气资源及环境和工程勘查，并具有高效、环保和低成本的优点，目前是近地表结构探测的重要勘探方法。
2、在对航空电磁数据反演解释工作提出了新要求下，传统的一维和二维反演解释不适合模拟复杂的地下结构，已无法满足资源探测和高精度构造研究的技术需求，因此必须发展更加适用于复杂三维模型的高精度电磁数据反演算法。
3、目前主流的航空电磁数据反演的方法是基于梯度类优化算法的正则化反演，需要通过构造目标函数求取模型解空间的梯度下降方向，然后经有限次迭代完成收敛并获得反演的最优解。各种梯度类优化算法，主要包括高斯-牛顿法、非线性共轭梯度法和有限内存的bfgs算法等被用于三维航空电磁数据反演。
4、公开号为cn115563791a的专利申请公开了一种基于压缩感知重构的大地电磁数据反演方法，包括：获取随机测点分布采集的大地电磁数据，并将大地电磁数据处理成电磁阻抗数据后，对电磁阻抗数据进行分类整理；采用基于曲波变换的压缩感知重构方式对整理后的每类电磁阻抗数据进行从随机测点分布到规则测点分布的数据重构，得到重构电磁阻抗数据；以重构电磁阻抗数据作为观测电磁数据进行大地电磁数据反演计算。
5、上述专利申请公开的反演方法主要用于大地电磁数据的反演，且在反演时采用的目标函数包括数据拟合项和全变差正则化项。当被用于航空电磁数据反演时，仅采用包含全变差正则化项和数据拟合项的目标函数会造成地下异常体边界分辨率损失严重，并不能满足航空电磁数据的高精度反演需求。
技术实现思路
1、鉴于上述，为了提高对异常体边界的分辨率并且获得清晰的背景电阻率数值，本发明的目的是提供一种基于混合正则化的航空电磁数据反演方法，剪切波系数正则化项和全变差正则化项以及数据拟合项构建目标函数，以提高航空电磁数据的反演精度，实现高精度反演，得到高精度反演结果。
2、为实现上述发明目的，本发明实施例提供的一种基于混合正则化的航空电磁数据反演方法，包括以下步骤：
3、获取航空电磁数据并进行预处理后作为观测电磁数据；
4、根据观测数据分布初始化用于反演的电阻率模型，同时设置剪切波系统参数；
5、依据基于剪切波系统参数生成的剪切波变换算子构建剪切波系数正则化项，并构建包括剪切波系数正则化项、全变差正则化项以及数据拟合项构建航空电磁三维混合正则化反演的目标函数；
6、依据目标函数、观测电磁数据以及电阻率模型进行航空电磁数据反演。
7、优选地，根据观测数据分布初始化用于反演的电阻率模型，包括：
8、根据观测数据分布确定反演数据信息，其中，反演数据信息包括测线个数、测点个数、位置、以及观测数据频率；
9、构建电阻率模型，并根据反演范围采用规则网格单元对电阻率模型进行几何剖分，设置水平和垂直方向上的网格单元个数、单元尺寸、以及每个单元的电阻率值。
10、优选地，所述设置剪切波系统参数，包括：
11、依据电阻率模型的网格规模设置剪切波系统参数，其中，剪切波系统由剪切波母函数组成，剪切波系统参数即为剪切波母函数参数，包括尺度矩阵a确定的分辨率、由剪切矩阵s确定的反向、以及由平移矩阵m确定的位置。
12、优选地，依据基于剪切波系统参数生成的剪切波变换算子表示ws，其具体并定义为集合
13、
14、其中，表示低通滤波算子的展开形式，ψ，表示三个方向的高通滤波算子的展开形式，c表示平移的采样常数，表示低通滤波算子，ψ(ψ；c)，表示三个方向的高通滤波算子；
15、
16、
17、
18、
19、
20、
21、mc＝diag(c1,c2,c2)
22、其中，a2j、为尺度矩阵a的组成分量，j表示尺度，j∈z，z表示整数集，sk、为剪切矩阵s的组成分量，k1和k2表示剪切参数，mc、为平移矩阵m的组成分量，c1和c2表示平移的采样常数，diag表示对角函数，k表示剪切参数的统称，m表示空间位置，ψj,k,m，表示三个方向的高通滤波算子的展开形式。
23、优选地，所述包括剪切波系数正则化项、全变差正则化项以及数据拟合项构建航空电磁三维混合正则化反演的目标函数f，表示为：
24、
25、其中，wd表示数据协方差矩阵，dobs和dpre表示观测电磁数据和预测电磁数据，表示数据拟合项，ws表示剪切波变换算子，λ||ws(m-mref)||1表示剪切波系数正则化项，wm表示差分算子，γ||wm(m-mref)||1表示全变差正则化项，实现相邻单元模型参数之间的约束，mref和m分别是反演参考模型和反演求解的电阻率模型，基于l2范数构建数据拟合项，采用l1范数约束剪切波系数正则化项和全变差正则化项，λ和γ是剪切波变换正则化项和全变差正则化项的平衡参数。
26、优选地，依据目标函数、观测电磁数据以及电阻率模型进行航空电磁数据反演，包括以下步骤：
27、(a)通过一次正演和一次伴随正演计算目标函数的梯度，并根据目标函数的梯度构建高斯牛顿反演方程；
28、(b)利用非线性共轭梯度算法对高斯牛顿反演方程进行求解得到模型更新量；
29、(c)基于模型更新量计算当前迭代的新电阻率模型，计算对新电阻率模型进行正演的预测电磁数据与观测电磁数据的拟合差；
30、(d)依据拟合差判定达到收敛或依据迭代次数判定达到收敛时，输出最终的电阻率模型，否则跳转执行步骤(a)。
31、优选地，第n次迭代的目标函数形式表示为：
32、
33、其中，表示第n次迭代的数据拟合项，u表示数据拟合项对应的变量，u＝wd(dobs-dn-1-jδm)，表示第n次迭代的剪切波系数正则化项，v表示剪切波系数正则化项对应的变量，v＝ws(mn-1+δm-mref)，表示第n次迭代的全变差正则化项，w为全变差正则化项对应的变量，w＝wm(mn-1+δm-mref)，λn和γn表示第n次迭代的权重系数，dn-1表示根据第n-次迭代的电导率模型mn-1计算的预测电磁数据，dobs表示观测电磁数据，δm表示第n次迭代的模型参数更新量，j表示雅克比矩阵，mref表示反演参考模型，wd表示数据协方差矩阵，ws表示剪切波变换算子，wm表示差分算子；
34、通过一次正演和一次伴随正演计算目标函数的梯度表示为：
35、
36、
37、
38、其中，符号表示求导，rd、rs、rm均表示对角矩阵。
39、优选地，根据目标函数的梯度构建高斯牛顿反演方程表示为：
40、
41、其中，符号t表示转置。
42、优选地，模型更新量计算当前迭代的新电阻率模型，表示为：
43、mn＝mn-1+sδmn
44、其中，mn表示第n次迭代的新电阻率模型，s是通过线性搜索确定的比例因子。
45、与现有技术相比，本发明具有的有益效果至少包括：
46、本发明实例提供的基于混合正则化的航空电磁反演方法中，通过基于交错网格的有限差分正演计算得到预测电磁数据；并根据预测电磁数据和观测电磁数据构建数据拟合项，结合剪切波正则化和全变差正则化，构建稀疏域系数正则化项和空间域模型正则化项，形成整体的反演目标函数；采用伴随正演技术实现雅可比矩阵计算和存储，得到目标函数的梯度，可进一步节省计算内存；同时根据高斯牛顿法实现反演模型的迭代更新，最终获得符合观测数据的地质电阻率模型。采用剪切波变换和全变分的混合正则化反演可以显著改善单一正则化反演方法的分辨率不足的问题，并在一定程度上加速反演收敛。
1.一种基于混合正则化的航空电磁数据反演方法，其特征在于，包括以下步骤：
2.根据权利要求1所述的基于混合正则化的航空电磁数据反演方法，其特征在于，根据观测数据分布初始化用于反演的电阻率模型，包括：
3.根据权利要求2所述的基于混合正则化的航空电磁数据反演方法，其特征在于，所述设置剪切波系统参数，包括：
4.根据权利要求2所述的基于混合正则化的航空电磁数据反演方法，其特征在于，依据基于剪切波系统参数生成的剪切波变换算子表示ws，其具体并定义为集合
5.根据权利要求1所述的基于混合正则化的航空电磁数据反演方法，其特征在于，所述包括剪切波系数正则化项、全变差正则化项以及数据拟合项构建航空电磁三维混合正则化反演的目标函数f，表示为：
6.根据权利要求1所述的基于混合正则化的航空电磁数据反演方法，其特征在于，依据目标函数、观测电磁数据以及电阻率模型进行航空电磁数据反演，包括以下步骤：
7.根据权利要求6所述的基于混合正则化的航空电磁数据反演方法，其特征在于，第n次迭代的目标函数形式表示为：
8.根据权利要求7所述的基于混合正则化的航空电磁数据反演方法，其特征在于，根据目标函数的梯度构建高斯牛顿反演方程表示为：
9.根据权利要求6所述的基于混合正则化的航空电磁数据反演方法，其特征在于，模型更新量计算当前迭代的新电阻率模型，表示为：
