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

权利要求书

1.一种频率域航空电磁法2.5维带地形反演方法，其特征在于：包括以下步骤，

1)定义目标函数，设置迭代次数为i＝0、拟合精度及最大迭代次数，输入初始模型及反演数据；

2)进行正演计算，解正演方程KF＝b得到二次磁场H_x和H_z；

3)计算拟合误差，如果达到设定精度或最大迭代次数，退出计算，否则继续；

4)用拟正演计算雅克比矩阵，得到模型更新步长；

5)更新模型参数，m^k+1＝m^k+Δm；

所述步骤1)中，定义的目标函数如式(1)所示，

Φ(m)=12||d(m)-dobs||2+β2||W(m-mref)||2---(1),]]>

式(1)中，d为正演模拟得到的数据向量，d_obs为观测数据向量，m为模型参数向量，m_ref为参考模型向量或先验信息模型向量，W为模型光滑度矩阵，β为正则化参数；

所述步骤1)中，假定初始模型为m_i，代入式(1)并用泰勒展开，对展开后的式(1)线性化，并忽略高阶项可得式(2)，

Φ(m)=12||(d+∂d∂mi·Δm)-dobs||2+β2||W(mi-mref)||2---(2),]]>

式(2)中，Δm为m-m_i；

对式(2)求导并令其等于零，得到式(3)所示的高斯牛顿法模型更新迭代公式，

式(3)中，H为近似海森矩阵，g为目标函数的梯度，J为雅克比矩阵或灵敏度矩阵，所述灵敏度矩阵的元素表示如(4)式，

所述步骤2)的正演计算中，假设时谐因子为e^iωt，将电磁场分解为一次场和二次场，基于二次场的双旋度电场方程表示为式(5)，

▽×▽×E_s+iωμ₀σE_s＝-iωμ₀ΔσE_p (5)，

式(5)中，E_s为二次电场，ω为角频率，μ₀为真空中的磁导率，σ为电导率，E_p为背景场，Δσ为总电导率与背景电导率之差，表示为Δσ＝σ-σ_p，σ_p为背景电导率；

所述步骤2)中，采用伽辽金加权余量法对式(5)计算，得式(6)，

R_k＝∫_ΩN_k·[▽×▽×E_s+iωμσE_s+iωμΔσE_p]dV＝0 (6)，由于式中Ω^e代表一个离散单元，N_k为离散单元数，将式(6)写成式(7)所示的离散形式，

Rk=Σe=1Ne[KeEse-iωμσMeEse-iωμΔσMeEpe]=0---(7);]]>

所述式(7)中K^e，M^e为单元刚度矩阵，表示如式(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ωμ₀)^-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)更新下一次迭代的模型，

m_i+1＝m_i+αΔm (17)，

式(17)中，Δm为模型更新向量，α为步长，取值范围为0<α≤1；

采用最速下降公式式(18)来选取合适的步长，

φ(m_i+1)＝φ(m_i+αΔm)≤φ(m_i)+c₁α▽φ(m_i)Δm (18)，

式中c₁为一常量，通常取值为10^-4，开始α取值为1，判断其是否满足式(18)，若满足则取该α值，更新模型，否则α减少为原来的1/2，重复上述步骤，直到满足为止；

所述步骤5)中，采用下列方法之一选择模型正则化因子，

①在整个反演过程中β固定为一个值；

②在每次迭代过程中逐渐减小β值。