CN106875334B - 页岩ct成像方法及装置 - Google Patents

页岩ct成像方法及装置 Download PDF

Info

Publication number
CN106875334B
CN106875334B CN201710165715.6A CN201710165715A CN106875334B CN 106875334 B CN106875334 B CN 106875334B CN 201710165715 A CN201710165715 A CN 201710165715A CN 106875334 B CN106875334 B CN 106875334B
Authority
CN
China
Prior art keywords
projection
random
vector
data
square
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201710165715.6A
Other languages
English (en)
Other versions
CN106875334A (zh
Inventor
王彦飞
唐巍
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Institute of Geology and Geophysics of CAS
Original Assignee
Institute of Geology and Geophysics of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Institute of Geology and Geophysics of CAS filed Critical Institute of Geology and Geophysics of CAS
Priority to CN201710165715.6A priority Critical patent/CN106875334B/zh
Publication of CN106875334A publication Critical patent/CN106875334A/zh
Application granted granted Critical
Publication of CN106875334B publication Critical patent/CN106875334B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00—Geometric image transformations in the plane of the image
    • G06T3/08—Projecting images onto non-planar surfaces, e.g. geodetic screens
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/70—Denoising; Smoothing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10081—Computed x-ray tomography [CT]

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

本发明提供一种页岩CT成像方法及装置，涉及微纳米CT成像技术领域。通过生成页岩的投影矩阵；依据所述投影矩阵生成随机列向量分布以及随机行向量分布；根据所述随机列向量分布从所述投影矩阵中获得一随机列向量，基于所述随机列向量对投影数据进行处理，得到降噪后的投影数据；根据所述随机行向量分布从所述投影矩阵中获得一行向量，基于所述随机行向量对所述降噪后的投影数据处理，最终得到CT图像。这样处理方法能有效的降低噪音的干扰，且使处理速度加快，准确性提高。对页岩内存储的页岩油气的研究具有重大意义。

Description

页岩CT成像方法及装置
技术领域
本发明涉及微纳米CT成像技术领域，具体而言，涉及一种页岩CT成像方法及装置。
背景技术
近年来，微纳米CT技术广泛的应用于地质、地球化学、地球物理等领域，并对这些领域的研究带来了积极的意义。以页岩为例，页岩内存储页岩油气，在能源稀缺的年代，对页岩内存储的页岩油气研究显得极为重要。但由于页岩油气在页岩内存储方式特殊，只有对于其内部的纳米级孔喉分布状况进行研究，才能有效获得页岩内页岩油气的分布状况。因此，对页岩进行微纳米CT成像对页岩内存储的页岩油气研究具有重大意义。
传统对页岩CT成像的方法是基于傅里叶中心切片定理进行滤波反投影。这是一种显式的，计算速度较快的图像重构方法。但这种方法会产生伪影，使用效果不佳。近年来，随着计算能力的逐年提高，人们用基于迭代的图像重构方法替代传统方法。但现有的基于迭代重构方法收敛速度较慢，且在有噪音的时候不一定能有效收敛到较好的结果。
发明内容
本发明的目的在于提供一种页岩CT成像方法，用以改善上述问题。
本发明的另一目的在于提供一种页岩CT成像装置，用以改善上述问题。
为了实现上述目的，本发明实施例采用的技术方案如下：
本发明实施例提供一种页岩CT成像方法，所述方法包括：生成页岩的投影矩阵；依据所述投影矩阵生成随机列向量分布以及随机行向量分布；根据所述随机列向量分布从所述投影矩阵中获得一随机列向量，基于所述随机列向量对投影数据进行处理，得到降噪后的投影数据；根据所述随机行向量分布从所述投影矩阵中获得一行向量，基于所述随机行向量对所述降噪后的投影数据处理，得到CT图像。
本发明实施例还提供一种页岩CT成像装置，所述装置包括：生成模块，用于生成页岩的投影矩阵；第一计算模块，用于依据所述投影矩阵生成随机列向量分布以及随机行向量分布；降噪模块，用于根据所述随机列向量分布从所述投影矩阵中获得一随机列向量，基于所述随机列向量对投影数据进行处理，得到降噪后的投影数据；迭代模块，用于根据所述随机行向量分布从所述投影矩阵中获得一行向量，基于所述随机行向量对所述降噪后的投影数据处理，得到CT图像。
与现有技术相比，本发明提供一种页岩CT成像方法及装置。通过根据所述随机列向量分布从所述投影矩阵中获得一随机列向量，基于所述随机列向量对投影数据进行处理，得到降噪后的投影数据；根据所述随机行向量分布从所述投影矩阵中获得一行向量，基于所述随机行向量对所述降噪后的投影数据处理，最终得到CT图像。这样处理方法能有效的降低噪音的干扰，且使处理速度加快，准确性提高。对页岩内存储的页岩油气的研究具有重大意义。
为使本发明的上述目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附附图，作详细说明如下。
附图说明
为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。
图1为本发明较佳实施例提供的客户端的方框示意图。
图2为本发明较佳实施例提供的页岩CT成像装置的功能模块示意图。
图3为图2中示出的第一计算模块的功能子模块示意图。
图4为图2中示出的降噪模块的功能子模块示意图。
图5为图2中示出的迭代模块的功能子模块示意图。
图6为本发明较佳实施例提供的页岩CT成像方法的流程图。
图7为图6中示出的生成随机列向量分布以及随机行向量分布的子步骤流程图。
图8为图6中示出的得到降噪后投影数据的子步骤流程图。
图9为图6中示出的得到CT图像的子步骤流程图。
图标：100-客户端；101-存储器；102-存储控制器；103-处理器；104-外设接口；105-显示单元；106-输入输出单元；200-页岩CT成像装置；201-获得模块；202-第二计算模块；203-离散化处理模块；204-生成模块；205-第一计算模块；2051-第一获取子模块；2052-运算子模块；206-降噪模块；2061-获取随机列向量子模块；2062-更新子模块；2063-第一比较子模块；2064-降噪子模块；207-迭代模块；2071-获取随机行向量子模块；2072-迭代子模块；2073-第二比较子模块；2074-第二获取子模块。
具体实施方式
下面将结合本发明实施例中附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。同时，在本发明的描述中，术语“第一”、“第二”等仅用于区分描述，而不能理解为指示或暗示相对重要性。
图1示出本发明较佳实施例提供的客户端100的方框示意图。所述客户端100可以是平板电脑、智能手机、个人数字助理(personal digital assistant，PDA)等。所述客户端100包括页岩CT成像装置200、存储器101、存储控制器102、处理器103、外设接口104、显示单元105、输入输出单元106。
所述存储器101、存储控制器102、处理器103、外设接口104、显示单元105、输入输出单元106各元件相互之间直接或间接地电性连接，以实现数据的传输或交互。例如，这些元件相互之间可通过一条或多条通讯总线或信号线实现电性连接。所述页岩CT成像装置200包括至少一个可以软件或固件(firmware)的形式存储于所述存储器101中或固化在所述客户端100的操作系统(operating system，OS)中的软件功能模块。所述处理器103用于执行存储器101中存储的可执行模块，例如所述页岩CT成像装置200包括的软件功能模块或计算机程序。
其中，存储器101可以是，但不限于，随机存取存储器(Random Access Memory，RAM)，只读存储器(Read Only Memory，ROM)，可编程只读存储器(Programmable Read-OnlyMemory，PROM)，可擦除只读存储器(Erasable Programmable Read-Only Memory，EPROM)，电可擦除只读存储器(Electric Erasable Programmable Read-Only Memory，EEPROM)等。其中，存储器101用于存储程序，所述处理器103在接收到执行指令后，执行所述程序，本发明任一实施例揭示的由过程定义的服务器所执行的方法可以应用于处理器103中，或者由处理器103实现。
处理器103可能是一种集成电路芯片，具有信号的处理能力。上述的处理器103可以是通用处理器，包括中央处理器(Central Processing Unit，简称CPU)、网络处理器(Network Processor，简称NP)等；还可以是数字信号处理器(DSP)、专用集成电路(ASIC)、现场可编程门阵列(FPGA)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件。可以实现或者执行本发明实施例中的公开的各方法、步骤及逻辑框图。通用处理器可以是微处理器或者该处理器103也可以是任何常规的处理器103等。
所述外设接口104将各种输入/输出装置耦合至处理器103以及存储器101。在一些实施例中，外设接口104，处理器103以及存储控制器102可以在单个芯片中实现。在其他一些实例中，他们可以分别由独立的芯片实现。
显示单元105在所述客户端100与用户之间提供一个交互界面(例如用户操作界面)或用于显示图像数据给用户参考。在本实施例中，所述显示单元105可以是液晶显示器或触控显示器。若为触控显示器，其可为支持单点和多点触控操作的电容式触控屏或电阻式触控屏等。支持单点和多点触控操作是指触控显示器能感应到来自该触控显示器上一个或多个位置处同时产生的触控操作，并将该感应到的触控操作交由处理器103进行计算和处理。
输入输出单元106用于提供给用户输入数据实现用户与所述客户端100的交互。所述输入输出单元106可以是，但不限于，鼠标和键盘等，所述键盘可以是虚拟键盘。
第一实施例
请参考图2，本发明较佳实施例提供的页岩CT成像装置200的功能模块示意图。页岩CT成像装置200包括：获得模块201、第二计算模块202、离散化处理模块203、生成模块204、第一计算模块205、降噪模块206以及迭代模块207。
获得模块201，用于获得对所述页岩不同投影角度的入射光强以及出射光强。用已知入射光强的入射光线照射作为目标物体的页岩，并按照入射光线入射角度依次获取反射光线的反射光强。例如，用感光平面获取反射光强。
第二计算模块202，用于依据亮场暗场数据对收集到的出射光强进行预处理。具体为，计算所述出射光强与所述入射光强的比值。以及对计算出的出射光强与入射光强的比值取对数。当入射光线沿着一个路径通过物体后强度的损失可以表示为如下公式：
I1/I0＝exp{-∫Lf(x)dx}，
上式中I0表示入射光强，I1表示出射光强，L表示射线路径，f(x)表示作为目标物体的页岩沿射线路径上对光强吸收的函数。光强吸收的函数f(x)为我们最终需要得出的CT成像的物理意义，因此我们在收集到出射光强后会进行以上预处理。
离散化处理模块203，对所述比值的对数离散化处理后，从而获取所述投影数据。因为无论是入射光强还是出射光强，其在物理上的意义均可视为连续函数，为了后续的处理，需将出射光强与入射光强的比值的对数离散化处理。获得后续需被处理的投影数据。
生成模块204，生成页岩的投影矩阵，并将其存储。
本实施例中，根据所有的参与生成投影数据的入射光线以及作为目标物体的页岩，按照预设的计算投影矩阵规则，生成投影矩阵。需要说明的是，为了方便说明，投影矩阵用A表示，投影矩阵中位于第i行第j列的元素用ai，j表示。具体地，预设的计算投影矩阵规则可以但不限于为，视入射光线为一条没有宽度的射线，将作为目标物体的页岩离散化成若干个像素；测算每一条入射光线与每一个像素点交叉部分的长度，并将第i条入射光线与第j个像素点交叉部分的长度存储为ai，j。
第一计算模块205，用于依据所述投影矩阵生成随机列向量分布以及随机行向量分布。
本实施例中，如图3所述，第一计算模块205包括第一获取子模块2051以及运算子模块2052。
第一获取子模块2051，可以用于根据投影矩阵，利用公式：
获取投影矩阵中每一列列向量的二范数的平方，其中代表投影矩阵中的第j列列向量的二范数的平方，i代表行数，j代表列数，ai，j表示投影矩阵中第i行第j列的元素，m代表投影矩阵的总行数；
第一获取子模块2051，还可以用于根据投影矩阵，利用公式：
获取投影矩阵中每一行行向量的二范数的平方，其中代表投影矩阵中的第i行行向量的二范数的平方，i代表行数，j代表列数，ai，j表示投影矩阵中第i行第j列的元素，n代表投影矩阵总的列数；
第一获取子模块2051，还可以用于根据投影矩阵，利用公式：
获取投影矩阵的Frobenius范数的平方，其中，表示投影矩阵的Frobenius范数的平方，ai，j表示投影矩阵中第i行第j列的元素，表示定义为，n代表投影矩阵总的列数，m代表投影矩阵的总行数。
运算子模块2052，可以用于依据列向量的二范数的平方以及投影矩阵的Frobenius范数的平方，利用公式：
经过运算获取所述随机列向量分布，其中，a(j)表示所述投影矩阵的第j列列向量，表示所述第j列列向量的二范数的平方，表示所述Frobenius范数的平方，P(Y＝j)表示随机列向量为所述第j列列向量的概率。
运算子模块2052，还可以用于还用于依据所述Frobenius范数的平方以及所述行向量的二范数的平方，获得所述随机行向量分布：
经过运算获取所述随机行向量分布，其中，a(i)表示所述投影矩阵的第i行行向量，表示所述第i行行向量的二范数的平方，表示所述Frobenius范数的平方，P(Y＝i)表示随机行向量为所述第i行行向量的概率。
降噪模块206，用于根据所述随机列向量分布从所述投影矩阵中获得一随机列向量，基于所述随机列向量对投影数据进行处理，得到降噪后的投影数据。降噪模块206可以降低投影数据中噪音的影响，提高对投影数据处理时的速度以及效果。
本实施例中，如图4所示，降噪模块206包括述获取随机列向量子模块2061、更新子模块2062、第一比较子模块2063以及降噪子模块2064。
获取随机列向量子模块2061，用于根据所述随机列向量分布从投影矩阵的列向量中随机抽取一个列向量作为随机列向量。具体为，获取随机列向量子模块2061执行按照随机列向量分布从投影矩阵中抽取一个列向量，任何一个列向量被抽到的概率都服从随机列向量分布。被抽到概率最大的列向量为后续处理中需要的最佳的列向量。需要说明的是本实施例通过抽取出随机向量的方式替代了传统的将投影矩阵通过把每一个列向量加入更新处理过程，这样大大提高了工作效率。
更新子模块2062，用于依据所述随机列向量以及所述随机列向量的二范数的平方，利用公式：
获取更新后投影数据，其中，K代表更新次数。ZK代表第K次更新后获得的所述更新后投影数据。ZK-1表示第K-1次更新后获得的所述更新后投影数据，K为不小于1的整数，且当K-1＝0时，ZK-1代表所述投影数据。也就是说当第一次进行更新时ZK-1为投影数据；ZK为对投影数据进行了一次更新后的更新后投影数据。a(Y)表示所述随机列向量。表示所述随机列向量的二范数的平方，<a(Y)，ZK-1>表示所述随机列向量与所述投影数据的内积。
第一比较子模块2063，用于比较所述更新次数是否小于预设的第一迭代次数。需要说明的是，为了更好地降低投影数据的噪声，需要多次更新，后续的降噪效果更好，因此预设一个第一迭代次数。第一迭代次数为最终需要更新的总次数。在其他实施例中，第一迭代次数也可以通过建立训练集获得，训练集通过将每一次降噪处理的更新次数以及通过该更新次数进行更新产生的效果存储，通过从训练集不断学习，使第一迭代次数靠近最优值。
本实施例中，当比较模块比较出当前的所述更新次数小于所述第一迭代次数时，获取随机列向量子模块2061将重复执行一次根据所述随机列向量分布从所述投影矩阵中获得的一随机列向量；更新子模块2062将依据更新获取的更新后的投影数据、获取随机列向量子模块2061再次从所述投影矩阵中获得的所述随机列向量以及该随机列向量的二范数的平方，执行获取更新后投影数据的步骤。同时将更新次数加一，第一比较子模块2063再继续执行将加一后的更新次数与第一迭代次数进行比较，直到更新次数等于第一迭代次数，并得到更新后投影数据。例如，第一迭代次数为三，当更新子模块2062执行完第二次更新后，第一比较子模块2063执行比较后发现更新次数二小于第一迭代次数三。获取随机列向量子模块2061再次从投影矩阵中抽取一个随机列向量。更新子模块2062，也依据重新抽取的随机列向量对第二次更新后的投影数据进行更新，并得到第三次更新后的投影数据；此时更新次数与第一迭代次数相等，因此第三次更新后得到的投影数据，为最终需要获得的更新后投影数据。
降噪子模块2064，用于当所述更新次数等于所述第一迭代次数时，依据所述更新后投影数据，利用公式：
B1＝B0-ZT1，
获取所述降噪后的投影数据，其中B1代表所述降噪后的投影数据，B0代表所述投影数据，ZT1代表当所述更新次数等于所述第一迭代次数时，所获得的所述更新后投影数据。完成对投影数据的降噪，以使后续的处理结果更加稳定。
迭代模块207，用于根据所述随机行向量分布从所述投影矩阵中获得一行向量，基于所述随机行向量对所述降噪后的投影数据处理，得到CT图像。
本实施例中，如图5所述，迭代模块207包括获取随机行向量子模块2071、迭代子模块2072、第二比较子模块2073以及第二获取子模块2074。
获取随机行向量子模块2071，用于根据所述随机行向量分布从投影矩阵的行向量中随机抽取一个行向量作为随机行向量。具体为，获取随机行向量子模块2071执行按照随机行向量分布从投影矩阵中抽取一个行向量，任何一个行向量被抽到的概率服从随机行向量分布。被抽到概率最大的行向量为后续处理中需要的最佳的行向量。需要说明的是本实施例通过抽取随机向量的方式替代了传统的将投影矩阵通过把每一个列向量加入更新处理过程，这样大大提高了工作效率。
迭代子模块2072，用于依据所述随机行向量、所述随机行向量的二范数的平方以及所述降噪后的投影数据，利用公式：
获取迭代图像，其中，M代表迭代次数，xM代表第M次迭代后获取的所述迭代图像。xM-1代表第M-1次迭代后的所述迭代图像，M为不小于1的整数，且当M-1＝0时，xM-1代表零向量。也就是说当第一次进行迭代时xM-1为零向量；xM为迭代子模块2072执行一次获取迭代图像步骤后获取到的迭代图像。a(X)表示所述随机行向量，表示所述随机行向量的二范数的平方，<a(X)，xM-1>表示所述随机行向量与所述M次迭代前的所述迭代图像的内积，X代表所述随机行向量在所述投影矩阵中的行数，BX代表所述降噪后的投影数据中位于第X行的元素，代表所述随机行向量的转置。
第二比较子模块2073，用于比较所述迭代次数是否小于预设的第二迭代次数。需要说明的是，为了迭代出的图像效果更加好，需要多次迭代，因此预设一个第二迭代次数。第二迭代次数为最终需要迭代的总次数。在其他实施例中，第二迭代次数也可以通过建立训练集获得，训练集通过将每一次迭代图像的迭代次数以及通过该迭代次数进行迭代产生的效果存储，通过从训练集不断学习，使第二迭代次数逐渐靠近最优值。
本实施例中，当比较模块比较出当前的所述迭代次数小于所述第二迭代次数时，获取随机行向量子模块2071将重复执行一次根据所述随机行向量分布从所述投影矩阵中获得的一随机行向量；迭代子模块2072将依据迭代后获取的迭代图像、获取随机行向量子模块2071再次从所述投影矩阵中获得的所述随机行向量以及该随机行向量的二范数的平方，执行获取迭代图像的步骤。同时将迭代次数加一，第二比较子模块2073再继续执行将加一后的迭代次数与第二迭代次数进行比较，直到迭代次数等于第二迭代次数，并获取迭代图像。例如，第二迭代次数为三，当迭代子模块2072执行完第二次迭代后，第二比较子模块2073执行比较后发现迭代次数二小于第二迭代次数三。获取随机行向量子模块2071再次执行从投影矩阵中抽取一个随机行向量；迭代子模块2072，也依据重新抽取的随机行向量对第二次迭代后获取的迭代图像进行更新，并得到第三次迭代后的迭代图像；此时迭代次数为三次，与第二迭代次数相等，停止继续迭代。
第二获取子模块2074，用于当所述迭代次数等于预设的第二迭代次数时，输出迭代图像，得到最终的CT图像。
第二实施例
请参照图6，图6为本发明较佳实施例提供的一种页岩CT成像方法的流程图。页岩CT成像方法包括：
步骤S101，获得入射光强以及出射光强。具体为，获得对所述页岩不同投影角度的入射光强以及出射光强。
在本发明实施例中，步骤S101可以由获得模块201执行。
步骤S102，计算比值。具体为，计算所述出射光强与所述入射光强的比值。
在本发明实施例中，步骤S102可以由第二计算模块202执行。
步骤S103，获取所述投影数据。具体为，对所述比值的对数离散化处理后，获取所述投影数据。
在本发明实施例中，步骤S103可以由离散化处理模块203执行。
步骤S104，生成页岩的投影矩阵。
在本发明实施例中，步骤S104可以由生成模块204执行。
步骤S105，依据所述投影矩阵生成随机列向量分布以及随机行向量分布。
在本发明实施例中，步骤S105可以由第一计算模块205执行。请参考图7，步骤S105包括一下子步骤：
子步骤S1051，获取所述投影矩阵的列向量的二范数的平方、行向量的二范数的平方以及所述投影矩阵的Frobenius范数的平方。
在本发明实施例中，步骤S1051可以由第一获取子模块2051执行。
子步骤S1052，获得所述随机列向量分布。具体为，依据所述Frobenius范数的平方以及所述列向量的二范数的平方，获得所述随机列向量分布：
其中，a(j)表示所述投影矩阵的第j列列向量，表示所述第j列列向量的二范数的平方，表示所述Frobenius范数的平方，P(Y＝j)表示随机列向量为所述第j列列向量的概率。
在本发明实施例中，步骤S1052可以由运算子模块2052执行。
子步骤S1053，获得所述随机行向量分布。具体为，依据所述Frobenius范数的平方以及所述行向量的二范数的平方，获得所述随机行向量分布：
其中，a(i)表示所述投影矩阵的第i行行向量，表示所述第i行行向量的二范数的平方，表示所述Frobenius范数的平方，P(Y＝i)表示随机行向量为所述第i行行向量的概率。
在本发明实施例中，步骤S1053可以由运算子模块2052执行。
步骤S106，得到降噪后投影数据。具体为根据所述随机列向量分布从所述投影矩阵中获得一随机列向量，基于所述随机列向量对投影数据进行处理，得到降噪后的投影数据。
在本发明实施例中，步骤S106可以由降噪模块206执行。请参考图8，步骤S106包括以下子步骤：
子步骤S1061，根据所述随机列向量分布从所述投影矩阵中获得一随机列向量。
在本发明实施例中，步骤S1061可以由获取随机列向量子模块2061执行。
子步骤S1062，获取更新后投影数据。具体为，依据所述随机列向量以及所述随机列向量的二范数的平方，利用公式：
获取更新后投影数据，其中，K代表更新次数，ZK代表第K次更新后获得的所述更新后投影数据，ZK-1表示第K-1次更新后获得的所述更新后投影数据，K为不小于1的整数，且当K-1＝0时，ZK-1代表所述投影数据，a(Y)表示所述随机列向量，表示所述随机列向量的二范数的平方，a(Y)，ZK-1表示所述随机列向量与所述投影数据的内积；
在本发明实施例中，步骤S1062可以由更新子模块2062执行。
子步骤S1063，更新次数小于第一迭代次数？具体为，比较所述更新次数是否小于预设的第一迭代次数。
在本发明实施例中，子步骤S1063可以由第一比较子模块2063执行。
当所述更新次数小于所述第一迭代次数时，流程进入子步骤S1061，再次获得的一随机列向量。并依据所述更新后投影数据、再次获得的所述随机列向量以及再次获得的所述随机列向量的二范数的平方，重复子步骤S1062，直至所述更新次数等于所述第一迭代次数。当所述更新次数等于所述第一迭代次数时，进入子步骤S1064。
子步骤S1064，获取所述降噪后的投影数据。具体为，依据所述更新后投影数据，利用公式：
B1＝B0-ZT1，
获取所述降噪后的投影数据，其中B1代表所述降噪后的投影数据，B0代表所述投影数据，ZT1代表当所述更新次数等于所述第一迭代次数时，所获得的所述更新后投影数据。
在本发明实施例中，子步骤S1064可以由降噪子模块2064执行。
步骤S107，得到CT图像。具体为，根据所述随机行向量分布从所述投影矩阵中获得一行向量，基于所述随机行向量对所述降噪后的投影数据处理，得到CT图像。
在本发明实施例中，步骤S107可以由迭代模块207执行。本实施例中，如图9所示，所述步骤S107包括以下子步骤：
子步骤S1071，根据所述随机行向量分布从所述投影矩阵中获得一随机行向量。
在本发明实施例中，子步骤S1071可以由获取随机行向量子模块2071执行。
子步骤S1072，获取迭代图像。具体为，依据所述随机行向量、所述随机行向量的二范数的平方以及所述降噪后的投影数据，利用公式：
获取迭代图像，其中，M代表迭代次数，xM代表第M次迭代后获取的所述迭代图像，xM-1代表M-1次迭代后的所述迭代图像，K为不小于1的整数，且当M-1＝0时，xM-1代表零向量，a(X)表示所述随机行向量，表示所述随机行向量的二范数的平方，<a(X)，xM-1>表示所述随机行向量与所述M次迭代前的所述迭代图像的内积，X代表所述随机行向量在所述投影矩阵中的行数，BX代表所述降噪后的投影数据中位于第X行的元素，代表所述随机行向量的转置。
在本发明实施例中，子步骤S1072可以由迭代子模块2072执行。
子步骤S1073，迭代次数小于所述第二迭代次数？具体为，比较所述迭代次数是否小于所述第二迭代次数。
在本发明实施例中，子步骤S1073可以由第二比较子模块2073执行。
当所述迭代次数小于所述第二迭代次数时，流程进入S1071，根据所述随机行向量分布从所述投影矩阵中再次获得的一随机行向量，并依据迭代后的所述迭代图像、再次获得的所述随机行向量以及再次获得的所述随机行向量的二范数的平方，重复所述获取迭代图像的步骤，直至所述迭代次数等于所述第二迭代次数。
当所述迭代次数等于所述第二迭代次数时，停止迭代，流程进入S1074。
子步骤S1074，得到所述CT图像。其中所述CT图像为所述迭代次数对应的所述迭代图像。
在本发明实施例中，子步骤S1074可以由第二获取子模块2074执行。
综上所述，本发明提供一种页岩CT成像方法及装置。通过生成页岩的投影矩阵；依据所述投影矩阵生成随机列向量分布以及随机行向量分布；根据所述随机列向量分布从所述投影矩阵中获得一随机列向量，基于所述随机列向量对投影数据进行处理，得到降噪后的投影数据；根据所述随机行向量分布从所述投影矩阵中获得一行向量，基于所述随机行向量对所述降噪后的投影数据处理，最终得到CT图像。这样处理方法能有效的较低噪音的干扰，且收敛速度快、效果好，对页岩内存储的页岩油气的研究具有重大意义。
在本申请所提供的几个实施例中，应该理解到，所揭露的装置和方法，也可以通过其它的方式实现。以上所描述的装置实施例仅仅是示意性的，例如，附图中的流程图和框图显示了根据本发明的多个实施例的装置、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段或代码的一部分，所述模块、程序段或代码的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。也应当注意，在有些作为替换的实现方式中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个连续的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或动作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。
另外，在本发明各个实施例中的各功能模块可以集成在一起形成一个独立的部分，也可以是各个模块单独存在，也可以两个或两个以上模块集成形成一个独立的部分。
所述功能如果以软件功能模块的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。
需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。
以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。

Claims (4)

1.一种页岩CT成像方法，其特征在于，所述方法包括：
生成页岩的投影矩阵；
获取所述投影矩阵的列向量的二范数的平方、行向量的二范数的平方以及所述投影矩阵的Frobenius范数的平方；
依据所述Frobenius范数的平方以及所述列向量的二范数的平方，获得随机列向量分布：
其中，a(j)表示所述投影矩阵的第j列列向量，表示所述第j列列向量的二范数的平方，表示所述Frobenius范数的平方，P(Y＝j)表示随机列向量为所述第j列列向量的概率；
依据所述Frobenius范数的平方以及所述行向量的二范数的平方，获得随机行向量分布：
其中，a(i)表示所述投影矩阵的第i行行向量，表示所述第i行行向量的二范数的平方，表示所述Frobenius范数的平方，P(Y＝i)表示随机行向量为所述第i行行向量的概率；
依据所述随机列向量以及所述随机列向量的二范数的平方，利用公式：
获取更新后投影数据，其中，K代表更新次数，ZK代表第K次更新后获得的所述更新后投影数据，ZK-1表示第K-1次更新后获得的所述更新后投影数据，K为不小于1的整数，且当K-1＝0时，ZK-1代表所述投影数据，a(Y)表示所述随机列向量，表示所述随机列向量的二范数的平方；在K-1＝0时，<a(Y)，ZK-1>表示所述随机列向量与所述投影数据的内积；K-1≠0时，<a(Y)，ZK-1>表示所述随机列向量与第K-1次更新后获得的所述更新后投影数据的内积；
比较所述更新次数是否小于预设的第一迭代次数；
当所述更新次数小于所述第一迭代次数时，根据所述随机列向量分布从所述投影矩阵中再次获得一随机列向量，并依据所述更新后投影数据、再次获得的所述随机列向量以及再次获得的所述随机列向量的二范数的平方，重复所述获取更新后投影数据的步骤，直至所述更新次数等于所述第一迭代次数；
当所述更新次数等于所述第一迭代次数时，依据所述更新后投影数据，利用公式：
B1＝B0-ZT1，
获取降噪后的投影数据，其中B1代表所述降噪后的投影数据，B0代表所述投影数据，ZT1代表当所述更新次数等于所述第一迭代次数时，所获得的所述更新后投影数据；
依据所述随机行向量、所述随机行向量的二范数的平方以及所述降噪后的投影数据，利用公式：
获取迭代图像，其中，M代表迭代次数，xM代表第M次迭代后获取的所述迭代图像，xM-1代表M-1次迭代后的所述迭代图像，M为不小于1的整数，且当M-1＝0时，xM-1代表零向量，a(X)表示所述随机行向量，表示所述随机行向量的二范数的平方，<a(X)，xM-1>表示所述随机行向量与所述M次迭代前的所述迭代图像的内积，X代表所述随机行向量在所述投影矩阵中的行数，BX代表所述降噪后的投影数据中位于第X行的元素，代表所述随机行向量的转置；
比较所述迭代次数是否小于预设的第二迭代次数；
当所述迭代次数小于所述第二迭代次数时，根据所述随机行向量分布从所述投影矩阵中再次获得一随机行向量，并依据迭代后的所述迭代图像、再次获得的所述随机行向量以及再次获得的所述随机行向量的二范数的平方，重复所述获取迭代图像的步骤，直至所述迭代次数等于所述第二迭代次数；
当所述迭代次数等于所述第二迭代次数时，停止迭代，得到CT图像，其中所述CT图像为所述迭代次数对应的所述迭代图像。
2.如权利要求1所述页岩CT成像方法，其特征在于，在获取投影矩阵的步骤之前，所述方法还包括：
获得对所述页岩不同投影角度的入射光强以及出射光强；
计算所述出射光强与所述入射光强的比值；
对所述比值的对数离散化处理后，获取所述投影数据。
3.一种页岩CT成像装置，其特征在于，所述装置包括：
生成模块，用于生成页岩的投影矩阵；
第一计算模块，用于：
获取所述投影矩阵的列向量的二范数的平方、行向量的二范数的平方以及所述投影矩阵的Frobenius范数的平方；
依据所述Frobenius范数的平方以及所述列向量的二范数的平方，获得随机列向量分布：
其中，a(j)表示所述投影矩阵的第j列列向量，表示所述第j列列向量的二范数的平方，表示所述Frobenius范数的平方，P(Y＝j)表示随机列向量为所述第j列列向量的概率；
依据所述Frobenius范数的平方以及所述行向量的二范数的平方，获得随机行向量分布：
其中，a(i)表示所述投影矩阵的第i行行向量，表示所述第i行行向量的二范数的平方，表示所述Frobenius范数的平方，P(Y＝i)表示随机行向量为所述第i行行向量的概率；
降噪模块，用于：
依据所述随机列向量以及所述随机列向量的二范数的平方，利用公式：
获取更新后投影数据，其中，K代表更新次数，ZK代表第K次更新后获得的所述更新后投影数据，ZK-1表示第K-1次更新后获得的所述更新后投影数据，K为不小于1的整数，且当K-1＝0时，ZK-1代表所述投影数据，a(Y)表示所述随机列向量，表示所述随机列向量的二范数的平方；在K-1＝0时，<a(Y)，ZK-1>表示所述随机列向量与所述投影数据的内积；K-1≠0时，<a(Y)，ZK-1>表示所述随机列向量与第K-1次更新后获得的所述更新后投影数据的内积；
比较所述更新次数是否小于预设的第一迭代次数；
当所述更新次数小于所述第一迭代次数时，根据所述随机列向量分布从所述投影矩阵中再次获得一随机列向量，并依据所述更新后投影数据、再次获得的所述随机列向量以及再次获得的所述随机列向量的二范数的平方，重复执行所述获取更新后投影数据的步骤，直至所述更新次数等于所述第一迭代次数；
当所述更新次数等于所述第一迭代次数时，依据所述更新后投影数据，利用公式：
B1＝B0-ZT1，
获取降噪后的投影数据，其中B1代表所述降噪后的投影数据，B0代表所述投影数据，ZT1代表当所述更新次数等于所述第一迭代次数时，所获得的所述更新后投影数据；
迭代模块，用于：
依据所述随机行向量、所述随机行向量的二范数的平方以及所述降噪后的投影数据，利用公式：
获取迭代图像，其中，M代表迭代次数，xM代表第M次迭代后获取的所述迭代图像，xM-1代表M-1次迭代后的所述迭代图像，M为不小于1的整数，且当M-1＝0时，xM-1代表零向量，a(X)表示所述随机行向量，表示所述随机行向量的二范数的平方，<a(X)，xM-1>表示所述随机行向量与所述M次迭代前的所述迭代图像的内积，X代表所述随机行向量在所述投影矩阵中的行数，BX代表所述降噪后的投影数据中位于第X行的元素，代表所述随机行向量的转置；
比较所述迭代次数是否小于预设的第二迭代次数；
当所述迭代次数小于所述第二迭代次数时，根据所述随机行向量分布从所述投影矩阵中再次获得一随机行向量，并依据迭代后获取的所述迭代图像、再次获得的所述随机行向量以及再次获得的所述随机行向量的二范数的平方，重复执行所述获取迭代图像的步骤，直至所述迭代次数等于所述第二迭代次数；
当所述迭代次数等于所述第二迭代次数时，停止迭代，得到CT图像，其中所述CT图像为所述迭代次数对应的所述迭代图像。
4.如权利要求3所述页岩CT成像装置，其特征在于，所述装置还包括：
获得模块，用于获得对所述页岩不同投影角度的入射光强以及出射光强；
第二计算模块，用于计算所述出射光强与所述入射光强的比值；
离散化处理模块，用于对所述比值的对数离散化处理后，获取所述投影数据。
CN201710165715.6A 2017-03-20 2017-03-20 页岩ct成像方法及装置 Active CN106875334B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710165715.6A CN106875334B (zh) 2017-03-20 2017-03-20 页岩ct成像方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710165715.6A CN106875334B (zh) 2017-03-20 2017-03-20 页岩ct成像方法及装置

Publications (2)

Publication Number Publication Date
CN106875334A CN106875334A (zh) 2017-06-20
CN106875334B true CN106875334B (zh) 2019-10-25

Family

ID=59173136

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710165715.6A Active CN106875334B (zh) 2017-03-20 2017-03-20 页岩ct成像方法及装置

Country Status (1)

Country Link
CN (1) CN106875334B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112233156B (zh) * 2020-10-14 2022-02-15 俐玛精密测量技术(苏州)有限公司 微纳米ct投影数据的中心切片对齐方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101073499A (zh) * 2007-04-20 2007-11-21 清华大学 一种基于支撑域的ct图像重建系统
CN103413280A (zh) * 2013-08-26 2013-11-27 南方医科大学 一种低剂量x射线ct图像重建方法
CN104240210A (zh) * 2014-07-21 2014-12-24 南京邮电大学 基于压缩感知的ct图像迭代重建方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101073499A (zh) * 2007-04-20 2007-11-21 清华大学 一种基于支撑域的ct图像重建系统
CN103413280A (zh) * 2013-08-26 2013-11-27 南方医科大学 一种低剂量x射线ct图像重建方法
CN104240210A (zh) * 2014-07-21 2014-12-24 南京邮电大学 基于压缩感知的ct图像迭代重建方法

Also Published As

Publication number Publication date
CN106875334A (zh) 2017-06-20

Similar Documents

Publication Publication Date Title
US20210004724A1 (en) Source identification by non-negative matrix factorization combined with semi-supervised clustering
Chamberland et al. Deep neural decoders for near term fault-tolerant experiments
Domingues et al. An adaptive multiresolution scheme with local time stepping for evolutionary PDEs
US11106790B2 (en) Dimensionality reduction of computer programs
Krishnan et al. Efficient preconditioning of laplacian matrices for computer graphics
Prato et al. Efficient deconvolution methods for astronomical imaging: algorithms and IDL-GPU codes
EP3218811B1 (en) Testing insecure computing environments using random data sets generated from characterizations of real data sets
CN106918838A (zh) 起伏地表条件下高斯束偏移成像方法及装置
Hutter et al. Improved hdrg decoders for qudit and non-abelian quantum error correction
Selig et al. Denoising, deconvolving, and decomposing photon observations-Derivation of the D3PO algorithm
Yu et al. Efficient euclidean projections onto the intersection of norm balls
CN106875334B (zh) 页岩ct成像方法及装置
Diener et al. Simulating neutron star mergers with the Lagrangian Numerical Relativity code SPHINCS_BSSN
Yuan et al. Adaptively sketched Bregman projection methods for linear systems
Altafini et al. An edge centrality measure based on the Kemeny constant
Uribe et al. Horseshoe priors for edge-preserving linear Bayesian inversion
Varadhan et al. Squared extrapolation methods (SQUAREM): A new class of simple and efficient numerical schemes for accelerating the convergence of the EM algorithm
Li et al. A learnable group-tube transform induced tensor nuclear norm and its application for tensor completion
Govindarajan et al. A Fast Algorithm for Computing Macaulay Null Spaces of Bivariate Polynomial Systems
Tursunov et al. Comparative variational studies of 0+ states in three-α models
Grant B-spline methods for radial Dirac equations
CN107589451B (zh) 地震全波形反演方法及装置
Camblong et al. Conformal tightness of holographic scaling in black hole thermodynamics
Pritchard et al. Towards Practical Large-Scale Randomized Iterative Least Squares Solvers through Uncertainty Quantification
Willwerscheid Empirical Bayes matrix factorization: methods and applications

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant