磁法数据多目标粒子群反演方法
【技术领域】
[0001] 本发明具体涉及一种用于磁法数据的多目标粒子群反演方法。
【背景技术】
[0002] 地球物理反演存在不稳定性和多解性，吉洪诺夫正则化方法[1'2]已被广泛接受并 应用于解决此类病态反问题。正则化反演的目标函数由数据拟合、模型约束和正则化因子 构成，引入模型约束和正则化因子可稳定反演过程，降低反演的多解性。其中，正则化因子 是数据拟合和模型约束之间的折衷系数，其取值直接影响到反演结果的好坏，诸多学者对 正则化参数选取问题进行了深入的研究。如何选取最合适的正则化因子，目前仍然是正则 化反演研究的热点与难点之一。
[0003] 从反演目标函数的求解方法来看，地球物理反演方法分为两类：局部搜索法和全 局搜索法。传统的局部搜索法以目标函数梯度信息为指导，线性迭代求解。该类方法稳定 性好、计算效率高，但过于依赖初始模型的选择，不易加入约束条件。蒙特卡罗、模拟退火、 遗传算法、人工神经网络、粒子群优化、蚁群优化、差分进化等经典的全局搜索法则不需要 目标函数梯度信息，采用启发式随机搜索策略直接搜索整个解空间，具有不依赖初始模型、 不易陷入局部最优解等优点。然而，目前该类方法更多地考虑数据拟合，而与正则化反演的 结合较少。由于全局搜索法比局部搜索法更耗时，其正则化因子的选取也更为困难。如何 既解决正则化因子选取难题，又克服初始模型依赖，是地球物理反演领域亟待解决的问题。
【发明内容】
[0004] 有鉴于此，有必要研究一种既能解决正则化因子选取难题，又克服初始模型依赖 的地球物理反演方法。
[0005] -种磁法数据多目标粒子群反演算法，其包括如下步骤：
[0006] S1、在磁化率上下界范围内，随机初始化各粒子的初始位置，各粒子初始速度为 〇 ；
[0007] S2、将每个粒子的位置看作模型向量，计算所述模型向量的理论响应，进而计算得 到每个模型向量的数据拟合和模型约束；
[0008] S3、初始化每个粒子的个体最优，选择与非受控目标函数向量对应的粒子构 成初始多目标反演解集；
[0009] S4、从初始多目标反演解集中选择一个全集最优LEADER ;
[0010] S5、根据全集最优LEADER和个体最优更新每个粒子的速度和位置；
[0011] S6、根据更新后的速度和位置计算粒子的模型向量的理论响应，进而得到新的模 型向量的数据拟合和模型约束；
[0012] S7、根据新的数据拟合和模型约束对每个粒子的个体最优进行更新，
[0013] S8、相应更新初始多目标反演解集；
[0014] S9、如果所述磁法数据多目标粒子群反演算法未收敛且未达到最大迭代次数，返 回步骤S4,进行下一轮迭代；
[0015] S10、直至所述磁法数据多目标粒子群反演算法收敛且达到最大迭代次数，则最终 得到最优多目标反演解集，同时输出最优多目标反演解集在目标函数空间的映射Pareto 前沿；
[0016] S11、根据最优多目标反演解集在目标函数空间的映射Pareto前沿，权衡数据拟 合和模型约束的相对重要程度，自多目标反演解集中选出一个最优解向量，从而得到最优 反演结果。
[0017] 本发明所述的磁法数据多目标粒子群反演算法，通过将正则化反演转换成多目标 最优化问题，先采用全局最优化方法同时求数据拟合和模型约束的多目标反演解集，再权 衡两者的相对重要程度，最后从反演解集中优选出最终反演结果，从而起到正则化因子的 作用，能同时解决正则化因子选取困难和初始模型依赖问题。与传统正则化反演方法只得 到一个反演结果不同，该算法能尽可能多地保留可行解，得到一个反演解集。地球物理研究 人员通过分析该解集，既可更加深入的理解反演的过程，又可灵活地从数据拟合和模型约 束两方面进行权衡与选择，得到优于正则化反演的结果。
【附图说明】
[0018] 图1为二维磁法正演示意图；
[0019] 图2为多目标反演示意图；
[0020] 图3为本发明所述的磁法数据多目标粒子群反演方法的流程图；
[0021 ] 图4为图3中步骤S4的子流程图；
[0022] 图5为图3中步骤S8的子流程图；
[0023] 图6为相邻两棱柱体模型剖面图；
[0024] 图7为相邻两个棱柱体模型无噪声数据正则化反演结果；
[0025] 图8是相邻两棱柱体模型无噪声多目标粒子群反演算法中反演解集的Pareto前 沿；
[0026] 图9是相邻两棱柱体模型无噪声多目标粒子群反演算法中最优粒子的收敛曲线；
[0027] 图10是相邻两棱柱体模型无噪声多目标粒子群反演算法的注重数据拟合的反演 结果；
[0028] 图11是相邻两棱柱体模型无噪声多目标粒子群反演算法的注重模型约束的反演 结果；
[0029] 图12是相邻两棱柱体模型无噪声多目标粒子群反演算法的最终反演结果。
[0030] 具体实施方法
[0031] 为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对 本发明进行进一步详细说明，应当理解，此处所述的具体实施案例仅仅用于解释本发明，并 不用于限定本发明。
[0032] 首先简要回顾经典的正则化反演方法，然后阐述本发明提出的磁法数据多目标粒 子群反演算法。
[0033] (1)现有技术的磁法数据正则化反演方法
[0034] 磁法数据二维物性正则化反演涉及到规则网格单元的正演计算，因此需对正演计 算进行讨论。将地下半空间剖分为规则排列的二维棱柱体单元，任意一个棱柱体单元如图 1所示。当磁化方向与地磁场方向一致时，磁异常计算公式如下：
[0035] F = 2kFe[sin (21) sin β In (r2r3/r4rl) + (cos2I + sin2 β - sin2I) (Φ1-Φ2-Φ3+Φ4)]. (I)
[0036] 其中，F是磁异常，Ff3是地磁场强度，k是磁化率，I是地磁倾角，β是棱柱体走向 与磁北的夹角，^，Φi分别是图中所示的距离和夹角。
[0037] 磁异常向量F可视为离散网格单元磁化率向量m与常量核矩阵A的乘积，因此， (1)式的矩阵形式如下：
[0038] F = Am. (2)
[0039] 其中，m = Qc1,…h，…kn)，Ici是第i个棱柱体的磁化率。
[0040] 反演是找到一个模型m，使其理论磁异常与实测磁异常拟合得最好，因此，数据拟 合目标函数定义如下 [2]:
[0041] f! (m) = I I d-Am | 1min, (3)
[0042] 其中，d是实测磁异常，Am是模型m的理论磁异常。
[0043] 为稳定反演过程，克服多解性，正则化反演引入模型约束f2(m)和正则化因子λ， 构成正则化反演目标函数如下[2]:
[0044] Ρλ (m) = f1 (m) + λ f2 (m) min. (4)
[0045] 最常见的模型约束是模型范数最小：
[0046] f2(m) = I |m| |2, (5)
[0047] 此时，正则化反演的矩阵形式目标函数为：
[0048] Pλ (m) = Π d-Am | |2+ λ | | m | |2 = (d-Am) T (d-Am) + λ mTm. (6)
[0049] 其正则化解为[1°]:
[0050] m = [ΑΤΑ+λ I] 1AtCL (7)
[0051] 公式（4) (6)表明，正则化反演的实质是将数据拟合和模型约束两个泛函最优化 问题转换为单目标最优化问题，其中，正则化因子决定数据拟合和模型约束的相对权重，其 取值直接影响到反演结果的好坏。
[0052] (2)本发明所述磁法数据多目标粒子群反演算法
[0053] 本发明所述磁法数据多目标粒子群反演算法包括如下步骤：
[0054] S1、在磁化率上下界范围内，随机初始化各粒子的初始位置，各粒子初始速度为 〇 ；
[0055] S2、将每个粒子的位置看作模型向量，根据磁异常计算公式计算所述模型向量的 理论响应，所述磁异常计算公式如下：
[0056] F = 2kFe[sin (21) sin β In (r2r3/r4rl) + (cos2I + sin2 β - sin2I) (Φ1-Φ2-Φ3+Φ4)]. (I)
[0057] 其中，F是磁异常向量，Ff3是地磁场强度，k是磁化率，I是地磁倾角，β是棱柱体 走向与磁北的夹角， ri，O1分别是图中所示的距离和夹角；
[0058] 并根据矩阵形式公式计算得到每个模型向量的数据拟合和模型约束，所述矩阵形 式公式为：
[0059] F = Am. (2)
[0060] 其中，A为常量核矩阵；m为离散网格单元磁化率向量；m= d · · · h，· · · kn)，1^是 第i个棱柱体的磁化率；
[0061] S3、初始化每个粒子的个体最优Phe<，选择与非受控目标函数向量对应的粒子构 成初始多目标反演解集；
[0062] S4、从初始多目标反演解集中选择一个全集最优LEADER ;
[0063] S5、根据全集最优LEADER和个体最优pliest吏新每个粒子的速度和位置；
[0064] S6、根据更新后的速度和位置，按照公式公式（1) (2)计算粒子的模型向量的理论 响应，进而按照数据拟合目标函数定义公式及最常见的模型约束是模型范数最小公式计算 得到新的模型向量的数据拟合和模型约束；
[0065] 所述数据拟合目标函数定义公式如下[2]:
[0066] f! (m) = I I d-Am | 1min, (3)
[0067] 其中，d是实测磁异常，Am是模型m的理论磁异常。
[0068] 所述最常见的模型约束是模型范数最小公式：
[0069] f2(m) = I |m| I2, (5)；
[0070] S7、根据新的数据拟合和模型约束对每个粒子的个体最优进行更新，
[0071] S8、相应更新初始多目标反演解集；
[0072] S9、如果所述磁法数据多目标粒子群反演算法未收敛且未达到最大迭代次数，返 回步骤S4,进行下一轮迭代；
[0073] S10、直至所述磁法数据多目标粒子群反演算法收敛且达到最大迭代次数，则最终 得到最优多目标反演解集，同时输出最优多目标反演解集在目标函数空间的映射Pareto 前沿；
[0074] S11、根据最优多目标反演解集在目标函数空间的映射Pareto前沿，权衡数据拟 合和模型约束的相对重要程度，自多目标反演解集中选出一个最优解向量，从而得到最优 反演结果。
[0075] 其中，所述步骤S4包括以下子步骤：
[0076] S41、将数据误差函数和模型约束函数张成的空间划分为IOX 10均匀网格；
[0077] S42、将REP中个粒子按其目标函数值对应到不同网格；
[0078] S43、给每个非空的网格赋予一个权值，该权值等于该网格包含的粒子数目的倒 数；
[0079] S44、以轮盘赌的方式从非空的网格中选择一个包含最少粒子的网格；
[0080] S45、从该网格中随机选择一个粒子作为LEADER。
[0081] 其中，所述步骤S8包括以下子步骤：
[0082] S81、检查所有粒子之间的控制和非控制关系