CN105785460B - 磁化率反演方法及装置 - Google Patents

Abstract

本发明提供了一种磁化率反演方法及装置，该方法包括：获取全张量磁梯度数据；根据全张量磁梯度数据建立Tikhonov正则化模型；基于CPU和GPU协同并行方式采用共轭梯度算法对正则化模型进行迭代求解，反演得到磁化率。本发明中，对模型和数据进行先验约束，在求解过程中，发展了一种快速收敛的预条件的凸集投影混合共轭梯度算法，并通过CPU与GPU联合计算来提高算法过程中大规模矩阵计算效率。二维和三维的拟合数据实例证明采用全张量磁梯度数据反演的磁化率结果相对由TMI数据计算所得的结果精确度更高。通过本发明中的磁化率反演方法及装置能够反演得到较为精确的磁化率，缓解相关技术中反演得到的磁化率不够精确的问题。

Description

磁化率反演方法及装置

技术领域

本发明涉及磁力勘探和地球物理勘探领域，尤其涉及一种磁化率反演方法及装置。

背景技术

磁测数据在地球物理勘探领域中有着广泛的应用。例如，在矿物勘探中，磁方法可以用于决定体参数；在石油勘探中，磁方法可以用来映射地下的沉积特性和缺点，用来控制沉积盆地的沉积背景。

物理参数的反演，例如计算磁化率和磁化强度，是使用磁场数据的主要科学问题。实际应用中，能够反演目标地理区域的磁化率，得到该目标地理区域的磁化率模型，通过该磁化率模型能够反映该目标地理区域内各个地理位置的磁化率，进而将磁化率符合预设要求的某地理位置选择出来，作为进一步的研究对象。相关技术中通常基于总磁场强度(TMI)数据反演求解磁化率。

然而，由于总磁场强度(TMI)数据包含的磁测数据不够全面，因此导致相关技术中反演得到的磁化率不够精确。

发明内容

有鉴于此，本发明提供了一种磁化率反演方法及装置，能够反演得到较为精确的磁化率，缓解相关技术中反演得到的磁化率不够精确的问题。

第一方面，本发明实施例提供了一种磁化率反演方法，包括：获取全张量磁梯度数据；根据所述全张量磁梯度数据建立Tikhonov正则化模型；基于CPU和GPU协同并行方式采用共轭梯度算法对所述正则化模型进行迭代求解，反演得到磁化率。

结合第一方面，本发明实施例提供了第一方面第一种可能的实施方式，其中，根据所述全张量磁梯度数据建立Tikhonov正则化模型，通过以下公式实现：

J^α(m)＝ρ²(Lm,d)+αΩ(m)

S_m＝W_mW_z

其中，J^α(m)表示Tikhonov正则化模型，ρ(Lm,d)表示定义在数据域的函数，Ω(m)表示定义在参数域的函数，α表示正则化参数，L表示离散化紧算子，m表示磁化率向量，d表示所述全张量磁梯度数据，S_d表示作用于数据的尺度算子，S_m表示作用于模型的尺度算子，k表示所述离散化紧算子L的列号，i表示所述离散化紧算子L的行号，diag表示对角化，W_m表示正则化模型的先验约束，W_z表示作用于模型深度的先验约束，表示大于0的常数。

结合第一方面，本发明实施例提供了第一方面第二种可能的实施方式，其中，采用共轭梯度算法对所述正则化模型进行迭代求解，反演得到磁化率，包括：通过以下公式对所述正则化模型进行迭代求解，反演得到磁化率：

m_k+1＝m_k+τ_kh_k,

τ_k＝arg_τ min J^α(m_k+τh_k),

其中，m表示磁化率向量，k表示迭代次数，τ表示步长，h表示在负梯度方向为初始方向的搜索方向，J表示所述正则化模型，α表示正则化参数，表示所述正则化模型的梯度，β^FR表示共轭方向参数变量。

结合第一方面，本发明实施例提供了第一方面第三种可能的实施方式，其中，采用共轭梯度算法对所述正则化模型进行迭代求解，反演得到磁化率，包括：采用预条件共轭梯度算法对所述正则化模型进行迭代求解，反演得到磁化率；其中，预条件共轭梯度算法通过以下公式实现：

m_k+1＝m_k+τ_kh_k,

τ_k＝arg min J^α(m_k+τh_k),

h_k＝P_g(m₀),if k＝0,

其中，m表示磁化率向量，k表示迭代次数，τ表示步长，h表示在负梯度方向为初始方向的搜索方向，J表示所述正则化模型，α表示正则化参数，P表示预条件矩阵，g表示所述正则化模型的梯度，表示预条件的共轭方向参数变量。

结合第一方面，本发明实施例提供了第一方面第四种可能的实施方式，其中，采用共轭梯度算法对所述正则化模型进行迭代求解，反演得到磁化率，包括：采用混合共轭梯度算法对所述正则化模型进行迭代求解，反演得到磁化率；其中，混合共轭梯度算法通过以下公式实现：

m_k+1＝m_k+τ_kh_k,

τ_k＝arg_τ min J^α(m_k+τh_k),

其中，m表示磁化率向量，k表示迭代次数，τ表示步长，h表示在负梯度方向为初始方向的搜索方向，J表示所述正则化模型，α表示正则化参数，表示所述正则化模型的梯度，β^hybrid表示共轭方向混合参数变量，β^HS表示HS型共轭方向参数变量，β^DY表示DY共轭方向参数变量，g表示所述正则化模型的梯度。

结合第一方面，本发明实施例提供了第一方面第五种可能的实施方式，其中，采用共轭梯度算法对所述正则化模型进行迭代求解，反演得到磁化率，包括：采用预条件混合共轭梯度算法对所述正则化模型进行迭代求解，反演得到磁化率；其中，预条件混合共轭梯度算法通过以下公式实现：

m_k+1＝m_k+τ_kh_k,

τ_k＝arg min J^α(m_k+τh_k),

