1.一种基于大地电磁和直流电阻率数据的联合反演方法,其特征在于:包括如下步骤:步骤S1:采集观测点的观测数据,所述观测数据包括大地电磁法观测的视电阻率和相位以及直流电阻率法观测的视电阻率;
步骤S2:根据观测数据确定反演区域,并对所述反演区域进行剖分得到初始网格,并设定所述初始网格对应的反演参数向量的初值;
步骤S3:基于反演参数向量的当前值进行高斯-牛顿反演迭代计算出模型改变量以及线性搜索步长,再根据模型改变量以及线性搜索步长计算出反演参数向量的迭代更新值;
其中,所述高斯-牛顿反演迭代对应的方程、线性搜索步长、以及迭代更新值的公式如下:mk+1=mk+λδmk
式中, Jk为第k次对应的雅克比矩阵J,Jk(MT)为观测数据中大地电磁法数据对第k次迭代对应的反演参数向量的偏导数矩阵,Jk(DC)为观测数据中直流电阻率法数据对第k次迭代对应的反演参数向量的偏导数矩阵,mk、mk+1分别为第k次、k+1次迭代对应的反演参数向量,k=1时,mk为反演参数向量的初值;Wd为数据方差矩阵,Wm为模型误差矩阵,mapr为先验模型,λ为线性搜索步长,δmk为第k次迭代对应的模型改变量,μ为正则化因子,δdk为第k次迭代对应的正演计算值与观测数据的差;
步骤S4:判断是否达到反演终止条件,若未达到,返回步骤S3进行下一次迭代;
若达到,判断是否达到渐进网格反演终止条件,若达到,反演结束;若未达到,对网格单元进行细化剖分并更新反演参数向量值,再返回步骤S3进行高斯-牛顿反演迭代,其中,基于细化后新增的网格单元个数更新反演参数向量中的元素个数,所述渐进网格反演终止条件为网格剖分终止条件。
2.根据权利要求1所述的方法,其特征在于:步骤S4中对网格单元再进行剖分时的剖分规则如下:按照如下公式计算当前每个网格单元对应的模型灵敏度,并对模型灵敏度超过预设阈值的网格单元进行标记,再对标记的网格单元进行剖分;
式中,Sj为第j个网格单元对应的模型灵敏度;Jij为当前雅克比矩阵J中第i个观测点第j个网格单元对应的元素,n为观测点的个数;
雅克比矩阵J如下所示:
式中, 分别表示为第1个、第2个、第n个观测点的观测数据,m1、m2、mM分别表示为第1个、第2个、第M个网格单元对应在反演参数向量中的值。
3.根据权利要求2所述的方法,其特征在于:判断是否达到所述渐进网格反演终止条件的过程为:计算模型灵敏度的方差,若当前方差与前一次的模型灵敏度的方差的变化值在预设变化范围,当前剖分合理,达到渐进网格反演终止条件;否则,未达到渐进网格反演终止条件,继续选择网格单元进行细化剖分。
4.根据权利要求1所述的方法,其特征在于:所述反演参数向量的元素个数与网格单元数相对应,所述初始网格对应的反演参数向量的初值获取过程如下:首先,获取反演区域的电阻率的下限和上限值;
其次,分别在所述下限和上限值的区间范围取值作为网格单元的电阻率值,再按照如下公式计算反演参数向量初值中各个元素值;
式中,mj为第j个网格单元对应在反演参数向量中的元素值,ρj为第j个网格单元对应的电阻率值,ρmin、ρmax分别为预设的反演区域的电阻率的下限和上限值。
5.根据权利要求1所述的方法,其特征在于:所述数据方差矩阵如下所示:式中,ε为确保分母不为零的一个任意正实数,χi为第i个观测点的观测数据方差,αi为第i个观测点对应的数据权重系数。
6.根据权利要求5所述的方法,其特征在于:数据权重系数αi的取值如下:式中,NMT为大地电磁法参与反演的数据个数、NDC为直流电阻率法参与反演的数据个数。
7.根据权利要求1所述的方法,其特征在于:步骤S4中所述反演终止条件为:当前迭代对应的数据均方差误差小于预设阈值或者迭代次数达到预设迭代总数;
RMS为数据均方根误差, 表示第i个观测点的观测数据, 表示当前反演迭代后再进行正演计算得到的第i个观测点对应的数据,n为观测点的个数。
8.根据权利要求1所述的方法,其特征在于:所述高斯-牛顿反演迭代对应的方程是依据目标函数最优化求解得到,所述目标函数如下:式中, 表示目标函数,dobs为反演数据向量;μ为正则化因子; 为数据误差函数, 为模型误差函数。
9.根据权利要求8所述的方法,其特征在于:所述模型误差函数选用最小结构模型约束,如下所示:式中,βs,βx,βy为模型误差三部分之间的比例系数,Ω为反演区域或网格区域,x,y为x.y坐标方向。
10.根据权利要求1所述的方法,其特征在于:随着反演迭代,所述正则化因子μ逐步减小。