CN105631912A - 页岩微米孔隙成像方法及装置 - Google Patents

Abstract

本发明提供了一种页岩微米孔隙成像方法及装置，涉及勘探地球物理技术。其中方法包括：采用同步辐射平行X射线束对页岩进行CT成像，得到页岩的成像数据；根据页岩的成像数据建立稀疏正则化模型；采用梯度下降算法对稀疏正则化模型进行求解，得到重构的页岩内部的微观图像。本发明能够对页岩内部进行微观成像，特别是对页岩内部的微米孔隙进行微观成像，得到精准的页岩内部微观结构图像，缓解相关技术中通过滤波反投影方法观测页岩内部微观结构，观测得到的图像具有噪声和重构伪影的问题。

Description

页岩微米孔隙成像方法及装置

技术领域

本发明涉及勘探地球物理技术，尤其涉及一种页岩微米孔隙成像方法及装置。

背景技术

页岩是沉积岩的一种，主要由黏土沉积经压力和温度所形成。页岩气是蕴藏于页岩层中可供开采的天然气资源，成分以甲烷为主，是一种清洁、高效的能源资源和化工原料。页岩气主要用于居民燃气、城市供热、发电、汽车燃料和化工生产等，用途广泛。中国的页岩气可采储量居世界第一。

尽管页岩气存在于页岩中非自然联通并且难以提取的小断裂中，但在过去的数十年中，页岩气勘探仍然发展迅猛。页岩的微观结构是页岩气勘探的关键问题，在勘探页岩气的过程中，必须对页岩的微观结构进行观测。

传统的基于表面的观测方法，如通过光学显微镜和扫描显微镜进行观测，是常用的能够提供物质微观结构信息的方法。然而，通过这种观测方法只能观测页岩的表面微观结构，无法观测页岩内部的3D细节信息，并且样品通常会受损，对于一些难以开采的页岩样本，并不适用。基于此，相关技术提供了一种X射线计算机断层扫描(CT)技术，并利用滤波反投影算法重构页岩内部微观结构的方法。

发明人在研究中发现，通过相关技术中的滤波反投影算法观测页岩内部微观结构，观测得到的图像具有噪声和重构伪影，具有图像不够精准的缺陷。

发明内容

有鉴于此，本发明提供了一种页岩微米孔隙成像方法及装置，能够对页岩内部进行微观成像，得到精准的页岩内部微观结构图像，缓解相关技术中通过滤波反投影方法观测页岩内部微观结构，观测得到的图像具有噪声和重构伪影的问题。

为达到上述目的，第一方面，本发明实施例提供了一种页岩微米孔隙成像方法，所述方法包括：采用同步辐射平行X射线束对页岩进行CT成像，得到所述页岩的成像数据；根据所述页岩的成像数据建立稀疏正则化模型；采用梯度下降算法对所述稀疏正则化模型进行求解，得到重构的所述页岩内部的微观图像。

结合第一方面，本发明实施例提供了第一方面第一种可能的实施方式，其中，通过以下公式根据所述页岩的成像数据建立稀疏正则化模型：

${\mathrm{minf}}^{\mathrm{\α}}\left(m\right):=\left|\right|Lm-d|{|}_{l_2}^2+\mathrm{\α}|\left|m\right|{|}_{l_1}$

其中，min表示最小化，f^α(m)表示所述稀疏正则化模型，m表示重构的所述页岩内部的微观图像，数学符号:＝表示定义为，L代表离散拉顿变换，d表示所述页岩的成像数据，α表示正则化因子，表示Lm-d的l₂范数，表示m的l₁范数。

结合第一方面，本发明实施例提供了第一方面第二种可能的实施方式，其中，采用梯度下降算法对所述稀疏正则化模型进行求解，得到重构的所述页岩内部的微观图像，包括：通过以下公式确定所述稀疏正则化模型的梯度；

g(m)≈L^T(Lm-d)+αγ(m)

$\mathrm{\γ}\left(m\right)={\left(\frac{m_1}{\sqrt{{\left(m_1\right)}^T{m}_1+\mathrm{\ϵ}}},\frac{m_2}{\sqrt{{\left(m_2\right)}^T{m}_2+\mathrm{\ϵ}}},\mathrm{...},\frac{m_n}{\sqrt{{\left(m_n\right)}^T{m}_n+\mathrm{\ϵ}}}\right)}^T$

其中，g(m)表示所述稀疏正则化模型的梯度，L代表离散拉顿变换，m表示重构的所述页岩内部的微观图像，d表示所述页岩的成像数据，L^T表示离散拉顿变换的转置，α表示正则化因子，ε表示大于零的常量，n表示m的维数，m_n表示m的第n个分量，(m_n)^T表示m_n的转置；

通过以下公式确定所述稀疏正则化模型的海森矩阵；

H(m)≈L^TL+αχ₃(m)

其中，H(m)表示所述稀疏正则化模型的海森矩阵，L代表离散拉顿变换，L^T表示离散拉顿变换的转置，α表示正则化因子，m表示重构的所述页岩内部的微观图像，ε表示大于零的常量，n表示m的维数，m_n表示m的第n个分量，(m_n)^T表示m_n的转置，i表示介于1和n之间的常整数；

根据梯度下降算法中的非单调梯度算法通过以下公式确定迭代步长；

${\mathrm{\μ}}_k^{Ray1}={\mathrm{\β}}_1{\mathrm{\μ}}_k^1+{\mathrm{\β}}_2{\mathrm{\μ}}_k^2$

${\mathrm{\μ}}_k^1=\frac{(g_{k-1},g_{k-1})}{(g_{k-1},H_k{g}_{k-1})},{\mathrm{\μ}}_k^2=\frac{(g_{k-1},H_k{g}_{k-1})}{(g_{k-1},H_k^T{g}_{k-1})}$

其中，表示所述迭代步长，k表示迭代次数，β₁和β₂均表示正实数，H表示所述海森矩阵，g表示所述梯度；

按照所述迭代步长和所述梯度采用梯度下降算法对所述稀疏正则化模型进行迭代求解，得到重构的所述页岩内部的微观图像。

结合第一方面第二种可能的实施方式，本发明实施例提供了第一方面第三种可能的实施方式，其中，在确定迭代步长之后，所述方法还包括：采用循环重用最陡下降步长法确定每次迭代时梯度下降算法的迭代步长，将确定后的所述迭代步长作为梯度下降算法的迭代步长；其中，所述循环重用最陡下降步长法通过以下公式实现：

${\mathrm{\μ}}_{Nl+i}={\mathrm{\μ}}_{Nl+1}^{Ray1}$

l＝0,1,2,...,m

i＝1,2,...,N