h_k＝P_g(m₀),if k＝0,

其中，m表示磁化率向量，k表示迭代次数，τ表示步长，h表示在负梯度方向为初始方向的搜索方向，J表示所述正则化模型，α表示正则化参数，P表示预条件矩阵，g表示所述正则化模型的梯度，β^hybrid表示共轭方向混合参数变量，β^HS表示HS型共轭方向参数变量，β^DY表示DY型共轭方向参数变量。

结合第一方面上述实施方式，本发明实施例提供了第一方面第六种可能的实施方式，其中，在基于CPU和GPU协同并行方式采用共轭梯度算法对所述正则化模型进行迭代求解之前，所述方法还包括：通过凸集投影技术限定所述正则化模型的有界域，其中，所述凸集投影技术通过以下公式实现：

m_k+1＝P_Π(m_k+τ_kh_k)

其中，m表示磁化率向量，k表示迭代次数，P_Π表示凸集投影算子，τ表示步长，h表示在负梯度方向为初始方向的搜索方向。

第二方面，本发明实施例提供了一种磁化率反演装置，包括：数据获取模块，用于获取全张量磁梯度数据；模型建立模块，用于根据所述全张量磁梯度数据建立Tikhonov正则化模型；磁化率反演模块，用于基于CPU和GPU协同并行方式采用共轭梯度算法对所述正则化模型进行迭代求解，反演得到磁化率。

结合第二方面，本发明实施例提供了第二方面第一种可能的实施方式，其中，所述磁化率反演模块用于：采用预条件混合共轭梯度算法对所述正则化模型进行迭代求解，反演得到磁化率；其中，预条件混合共轭梯度算法通过以下公式实现：

m_k+1＝m_k+τ_kh_k,

τ_k＝arg min J^α(m_k+τh_k),

h_k＝P_g(m₀),if k＝0,

其中，m表示磁化率向量，k表示迭代次数，τ表示步长，h表示在负梯度方向为初始方向的搜索方向，J表示所述正则化模型，α表示正则化参数，P表示预条件矩阵，g表示所述正则化模型的梯度，β^hybrid表示共轭方向混合参数变量，β^HS表示HS型共轭方向参数变量，β^DY表示DY型共轭方向参数变量。

结合第二方面上述实施方式，本发明实施例提供了第二方面第二种可能的实施方式，其中，所述装置还包括：有界域限定模块，用于通过凸集投影技术限定所述正则化模型的有界域，其中，所述凸集投影技术通过以下公式实现：

m_k+1＝P_Π(m_k+τ_kh_k)

其中，m表示磁化率向量，k表示迭代次数，P_Π表示凸集投影算子，τ表示步长，h表示在负梯度方向为初始方向的搜索方向。

本发明实施例中的磁化率反演方法及装置，首先获取全张量磁梯度数据，然后根据全张量磁梯度数据建立Tikhonov正则化模型，最后基于CPU和GPU协同并行方式采用共轭梯度算法对正则化模型进行迭代求解，反演得到磁化率。由于全张量磁梯度数据与基于总磁场强度(TMI)数据相比，包含的物理信息更加全面，因此通过本实施例中的磁化率反演方法及装置，基于全张量磁梯度数据求解磁化率，能够反演得到较精确的磁化率，缓解相关技术中反演得到的磁化率不够精确的问题。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。

图1示出本发明实施例所提供的磁化率反演方法的流程示意图；

图2示出本发明实施例所提供的利用CPU和GPU并行协同求解磁化率的方法的流程示意图；

图3示出本发明实施例所提供的磁化率反演装置的结构示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

考虑到相关技术中通常基于总磁场强度(TMI)数据反演求解磁化率，由于总磁场强度(TMI)数据包含的磁测数据不够全面，因此导致相关技术中反演得到的磁化率不够精确，本发明提供了一种磁化率反演方法及装置，下面结合实施例进行具体描述。

如图1所示，本发明实施例提供了一种磁化率反演方法，该方法包括以下步骤：

步骤S102，获取全张量磁梯度数据；

步骤S104，根据全张量磁梯度数据建立Tikhonov正则化模型；

步骤S106，基于CPU(Central Processing Unit，中央处理器)和GPU(GraphicsProcessing Unit，图形处理器)协同并行方式采用共轭梯度算法对正则化模型进行迭代求解，反演得到磁化率。

本发明实施例中的磁化率反演方法，首先获取全张量磁梯度数据，然后根据全张量磁梯度数据建立Tikhonov正则化模型，最后基于CPU和GPU协同并行方式采用共轭梯度算法对正则化模型进行迭代求解，反演得到磁化率。由于全张量磁梯度数据与基于总磁场强度(TMI)数据相比，包含的物理信息更加全面，因此通过本实施例中的磁化率反演方法，基于全张量磁梯度数据求解磁化率，能够反演得到精确的磁化率，缓解相关技术中反演得到的磁化率不够精确的问题。

本领域技术人员能够理解，对正则化模型进行求解，反演磁化率，实际解决的就是正则化模型最小化的问题，满足正则化模型的值最小的磁化率，就是最终反演得到的磁化率。

步骤S102中，能够通过相关技术中的设备、器材以及多种方法获取到目标地理区域的全张量磁梯度数据，这里不再赘述。

考虑到介质的磁感应强度与磁化率的关系公式，本实施例中，将磁化率反演问题归结为正则化模型的求解问题。上述步骤S104中，根据全张量磁梯度数据建立Tikhonov正则化模型，通过以下公式(1)至公式(7)实现：

J^α(m)＝ρ²(Lm,d)+αΩ(m) (1)

S_m＝W_mW_z (5)

