[发明专利]一种频率域航空电磁法2.5维带地形反演方法

摘要

本发明公开了一种频率域航空电磁法2.5维带地形反演方法，包括以下步骤1)定义目标函数，设置迭代次数为i＝0、拟合精度及最大迭代次数，输入初始模型及反演数据；2)进行正演计算，解正演方程KF＝b得到二次磁场Hx和Hz；3)计算拟合误差，如果达到设定精度或最大迭代次数，退出计算，否则继续；4)用拟正演计算雅克比矩阵，得到模型更新步长；5)更新模型参数，mk+1＝mk+Δm。本发明考虑带地形反演，消除了地形的影响，实现了起伏地表条件下的频率域航空电磁法2.5维反演，为山区航空电磁数据处理与解释提供了一种行之有效的计算方法。

主权项

一种频率域航空电磁法2.5维带地形反演方法，其特征在于：包括以下步骤，1)定义目标函数，设置迭代次数为i＝0、拟合精度及最大迭代次数，输入初始模型及反演数据；2)进行正演计算，解正演方程KF＝b得到二次磁场Hx和Hz；3)计算拟合误差，如果达到设定精度或最大迭代次数，退出计算，否则继续；4)用拟正演计算雅克比矩阵，得到模型更新步长；5)更新模型参数，mk+1＝mk+Δm；所述步骤1)中，定义的目标函数如式(1)所示，Φ(m)=12||d(m)-dobs||2+β2||W(m-mref)||2---(1),]]>式(1)中，d为正演模拟得到的数据向量，dobs为观测数据向量，m为模型参数向量，mref为参考模型向量或先验信息模型向量，W为模型光滑度矩阵，β为正则化参数；所述步骤1)中，假定初始模型为mi，代入式(1)并用泰勒展开，对展开后的式(1)线性化，并忽略高阶项可得式(2)，Φ(m)=12||(d+∂d∂mi·Δm)-dobs||2+β2||W(mi-mref)||2---(2),]]>式(2)中，Δm为m‑mi；对式(2)求导并令其等于零，得到式(3)所示的高斯牛顿法模型更新迭代公式，式(3)中，H为近似海森矩阵，g为目标函数的梯度，J为雅克比矩阵或灵敏度矩阵，所述灵敏度矩阵的元素表示如(4)式，所述步骤2)的正演计算中，假设时谐因子为eiωt，将电磁场分解为一次场和二次场，基于二次场的双旋度电场方程表示为式(5)，▽×▽×Es+iωμ0σEs＝‑iωμ0ΔσEp   (5)，式(5)中，Es为二次电场，ω为角频率，μ0为真空中的磁导率，σ为电导率，Ep为背景场，Δσ为总电导率与背景电导率之差，表示为Δσ＝σ‑σp，σp为背景电导率；所述步骤2)中，采用伽辽金加权余量法对式(5)计算，得式(6)，Rk＝∫ΩNk·[▽×▽×Es+iωμσEs+iωμΔσEp]dV＝0   (6)，由于式中Ωe代表一个离散单元，Nk为离散单元数，将式(6)写成式(7)所示的离散形式，Rk=Σe=1Ne[KeEse-iωμσMeEse-iωμΔσMeEpe]=0---(7);]]>所述式(7)中Ke，Me为单元刚度矩阵，表示如式(8)、式(9)所示，Kkle=∫Ωe(▿×Nke)·(▿×Nle)dV---(8),]]>Mkle=∫ΩeNke·NledV---(9),]]>式(8)、(9)中为矢量基函数，式(8)和式(9)采用27点高斯积分进行计算；将单元刚度矩阵分配到全局刚度矩阵，得到式(10)所示的大型线性方程组，KF＝b   (10)，式(10)中，K为稀疏复对称矩阵，F为场值，b为源项；采用简单的狄利克雷边界条件，认为在离异常体足够远的边界处二次异常场已经衰减为零，即在边界处满足式(11)所示，n×E|∂Ω=0---(11),]]>求解得到电场后，磁场用法拉第定律求得，即式(12)所示，H＝(‑iωμ0)‑1▽×E   (12)；所述步骤4)中，采用求解伴随方程的方法来求解雅克比矩阵，其中磁场的计算公式为式(13)所示，H=(-iωμ0)-1Σi=112▿×NieEi---(13),]]>磁场对模型参数的导数为式(14)所示，∂H∂m=(-iωμ0)-1Σi=112▿×Ni∂Ei∂m---(14),]]>磁场对模型参数的导数可以转换成电场对模型参数的导数，由于正演最后得到式(15)所示的大型复数线性方程组，KF＝b   (15)，对式(14)两边同时对m求导数，得式(16)，K∂F∂m=-∂K∂mF+∂b∂m---(16),]]>通过求解式(16)得到电场对模型参数的导数；所述步骤5)中，通过求解线性方程组得到模型参数更新量，即采用伴随方程的方法来求解雅克比矩阵，进而对模型更新方程组求解；所述步骤5)中，采用式(17)更新下一次迭代的模型，mi+1＝mi+αΔm   (17)，式(17)中，Δm为模型更新向量，α为步长，取值范围为0<α≤1；采用最速下降公式式(18)来选取合适的步长，φ(mi+1)＝φ(mi+αΔm)≤φ(mi)+c1α▽φ(mi)Δm   (18)，式中c1为一常量，通常取值为10‑4，开始α取值为1，判断其是否满足式(18)，若满足则取该α值，更新模型，否则α减少为原来的1/2，重复上述步骤，直到满足为止；所述步骤5)中，采用下列方法之一选择模型正则化因子，①在整个反演过程中β固定为一个值；②在每次迭代过程中逐渐减小β值。