其中，μ_Nl+i表示循环重用后的迭代步长，N表示迭代周期，l表示当前迭代周期的序号，i当前迭代周期内的当前迭代的序号，表示循环重用前的迭代步长。

结合第一方面第二种或第三种可能的实施方式，本发明实施例提供了第一方面第四种可能的实施方式，其中，按照所述迭代步长和所述梯度采用梯度下降算法对所述稀疏正则化模型进行迭代求解，得到重构的所述页岩内部的微观图像，包括：

通过迭代公式m_k+1＝m_k-μ_kg(m_k)计算重构的所述页岩内部的微观图像m_k+1；其中，m表示重构的所述页岩内部的微观图像，k表示迭代次数，μ_k表示所述迭代步长，g表示所述梯度；

根据重构的所述页岩内部的微观图像m_k+1计算梯度g(m_k+1)，当满足公式||g(m_k+1)||≤ε时，迭代停止，将当前计算得到的m_k+1作为重构的所述页岩内部的微观图像，否则，继续迭代；其中，ε表示大于零的常量。

结合第一方面第二种或第三种可能的实施方式，本发明实施例提供了第一方面第五种可能的实施方式，其中，按照所述迭代步长和所述梯度采用梯度下降算法对所述稀疏正则化模型进行迭代求解，得到重构的所述页岩内部的微观图像，包括：

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

Π＝{m∈Rⁿ:0≤m≤∞}

P_Π(m)＝χ_Π(m)