公式(1)至公式(7)中，J^α(m)表示Tikhonov正则化模型，ρ(Lm,d)表示定义在数据域的函数，Ω(m)表示定义在参数域的函数，α表示正则化参数，L表示离散化紧算子，m表示磁化率向量，d表示全张量磁梯度数据，S_d表示作用于数据的尺度算子，S_m表示作用于模型的尺度算子，k表示离散化紧算子L的列号，i表示离散化紧算子L的行号，diag(.)表示对角化，W_m表示正则化模型的先验约束，W_z表示作用于模型深度的先验约束，由函数w(z)在深度方向离散生成，表示大于0的常数，w(z)表示基于深度的权重函数，z表示埋深，η和z₀均表示非负常量。

上述公式(1)至公式(7)中，正则化参数α>0，优选α＝0.001。离散化紧算子L是M×N的矩阵，磁化率向量m的长度为N，全张量磁梯度数据d的长度为M。在确定作用于模型的尺度算子S_m的过程中，表示将向量m的每个元素插入公式中得到的对角形式的新向量，其中表示大于0的较小的常数，优选取1.0×10^-10。基于深度的权重函数中，η优选2，z₀优选0。基于深度的权重函数可以使∫_V(w(z)m(x,y,z))²dv达到最小。基于深度的权重函数对于磁化率反演至关重要，能够避免磁化率分布集中在表面的情况。在确定作用于数据的尺度算子S_d的过程中，表示将离散化紧算子L的每一行2范数平方后进行对角化形成的新矩阵，这种选择方式可以减小反演计算对观测数据(即测量到的全张量磁梯度数据d)的依赖性，得到一个稳定的算法。

本领域技术人员能够理解，由于W_m和W_z表示正则化模型的先验约束，因此通过上述步骤S104得到的正则化模型与先验约束正则化相关。

本实施例中，通过建立正则化模型，将物理问题转变成数学问题，只需要对正则化模型的最小化问题进行求解，就能够反演得到磁化率。

一种优选的实施方式中，上述步骤S106中，采用共轭梯度算法对正则化模型进行迭代求解，反演得到磁化率，包括：通过以下公式(8)对正则化模型进行迭代求解，反演得到磁化率：

公式(8)中，m表示磁化率向量，k表示迭代次数，τ表示步长，h表示在负梯度方向为初始方向的搜索方向，J表示正则化模型，α表示正则化参数，表示正则化模型的梯度，β^FR表示共轭方向参数变量。需要说明的是，β^FR同样表示Flectcher-Reeves公式。

上述公式(8)示出了采用基础共轭梯度算法求解正则化模型，反演得到磁化率的过程。上述公式(8)中，正则化模型的梯度可以通过公式得到，其中，g(m)表示正则化模型的梯度，α表示正则化参数，m表示磁化率向量，ρ表示定义在数据域的函数，Ω表示定义在参数域的函数。

考虑到通过基础共轭梯度算法求解正则化模型，由于全张量磁梯度数据包含的数据量大，导致计算量大，可能存在算法收敛速度过慢的问题，另一个优选的实施例中，上述步骤S106中，采用共轭梯度算法对正则化模型进行迭代求解，反演得到磁化率，包括：采用预条件共轭梯度算法对正则化模型进行迭代求解，反演得到磁化率；其中，预条件共轭梯度算法通过以下公式(9)实现：

公式(9)中，m表示磁化率向量，k表示迭代次数，τ表示步长，h表示在负梯度方向为初始方向的搜索方向，J表示正则化模型，α表示正则化参数，P表示预条件矩阵，g表示正则化模型的梯度，表示预条件的共轭方向参数变量。

具体地，在基础共轭梯度算法的技术公式上增加预条件技术，定义一个新的变量ν使m＝Sv，能够由公式(8)推导出公式(9)，将公式(9)称为预条件共轭梯度算法，其中P是预条件矩阵并有P＝SS^T。对于适当的正则化参数α>0，可以采用矩阵的对角形式来确定矩阵P。

通过如公式(9)所示的预条件共轭梯度算法求解正则化模型，具有算法收敛速度快的效果，能够提高求解速度。

考虑到通过基础共轭梯度算法求解正则化模型，可能存在步长小、以至于引起欠定反演问题收敛所需的代次数多、导致求解时间过长的问题，另一个优选的实施例中，上述步骤S106中，采用共轭梯度算法对正则化模型进行迭代求解，反演得到磁化率，包括：采用混合共轭梯度算法对正则化模型进行迭代求解，反演得到磁化率；其中，混合共轭梯度算法通过以下公式(10)至公式(13)实现：

公式(10)至公式(13)中，m表示磁化率向量，k表示迭代次数，τ表示步长，h表示在负梯度方向为初始方向的搜索方向，J表示正则化模型，α表示正则化参数，表示正则化模型的梯度，β^hybrid表示共轭方向混合参数变量，β^HS表示HS型(Hestenes-Stiefel)共轭方向参数变量，β^DY表示DY型(Dai-Yuan)共轭方向参数变量，g表示正则化模型的梯度。

通过如公式(10)至公式(13)所示的混合共轭梯度算法，能够在求解正则化模型的过程中避免小步长的出现，尤其防止病态问题中，小步长可能会引起误差累计并导致迭代步数增加才能收敛的问题，从而加快求解速度。

在另一个优选的实施方式中，为了加快共轭梯度算法的收敛速度，并减少求解过程中的迭代次数，上述步骤S106中，采用共轭梯度算法对正则化模型进行迭代求解，反演得到磁化率，包括：采用预条件混合共轭梯度算法对正则化模型进行迭代求解，反演得到磁化率；其中，预条件混合共轭梯度算法通过以下公式(14)至公式(17)实现：

公式(14)至公式(17)中，m表示磁化率向量，k表示迭代次数，τ表示步长，h表示在负梯度方向为初始方向的搜索方向，J表示正则化模型，α表示正则化参数，P表示预条件矩阵，g表示正则化模型的梯度，β^hybrid表示共轭方向混合参数变量，β^HS表示HS型(Hestenes-Stiefel)共轭方向参数变量，β^DY表示DY型(Dai-Yuan)共轭方向参数变量。

