本发明涉及多孔介质孔隙结构特征定量表征，具体涉及一种基于树形拓扑结构骨架的多孔介质迂曲度计算方法。
背景技术：
1、多孔介质材料内部的微观孔隙结构特征对其各类宏观物理特性均有着直接且十分重要的影响，如抗拉/压强度、导电性及渗透特性等。岩土体等天然多孔介质材料内部往往分布有复杂的孔隙结构，准确且定量的提取表征各类孔隙结构特征参数对于高精度数字岩心模型的构建以及后续孔隙尺度模拟至关重要。
2、用于多孔介质孔隙结构特征表征的参数指标主要有孔隙度、连通度、孔径分布和迂曲度等。其中，对于孔隙度、连通度及孔径分布等较低阶孔隙结构特征参数而言，现在已经有十分成熟的评估方法和计算提取软件；而对于迂曲度这类较高阶孔隙结构特征参数的计算，目前仅存在少数几种计算方法，它们或基于简化的孔隙网络模型，或使用修剪后含有闭合环的孔隙骨架，只能计算得到某段或者少量几条路径的孔隙迂曲度，而无法对多孔介质孔隙结构的迂曲度特征进行一个全面且准确的提取评估。此外，现有的每种方法都有其自身相应的使用限制条件和适用对象，也即具有较弱的通用扩展性，从而导致对于同一多孔介质材料孔隙结构迂曲度的计算，不同方法的提取结果存在较大差异。
3、综上，可以看出对于孔隙结构迂曲度的计算还尚未达成统一共识，这导致无法对不同方法计算得到的孔隙结构迂曲度结果进行可靠且有效的对比和评估，并势必会对后续基于该计算结果的多孔介质材料物理过程的模拟仿真及物性参数的预测计算产生不利影响。
技术实现思路
1、针对上述问题，本发明公开了一种基于树形拓扑结构骨架的多孔介质迂曲度计算方法，可以准确且全面的评估多孔介质材料内部高阶孔隙结构特征，实现多孔介质材料内部孔隙结构迂曲度特征的定量化提取表征，进而提高后续物理过程模拟仿真和物性参数预测计算结果的精度和可靠度。
2、根据本发明的目的提出的一种基于树形拓扑结构骨架的多孔介质迂曲度计算方法，包括以下步骤：
3、步骤一：通过成像技术获取多孔介质材料内部结构图像，并对其进行预处理；
4、步骤二：利用阈值自动化确定算法，计算最优像素灰度阈值，完成图像中孔隙结构的提取分割，生成只包含孔隙结构与固体基质相的二值图像，并对图像中存在的与物理事实相悖的孤立固体基质相进行滤除优化；
5、步骤三：基于步骤二分割优化后的二值图像，将其直接导入avizo一体式图像分析平台，使用内置centerlinetree功能模块提取与真实孔隙结构相对应的树形拓扑结构骨架，并通过调整模块中的slope和zeroval参数来对所生成骨架的具体细节进行控制；所提取的树形拓扑结构骨架以圆点和曲线段表示，圆点代表了表征孔隙结构拓扑特征的两类关键位置节点：分支点和端点，曲线段则代表了两关键位置点间的真实连通路径；树形拓扑结构骨架中任意两端点间存在且仅存在一条唯一确定的路径能够将两端点进行连通；
6、步骤四：使用avizo一体式图像分析平台内置centerlinetree功能模块提取能够刻画树形拓扑结构骨架的四类基本数据，包括每一关键位置节点的空间位置坐标、每一关键位置节点的配位数也即与该节点相连接的曲线段的数目、每一曲线段的真实长度以及组成每一曲线段的两个关键位置节点的编号，并对基本数据进行进一步的分离提取和存储形式调整处理；
7、步骤五：确定树形拓扑结构骨架端点-端点间连通路径并进行长度计算；
8、步骤六：结合步骤五计算的树形拓扑结构骨架端点-端点间连通路径的真实长度以及两端点间的最短长度计算迂曲度；
9、步骤七：结合步骤六的计算结果，进一步筛选出所有能够形成有效渗透流动通道的连通路径的迂曲度，并对其进行统计分析，以评估该多孔介质孔隙结构的渗透特性。
10、优选的，成像技术包括扫描电子显微镜(sem)、透射电子显微镜(tem)、聚焦离子束显微镜(fib-sem)、核磁共振成像(nmri)及计算机断层扫描(ct)，直接获取多孔介质材料内部结构的二维或三维图像。
11、优选的，步骤一中，图像预处理包括对其进行噪点滤除、整体亮度调整及不同相区域间对比度增强，以改善原始扫描图像质量。
12、优选的，步骤四中，对树形拓扑结构骨架的四类基本数据进行进一步的分离提取和存储形式调整处理的方法，包括以下步骤：
13、s4.1、对所有关键位置节点按照配位数从小到大进行排序，筛选出配位数为1的所有位置节点，也即骨架中的全部端点，并将其编号信息单独存储于一维数组end[:]中，例如end(m)＝i表示骨架中第m个端点的编号为i；
14、s4.2、对组成骨架的曲线段进行遍历，并将其真实长度存储于二维数组len[:,:]中，其中具体存储位置由组成该曲线段的关键位置节点编号来进行确定，例如对于由关键位置节点编号i和j所组成的真实长度为l的曲线段而言，令len[i,j]＝l，此外考虑到关键位置节点编号的先后顺序对长度值索引的影响，同时令len[j,i]＝l，也即用于存储曲线段真实长度的len[:,:]为对称数组len[i,j]＝len[j,i]；
15、s4.3、遍历骨架中的关键位置节点，以节点编号作为判断基准，逐一搜索对比曲线段的两关键位置节点编号，以找到与遍历节点直接相连接的其他关键位置节点编号信息，并将其依次顺序存储于二维数组con[:,:]中，例如con[m,1:n]＝[k1,k2,…,kn]表示骨架中第m个关键位置节点分别与编号为k1,k2,…,kn的n个关键位置节点相连接，其中n等于该关键位置节点的配位数。
16、优选的，步骤五中，树形拓扑结构骨架端点-端点间连通路径确定及长度计算方法，包括以下步骤：
17、s5.1、建立一维动态序列数组seq[:]，用于存储连接任意两端点的全部关键位置节点编号信息，其中数组首尾元素为端点的编号，中间元素为连接两端点的分叉点的编号，由于不同端点间连通路径的差异，故该数组的大小是动态可变的；
18、s5.2、对存储骨架中全部端点编号信息的一维数组end[:]元素逐一进行遍历，并将端点元素编号end[m]取出放入一维动态序列数组seq[:]的首位，即seq[1]＝end[m]；
19、s5.3、以seq[1]作为查询索引值，通过二维数组con[:,:]提取与端点seq[1]直接相连接的关键位置节点的编号信息，并将该编号顺序存入序列数组seq[:]的次位，为便于后续的引用说明，我们选择使用o来对次位进行表示，即seq[o]＝con[seq[1],1]，此处o＝2；
20、s5.4、检查判断seq[o]的配位数信息，若其等于1，则表明该关键位置节点为端点，也即完成一条端点-端点间路径的确定，随后对组成该路径的所有曲线段长度进行累加求和，以得到两端点间的真实长度；紧接着将一维动态序列数组seq[:]中的元素全部清空，并直接跳转到s5.2开始遍历下一端点编号；
21、s5.5、若其大于1，则表明该关键位置节点为分叉点，以seq[o]作为查询索引值，再次通过二维数组con[:,:]提取与该分叉点直接相连接的关键位置节点的编号信息，也即得到一维数组con[seq[o],:]；此时，对该数组中所包含的关键位置节点编号信息进行循环遍历，当遍历节点编号值异于动态序列数组中查询索引值seq[o]的前一位时，也即seq[o-1]，将该遍历到的节点编号继续顺序存入序列数组seq[:]的次位，也即seq[o+1]，同时将该遍历到的节点编号从一维数组con[seq[o],:]中删除；随后，直接跳转到s5.4，并将s5.4中的seq[o]替换为seq[o+1]；
22、s5.6、不断迭代循环上述步骤，直至存储骨架中全部端点编号信息的一维数组end[:]元素全部遍历完毕，从而最终完成骨架中任意两端点间连通路径的确定。
23、优选的，在s5.4中，当完成一条端点-端点路径的确定之后，需顺次完成该两端点间真实路径长度的计算，具体计算步骤为：首先，获取一维动态序列数组seq[:]中节点元素的总数目p；其次，由首位开始从数组seq[:]中按顺序依次取出两个相邻元素节点seq[q]和seq[q+1]；再次，以seq[q]和seq[q+1]作为位置索引，从存储曲线段长度的二维数组len[:,:]中，提取得到由上述两节点所组成的曲线段的长度len[seq[q],seq[q+1]]；然后，不断重复循环相邻元素节点的取出和曲线段长度的提取，并将每次提取得到的曲线段长度值进行累加，共计循环累加p-1次；最后，得到两端点间由p-1条曲线段所组成的连通路径的总长度lt，并将其存储于二维数组t_len[:,:]中，其中具体存储位置由组成该连通路径的两个端点编号来进行确定，也即t_len[i,j]＝lt，其中i和j分别为一维动态序列数组seq[:]中首位和末位端点编号。
24、优选的，步骤六中，迂曲度的计算方法，包括以下步骤：
25、s6.1、对存储全部端点编号信息的一维数组end[:]进行内外两层的嵌套遍历，为避免重复计算，内层遍历的起始位置在外层遍历的基础上进行后移一位处理，例如对于由r个端点元素所组成的一维数组end[:]，当外层遍历至第m个端点时(1≤m≤r-1)，内层遍历的起始位置为m+1；
26、s6.2、使用内外层所遍历到的端点编号信息作为位置索引，通过已有的数据集得到两端点的空间位置坐标，并直接计算其相应的欧几里得长度；
27、s6.3、再次使用内外层所遍历到的端点编号信息作为位置索引，通过存储两端点真实长度的二维数组t_len[:,:]提取其真实长度，并除以欧几里得长度，从而计算得到该连通路径的最终迂曲度数据；
28、s6.4、不断循环上述过程，直至完成端点-端点间所有连通路径迂曲度的计算，通过对该计算结果进行统计分析，定量评估多孔介质孔隙结构空间分布的复杂程度。
29、优选的，步骤七中，筛选所有能够形成有效渗透流动通道的连通路径的迂曲度的方法，包括以下步骤：
30、s7.1、选定需要分析研究的渗透方向，以x方向为例；
31、s7.2、通过已有的端点坐标数据集，分别筛选出x＝0和x＝lx的端点，其中lx为多孔介质孔隙结构在x方向上的长度；
32、s7.3、以上述筛选出的端点为索引，从s6所得到的迂曲度计算结果中进一步提取得到由横坐标为x＝0和x＝lx的端点所组成的有效连通路径的迂曲度结果，并对其进行统计分析，以评估该多孔介质孔隙结构的渗透特性。
33、与现有技术相比，本发明公开的一种基于树形拓扑结构骨架的多孔介质迂曲度计算方法的优点是：
34、本发明综合考虑孔隙结构的空间拓扑结构特征以及真实流体在其中的流动传输过程，提出了一种基于树形拓扑结构骨架，能够用于计算各类多孔介质材料孔隙结构迂曲度的通用方法，可以准确且全面的评估多孔介质材料内部高阶孔隙结构特征，实现多孔介质材料内部孔隙结构迂曲度特征的定量化提取表征，进而提高后续物理过程模拟仿真和物性参数预测计算结果的精度和可靠度。
1.一种基于树形拓扑结构骨架的多孔介质迂曲度计算方法，其特征在于，包括以下步骤：
2.根据权利要求1所述的一种基于树形拓扑结构骨架的多孔介质迂曲度计算方法，其特征在于，成像技术包括扫描电子显微镜(sem)、透射电子显微镜(tem)、聚焦离子束显微镜(fib-sem)、核磁共振成像(nmri)及计算机断层扫描(ct)，直接获取多孔介质材料内部结构的二维或三维图像。
3.根据权利要求1所述的一种基于树形拓扑结构骨架的多孔介质迂曲度计算方法，其特征在于，步骤一中，图像预处理包括对其进行噪点滤除、整体亮度调整及不同相区域间对比度增强，以改善原始扫描图像质量。
4.根据权利要求1所述的一种基于树形拓扑结构骨架的多孔介质迂曲度计算方法，其特征在于，步骤四中，对树形拓扑结构骨架的四类基本数据进行进一步的分离提取和存储形式调整处理的方法，包括以下步骤：
5.根据权利要求1所述的一种基于树形拓扑结构骨架的多孔介质迂曲度计算方法，其特征在于，步骤五中，树形拓扑结构骨架端点-端点间连通路径确定及长度计算方法，包括以下步骤：
6.根据权利要求5所述的一种基于树形拓扑结构骨架的多孔介质迂曲度计算方法，其特征在于，在s5.4中，当完成一条端点-端点路径的确定之后，需顺次完成该两端点间真实路径长度的计算，具体计算步骤为：首先，获取一维动态序列数组seq[:]中节点元素的总数目p；其次，由首位开始从数组seq[:]中按顺序依次取出两个相邻元素节点seq[q]和seq[q+1]；再次，以seq[q]和seq[q+1]作为位置索引，从存储曲线段长度的二维数组len[:,:]中，提取得到由上述两节点所组成的曲线段的长度len[seq[q],seq[q+1]]；然后，不断重复循环相邻元素节点的取出和曲线段长度的提取，并将每次提取得到的曲线段长度值进行累加，共计循环累加p-1次；最后，得到两端点间由p-1条曲线段所组成的连通路径的总长度lt，并将其存储于二维数组t_len[:,:]中，其中具体存储位置由组成该连通路径的两个端点编号来进行确定，也即t_len[i,j]＝lt，其中i和j分别为一维动态序列数组seq[:]中首位和末位端点编号。
7.根据权利要求6所述的一种基于树形拓扑结构骨架的多孔介质迂曲度计算方法，其特征在于，步骤六中，迂曲度的计算方法，包括以下步骤：
8.根据权利要求1所述的一种基于树形拓扑结构骨架的多孔介质迂曲度计算方法，其特征在于，步骤七中，筛选所有能够形成有效渗透流动通道的连通路径的迂曲度的方法，包括以下步骤：