${\&lsqb;\left(P_{\mathrm{\&Pi;}}m\right)\&rsqb;}_i=max(m_i,0)=\left\{\begin{array}{cc}{m}_i,& m_i\&GreaterEqual;0\\ 0,& m_i<0\end{array}\right.$

其中，m表示重构的所述页岩内部的微观图像，Rⁿ表示n维实数空间，P_Π表示作用于可行域Π的正交投影算子，Π表示所述稀疏正则化模型的可行域，χ_Π(m)表示可行域Π的特征函数，χ_Π(m)投影到可行域Π外全为0的所有函数的子空间，[(P_Π(m))]_i表示P_Π(m)的第i个元素，m_i表示m的第i个元素；

通过以下迭代公式计算重构的所述页岩内部的微观图像m_k+1；

m_k+1＝P_Π(m_k-μ_kg(m_k))

其中，m表示重构的所述页岩内部的微观图像，k表示迭代次数，P_Π表示作用于可行域Π的正交投影算子，μ_k表示所述迭代步长，g(m_k)表示所述稀疏正则化模型的梯度；

根据重构的所述页岩内部的微观图像m_k+1计算梯度g(m_k+1)，当满足公式||g(m_k+1)||≤ε时，迭代停止，将当前计算得到的m_k+1作为重构的所述页岩内部的微观图像，否则，继续迭代；其中，ε表示大于零的常量。

对应地，第二方面，本发明实施例提供了一种页岩微米孔隙成像装置，所述装置包括：成像模块，用于采用同步辐射平行X射线束对页岩进行CT成像，得到所述页岩的成像数据；模型建立模块，用于根据所述页岩的成像数据建立稀疏正则化模型；求解模块，用于采用梯度下降算法对所述稀疏正则化模型进行求解，得到重构的所述页岩内部的微观图像。

结合第二方面，本发明实施例提供了第二方面第一种可能的实施方式，其中，所述模型建立模块用于：通过以下公式根据所述页岩的成像数据建立稀疏正则化模型：

${\mathrm{minf}}^{\mathrm{\&alpha;}}\left(m\right):=\left|\right|Lm-d|{|}_{l_2}^2+\mathrm{\&alpha;}|\left|m\right|{|}_{l_1}$

其中，min表示最小化，f^α(m)表示所述稀疏正则化模型，m表示重构的所述页岩内部的微观图像，数学符号:＝表示定义为，L代表离散拉顿变换，d表示所述页岩的成像数据，α表示正则化因子，表示Lm-d的l₂范数，表示m的l₁范数。

结合第二方面，本发明实施例提供了第二方面第二种可能的实施方式，其中，所述求解模块包括：梯度确定单元，用于通过以下公式确定所述稀疏正则化模型的梯度；

g(m)≈L^T(Lm-d)+αγ(m)

$\mathrm{\&gamma;}\left(m\right)={\left(\frac{m_1}{\sqrt{{\left(m_1\right)}^T{m}_1+\mathrm{\&epsiv;}}},\frac{m_2}{\sqrt{{\left(m_2\right)}^T{m}_2+\mathrm{\&epsiv;}}},\mathrm{...},\frac{m_n}{\sqrt{{\left(m_n\right)}^T{m}_n+\mathrm{\&epsiv;}}}\right)}^T$

其中，g(m)表示所述稀疏正则化模型的梯度，L代表离散拉顿变换，m表示重构的所述页岩内部的微观图像，d表示所述页岩的成像数据，L^T表示离散拉顿变换的转置，α表示正则化因子，ε表示大于零的常量，n表示m的维数，m_n表示m的第n个分量，(m_n)^T表示m_n的转置；

海森矩阵确定单元，用于通过以下公式确定所述稀疏正则化模型的海森矩阵；

H(m)≈L^TL+αχ₃(m)

其中，H(m)表示所述稀疏正则化模型的海森矩阵，L代表离散拉顿变换，L^T表示离散拉顿变换的转置，α表示正则化因子，m表示重构的所述页岩内部的微观图像，ε表示大于零的常量，n表示m的维数，m_n表示m的第n个分量，(m_n)^T表示m_n的转置，i表示介于1和n之间的常整数；

迭代步长确定单元，用于根据梯度下降算法中的非单调梯度算法通过以下公式确定迭代步长；

${\mathrm{\&mu;}}_k^{Ray1}={\mathrm{\&beta;}}_1{\mathrm{\&mu;}}_k^1+{\mathrm{\&beta;}}_2{\mathrm{\&mu;}}_k^2$

${\mathrm{\&mu;}}_k^1=\frac{(g_{k-1},g_{k-1})}{(g_{k-1},H_k{g}_{k-1})},{\mathrm{\&mu;}}_k^2=\frac{(g_{k-1},H_k{g}_{k-1})}{(g_{k-1},H_k^T{g}_{k-1})}$

其中，表示所述迭代步长，k表示迭代次数，β₁和β₂均表示正实数，H表示所述海森矩阵，g表示所述梯度；

迭代求解单元，用于按照所述迭代步长和所述梯度采用梯度下降算法对所述稀疏正则化模型进行迭代求解，得到重构的所述页岩内部的微观图像。

结合第二方面第二种可能的实施方式，本发明实施例提供了第二方面第三种可能的实施方式，其中，所述求解模块还包括：循环重用单元，用于采用循环重用最陡下降步长法确定每次迭代时梯度下降算法的迭代步长，将确定后的所述迭代步长作为梯度下降算法的迭代步长；其中，所述循环重用最陡下降步长法通过以下公式实现：

${\mathrm{\&mu;}}_{Nl+i}={\mathrm{\&mu;}}_{Nl+1}^{Ray1}$

l＝0,1,2,...,m

i＝1,2,...,N

其中，μ_Nl+i表示循环重用后的迭代步长，N表示迭代周期，l表示当前迭代周期的序号，i当前迭代周期内的当前迭代的序号，表示循环重用前的迭代步长。

结合第二方面第二种或第三种可能的实施方式，本发明实施例提供了第二方面第四种可能的实施方式，其中，所述迭代求解单元包括：

第一重构子单元，用于通过迭代公式m_k+1＝m_k-μ_kg(m_k)计算重构的所述页岩内部的微观图像m_k+1；其中，m表示重构的所述页岩内部的微观图像，k表示迭代次数，μ_k表示所述迭代步长，g(m_k)表示所述梯度；

第一迭代子单元，用于根据重构的所述页岩内部的微观图像m_k+1计算梯度g(m_k+1)，当满足公式||g(m_k+1)||≤ε时，迭代停止，将当前计算得到的m_k+1作为重构的所述页岩内部的微观图像，否则，继续迭代；其中，ε表示大于零的常量。

结合第二方面第二种或第三种可能的实施方式，本发明实施例提供了第二方面第五种可能的实施方式，其中，所述迭代求解单元包括：

凸集投影子单元，用于通过凸集投影技术限定所述稀疏正则化模型的可行域，其中，所述凸集投影技术通过以下公式实现：

Π＝{m∈Rⁿ:0≤m≤∞}

P_Π(m)＝χ_Π(m)

${\&lsqb;\left(P_{\mathrm{\&Pi;}}m\right)\&rsqb;}_i=max(m_i,0)=\left\{\begin{array}{cc}{m}_i,& m_i\&GreaterEqual;0\\ 0,& m_i<0\end{array}\right.$

其中，m表示重构的所述页岩内部的微观图像，Rⁿ表示n维实数空间，P_Π表示作用于可行域Π的正交投影算子，Π表示所述稀疏正则化模型的可行域，χ_Π(m)表示可行域Π的特征函数，χ_Π(m)投影到可行域Π外全为0的所有函数的子空间，[(P_Π(m))]_i表示P_Π(m)的第i个元素，m_i表示m的第i个元素；

第二重构子单元，用于通过以下迭代公式计算重构的所述页岩内部的微观图像m_k+1；

m_k+1＝P_Π(m_k-μ_kg(m_k))

其中，m表示重构的所述页岩内部的微观图像，k表示迭代次数，P_Π表示作用于可行域Π的正交投影算子，μ_k表示所述迭代步长，g(m_k)表示所述稀疏正则化模型的梯度；

第二迭代子单元，用于根据重构的所述页岩内部的微观图像m_k+1计算梯度g(m_k+1)，当满足公式||g(m_k+1)||≤ε时，迭代停止，将当前计算得到的m_k+1作为重构的所述页岩内部的微观图像，否则，继续迭代；其中，ε表示大于零的常量。

本发明实施例中，首先采用同步辐射平行X射线束对页岩进行CT成像，得到页岩的成像数据，然后根据页岩的成像数据建立稀疏正则化模型，最后采用梯度下降算法对稀疏正则化模型进行求解，得到重构的页岩内部的微观图像。通过本实施例中的页岩微米孔隙成像方法及装置，能够对页岩内部进行微观成像，特别是对页岩内部的微米孔隙进行成像，得到精准的页岩内部微观结构图像，缓解相关技术中通过滤波反投影方法观测页岩内部微观结构，观测得到的图像具有噪声和重构伪影的问题。

附图说明

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

图1示出本发明实施例所提供的页岩微米孔隙成像方法的一种流程示意图；

图2示出本发明实施例所提供的多能量X射线计算机断层扫描原理示意图；

图3示出本发明实施例所提供的页岩微米孔隙成像方法的另一种流程示意图；

图4示出本发明实施例所提供的页岩微米孔隙成像装置的一种结构示意图。

具体实施方式

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

考虑到通过相关技术中的滤波反投影算法观测页岩内部微观结构，观测得到的图像具有噪声和重构伪影，具有图像不够精准的缺陷，基于此，本发明提供了一种页岩微米孔隙成像方法及装置，下面结合实施例进行详细描述。

图1示出本发明实施例所提供的页岩微米孔隙成像方法的一种流程示意图。如图1所示，本实施例提供的页岩微米孔隙成像方法包括以下步骤：

步骤S102，采用同步辐射平行X射线束对页岩进行CT(ComputedTomography，电子计算机断层扫描)成像，得到页岩的成像数据；

步骤S104，根据页岩的成像数据建立稀疏正则化模型；

步骤S106，采用梯度下降算法对稀疏正则化模型进行求解，得到重构的页岩内部的微观图像。

本发明实施例中，首先采用同步辐射平行X射线束对页岩进行CT成像，得到页岩的成像数据，然后根据页岩的成像数据建立稀疏正则化模型，最后采用梯度下降算法对稀疏正则化模型进行求解，得到重构的页岩内部的微观图像。通过本实施例中的方法，能够对页岩内部进行微观成像，特别是对页岩内部的微米孔隙进行成像，得到精准的页岩内部微观结构图像，缓解相关技术中通过滤波反投影方法观测页岩内部微观结构，观测得到的图像具有噪声和重构伪影的问题。

上述步骤S102中，将同步辐射X射线断层扫描技术，即CT技术，用于页岩气微观结构成像问题。采用同步辐射源产生能量可调、单色、自然平行的X射线光束。具体地，采用上海光源同步辐射装置BL13W1来获取页岩样品的成像数据。

图2示出本发明实施例所提供的多能量X射线计算机断层扫描原理图。对于计算机断层扫描问题，采用平行X射线束。多能量X射线微型断层扫描过程通过上海光源同步辐射装置(SSRF)BL13W1实现。储存环能量在3.5GeV，电流230mA准恒流模式下运行。BL13W1的组装如图2所示。X射线源产生自一个16极正弦式轨迹激励磁铁，磁感应强度1.9T，磁周期长度为14cm。X射线源是同步加速器，例如，一束电子在储存环中被加速，当它通过一个弯曲的磁系统时释放出高频X射线辐射。通过这种方式产生的X射线具有较宽的频谱。为了使频谱变窄，可以采用Si(111)双晶单色仪来单色化白光。经过单色仪处理以后，此光束通过一个绕垂直轴旋转的样本并使其被闪烁基数器吸收，并发出可见光。对于现代的闪烁基数器来说，发射光线的强度几乎与吸收的X射线的强度成线性正比。

本实施例中，圆柱形页岩样本被安置在一个旋转台上，并以0.167°的增量步长连续旋转180°。X射线能量分别为20keV和30keV。样本到探测器的距离为10cm，曝光时间为4s。传播的X射线被一层薄薄的闪烁体探测器屏幕吸收并将X射线转换成一种特定波长的可见光。一个电荷耦合器件(CCD)探测器通过一个10倍物镜，可以用于获取图像，本地像素尺寸为0.74微米×0.74微米。

本实施例中，将稀疏正则化模型用于解决页岩稀疏恢复问题。上述步骤S104中，通过以下公式(1)根据页岩的成像数据建立稀疏正则化模型：

${\mathrm{minf}}^{\mathrm{\&alpha;}}\left(m\right):=\left|\right|Lm-d|{|}_{l_2}^2+\mathrm{\&alpha;}|\left|m\right|{|}_{l_1}---\left(1\right)$

其中，min表示最小化，f^α(m)表示建立的稀疏正则化模型，m表示重构的页岩内部的微观图像，数学符号:＝表示定义为，L代表离散拉顿变换，d表示页岩的成像数据，α表示正则化因子，表示Lm-d的l₂范数，表示m的l₁范数。

公式(1)中，m为最终要求解得到的量，d为通过CT成像技术对页岩进行观测得到的量。正则化因子α能够由用户按照需要设定。通过上述公式(1)能够根据页岩的成像数据建立稀疏正则化模型，以便通过稀疏正则化模型求解得到重构的页岩内部的微观图像，即最终要得到的页岩内部图像。

为了得到精准的页岩内部的微观图像，上述步骤S106中，采用梯度下降算法对稀疏正则化模型进行求解，得到重构的页岩内部的微观图像，具体包括：

过程11，通过以下公式(2)和(3)确定稀疏正则化模型的梯度；

g(m)≈L^T(Lm-d)+αγ(m)(2)

$\mathrm{\&gamma;}\left(m\right)={\left(\frac{m_1}{\sqrt{{\left(m_1\right)}^T{m}_1+\mathrm{\&epsiv;}}},\frac{m_2}{\sqrt{{\left(m_2\right)}^T{m}_2+\mathrm{\&epsiv;}}},\mathrm{...},\frac{m_n}{\sqrt{{\left(m_n\right)}^T{m}_n+\mathrm{\&epsiv;}}}\right)}^T---\left(3\right)$

其中，g(m)表示确定的稀疏正则化模型的梯度，L代表离散拉顿变换，m表示重构的页岩内部的微观图像，d表示页岩的成像数据，L^T表示离散拉顿变换的转置，α表示正则化因子，ε表示大于零的常量，n表示m的维数，m_n表示m的第n个分量，(m_n)^T表示m_n的转置。

过程11中，m为最终要求解得到的量，d为通过CT成像技术对页岩进行观测得到的量。m的维数n、m_n、(m_n)^T为与m有关的量。离散拉顿变换L、离散拉顿变换的转置L^T、正则化因子α和大于零的常量ε均为已知量，其中，正则化因子α和大于零的常量ε能够由用户设定。

过程22，通过以下公式(4)和(5)确定稀疏正则化模型的海森(Hessian)矩阵；

H(m)≈L^TL+αχ₃(m)(4)

其中，H(m)表示确定的稀疏正则化模型的海森矩阵，L代表离散拉顿变换，L^T表示离散拉顿变换的转置，α表示正则化因子，m表示重构的页岩内部的微观图像，ε表示大于零的常量，n表示m的维数，m_n表示m的第n个分量，(m_n)^T表示m_n的转置，i表示介于1和n之间的常整数。

过程22中，m为最终要求解得到的量。m的维数n、m_n、(m_n)^T为与m有关的量。离散拉顿变换L、离散拉顿变换的转置L^T、正则化因子α和大于零的常量ε均为已知量，其中，正则化因子α和大于零的常量ε能够由用户设定。

过程33，根据梯度下降算法中的非单调梯度算法通过以下公式(6)和(7)确定迭代步长；

${\mathrm{\&mu;}}_k^{Ray1}={\mathrm{\&beta;}}_1{\mathrm{\&mu;}}_k^1+{\mathrm{\&beta;}}_2{\mathrm{\&mu;}}_k^2---\left(6\right)$

${\mathrm{\&mu;}}_k^1=\frac{(g_{k-1},g_{k-1})}{(g_{k-1},H_k{g}_{k-1})},{\mathrm{\&mu;}}_k^2=\frac{(g_{k-1},H_k{g}_{k-1})}{(g_{k-1},H_k^T{g}_{k-1})}---\left(7\right)$

其中，表示最终确定的迭代步长，k表示迭代次数，β₁和β₂均表示正实数，H表示上述海森矩阵，g表示上述梯度。

过程33中，β₁和β₂能够由于用户设定，通常选择β₂＝1-β₁并β₂≤β₁。典型值可以选取β₁＝0.6和β₂＝0.4。

过程44，按照上述迭代步长和上述梯度采用梯度下降算法对稀疏正则化模型进行迭代求解，得到重构的页岩内部的微观图像。

具体地，过程44中，按照上述迭代步长和上述梯度采用梯度下降算法对稀疏正则化模型进行迭代求解，得到重构的页岩内部的微观图像，包括：

(a)通过迭代公式m_k+1＝m_k-μ_kg(m_k)计算重构的页岩内部的微观图像m_k+1；其中，m表示重构的页岩内部的微观图像，k表示迭代次数，μ_k表示迭代步长，g(m_k)表示梯度；

(b)根据重构的页岩内部的微观图像m_k+1计算梯度g(m_k+1)，当满足公式||g(m_k+1)||≤ε时，迭代停止，将当前计算得到的m_k+1作为重构的页岩内部的微观图像，否则，继续迭代；其中，ε表示大于零的常量。

通过上述步骤S106，能够利用梯度下降算法对稀疏正则化模型进行迭代求解，得到精准的页岩内部微观图像，特别是对页岩内部的微米孔隙进行成像，克服相关技术在处理页岩微米孔隙成像方面分辨率、噪声压制以及计算速率的不足。

结合上述步骤S102至步骤S106，对稀疏正则化模型进行迭代求解包括以下过程：

1)获取页岩的成像数据；

2)根据页岩的成像数据通过公式(1)建立稀疏正则化模型；

3)设定m的初始值m₀，ε＞0，令k:＝1，通过公式(2)和(3)计算梯度g(m_k)；

4)判断||g(m_k)||≤ε是否成立，如果成立，执行5)，否则执行6)；

5)迭代停止，确定此时的m为最终得到的结果；

6)通过公式(4)和(5)计算海森矩阵；

7)通过公式(6)和(7)利用梯度和海森矩阵计算迭代步长；

8)通过公式m_k+1＝m_k-μ_kg(m_k)计算m_k+1；

9)计算g(m_k+1)，令k:＝k+1，并跳转至4)。

循环执行上述1)至9)，直至求解得到m。

考虑到在梯度下降算法中选择合适的迭代步长，能够减少迭代次数，提到求解速度，本实施例中，在上述步骤S106过程33之后，还包括以下过程34：

过程34，采用循环重用最陡下降步长法确定每次迭代时梯度下降算法的迭代步长，将确定后的迭代步长作为梯度下降算法的迭代步长；

其中，循环重用最陡下降步长法通过以下公式实现：

${\mathrm{\&mu;}}_{Nl+i}={\mathrm{\&mu;}}_{Nl+1}^{Ray1}---\left(8\right)$

l＝0,1,2,...,m

i＝1,2,...,N

其中，μ_Nl+i表示循环重用后的迭代步长，N表示迭代周期，l表示当前迭代周期的序号，i当前迭代周期内的当前迭代的序号，表示循环重用前的迭代步长。