通过如公式(14)至公式(17)所示的预条件混合共轭梯度算法求解正则化模型，能够加快共轭梯度算法的收敛速度，并减少求解过程中的迭代次数，达到快速反演磁化率的效果。

本发明实施例提供的方法流程中，为了使反演得到的磁化率有界，保证反演得到的磁化率的值符合物理意义，得到高效的共轭梯度算法，在基于CPU和GPU协同并行方式采用共轭梯度算法对正则化模型进行迭代求解之前，还可以包括以下动作：通过凸集投影技术限定正则化模型的有界域，其中，凸集投影技术通过以下公式(18)至公式(19)实现：

m_k+1＝P_Π(m_k+τ_kh_k) (18)

公式(18)至公式(19)中，m表示磁化率向量，k表示迭代次数，P_Π表示凸集投影算子，τ表示步长，h表示在负梯度方向为初始方向的搜索方向。

具体地，对于计算磁化率问题，假设磁化率m是有界的。在本实施例中，只考虑铁磁性材料，模型磁化率值为正。定义正则化模型J^α(m)＝ρ²(Lm,d)+αΩ(m)的可行域为

对于迭代方法，一旦出现一个迭代点超出可行域Π，就会通过凸集投影算子P_Π投影回可行域中。在有界域Π的凸集投影算子可以通过等式(P_Πx)(t)＝χ_Π(t)x(t)决定，其中，χ_Ω(t)是有界域Π的特征方程。凸集投影算子P_Π作用于在Π域外全为0的所有函数构成的子空间。P_Π(x)的第i个构成为假设当前循环中m_k可用，则下次循环可以通过公式m_k+1＝P_Π(m_k+τ_kh_k)获得，其中，h_k是搜索方向，τ_k是由前面提到的非确定线性搜索计算得到的步长。

通过如公式(18)至公式(19)所示的凸集投影技术限定正则化模型的有界域，能够增加正则化模型的解的物理约束，采用基于投影到凸集的梯度投影算法来解决正则化问题，能够生成高效的共轭梯度算法。

另外，本实施例中，为了保证正则化模型J^α(m)＝ρ²(Lm,d)+αΩ(m)是递减的，在基于CPU和GPU协同并行方式采用共轭梯度算法对正则化模型进行迭代求解之前，还可以包括以下动作：确定共轭梯度算法的迭代步长的搜索条件，将搜索条件作为求解正则化模型的约束条件；其中，搜索条件通过以下公式(20)至公式(21)确定：

公式(20)至公式(21)中，表示目标函数，τ表示步长，k表示迭代次数，γ₁和γ₂均表示正常数，表示的一阶导数，J表示正则化模型，α表示正则化参数，m表示磁化率向量，h表示在负梯度方向为初始方向的搜索方向。

需要说明的是，γ₁和γ₂是两个正常数并满足γ₁<γ₂<1。优选γ₁＝0.4，γ₂＝0.6。由于φ(τ)是下有界的并因此，由中值定理可以得出，存在一个参数τ使公式成立。

本实施例中，基于CPU和GPU协同并行的方式求解磁化率，利用GPU处理大规模数值计算问题，能够加快运算速度。利用本实施例中的方法，将TMI数据、2D和3D的全张量磁梯度数据都用于进行对比，结合大量实验表明全张量磁梯度数据可以得到更高精度的反演结果，并且利用CPU\GPU并行方式可以显著提高计算效率。

本实施例中，上述步骤S106中，基于CPU和GPU协同并行方式采用共轭梯度算法对正则化模型进行迭代求解，反演得到磁化率，包括以下过程1)、2)、3)、4)：

1)给出初始磁化率向量m₀，预条件矩阵P，输入磁梯度张量数据d，设置参数γ₁、γ₂、α，停止迭代求解的阈值ε>0，并设置迭代索引k:＝0；

2)计算g_k＝g(m_k)；如果||g_k||≤ε，停止计算；否则，运行步骤3)；其中，ε表示停止迭代求解的阈值；

3)计算m_k+1＝m_k+τ_kh_k，其中τ_k、h_k由公式(8)、公式(9)、公式(10)或公式(14)决定；

4)设置k:＝k+1，转到过程2)。

能够理解，上述过程1)2)3)4)中，迭代求解的停止标准由公式||g_k||≤ε决定，当满足公式||g_k||≤ε时，确定求解得到磁化率，迭代停止。

本实施例中的方法在实际执行时，可以将正则化模型的建立、正则化模型的求解过程以程序的形式输入计算机，操作人员只需在计算机中输入初始磁化率向量m₀，预条件矩阵P，磁梯度张量数据d，设置γ₁＝0.4，γ₂＝0.6，α＝0.001，ε>0，并设置迭代索引k＝0，然后由计算机自动执行上述方法，求解得到磁化率。

本发明实施例还提供了一种如图2所示的利用CPU和GPU并行协同求解磁化率的方法，如图2所示，该方法包括以下步骤：

步骤S202，在CPU中初始化算法：给出初始磁化率向量m₀，预条件矩阵P，输入磁梯度张量数据d，设置γ₁＝0.4，γ₂＝0.6，α＝0.001，ε>0，并设置迭代索引k＝0；

步骤S204，将初始化得到的初始磁化率向量，预条件矩阵，磁梯度张量数据复制到GPU全局存储器中；

步骤S206，将初始化数据从GPU全局存储器复制到GPU共享存储器中；

步骤S208，开始共轭梯度算法，在GPU端进行算法中大规模矩阵向量的计算、更新操作；

步骤S210，将需要进行逻辑运算的标量数据传回CPU；

步骤S212，判断是否达到迭代停止标准(如达到规定的迭代次数或误差精度)，若达到停止标准转到步骤S214，否则转到步骤S208，继续迭代操作；

步骤S214，将计算得到的最优化结果复制到CPU中，并输出结果。

与相关技术相比，本发明实施例具有以下优点：

