[发明专利]一种基于数据空间快速计算的联合反演方法

权利要求书

1.一种基于数据空间快速计算的联合反演方法，其特征在于，包括如下步骤：

S1：建立重磁电震联合反演目标函数，所述联合反演目标函数如式(1)：

约束条件：τ(m)＝0

其中：

C_d＝diag[C_d1,C_d2,C_d3,C_d4]，C_m＝diag[C_m1,C_m2,C_m3,C_m4]，

α＝[α₁,α₂,α₃,α₄]，

其中，m为模型参数，m₁,m₂,m₃,m₄分别为电阻率、速度、密度和磁化率模型参数；m₀为先验模型参数，m₀₁,m₀₂,m₀₃,m₀₄分别为电阻率、速度、密度和磁化率先验模型参数；d为观测数据，d₁,d₂,d₃,d₄分别为视电阻率、地震走时、重力异常和磁异常；C_d为观测数据d的数据协方差矩阵，C_d1,C_d2,C_d3,C_d4分别为视电阻率、地震走时、重力异常和磁异常的数据协方差矩阵；C_m为模型参数m的模型协方差矩阵，C_m1,C_m2,C_m3,C_m4分别为电阻率、速度、密度和磁化率的模型协方差矩阵；α为阻尼因子，α₁,α₂,α₃,α₄分别为大地电磁、地震初至波走时、重力和磁法的阻尼因子，f(m)表示正演响应，f₁(m),f₂(m),f₃(m),f₄(m)分别为大地电磁、地震初至波走时、重力和磁法的正演响应，为梯度，τ为交叉梯度，τ_ij为大地电磁、地震初至波走时、重力和磁法两两之间的交叉梯度；为数据拟合项、为模型平滑约束项；T和-1分别为矩阵转置和矩阵求逆；

S2：建立电阻率、速度、密度和磁化率的初始模型：

在笛卡尔坐标系下沿x，z二个坐标轴分别将初始模型空间划分成为N_x，N_z个小的矩形网格单元，间距为Δx(i)(i＝1,...,N_x)，Δz(i)(i＝1,...,N_z)，初始模型均采用均匀半空间模型；

S3：正演计算和雅克比矩阵的求取：

采用Wanamaker提出的三角形剖分有限元算法进行大地电磁正演计算；采用了Singh提出的对任意多边形棱柱进行重磁异常正演计算；采用Colin Zelt改进的有限差分方法求解Eikonal方程进行地震初至波走时正演计算；

通过公式(3)求取雅克比矩阵：

其中，A为正演响应f(m)的雅克比矩阵；

S4：计算交叉梯度偏导数和拉格朗日算子：

通过公式(4)计算交叉梯度偏导数

其中，B为交叉梯度函数τ的导数，τ_ij为大地电磁、地震初至波走时、重力和磁法两两之间的交叉梯度，i,j＝1,2,3,4；

拉格朗日算子计算公式如下：

首先将目标函数中的交叉梯度约束项通过拉格朗日算子法加入到目标函数中，如式(5)：

其中，Ψ为加入约束条件的目标函数，Γ为拉格朗日算子，

然后，对正演响应f(m)和交叉梯度约束条件τ(m)进行泰勒级数展开，如式(6)：

将式(6)代入到式(5)中，得式(7)：

其中，

对公式(7)求极值推导出模型改变量的表达式如式(8)：

Δm＝N^-1·n-N^-1·B^T·Γ (8)

其中：

将公式(8)带入到公式(1)中的交叉梯度约束项中，求取拉格朗日算子Γ，如式(11)：

S5：计算获得模型空间迭代模型：

将式(11)得到的拉格朗日算子反代入到式(8)中，得到模型空间模型改变量式(12)：

Δm_m＝N^-1·n-N^-1·B^T·Γ (12)

之后，获得模型空间迭代模型如式(13)：

m_m＝m₀+Δm_m (13)

S6：计算获得数据空间迭代模型：

对模型空间模型改变量式(12)进行恒等式变换

将式(14)和式(15)带入到式(8)中，得到数据空间模型改变量

Δm_s＝N^-1·n-N^-1·B^T·Γ (16)

之后，获得数据空间迭代模型如式(17)：

m_s＝m₀+Δm_s (17)

S7：对获得的数据空间迭代模型进行正演计算，并与观测数据进行拟合差求取，如果拟合差增大、达到最大迭代次数或者拟合差达到阈值，则停止循环反演，输出该数据空间联合反演模型，否则执行S3，

其中，拟合差通过式(18)计算获得

其中，N为观测数据个数。