具体地，N表示迭代周期，如在求解过程中，每迭代5次为一个迭代周期，则N为5。l表示当前迭代周期的序号，如当前处于第一轮迭代周期内，则l取0，当前处于第二轮迭代周期内，则l取1。i表示当前迭代周期内的当前迭代的序号，如i为3，则表示当前迭代周期内的第3次迭代。

当N为5、l取1、i为3时，公式(8)变为也就是说，当前第二轮迭代周期内，第三次迭代的步长与第一次迭代的步长相等。观察公式(8)可知，等式左右两侧均包含迭代周期参量Nl，不同的是，等式左侧表示每个迭代周期内的第i个迭代步长，等式右侧表示每个迭代周期内的第一个迭代步长，也就是说，通过公式(8)，当前迭代周期内每个步长均与第一个步长相等。其中，迭代周期参量Nl表示进入当前迭代周期之前已完成的迭代次数。通过如公式(8)示出的循环重用最陡下降步长法能够减少迭代次数，提到求解速度。

考虑到最终求解结果m应当是有界的，本实施例中，上述步骤S106过程44还能够通过以下过程实现：

(a)通过凸集投影技术限定稀疏正则化模型的可行域，其中，凸集投影技术通过以下公式实现：

Π＝{m∈Rⁿ:0≤m≤∞}(9)

P_Π(m)＝χ_Π(m)(10)