使用全张量磁梯度数据来反演感兴趣的参数是相对较新的研究，关键原因之一是全张量磁梯度数据很难获得。本发明实施例中，将Tikhonov正则化模型应用于反演物理参数。为了实现这个目的，考虑对模型和数据进行先验约束。在求解过程中，发展了一种快速收敛的预条件的凸集投影混合共轭梯度算法，并通过CPU与GPU联合计算来提高算法过程中大规模矩阵计算效率。二维和三维的拟合数据实例证明采用全张量磁梯度数据反演的磁化率结果相对由TMI数据计算所得的结果精确度更高。

对应上述提供的磁化率反演方法，本发明实施例还提供了一种磁化率反演装置，该磁化率反演装置用于执行上述磁化率反演方法，因此适用于上述方法的描述同样适用于该磁化率反演装置。

如图3所示，本实施例中的磁化率反演装置包括：

数据获取模块31，用于获取全张量磁梯度数据；

模型建立模块32，用于根据全张量磁梯度数据建立Tikhonov正则化模型；

磁化率反演模块33，用于基于CPU和GPU协同并行方式采用共轭梯度算法对正则化模型进行迭代求解，反演得到磁化率。

本发明实施例中的磁化率反演装置，首先获取全张量磁梯度数据，然后根据全张量磁梯度数据建立Tikhonov正则化模型，最后基于CPU和GPU协同并行方式采用共轭梯度算法对正则化模型进行迭代求解，反演得到磁化率。由于全张量磁梯度数据与基于总磁场强度(TMI)数据相比，包含的物理信息更加全面，因此通过本实施例中的磁化率反演装置，基于全张量磁梯度数据求解磁化率，能够反演得到精确的磁化率，缓解相关技术中反演得到的磁化率不够精确的问题。

考虑到介质的磁感应强度与磁化率的关系公式，上述装置中，模型建立模块32用于：通过以下公式建立Tikhonov正则化模型：

J^α(m)＝ρ²(Lm,d)+αΩ(m)

S_m＝W_mW_z

其中，J^α(m)表示Tikhonov正则化模型，ρ(Lm,d)表示定义在数据域的函数，Ω(m)表示定义在参数域的函数，α表示正则化参数，L表示离散化紧算子，m表示磁化率向量，d表示全张量磁梯度数据，S_d表示作用于数据的尺度算子，S_m表示作用于模型的尺度算子，k表示离散化紧算子L的列号，i表示离散化紧算子L的行号，diag(.)表示对角化，W_m表示正则化模型的先验约束，W_z表示作用于模型深度的先验约束，由函数w(z)在深度方向离散生成，表示大于0的常数，w(z)表示基于深度的权重函数，z表示埋深，η和z₀均表示非负常量。

一种优选的实施方式中，上述装置中，磁化率反演模块33用于：通过以下公式对正则化模型进行迭代求解，反演得到磁化率：

m_k+1＝m_k+τ_kh_k,

τ_k＝arg_τ min J^α(m_k+τh_k),

其中，m表示磁化率向量，k表示迭代次数，τ表示步长，h表示在负梯度方向为初始方向的搜索方向，J表示正则化模型，α表示正则化参数，表示正则化模型的梯度，β^FR表示共轭方向参数变量。

考虑到通过基础共轭梯度算法求解正则化模型，由于全张量磁梯度数据包含的数据量大，导致计算量大，可能存在算法收敛速度过慢的问题，另一个优选的实施例中，磁化率反演模块33用于：采用预条件共轭梯度算法对正则化模型进行迭代求解，反演得到磁化率；其中，预条件共轭梯度算法通过以下公式实现：

m_k+1＝m_k+τ_kh_k,

τ_k＝arg min J^α(m_k+τh_k),

h_k＝P_g(m₀),if k＝0,

其中，m表示磁化率向量，k表示迭代次数，τ表示步长，h表示在负梯度方向为初始方向的搜索方向，J表示正则化模型，α表示正则化参数，P表示预条件矩阵，g表示正则化模型的梯度，表示预条件的共轭方向参数变量。

通过预条件共轭梯度算法求解正则化模型，具有算法收敛速度快的效果，能够提高求解速度。

考虑到通过基础共轭梯度算法求解正则化模型，可能存在步长小、以至于引起欠定反演问题收敛所需的代次数多、导致求解时间过长的问题，另一个优选的实施例中，磁化率反演模块33用于：采用混合共轭梯度算法对正则化模型进行迭代求解，反演得到磁化率；其中，混合共轭梯度算法通过以下公式实现：

m_k+1＝m_k+τ_kh_k,

τ_k＝arg_τ min J^α(m_k+τh_k),

其中，m表示磁化率向量，k表示迭代次数，τ表示步长，h表示在负梯度方向为初始方向的搜索方向，J表示正则化模型，α表示正则化参数，表示正则化模型的梯度，β^hybrid表示共轭方向混合参数变量，β^HS表示HS型(Hestenes-Stiefel)共轭方向参数变量，β^DY表示DY型(Dai-Yuan)共轭方向参数变量，g表示正则化模型的梯度。

通过混合共轭梯度算法，能够在求解正则化模型的过程中避免小步长的出现，尤其防止病态问题中，小步长可能会引起误差累计并导致迭代步数增加才能收敛的问题，从而加快求解速度。

在另一个优选的实施方式中，为了加快共轭梯度算法的收敛速度，并减少求解过程中的迭代次数，磁化率反演模块33用于：采用预条件混合共轭梯度算法对正则化模型进行迭代求解，反演得到磁化率；其中，预条件混合共轭梯度算法通过以下公式实现：

m_k+1＝m_k+τ_kh_k,

τ_k＝arg min J^α(m_k+τh_k),

h_k＝P_g(m₀),if k＝0,

其中，m表示磁化率向量，k表示迭代次数，τ表示步长，h表示在负梯度方向为初始方向的搜索方向，J表示正则化模型，α表示正则化参数，P表示预条件矩阵，g表示正则化模型的梯度，β^hybrid表示共轭方向混合参数变量，β^HS表示HS型(Hestenes-Stiefel)共轭方向参数变量，β^DY表示DY型(Dai-Yuan)共轭方向参数变量。

通过预条件混合共轭梯度算法求解正则化模型，能够加快共轭梯度算法的收敛速度，并减少求解过程中的迭代次数，达到快速反演磁化率的效果。

为了使反演得到的磁化率有界，保证反演得到的磁化率的值符合物理意义，得到高效的共轭梯度算法，上述装置还包括：有界域限定模块，用于通过凸集投影技术限定正则化模型的有界域，其中，凸集投影技术通过以下公式实现：

m_k+1＝P_Π(m_k+τ_kh_k)

其中，m表示磁化率向量，k表示迭代次数，P_Π表示凸集投影算子，τ表示步长，h表示在负梯度方向为初始方向的搜索方向。

通过凸集投影技术限定正则化模型的有界域，能够增加正则化模型的解的物理约束，采用基于投影到凸集的梯度投影算法来解决正则化问题，能够生成高效的共轭梯度算法。

与相关技术相比，本发明实施例具有以下优点：

使用全张量磁梯度数据来反演感兴趣的参数是相对较新的研究，关键原因之一是全张量磁梯度数据很难获得。本发明实施例中，将Tikhonov正则化模型应用于反演物理参数。为了实现这个目的，考虑对模型和数据进行先验约束。在求解过程中，发展了一种快速收敛的预条件的凸集投影共轭梯度算法，并通过CPU与GPU联合计算来提高算法过程中大规模矩阵计算效率。二维和三维的拟合数据实例证明采用全张量磁梯度数据反演的磁化率结果相对由TMI数据计算所得的结果精确度更高。

本发明实施例所提供的磁化率反演装置可以为设备上的特定硬件或者安装于设备上的软件或固件等。本发明实施例所提供的装置，其实现原理及产生的技术效果和前述方法实施例相同，为简要描述，装置实施例部分未提及之处，可参考前述方法实施例中相应内容。所属领域的技术人员可以清楚地了解到，为描述的方便和简洁，前述描述的系统、装置和单元的具体工作过程，均可以参考上述方法实施例中的对应过程，在此不再赘述。

在本发明所提供的实施例中，应该理解到，所揭露装置和方法，可以通过其它的方式实现。以上所描述的装置实施例仅仅是示意性的，例如，所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，又例如，多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些通信接口，装置或单元的间接耦合或通信连接，可以是电性，机械或其它的形式。

所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。

另外，在本发明提供的实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。

所述功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。

应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释，此外，术语“第一”、“第二”、“第三”等仅用于区分描述，而不能理解为指示或暗示相对重要性。

最后应说明的是：以上所述实施例，仅为本发明的具体实施方式，用以说明本发明的技术方案，而非对其限制，本发明的保护范围并不局限于此，尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，其依然可以对前述实施例所记载的技术方案进行修改或可轻易想到变化，或者对其中部分技术特征进行等同替换；而这些修改、变化或者替换，并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围。都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。

Claims (8)

1.一种磁化率反演方法，其特征在于，包括：

获取全张量磁梯度数据；

根据所述全张量磁梯度数据建立Tikhonov正则化模型；

基于CPU和GPU协同并行方式采用共轭梯度算法对所述正则化模型进行迭代求解，反演得到磁化率；

其中，根据所述全张量磁梯度数据建立Tikhonov正则化模型，通过以下公式实现：

J^α(m)＝ρ²(Lm,d)+αΩ(m)

$\mathrm{\&rho;}(Lm,d)=\frac{1}{2}\left|\right|S_d(Lm-d)|{|}_2^2$

$\mathrm{\&Omega;}\left(m\right)=\frac{1}{2}\left|\right|S_{m}m|{|}_2^2$

$S_d=diag(1/\underset{k}{\&Sigma;}{\left(L_{ik}\right)}^2)$

S_m＝W_mW_z

其中，J^α(m)表示Tikhonov正则化模型，ρ(Lm,d)表示定义在数据域的函数，Ω(m)表示定义在参数域的函数，α表示正则化参数，L表示离散化紧算子，m表示磁化率向量，d表示所述全张量磁梯度数据，S_d表示作用于数据的尺度算子，S_m表示作用于模型的尺度算子，k表示所述离散化紧算子L的列号，i表示所述离散化紧算子L的行号，diag表示对角化，W_m表示正则化模型的先验约束，W_z表示作用于模型深度的先验约束，表示大于0的常数；

其中，采用共轭梯度算法对所述正则化模型进行迭代求解，反演得到磁化率，包括：

通过以下公式对所述正则化模型进行迭代求解，反演得到磁化率：

m_k+1＝m_k+τ_kh_k,

τ_k＝arg_τminJ^α(m_k+τh_k),

$h_k=-\&dtri;J^{\mathrm{\&alpha;}}\left(m_k\right)+{\mathrm{\&beta;}}_{k-1}^{FR}{h}_{k-1},$

$\begin{array}{ccc}{h}_k=-\&dtri;J^{\mathrm{\&alpha;}}\left(m_k\right),& if& k=0\end{array},$

${\mathrm{\&beta;}}_{k-1}^{FR}=\left|\right|\&dtri;J^{\mathrm{\&alpha;}}\left(m_k\right)|{|}^2/||\&dtri;J^{\mathrm{\&alpha;}}\left(m_{k-1}\right)|{|}^2$

其中，m表示磁化率向量，k表示迭代次数，τ表示步长，h表示在负梯度方向为初始方向的搜索方向，J表示所述正则化模型，α表示正则化参数，表示所述正则化模型的梯度，β^FR表示共轭方向参数变量。

2.根据权利要求1所述的方法，其特征在于，采用共轭梯度算法对所述正则化模型进行迭代求解，反演得到磁化率，包括：

采用预条件共轭梯度算法对所述正则化模型进行迭代求解，反演得到磁化率；

其中，预条件共轭梯度算法通过以下公式实现：

m_k+1＝m_k+τ_kh_k,

τ_k＝argminJ^α(m_k+τh_k),

$h_k=P_g\left(m_k\right)+{\widetilde{\mathrm{\&beta;}}}_{k-1}{h}_{k-1},$

h_k＝P_g(m₀),if k＝0,

${\widetilde{\mathrm{\&beta;}}}_{k-1}=\left|\right|P^{1/2}g\left(m_k\right)|{|}^2/|\left|P^{1/2}g\left(m_{k-1}\right)\right|{|}^2$

其中，m表示磁化率向量，k表示迭代次数，τ表示步长，h表示在负梯度方向为初始方向的搜索方向，J表示所述正则化模型，α表示正则化参数，P表示预条件矩阵，g表示所述正则化模型的梯度，表示预条件的共轭方向参数变量。

3.根据权利要求1所述的方法，其特征在于，采用共轭梯度算法对所述正则化模型进行迭代求解，反演得到磁化率，包括：

采用混合共轭梯度算法对所述正则化模型进行迭代求解，反演得到磁化率；

其中，混合共轭梯度算法通过以下公式实现：

m_k+1＝m_k+τ_kh_k,

τ_k＝arg_τminJ^α(m_k+τh_k),

$h_k=-\&dtri;J^{\mathrm{\&alpha;}}\left(m_k\right)+{\mathrm{\&beta;}}_{k-1}^{hybrid}{h}_{k-1},$

$\begin{array}{ccc}{h}_k=-\&dtri;J^{\mathrm{\&alpha;}}\left(m_k\right),& if& k=0\end{array},$