${\&lsqb;\left(P_{\mathrm{\&Pi;}}m\right)\&rsqb;}_i=max(m_i,0)=\left\{\begin{array}{cc}{m}_i,& m_i\&GreaterEqual;0\\ 0,& m_i<0\end{array}\right.---\left(11\right)$

其中，m表示重构的页岩内部的微观图像，Rⁿ表示n维实数空间，P_Π表示作用于可行域Π的正交投影算子，Π表示稀疏正则化模型的可行域，χ_Π(m)表示可行域Π的特征函数，χ_Π(m)投影到可行域Π外全为0的所有函数的子空间，[(P_Π(m))]_i表示P_Π(m)的第i个元素，m_i表示m的第i个元素；

(b)通过以下迭代公式计算重构的页岩内部的微观图像m_k+1；

m_k+1＝P_Π(m_k-μ_kg(m_k))(12)

其中，m表示重构的页岩内部的微观图像，k表示迭代次数，P_Π表示作用于可行域Π的正交投影算子，μ_k表示迭代步长，g(m_k)表示稀疏正则化模型的梯度；

(c)根据重构的页岩内部的微观图像m_k+1计算梯度g(m_k+1)，当满足公式||g(m_k+1)||≤ε时，迭代停止，将当前计算得到的m_k+1作为重构的页岩内部的微观图像，否则，继续迭代；其中，ε表示大于零的常量。

通过上述过程(a)(b)(c)能够限定稀疏正则化模型的可行域，从而加快求解速度，得到符合物理意义的重构的页岩内部的微观图像m。

上述过程34和过程(a)(b)(c)能够结合使用，从而充分加速求解速度，提到求解效率。当上述过程34和(a)(b)(c)结合使用时，如图3所示，对稀疏正则化模型进行迭代求解包括以下过程：

步骤S301，获取页岩的成像数据；

步骤S302，根据页岩的成像数据通过公式(1)建立稀疏正则化模型；

步骤S303，设定m的初始值m₀，ε＞0，令k:＝1，通过公式(2)和(3)计算梯度g(m_k)；

步骤S304，判断||g(m_k)||≤ε是否成立，如果成立，执行步骤S305，否则执行步骤S306；

步骤S305，迭代停止，确定此时的m为最终得到的结果；

步骤S306，通过公式(4)和(5)计算海森矩阵；

步骤S307，通过公式(6)和(7)利用梯度和海森矩阵计算迭代步长；

步骤S308，通过公式(8)采用循环重用最陡下降步长法，根据步骤S307中计算出的步长重新确定迭代步长：

步骤S309，通过公式(9)至(11)通过凸集投影技术限定稀疏正则化模型的可行域；

步骤S310，按照公式m_k+1＝P_Π(m_k-μ_kg(m_k))更新重构的页岩内部的微观图像m，得到m_k+1，并计算梯度g(m_k+1)，令k:＝k+1跳转至步骤S304。

需要说明的是，上述与μ_k均表示迭代步长，因不同公式所以采用不同写法，二者指代同一物理参量。

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

根据本发明实施例提供的技术方案，本发明实施例中的方法可以稳定的重建页岩微观结构，特别是对页岩内部的微米孔隙进行成像，提供一系列高质量的CT成像数据以便用于进一步的处理，例如，矿物组分分析、图像分割等等。本实施例与商业软件相比，重构结构可以消除可见的环状噪声(这些噪声很可能为矿物组分分析带来错误的解释，所以必须消除)。本实施例能够得到更高的成像分辨率，能够更大程度的消除噪声，能够更准确的区分高、低密度的矿物成分，提高对矿物组分百分比估计的精度。由此可以看出，本发明实施例提出的方法可以很好的替代商用软件常用的滤波反投影方法用以CT图像的重构。

进一步地，本实施例还提供了一种页岩微米孔隙成像装置，用于执行上述方法，如图4所示，该装置包括：

成像模块41，用于采用同步辐射平行X射线束对页岩进行CT成像，得到页岩的成像数据；

模型建立模块42，用于根据页岩的成像数据建立稀疏正则化模型；

求解模块43，用于采用梯度下降算法对稀疏正则化模型进行求解，得到重构的页岩内部的微观图像。

本发明实施例中，首先采用同步辐射平行X射线束对页岩进行CT成像，得到页岩的成像数据，然后根据页岩的成像数据建立稀疏正则化模型，最后采用梯度下降算法对稀疏正则化模型进行求解，得到重构的页岩内部的微观图像。通过本实施例中的装置，能够对页岩内部进行微观成像，特别是对页岩内部的微米孔隙进行成像，得到精准的页岩内部微观结构图像，缓解相关技术中通过滤波反投影方法观测页岩内部微观结构，观测得到的图像具有噪声和重构伪影的问题。

本实施例中，将稀疏正则化模型用于解决页岩稀疏恢复问题。上述模型建立模块42用于：通过以下公式根据页岩的成像数据建立稀疏正则化模型：

${\mathrm{minf}}^{\mathrm{\&alpha;}}\left(m\right):=\left|\right|Lm-d|{|}_{l_2}^2+\mathrm{\&alpha;}|\left|m\right|{|}_{l_1}$

其中，min表示最小化，f^α(m)表示稀疏正则化模型，m表示重构的页岩内部的微观图像，数学符号:＝表示定义为，L代表离散拉顿变换，d表示页岩的成像数据，α表示正则化因子，表示Lm-d的l₂范数，表示m的l₁范数。

通过上述模型建立模块42能够根据页岩的成像数据建立稀疏正则化模型，以便通过稀疏正则化模型求解得到重构的页岩内部的微观图像，重构的页岩内部的微观图像即为最终要得到的页岩内部图像。

为了得到精准的页岩内部的微观图像，上述求解模块43包括：梯度确定单元，用于通过以下公式确定稀疏正则化模型的梯度；

g(m)≈L^T(Lm-d)+αγ(m)

$\mathrm{\&gamma;}\left(m\right)={\left(\frac{m_1}{\sqrt{{\left(m_1\right)}^T{m}_1+\mathrm{\&epsiv;}}},\frac{m_2}{\sqrt{{\left(m_2\right)}^T{m}_2+\mathrm{\&epsiv;}}},\mathrm{...},\frac{m_n}{\sqrt{{\left(m_n\right)}^T{m}_n+\mathrm{\&epsiv;}}}\right)}^T$

其中，g(m)表示稀疏正则化模型的梯度，L代表离散拉顿变换，m表示重构的页岩内部的微观图像，d表示页岩的成像数据，L^T表示离散拉顿变换的转置，α表示正则化因子，ε表示大于零的常量，n表示m的维数，m_n表示m的第n个分量，(m_n)^T表示m_n的转置；

海森矩阵确定单元，用于通过以下公式确定稀疏正则化模型的海森矩阵；

H(m)≈L^TL+αχ₃(m)

其中，H(m)表示稀疏正则化模型的海森矩阵，L代表离散拉顿变换，L^T表示离散拉顿变换的转置，α表示正则化因子，m表示重构的页岩内部的微观图像，ε表示大于零的常量，n表示m的维数，m_n表示m的第n个分量，(m_n)^T表示m_n的转置，i表示介于1和n之间的常整数；

迭代步长确定单元，用于根据梯度下降算法中的非单调梯度算法通过以下公式确定迭代步长；

${\mathrm{\&mu;}}_k^{Ray1}={\mathrm{\&beta;}}_1{\mathrm{\&mu;}}_k^1+{\mathrm{\&beta;}}_2{\mathrm{\&mu;}}_k^2$

${\mathrm{\&mu;}}_k^1=\frac{(g_{k-1},g_{k-1})}{(g_{k-1},H_k{g}_{k-1})},{\mathrm{\&mu;}}_k^2=\frac{(g_{k-1},H_k{g}_{k-1})}{(g_{k-1},H_k^T{g}_{k-1})}$

其中，表示迭代步长，k表示迭代次数，β₁和β₂均表示正实数，H表示海森矩阵，g表示梯度；

迭代求解单元，用于按照迭代步长和梯度采用梯度下降算法对稀疏正则化模型进行迭代求解，得到重构的页岩内部的微观图像。

该迭代求解单元能够通过以下子单元实现：

第一重构子单元，用于通过迭代公式m_k+1＝m_k-μ_kg(m_k)计算重构的页岩内部的微观图像m_k+1；其中，m表示重构的页岩内部的微观图像，k表示迭代次数，μ_k表示迭代步长，g(m_k)表示梯度；

第一迭代子单元，用于根据重构的页岩内部的微观图像m_k+1计算梯度g(m_k+1)，当满足公式||g(m_k+1)||≤ε时，迭代停止，将当前计算得到的m_k+1作为重构的页岩内部的微观图像，否则，继续迭代；其中，ε表示大于零的常量。

通过上述求解模块43，能够利用梯度下降算法对稀疏正则化模型进行迭代求解，得到精准的页岩内部微观图像，特别是对页岩内部的微米孔隙进行成像，克服相关技术在处理页岩微米孔隙成像方面分辨率、噪声压制以及计算速率的不足。

考虑到在梯度下降算法中选择合适的迭代步长，能够减少迭代次数，提到求解速度，本实施例中，求解模块43还包括：循环重用单元，用于采用循环重用最陡下降步长法确定每次迭代时梯度下降算法的迭代步长，将确定后的迭代步长作为梯度下降算法的迭代步长；其中，循环重用最陡下降步长法通过以下公式实现：

${\mathrm{\&mu;}}_{Nl+i}={\mathrm{\&mu;}}_{Nl+1}^{Ray1}$

l＝0,1,2,...,m

i＝1,2,...,N

其中，μ_Nl+i表示循环重用后的迭代步长，N表示迭代周期，l表示当前迭代周期的序号，i当前迭代周期内的当前迭代的序号，表示循环重用前的迭代步长。通过循环重用单元能够减少迭代次数，提到求解速度。

考虑到最终求解结果m应当是有界的，上述迭代求解单元包括：

凸集投影子单元，用于通过凸集投影技术限定稀疏正则化模型的可行域，其中，凸集投影技术通过以下公式实现：

Π＝{m∈Rⁿ:0≤m≤∞}

P_Π(m)＝χ_Π(m)

${\&lsqb;\left(P_{\mathrm{\&Pi;}}m\right)\&rsqb;}_i=max(m_i,0)=\left\{\begin{array}{cc}{m}_i,& m_i\&GreaterEqual;0\\ 0,& m_i<0\end{array}\right.$

其中，m表示重构的页岩内部的微观图像，Rⁿ表示n维实数空间，P_Π表示作用于可行域Π的正交投影算子，Π表稀疏正则化模型的可行域，χ_Π(m)表示可行域Π的特征函数，χ_Π(m)投影到可行域Π外全为0的所有函数的子空间，[(P_Π(m))]_i表示P_Π(m)的第i个元素，m_i表示m的第i个元素；

第二重构子单元，用于通过以下迭代公式计算重构的页岩内部的微观图像m_k+1；

m_k+1＝P_Π(m_k-μ_kg(m_k))

其中，m表示重构的页岩内部的微观图像，k表示迭代次数，P_Π表示作用于可行域Π的正交投影算子，μ_k表示迭代步长，g(m_k)表示稀疏正则化模型的梯度；

第二迭代子单元，用于根据重构的页岩内部的微观图像m_k+1计算梯度g(m_k+1)，当满足公式||g(m_k+1)||≤ε时，迭代停止，将当前计算得到的m_k+1作为重构的页岩内部的微观图像，否则，继续迭代；其中，ε表示大于零的常量。

通过上述迭代求解单元能够限定稀疏正则化模型的可行域，从而加快求解速度，得到符合物理意义的重构的页岩内部的微观图像。

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

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

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

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

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

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

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

Claims (10)

1.页岩微米孔隙成像方法，其特征在于，所述方法包括：

采用同步辐射平行X射线束对页岩进行CT成像，得到所述页岩的成像数据；

根据所述页岩的成像数据建立稀疏正则化模型；

采用梯度下降算法对所述稀疏正则化模型进行求解，得到重构的所述页岩内部的微观图像。

2.根据权利要求1所述的方法，其特征在于，通过以下公式根据所述页岩的成像数据建立稀疏正则化模型：

$\min f^{\mathrm{\&alpha;}}\left(m\right):={\left|\right|\mathrm{Lm}-d\left|\right|}_{l_2}^2+\mathrm{\&alpha;}{\left|\right|m\left|\right|}_{l_1}$

其中，min表示最小化，f^α(m)表示所述稀疏正则化模型，m表示重构的所述页岩内部的微观图像，数学符号:＝表示定义为，L代表离散拉顿变换，d表示所述页岩的成像数据，α表示正则化因子，表示Lm-d的l₂范数，表示m的l₁范数。

3.根据权利要求1所述的方法，其特征在于，采用梯度下降算法对所述稀疏正则化模型进行求解，得到重构的所述页岩内部的微观图像，包括：

通过以下公式确定所述稀疏正则化模型的梯度；

g(m)≈L^T(Lm-d)+αγ(m)

$\mathrm{\&gamma;}\left(m\right)={\left(\frac{m_1}{\sqrt{{\left(m_1\right)}^T{m}_1+\mathrm{\&epsiv;}}},\frac{m_2}{\sqrt{{\left(m_2\right)}^T{m}_2+\mathrm{\&epsiv;}}},...,\frac{m_n}{\sqrt{{\left(m_n\right)}^T{m}_n+\mathrm{\&epsiv;}}}\right)}^T$

其中，g(m)表示所述稀疏正则化模型的梯度，L代表离散拉顿变换，m表示重构的所述页岩内部的微观图像，d表示所述页岩的成像数据，L^T表示离散拉顿变换的转置，α表示正则化因子，ε表示大于零的常量，n表示m的维数，m_n表示m的第n个分量，(m_n)^T表示m_n的转置；

通过以下公式确定所述稀疏正则化模型的海森矩阵；

H(m)≈L^TL+αχ₃(m)

其中，H(m)表示所述稀疏正则化模型的海森矩阵，L代表离散拉顿变换，L^T表示离散拉顿变换的转置，α表示正则化因子，m表示重构的所述页岩内部的微观图像，ε表示大于零的常量，n表示m的维数，m_n表示m的第n个分量，(m_n)^T表示m_n的转置，i表示介于1和n之间的常整数；

根据梯度下降算法中的非单调梯度算法通过以下公式确定迭代步长；

${\mathrm{\&mu;}}_k^{Ray1}={\mathrm{\&beta;}}_1{\mathrm{\&mu;}}_k^1+{\mathrm{\&beta;}}_2{\mathrm{\&mu;}}_k^2$

${\mathrm{\&mu;}}_k^1=\frac{(g_{k-1},g_{k-1})}{(g_{k-1},H_k{g}_{k-1})},{\mathrm{\&mu;}}_k^2=\frac{(g_{k-1},H_k{g}_{k-1})}{(g_{k-1},H_k^T{g}_{k-1})}$

其中，表示所述迭代步长，k表示迭代次数，β₁和β₂均表示正实数，H表示所述海森矩阵，g表示所述梯度；

按照所述迭代步长和所述梯度采用梯度下降算法对所述稀疏正则化模型进行迭代求解，得到重构的所述页岩内部的微观图像。

4.根据权利要求3所述的方法，其特征在于，在确定迭代步长之后，所述方法还包括：

采用循环重用最陡下降步长法确定每次迭代时梯度下降算法的迭代步长，将确定后的所述迭代步长作为梯度下降算法的迭代步长；

其中，所述循环重用最陡下降步长法通过以下公式实现：

${\mathrm{\&mu;}}_{Nl+i}={\mathrm{\&mu;}}_{Nl+1}^{Ray1}$

l＝0,1,2,...,m

i＝1,2,...,N

其中，μ_Nl+i表示循环重用后的迭代步长，N表示迭代周期，l表示当前迭代周期的序号，i当前迭代周期内的当前迭代的序号，表示循环重用前的迭代步长。

5.根据权利要求3或4所述的方法，其特征在于，按照所述迭代步长和所述梯度采用梯度下降算法对所述稀疏正则化模型进行迭代求解，得到重构的所述页岩内部的微观图像，包括：

通过迭代公式m_k+1＝m_k-μ_kg(m_k)计算重构的所述页岩内部的微观图像m_k+1；其中，m表示重构的所述页岩内部的微观图像，k表示迭代次数，μ_k表示所述迭代步长，g(m_k)表示所述梯度；

根据重构的所述页岩内部的微观图像m_k+1计算梯度g(m_k+1)，当满足公式||g(m_k+1)||≤ε时，迭代停止，将当前计算得到的m_k+1作为重构的所述页岩内部的微观图像，否则，继续迭代；其中，ε表示大于零的常量。

6.根据权利要求3或4所述的方法，其特征在于，按照所述迭代步长和所述梯度采用梯度下降算法对所述稀疏正则化模型进行迭代求解，得到重构的所述页岩内部的微观图像，包括：

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

Π＝{m∈Rⁿ:0≤m≤∞}

P_Π(m)＝χ_Π(m)

${\&lsqb;\left(P_{\mathrm{\&Pi;}}\right(m\left)\right)\&rsqb;}_i=max(m_i,0)=\left\{\begin{array}{cc}{m}_i,& m_i\&GreaterEqual;0\\ 0,& m_i<0\end{array}\right.$

其中，m表示重构的所述页岩内部的微观图像，Rⁿ表示n维实数空间，P_Π表示作用于可行域Π的正交投影算子，Π表示所述稀疏正则化模型的可行域，χ_Π(m)表示可行域Π的特征函数，χ_Π(m)投影到可行域Π外全为0的所有函数的子空间，[(P_Π(m))]_i表示P_Π(m)的第i个元素，m_i表示m的第i个元素；

通过以下迭代公式计算重构的所述页岩内部的微观图像m_k+1；

m_k+1＝P_Π(m_k-μ_kg(m_k))

其中，m表示重构的所述页岩内部的微观图像，k表示迭代次数，P_Π表示作用于可行域Π的正交投影算子，μ_k表示所述迭代步长，g(m_k)表示所述稀疏正则化模型的梯度；

根据重构的所述页岩内部的微观图像m_k+1计算梯度g(m_k+1)，当满足公式||g(m_k+1)||≤ε时，迭代停止，将当前计算得到的m_k+1作为重构的所述页岩内部的微观图像，否则，继续迭代；其中，ε表示大于零的常量。

7.页岩微米孔隙成像装置，其特征在于，所述装置包括：

成像模块，用于采用同步辐射平行X射线束对页岩进行CT成像，得到所述页岩的成像数据；

模型建立模块，用于根据所述页岩的成像数据建立稀疏正则化模型；

求解模块，用于采用梯度下降算法对所述稀疏正则化模型进行求解，得到重构的所述页岩内部的微观图像。

8.根据权利要求7所述的装置，其特征在于，所述求解模块包括：

梯度确定单元，用于通过以下公式确定所述稀疏正则化模型的梯度；

g(m)≈L^T(Lm-d)+αγ(m)

$\mathrm{\&gamma;}\left(m\right)={\left(\frac{m_1}{\sqrt{{\left(m_1\right)}^T{m}_1+\mathrm{\&epsiv;}}},\frac{m_2}{\sqrt{{\left(m_2\right)}^T{m}_2+\mathrm{\&epsiv;}}},...,\frac{m_n}{\sqrt{{\left(m_n\right)}^T{m}_n+\mathrm{\&epsiv;}}}\right)}^T$

其中，g(m)表示所述稀疏正则化模型的梯度，L代表离散拉顿变换，m表示重构的所述页岩内部的微观图像，d表示所述页岩的成像数据，L^T表示离散拉顿变换的转置，α表示正则化因子，ε表示大于零的常量，n表示m的维数，m_n表示m的第n个分量，(m_n)^T表示m_n的转置；

海森矩阵确定单元，用于通过以下公式确定所述稀疏正则化模型的海森矩阵；

H(m)≈L^TL+αχ₃(m)

其中，H(m)表示所述稀疏正则化模型的海森矩阵，L代表离散拉顿变换，L^T表示离散拉顿变换的转置，α表示正则化因子，m表示重构的所述页岩内部的微观图像，ε表示大于零的常量，n表示m的维数，m_n表示m的第n个分量，(m_n)^T表示m_n的转置，i表示介于1和n之间的常整数；

迭代步长确定单元，用于根据梯度下降算法中的非单调梯度算法通过以下公式确定迭代步长；

${\mathrm{\&mu;}}_k^{Ray1}={\mathrm{\&beta;}}_1{\mathrm{\&mu;}}_k^1+{\mathrm{\&beta;}}_2{\mathrm{\&mu;}}_k^2$

${\mathrm{\&mu;}}_k^1=\frac{(g_{k-1},g_{k-1})}{(g_{k-1},H_k{g}_{k-1})},{\mathrm{\&mu;}}_k^2=\frac{(g_{k-1},H_k{g}_{k-1})}{(g_{k-1},H_k^T{g}_{k-1})}$

其中，表示所述迭代步长，k表示迭代次数，β₁和β₂均表示正实数，H表示所述海森矩阵，g表示所述梯度；

迭代求解单元，用于按照所述迭代步长和所述梯度采用梯度下降算法对所述稀疏正则化模型进行迭代求解，得到重构的所述页岩内部的微观图像。

9.根据权利要求8所述的装置，其特征在于，所述求解模块还包括：

循环重用单元，用于采用循环重用最陡下降步长法确定每次迭代时梯度下降算法的迭代步长，将确定后的所述迭代步长作为梯度下降算法的迭代步长；

其中，所述循环重用最陡下降步长法通过以下公式实现：

${\mathrm{\&mu;}}_{Nl+i}={\mathrm{\&mu;}}_{Nl+1}^{Ray1}$

l＝0,1,2,...,m

i＝1,2,...,N

其中，μ_Nl+i表示循环重用后的迭代步长，N表示迭代周期，l表示当前迭代周期的序号，i当前迭代周期内的当前迭代的序号，表示循环重用前的迭代步长。

10.根据权利要求8或9所述的装置，其特征在于，所述迭代求解单元包括：

凸集投影子单元，用于通过凸集投影技术限定所述稀疏正则化模型的可行域，其中，所述凸集投影技术通过以下公式实现：

Π＝{m∈Rⁿ:0≤m≤∞}

P_Π(m)＝χ_Π(m)

${\&lsqb;\left(P_{\mathrm{\&Pi;}}m\right)\&rsqb;}_i=max(m_i,0)=\left\{\begin{array}{cc}{m}_i,& m_i\&GreaterEqual;0\\ 0,& m_i<0\end{array}\right.$

其中，m表示重构的所述页岩内部的微观图像，Rⁿ表示n维实数空间，P_Π表示作用于可行域Π的正交投影算子，Π表示所述稀疏正则化模型的可行域，χ_Π(m)表示可行域Π的特征函数，χ_Π(m)投影到可行域Π外全为0的所有函数的子空间，[(P_Π(m))]_i表示P_Π(m)的第i个元素，m_i表示m的第i个元素；

第二重构子单元，用于通过以下迭代公式计算重构的所述页岩内部的微观图像m_k+1；

m_k+1＝P_Π(m_k-μ_kg(m_k))

其中，m表示重构的所述页岩内部的微观图像，k表示迭代次数，P_Π表示作用于可行域Π的正交投影算子，μ_k表示所述迭代步长，g(m_k)表示所述稀疏正则化模型的梯度；

第二迭代子单元，用于根据重构的所述页岩内部的微观图像m_k+1计算梯度g(m_k+1)，当满足公式||g(m_k+1)||≤ε时，迭代停止，将当前计算得到的m_k+1作为重构的所述页岩内部的微观图像，否则，继续迭代；其中，ε表示大于零的常量。