${\mathrm{\&beta;}}_k^{hybrid}=\left\{\begin{array}{ccc}0,& if& {\mathrm{\&beta;}}_k^{HS}<0\\ {\mathrm{\&beta;}}_k^{HS},& if& 0\&le;{\mathrm{\&beta;}}_k^{HS}<{\mathrm{\&beta;}}_k^{DY}\\ {\mathrm{\&beta;}}_k^{DY},& if& {\mathrm{\&beta;}}_k^{HS}>{\mathrm{\&beta;}}_k^{DY}\end{array}\right.$

${\mathrm{\&beta;}}_k^{HS}=\frac{g{\left(m_k\right)}^T\left(g\right(m_k)-g(m_{k-1}\left)\right)}{h_{k-1}^T\left(g\right(m_k)-g(m_{k-1}\left)\right)}$

${\mathrm{\&beta;}}_k^{DY}=\frac{g{\left(m_k\right)}^{T}g\left(m_k\right)}{h_{k-1}^T\left(g\right(m_k)-g(m_{k-1}\left)\right)}$

其中，m表示磁化率向量，k表示迭代次数，τ表示步长，h表示在负梯度方向为初始方向的搜索方向，J表示所述正则化模型，α表示正则化参数，表示所述正则化模型的梯度，β^hybrid表示共轭方向混合参数变量，β^HS表示HS型共轭方向参数变量，β^DY表示DY型共轭方向参数变量，g表示所述正则化模型的梯度。

4.根据权利要求1所述的方法，其特征在于，采用共轭梯度算法对所述正则化模型进行迭代求解，反演得到磁化率，包括：

采用预条件混合共轭梯度算法对所述正则化模型进行迭代求解，反演得到磁化率；

其中，预条件混合共轭梯度算法通过以下公式实现：

m_k+1＝m_k+τ_kh_k,

τ_k＝argminJ^α(m_k+τh_k),

$h_k=P_g\left(m_k\right)+{\mathrm{\&beta;}}_{k-1}^{hybrid}{h}_{k-1},$

h_k＝P_g(m₀),if k＝0,

${\mathrm{\&beta;}}_k^{hybrid}=\left\{\begin{array}{ccc}0,& if& {\mathrm{\&beta;}}_k^{HS}<0\\ {\mathrm{\&beta;}}_k^{HS},& if& 0\&le;{\mathrm{\&beta;}}_k^{HS}<{\mathrm{\&beta;}}_k^{DY}\\ {\mathrm{\&beta;}}_k^{DY},& if& {\mathrm{\&beta;}}_k^{HS}>{\mathrm{\&beta;}}_k^{DY}\end{array}\right.$

${\mathrm{\&beta;}}_k^{HS}=\frac{g{\left(m_k\right)}^{T}P\left(g\right(m_k)-g(m_{k-1}\left)\right)}{h_{k-1}^{T}P\left(g\right(m_k)-g(m_{k-1}\left)\right)}$

${\mathrm{\&beta;}}_k^{DY}=\frac{g{\left(m_k\right)}^{T}Pg\left(m_k\right)}{h_{k-1}^{T}P\left(g\right(m_k)-g(m_{k-1}\left)\right)}$

其中，m表示磁化率向量，k表示迭代次数，τ表示步长，h表示在负梯度方向为初始方向的搜索方向，J表示所述正则化模型，α表示正则化参数，P表示预条件矩阵，g表示所述正则化模型的梯度，β^hybrid表示共轭方向混合参数变量，β^HS表示HS型共轭方向参数变量，β^DY表示DY型共轭方向参数变量。

5.根据权利要求1至4任一项所述的方法，其特征在于，在基于CPU和GPU协同并行方式采用共轭梯度算法对所述正则化模型进行迭代求解之前，所述方法还包括：

通过凸集投影技术限定所述正则化模型的有界域，其中，所述凸集投影技术通过以下公式实现：

m_k+1＝P_Π(m_k+τ_kh_k)

$P_{\&Pi;}(m_k+{\mathrm{\&tau;}}_k{h}_k)=\left\{\begin{array}{cc}(m_k+{\mathrm{\&tau;}}_k{h}_k),& (m_k+{\mathrm{\&tau;}}_k{h}_k)\&GreaterEqual;0\\ 0,& (m_k+{\mathrm{\&tau;}}_k{h}_k)<0\end{array}\right.$

其中，m表示磁化率向量，k表示迭代次数，P_Π表示凸集投影算子，τ表示步长，h表示在负梯度方向为初始方向的搜索方向。

6.一种磁化率反演装置，其特征在于，包括：

数据获取模块，用于获取全张量磁梯度数据；

模型建立模块，用于根据所述全张量磁梯度数据建立Tikhonov正则化模型；

磁化率反演模块，用于基于CPU和GPU协同并行方式采用共轭梯度算法对所述正则化模型进行迭代求解，反演得到磁化率；

其中，所述数据获取模块、所述模型建立模块、所述磁化率反演模块依次连接；

所述模型建立模块用于：通过以下公式建立Tikhonov正则化模型：

J^α(m)＝ρ²(Lm,d)+αΩ(m)

$\mathrm{\&rho;}(Lm,d)=\frac{1}{2}\left|\right|S_d(Lm-d)|{|}_2^2$

$\mathrm{\&Omega;}\left(m\right)=\frac{1}{2}\left|\right|S_{m}m|{|}_2^2$

$S_d=diag(1/\underset{k}{\&Sigma;}{\left(L_{ik}\right)}^2)$

S_m＝W_mW_z

其中，J^α(m)表示Tikhonov正则化模型，ρ(Lm,d)表示定义在数据域的函数，Ω(m)表示定义在参数域的函数，α表示正则化参数，L表示离散化紧算子，m表示磁化率向量，d表示所述全张量磁梯度数据，S_d表示作用于数据的尺度算子，S_m表示作用于模型的尺度算子，k表示所述离散化紧算子L的列号，i表示所述离散化紧算子L的行号，diag表示对角化，W_m表示正则化模型的先验约束，W_z表示作用于模型深度的先验约束，表示大于0的常数；

所述磁化率反演模块用于：通过以下公式对正则化模型进行迭代求解，反演得到磁化率：

m_k+1＝m_k+τ_kh_k,

τ_k＝arg_τminJ^α(m_k+τh_k),

$h_k=-\&dtri;J^{\mathrm{\&alpha;}}\left(m_k\right)+{\mathrm{\&beta;}}_{k-1}^{FR}{h}_{k-1},$

$\begin{array}{ccc}{h}_k=-\&dtri;J^{\mathrm{\&alpha;}}\left(m_k\right),& if& k=0\end{array},$

${\mathrm{\&beta;}}_{k-1}^{FR}=\left|\right|\&dtri;J^{\mathrm{\&alpha;}}\left(m_k\right)|{|}^2/||\&dtri;J^{\mathrm{\&alpha;}}\left(m_{k-1}\right)|{|}^2$

其中，m表示磁化率向量，k表示迭代次数，τ表示步长，h表示在负梯度方向为初始方向的搜索方向，J表示所述正则化模型，α表示正则化参数，表示所述正则化模型的梯度，β^FR表示共轭方向参数变量。

7.根据权利要求6所述的装置，其特征在于，所述磁化率反演模块用于：

采用预条件混合共轭梯度算法对所述正则化模型进行迭代求解，反演得到磁化率；

其中，预条件混合共轭梯度算法通过以下公式实现：

m_k+1＝m_k+τ_kh_k,

τ_k＝argminJ^α(m_k+τh_k),

$h_k=P_g\left(m_k\right)+{\mathrm{\&beta;}}_{k-1}^{hybrid}{h}_{k-1},$

h_k＝P_g(m₀),if k＝0,

${\mathrm{\&beta;}}_k^{hybrid}=\left\{\begin{array}{ccc}0,& if& {\mathrm{\&beta;}}_k^{HS}<0\\ {\mathrm{\&beta;}}_k^{HS},& if& 0\&le;{\mathrm{\&beta;}}_k^{HS}<{\mathrm{\&beta;}}_k^{DY}\\ {\mathrm{\&beta;}}_k^{DY},& if& {\mathrm{\&beta;}}_k^{HS}>{\mathrm{\&beta;}}_k^{DY}\end{array}\right.$

${\mathrm{\&beta;}}_k^{HS}=\frac{g{\left(m_k\right)}^T\left(g\right(m_k)-g(m_{k-1}\left)\right)}{h_{k-1}^T\left(g\right(m_k)-g(m_{k-1}\left)\right)}$

${\mathrm{\&beta;}}_k^{DY}=\frac{g{\left(m_k\right)}^{T}g\left(m_k\right)}{h_{k-1}^T\left(g\right(m_k)-g(m_{k-1}\left)\right)}$

其中，m表示磁化率向量，k表示迭代次数，τ表示步长，h表示在负梯度方向为初始方向的搜索方向，J表示所述正则化模型，α表示正则化参数，P表示预条件矩阵，g表示所述正则化模型的梯度，β^hybrid表示共轭方向混合参数变量，β^HS表示HS型共轭方向参数变量，β^DY表示DY型共轭方向参数变量。

8.根据权利要求6或7所述的装置，其特征在于，所述装置还包括：

有界域限定模块，用于通过凸集投影技术限定所述正则化模型的有界域，其中，所述凸集投影技术通过以下公式实现：

m_k+1＝P_Π(m_k+τ_kh_k)

$P_{\&Pi;}(m_k+{\mathrm{\&tau;}}_k{h}_k)=\left\{\begin{array}{cc}(m_k+{\mathrm{\&tau;}}_k{h}_k),& (m_k+{\mathrm{\&tau;}}_k{h}_k)\&GreaterEqual;0\\ 0,& (m_k+{\mathrm{\&tau;}}_k{h}_k)<0\end{array}\right.$

其中，m表示磁化率向量，k表示迭代次数，P_Π表示凸集投影算子，τ表示步长，h表示在负梯度方向为初始方向的搜索方向。