基于三维激光扫描技术的裂缝随机模拟：以库车坳陷库车河剖面为例

Stochastic Fracture Simulation Based on Three-Dimensional Laser Scan Technology: a Case Study of Kuqa River Outcrop in Kuqa Depression

摘要

库车坳陷是塔里木盆地重要的天然气勘探开发目标区，其大北-克深地区含气层系下白垩统埋深超过7 500 m，岩性致密孔隙度平均为4.8%，但仍然获得高产工业气流，岩层内所发育的裂缝成为天然气高产的核心要素. 因此利用合理技术对裂缝数据进行有效表达是明确深层致密砂岩裂缝发育特征及空间分布规律的关键，是有效开发裂缝型气藏前提. 利用三维激光扫描能够获得大数据的优势，运用点数据拼接、去噪、切割处理方法，直接识别点云数据中的有效裂缝信息，并结合野外剖面实测获得的裂缝数据对所识别的裂缝信息进行修正处理，建立目标区域裂缝信息地质知识库. 确定目标区域裂缝发育的空间范围，裂缝密度发育规律，裂缝方位分布状态，为后期随机模拟确定约束条件. 在三维激光扫描点云处理数据的约束下，应用随机模拟技术的fisher函数能够较好的表现裂缝发育区域单条裂缝之间方位偏转现象，体现自然裂缝的随机发育特征. 利用裂缝开度信息，应用Oda计算方法对对每条裂缝的渗流能力进行了计算分析，定量确定裂缝发育区对围岩流体渗流影响的大小.库车坳陷库车河剖面计算结果显示：右侧裂缝簇导致的渗流带宽度可达3.1 m，延伸长度平均为2.6 m；中间区域所展现裂缝簇渗流带宽度1.2 m，延伸长度平均为1.4 m；左侧裂缝簇渗流带宽度约0.9 m，延伸长度仅仅为0.33 m左右，这与野外数据吻合. 与人工对比显示：数据处理面积是人工测量的25倍，所消费的时间仅仅是人工处理的1/10，大大提高了工作效率.由此可见裂缝簇内裂缝密集程度越高，对流体渗流影响就越大，流体在裂缝簇渗流扩展的宽度，影响的范围就越广，较好的表达了库车深层裂缝致密气藏的特征. 说明基于激光扫描技术的裂缝随机模拟技术能够真实的表达裂缝三维空间展布状态，以及对流体渗流的贡献.

关键词

裂缝 / 三维激光扫描 / 随机模拟 / 库车坳陷 / 定量 / 石油地质

Key words

fractures / three-dimensional laser scanning / stochastic modeling / Kuqa Depression / quantitative analysis / petroleum goelogy

由于埋藏深度大，岩层所经历的成岩作用较强，导致岩层相对较致密（丁文龙等，2015；巩磊等，2017；朱筱敏等，2018）. 而致密岩层中所蕴藏的油气资源在我国广为分布（贾承造等，2012），致密油气勘探进一步深入化研究是中国后继油气资源拓展的核心主题（林良彪等，2021）. 研究表明裂缝是改善致密岩层质量的重要因素（Li et al.，2018），裂缝已成为致密油气藏研究的一个重要关切点. 大量学者从裂缝的识别（岩心识别、露头识别）（吕文雅等，2016；鲍怡晨等，2021）、裂缝的探测（常规测井、成像测井、地震资料处理）（邓少贵，2009；黄玉越等，2022）、裂缝模拟（构造应力场计算、三维裂缝空间展布随机模拟）进行了大量富有成效的研究（田宜平，2011；许俊闪，2022）. 其中野外露头最为直观的展现了裂缝全方位的信息，受到大量裂缝研究者的喜爱. 虽然露头承载了大量的裂缝，但由于露头自然条件的限制，很多裂缝信息无法获得，这很大程度上制约了目标体裂缝知识库的建立. 近年来，国内外石油工作者尝试利用三维激光扫描技术获取露头裂缝信息（任宇鹏，2019），建立目标体裂缝信息知识库，较好展现了地质体的裂缝信息，对传统人工测量所获得裂缝信息进行了修正和补充，并较好的应用于油气田的勘探开发研究工作. 虽然三维激光扫描技术存在数据量大，地质信息丰富的特点，能够全方位的识别裂缝优点，但由于分辨率和野外条件的限制，在裂缝表征上仍以下几个明显劣势：（1）由方形面来展示直接识别的裂缝与真实状态存在一定的误差；（2）由于观测的局限，岩芯、测井、地震很难表示出小间距裂缝簇，往往用一个方形面来表示裂缝簇，但实际上一个方形面很难代表小间距的裂缝簇，裂缝簇对流体渗流的影响更难定量表达清楚（图1），如库车河剖面间距2~12 cm裂缝簇，地震和测井是很难表征；（3）野外实际观测裂缝发育，经常表现为多组（由于受不同方向、不同期次构造应力场作用所致），例如，米斯布拉克剖面表现为3组剪切裂缝set1、set2、set3，其中set1、set3互为切割关系，为共轭剪切节理，set2被set1所限制，为晚期剪切节理. 因此对于多方向多期次裂缝，单纯的用方形面表示裂缝组系存在一定的局限性（图2）. 而裂缝随机模拟技术，能够很好的展现裂缝多组系、簇状发育的三维空间展布特征，但需要大量的地质信息建立裂缝知识库，进而选择合理的概率分布函数进行计算分析. 将两种技术耦合：三维激光扫描提供丰富地质信息，建立裂缝知识库；在建立的裂缝知识库上利用随机模拟技术处理知识库的数据，获得符合实际地质状况的裂缝空间模型；进而来定量表征复杂的裂缝信息. 这对裂缝定量表征研究方面具有较好的实用价值和应用前景.

库车坳陷油气资源十分丰富，是塔里木盆地重要的油气产区. 库车坳陷克拉苏构造带新发现了大北14、大北12、博孜12、博孜9等多个气藏，并获得了一系列的突破（杨海军等，2018），是继克拉-克深气田之后的又一个万亿方大气区，研究显示该构造带是超深层的裂缝性致密砂岩气藏. 钻井显示克拉苏构造带下白垩统巴什基奇克组和巴西盖组致密砂岩储层发育大量的裂缝（史超群等，2021；王珂等，2022），裂缝发育影响了气藏分布，进而塔里木油田科技工作者进行地震、地质、工程联合攻关，明确了大北气田巴什基奇克组裂缝发育的地质规律，解决了裂缝对储层性能的改善程度、裂缝对油气运移的影响、裂缝与油气藏产能关系等一系列问题. 但受裂缝属性缝数据来源的不完备的影响，导致目前对深层致密砂岩裂缝定量刻画、空间展布的三维雕刻以及裂缝三维空间知识库的建立仍然存在不足，特别是多条裂缝密集发育层段，裂缝信息损失量更大，裂缝数据不充足的短板暴露的更为明显，严重制约裂缝性致密砂岩气藏勘探开发一体化进程. 本文以库车河剖面为研究对象，基于三维激光扫描仪采集的裂缝点云信息为数据源，配合人工实测数据的校正，经过统计分析建立符合复杂山前构造带的裂缝信息库，借助三维随机模拟技术重构研究对象的三维裂缝模型，为地下储层裂缝预测提供保障.

1 地质背景

库车坳陷位于塔里木盆地北部，经历了多期构造运动（卢华复等，2017），与南天山造山带以大型逆冲褶皱及一系列逆冲断层相接，主要受燕山期以来构造活动的影响（贾承造，1999），尤其是新近纪-第四纪受印度板块和欧亚板块碰撞远程效应，整个研究区陆内前陆冲断活化，冲断构造叠加改造最为强烈，形成了一个在被动大陆边缘基础之上发育起来的中新生代叠合前陆盆地（何登发等，2009）. 中新生代地层完整，厚度变化较大，由北部的克拉苏-依奇克里克构造带向南急剧减薄（李勇等，2017）. 库车坳陷大北-博孜地区天然气主要来源于三叠系（黄山街组、塔里奇克组）和侏罗系（阳霞组、克孜勒努尔组、恰克马克组）烃源岩，主要含气层段位于下白垩统巴什基奇克组和巴西盖组致密砂岩储层.

自新生代以来，库车坳陷经历构造平静期、萌发期、发展期、顶峰期. 构造活动平静期（65~36 Ma），构造活动较弱，岩层的缩短率、上覆岩层的沉积速率均表现为相对低值特征. 由于整体呈现较弱的挤压环境，因此库车坳陷没有发生较大规模的断层及变形作用，沉降有限，沉积物可容空间较小，基本处于构造稳定期；构造活动萌芽期（36~13 Ma），岩层的缩短率开始增加，发育于南天山的逆冲断层沿着不同深度的滑脱面向库车坳陷内部传播，并在库车坳陷北部发育褶皱及断层，库车坳陷所堆积沉积厚度开始大量增加. 由于挤压隆升，库车坳陷北部也由湖相演变为三角洲；构造活动发展期（13~6.5 Ma），南天山强烈的隆升变形，发育南天山逆冲断裂持续向库车坳陷扩展，岩层缩短变形开始加速，坳陷沉降明显加大，发育一系列断层相关褶皱，库车坳陷北部也由三角洲转变为辫状河流沉积；构造活动顶峰期（6.5 Ma至今），南天山快速隆升，岩层缩短率增大为发展期的两倍，并通过断层传播和皱褶作用控制坳陷内部的沉积与构造，该时期以粗粒辫状河沉积为主，岩性主要为砾岩-粗砂.

多期次构造挤压活动导致库车坳陷裂缝大量发育，露头、岩心、成像测井等资料均已得到验证，如大北202井下白垩统巴什基奇克组裂缝密度为1.29 条/m（图3）. 为了获得裂缝三维空间建模的关键属性参数，我们应用三维激光扫描仪详细采集库车河剖面的裂缝数据，并进行了数据处理分析，为地下储层裂缝建模提供复合库车坳陷地质实际的关键参数奠定扎实基础（图4）.

2 三维激光扫描数据处理技术

随着大数据、云计算技术不断发展，野外露头裂缝测量方法也逐渐多样化. 利用传统人工测量裂缝数据信息的弊端已经被逐渐放大：效率低，获得数据量少. 为此传统人工测量露头裂缝数据与三维激光扫描裂缝信息互相补充，能够极大丰富裂缝资料信息库，同时也便于室内数据处理和补充（Cao et al.，2017）.

2.1　三维激光扫描数据处理流程

三维激光扫描仪体积小，可以采集工作人员无法达到区域的裂缝数据，并快速获取三维空间形态或者高精度的地貌数字模型. 我们在野外采集裂缝数据之前会对研究对象进行详细的分析，明确其沉积岩性特征以及所经历的构造变形，依据岩石破裂强度理论了解裂缝典型发育的部位. 结合野外踏勘，初步筛选多个合适的扫描点，并对这些扫描点的岩层出露程度，裂缝发育状态以及色彩的匹配度，选择典型剖面进行裂缝数据采集（曾庆鲁等，2017）. 现场扫描过程中要对不同机位的扫描数据与露头特征和人工测量的标定数据进行反复对比，选择最优3个机位进行精细的数据采集. 获得目标体点云数据信息，为获得更为详尽的裂缝信息需要将3个机位所采集的数据进行坐标转换，建立一个统一的坐标系统，然后将3个机位数据进行镶嵌拼合（点云数据进行配准），并删除冗余的计算数据，然后与数据的RGB信息进行融合处理，并导出数据，最终获得目标露头的点云数据文件（包含X、Y、Z、R、G、B）. 点云数据是离散数据源，并不能提供目标对象的连续性信息，则需要我们根据实际所拍摄的露头照片进行插值处理，并从中提取裂缝信息，并依据随机建模理论建立裂缝三维空间模型（图5）. 本次研究所应用的是奥地利Rigel公司的Riegl VZ1000，其扫描精度为5 mm（也就是100 m范围内，间距5 mm返回一个点数据），扫描距离1.5~1 400 m.

2.2　裂缝位置确定

由于三维激光扫描技术广泛应用，大量学者对三维激光点云数据的信息提取进行探讨（Bellian et al.，2005；Lato et al.，2010；任宇鹏，2019）. 由于原始点云数据量大，原始点云图中无法直接将岩层中裂缝提取出来，需要将点云数据依据我们的研究区需要进行裁剪. 对裁剪后的点云数据进行三维曲面重构，并按照高程值大小进行相应的灰度量化，生成与灰度图像类似的数字模型. 对数字模型进行遍历，得到Z值中最大值Z _max和最小值Z _min . 并对目标区域点云数据依据下式进行灰度量化，获得每个点数据的灰度值大小.

G_i = \frac{Z_i-Z_{\mathrm{m}\mathrm{i}\mathrm{n}}}{Z_{\mathrm{m}\mathrm{a}\mathrm{x}}-Z_{\mathrm{m}\mathrm{i}\mathrm{n}}} \times \mathrm{255}，(1)

式中： G_i表示第i个点数据灰度量化值； Z_i表示第i个点数据的Z值. 然后依据Marr提出的LOG算子边缘提取思路（Marr et al.，1980）：对灰度量化值进行最佳平滑处理，最大程度上抑制噪声，再进行边缘计算. 平滑处理选择二维高斯函数G（r，σ）：

G \left(r,\sigma \right) = \left(\frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }{\sigma }^{\mathrm{2}}}\right) \mathrm{e} \mathrm{x} \mathrm{p} ( - \frac{r^{\mathrm{2}}}{\mathrm{2}{\sigma }^{\mathrm{2}}} ),(2)

式中：σ为尺度因子，r表示数据体f（x，y）中点数据（x，y）为中心的收寻半径，该数据受尺度因子控制，用 G \left(r,\sigma \right)对目标点位灰度值进行卷积运算，得到平滑后的数据体I（x，y）=G（r，σ）×f（x，y）. 然后利用拉普拉斯算子（▽²），对平滑后的数据体I（x，y）计算方向导数得到数据体g（x，y）.

g（x，y）=▽² I（x，y）=▽²［G（r，σ）×f（x，y）］=［▽² G（r，σ）］×f（x，y），（3）

式中：▽² G（r，σ）即为LOG算子：▽² G（r，σ）=（-1/πσ ²）（1-r ²/2σ ²）exp（-r ²/2σ ²）. 求取g（x，y）的所有过零点轨迹即可得数据体f（x，y）的边缘.

我们知道裂缝是岩层破裂的损伤，露头表现较深较窄的特征，剖面上裂缝内基本无反射光. 因此无法通过反射光获得裂缝反射点的三维坐标，导致裂缝附近较小范围内点云密度骤减，导致裂缝边缘形成阶跃边缘，这种阶跃边缘是由z值的突变所形成的. 这种阶跃边界实际上就是裂缝所发育的位置. 由此我们在研究区的灰度数据体中寻找边缘梯度方向上梯度值急剧变化的点，这些梯度值急剧变化的连续点就是我们所阐释的裂缝体，最后根据实际观测结果进行修正.

2.3　裂缝产状数据提取

针对点云数据进行裂缝产状的提取，必须依据扫描点云数据进行扫描露头剖面进行三维虚拟重构. 目前对结构面进行虚拟重构较常用的方法是最小二乘法和克里金法. 最小二乘法是特征值在满足条件a ²+b ²+c ²=1的条件下，依据平面方程得到平面参数，在数据点较多的情况下能够获得最优的平面估计. 但相对三维激光扫描所获得的数据点，数据量太大，针对200 m²的剖面动则好几千万个点源数据，计算极为繁琐，同时获得的结构面异常值较多，很难获得自动平滑性较好的结构面. 而在地质研究应用较多的克里金法，则在数据计算上和结构面的平滑处理上要理想很多. 这主要归因于克里金法可以根据创建的变异函数和协方差函数来估算数据的空间自相关，并进行插值. 本次研究是利用克里金法进行插值计算，来提取裂缝产状信息.

裂缝产状计算分析是建立在剖面结构面的建立和人工实测定位的基础上. 三维激光扫描仪点云数据依据工程坐标系导出的，在获得点云数据之前已经将Y轴设置为正北，X轴设置为正东方向，以机位点为坐标原点，以XOY面为水平面，依照右手法则确定Z轴正方向（Cao，2017）. 可以通过数据转换建立真实的三维空间野外露头模型，然后利用平面上的3点来确定裂缝面的产状. 利用克里金方法确定裂缝面以后，确定裂缝面的法向量为FN（a，b，c）. 根据裂缝面法向量的3个坐标值来判定裂缝面的方位，定量分析裂缝面的倾向和倾角. 假定Y轴正方向向量为n（0，1，0），α为裂缝面法向量在XOY面上的投影向量与Y轴正方向之间的夹角： \mathrm{\alpha } = \mathrm{c} \mathrm{o} {\mathrm{s}}^{-\mathrm{1}} ( \frac{FN\times n}{\sqrt[]{F{N}^{\mathrm{2}}\times n^{\mathrm{2}}}} ) . 若XOY平面法向量为m（0，0，1），β为裂缝面法向量与裂缝面法向量在XOY面上的投影之间的夹角： \mathrm{\beta } = \mathrm{c} \mathrm{o} {\mathrm{s}}^{-\mathrm{1}} ( \frac{FN\times m}{\sqrt[]{F{N}^{\mathrm{2}}\times m^{\mathrm{2}}}} ) . 通过α、β的分析即可获得裂缝面的产状.

3 裂缝随机模拟技术

裂缝随机模拟是在大量裂缝信息知识库的基础上，利用数据统计分析原理对裂缝属性数据进行定量表征的方法.

3.1　裂缝属性分析

完整细致的描述单个裂缝是构建裂缝网络模型，但这些裂缝具有不同形状、坐标、尺寸、方位、开度、及所附带的基质块等属性. 因此模型中要全面展现裂缝信息就需要：空间分析裂缝位置及密度，表现裂缝空间分布特征；统计分析裂缝大小及方位特征，弄清地下裂缝规模，主要的发育方向；为裂缝预测建立合理的概率函数.

3.1.1　裂缝走向和倾角分析

在裂缝空间状态的表征中，裂缝产状描述是最为重要的一环，裂缝方向是控制流体渗流核心参数. 如何符合实际的描述裂缝方位，并将其用数值方法表达出来显的尤为重要. 野外露头显示天然裂缝随着岩石矿物、颗粒、结构等差异会对裂缝扩展方向有一定程度的影响. 野外裂缝观测发现，沿着某一方向发育的裂缝集，裂缝并不是稳定沿着一个方向发展，裂缝之间会存在一定的角度偏转. 例如库车坳陷米斯布拉克露头上所发育的剪切裂缝set3，相邻裂缝之间存在5°左右的差异，单条裂缝之间的最大角度差约30°（图2）. 如何定量描述裂缝之间的角度偏差？Kemeny et al.（2003）提出Fisher分布能很好的描述这一地质现象，裂缝产状概率密度函数能较好的表征这一角度差.

f \left({\phi }^{,},{\theta }^{,}\right) = \frac{k\mathrm{s}\mathrm{i}\mathrm{n}{\phi }^{,}{\mathrm{e}}^{k\mathrm{c}\mathrm{o}\mathrm{s}{\theta }^{,}}}{\mathrm{2}\mathrm{\pi }({\mathrm{e}}^k-\mathrm{1})} , \mathrm{0} \le {\theta }^{,} \le \mathrm{2} \mathrm{\pi } , （4）

式中： {\phi }^{,}为裂缝走向；k为Fisher常量，k值越大，表示数据聚集得越集中； {\theta }^{,}裂缝集平均向量的角偏差.

3.1.2　裂缝密度分析

裂缝密度直接体现裂缝与目标地质体的比例关系，是计算裂缝性储层非均质性最重要的参数之一. 它不仅展现了地质体不同部位裂缝发育的差异性，同时定量反映地质体裂缝密集程度以及裂缝在地质体内的随机、聚类等分布模式. 目前研究人员计算裂缝密度时将裂缝大小考虑进来，通过测量单位长度裂缝条数、单位面积裂缝长度、单位体积裂缝面积来量化裂缝分布的差异程度. 现阶段从岩心、露头、测井资料，利用面积取样获得裂缝密度的方法耗时长、人为统计误差大. 激光扫描数据量丰富、裂缝提取速度快、计算准确的特点能够很好的规避这些问题.

裂缝密度空间分布模型需要根据采集的裂缝密度数据建立数学模型进行计算. 依托已知的资料数据统计分析，构建准确的裂缝密度分布的概率函数. 大量的空间数据分析表明，裂缝分布完全符合空间随机分布（泊松过程模拟）、规则分布（均匀分布模拟）和聚类分布（Dowd et al.，2007；唐永等，2010）（图6）. 其实这3种空间分布模型均有各自使用范畴：（1）随机分布主要用于表征多性质、多方向、强度高的构造作用力对非均质性较强的岩层作用所形成的裂缝，由于多期次、多性质的裂缝叠加复合导致该类裂缝变异性较强，很难用某一个或两个规律总结归纳裂缝发育特征和形成的力学机理，多发育构造变形的较强区域（逆冲构造的根带）或者不同方向构造交汇部位；（2）规则分布主要用于表征单一期次、单一方向、强度中等的构造作用力对相对均质岩层作用所形成的裂缝，由于均质岩层受到的构造作用相对较弱，形成的裂缝发育规模，裂缝之间的间距，裂缝空间方位均较规则，可以较好的用一个或两个规律总结归纳裂缝发育特征和形成的力学机理，多发育于构造变形相对较弱区域（逆冲构造的前锋带）；（3）聚类分布，由于岩层非均质性或构造应力作用的非均质性，导致裂缝在某些区域不发育，在某些区域大量发育，很难用一个模式将裂缝表达完整，一般用嵌套的概率函数来表达这一现象，在构造作用相对较强，构造作用期次多的区域较为常见（董少群，2018）.

3.2　裂缝模型

如何将裂缝与基质较好的耦合在一起？或者说以更好的方式去展现裂缝属性及其裂缝之间的空间交切关系. 目前应用较多的有等效连续模型、离散裂缝网络模型以及复合模型. 其中离散裂缝网络模型最近发展非常迅猛，应用程度逐渐增高，其中最为成熟的Baecher模型是通过示性过程和点过程共同确定目标体中的各条裂缝的空间展布（Baecher，1983； Ivanova et al.，2014）. Baecher模型中默认自然界裂缝中心位置符合随机分布规律，每条裂缝是彼此独立的，单条裂缝所发育的位置与周围裂缝发育程度没有任何关系；裂缝大小的确定是依据裂缝样本数据，通过对裂缝样本进行处理分析，获得裂缝半径的平均数，在平均数的基础上利用对数正态分布对三维空间裂缝大小进行赋值；裂缝大小与裂缝中心位置在统计分析上没有任何关系. 其实我们野外观测发现裂缝的三维空间分布不是无条件随机分布，而是受到多个的聚类特征（优势方位）的约束，比如受断层、褶皱的发育部位、沉积岩性的非均质性控制. 为此研究者对Baecher模型进行改进，用已知的概率分布函数来确定多个裂缝集的分布，随机分布则是在每个裂缝集控制范围内随机生成单条裂缝，并且单条裂缝之间考虑截断与切割.

3.2.1　裂缝集的生成

自然界裂缝多呈现聚集状态. 利用指数分布来描述裂缝集位置及其聚集分布，将随机行走过程带入裂缝集分布表达中，用于确定裂缝的位置. 若多个裂缝集发育于目标研究区，可以用下面概率密度函数分布来表示（Karim et al.，2004；Casini et al.，2016）.

p（x_i）= λ>0，（5）

式中：λ为研究区裂缝集平均密度；x_i 表示第i个裂缝集范围.

3.2.2　单条裂缝的随机生成——裂缝集内的单条裂缝

在生成聚类裂缝的范围下，对每个裂缝集应用Lehmer（1938）提出的同余法. 随机选择一个点，作为初始种子点，计算其它点与初始种子点的距离（d ₁，d ₂，d ₃.....d_n ），并选择任意一个正整数m作为检验数据，若任意两个距离d_x，d_y 分别与正整数m相除所得余数相等，则称d_x，d_y 关于m同余，记为d_x =d_y （mod m）. 那么同余随机位置可以表达为（党发宁，2007）：

d_n =（a×d_{n -1}+c）（mod m），（6）

其中：d_n 是任意位置与初始种子点的距离；a为因子常量2¹⁶+1（32位计算机）；c为常量（0.5+ \sqrt[]{\mathrm{3}}）/m（32位计算机）；m为任意正整数

4 库车河剖面裂缝识别及三维空间展布

针对复杂裂缝系统，空间分布非均质性强的特点，在野外裂缝实测基础上，以三维激光扫描数据处理为突破口，建立符合地质实际的三维裂缝模型. 本次库车河剖面裂缝露头扫描点位于南天山造山带南坡的比尤勒包谷孜鼻状背斜的南翼，产状334°/41°. 露头主要发育下白垩统舒善河组（紫红色、灰紫色粉砂质泥岩夹薄层粉砂岩）、巴西盖组（灰褐色粉砂岩、细砂岩，内夹薄层泥质粉砂岩和泥岩）、巴什基奇克组（紫红色泥岩、粉砂质泥岩为主），其中扫描的核心区为巴西盖组相对均质一些的细砂岩层段（图4）.

4.1　库车河剖面裂缝特征

库车河剖面位于库车坳陷北部单斜带，由于挤压构造作用形成了一系列的构造裂缝，其中95%以上均为剪切裂缝（图7）. 露头剖面上裂缝分布具有明显的聚类特征，裂缝产状以中角度的斜裂缝为主，裂缝沿岩层平面上呈现平行线状分布，裂缝面平直，未见充填. 均发育于砂岩层内，泥岩较少见. 现场实测裂缝延伸长度介于0.3~12.7 m，平均1.08 m，具有双峰状分布特征，延伸长度较小的裂缝介于0.3~1.7 m，平均0.57 m；延伸长度较大的裂缝介于1.5~7.7 m，平均5.24 m. 裂缝间距范围为0.17~2.64 m，平均为0.69 m. 通过对80个裂缝范围统计分析：裂缝方位主要表现为NNW向，少量的呈现NW向（图8），倾角分布呈现较大的差异，高角度裂缝（>70°）约占20%，中角度裂缝（25°~70°）约占60%，低角度裂缝约占20%（图9）.

4.2　激光扫描数据处理

根据激光扫描仪所需要的采光、最小分辨率的限制，选择库车河巴西盖组灰褐色砂岩为目标进行数据采集（图10）. 获得三维点云数据，原始点云数据共二千五百万个节点. 同时由于目标剖面面积大，激光扫描仪旋转角度有限，整个剖面分为4个机位进行扫描. 为重构三维空间剖面必须对原始点云数据进行拼接、去噪、分割处理，进而获得目标清晰、完整的结构数据.

4.2.1　原始点云数据的拼接

三维激光扫描仪受限于扫描距离、扫描分辨率、扫描视角的限制，需要对多次测量的点云数据进行配准，将不同点位的点云数据归并到大地坐标系统中. 点云数据的拼接过程就是对原始扫描点云数据进行旋转、平移等坐标体系的配准. 现阶段点云拼接主要有多种方法：基于测量点的拼接、基于标靶的拼接、基于几何信息的拼接、混合拼接等（王智，2020）. 本次研究是基于几何信息的拼接方法对数据体进行计算分析，根据扫描对象曲率、轮廓线、形态异常区域等明显的形状特征参数进行拼接（图11）. 库车河实际野外地质剖面显示3个“凹槽”与三维点云构面重现为“凹槽”区位置吻合，同时砂岩抗风化能力较强的“凸出”区域与三维点云构面也表现为“凸出”区，较好的吻合（图12）. 为了保证拼接正确性，对裂缝发育密集区宽度也进行了对比分析，吻合程度达98%以上.

4.2.2　点云数据的去噪

由于物体表面反光程度差异、激光的散射性、浮土及植被的干扰等影响，三维激光扫描获取的数据，会造成大量的噪声. 这些噪声点不仅影响数据质量，造成点云数据处理分析出现偏差，而且大大增加无效的点云数据量，使无效的数据计算时间成倍增加. 为了精准聚焦扫描对象，降低处理数据无效耗时，需要对扫描得到的原始数据进行噪声剔除. 噪声主要有3类：（1）孤立点：小而密集且远离目标点云集中区域的点数据；（2）冗余点：由于反射的扰动造成多点位置重叠和极小的位置偏移，造成无效点数据的形成；（3）漂移点：明显远离点数据主体区域，少量散乱、稀疏漂浮于目标区之外（图13）. 对于这些噪声数据，利用点云处理软件Riscan Pro进行手动、距离过滤等方式进行剔除是最简单、直接的方式，对于结构面比较有规律的区域（稳定的斜坡区、褶皱区等）可以通过高斯滤波法、最小二乘滤波法进行去噪确保后续曲面拟合、网格划分的准确性.

4.2.3　点云数据的切割与结构面重构

由于野外剖面范围较大，三维激光扫描仪采集的点数据量非常大，并且剖面曲面非常不规则，导致全局曲面拟合会产生复杂的高次插值多项式. 这不仅仅增加了计算负荷，大大降低点云数据的运算效率，而且计算结果产生误差的可能性大大增高. 为了增加裂缝信息提取的准确度，需要对点云数据进行切割处理，将研究目标聚焦某个范围裂缝发育密集层段. 本次研究利用手动选择的方式完成点云切割处理，从而获得厚约7.7 m、长度27.1 m范围内的细砂岩层段局部点云数据，经切割处理后的点数据共计199 206个. 依据Delaunay三角剖分法对点云数据结构面的重建（冯驰，2019）（图14）

4.3　裂缝面识别及三维空间展布

应用点云所切割数据重构的结构面对裂缝数据进行了提取，识别出52条产状各异，延伸长度大小不均的裂缝，裂缝所在的位置和空间状态用正方形表示，正方形的产状表示了识别出的裂缝产状，相邻之间裂缝用不同颜色表示，裂缝的延伸长度越大，正方形面积越大. 从识别的裂缝来看，由左至右裂缝密度逐渐增加，单条裂缝延伸长度也是成倍增长，右边裂缝密度可达1.2 m/m²，单条裂缝延伸长度最大5.2 m，平均2.61 m；左边裂缝密度仅为0.16 m/m²，单条裂缝延伸长度最大2.3 m，平均0.33 m（图15）. 从左至右裂缝产状也有明显的变化，由中低角度向中高角度转变，可能是横向上岩性非均质性变化影响所致. 所获得裂缝与实际测试所获得裂缝吻合程度较高.

5 确定性约束的裂缝随机模拟

裂缝随机模拟能够很形象表达裂缝信息，能够很好表现自然界中裂缝不稳定性的特征. 但一旦没有大量的裂缝知识库作为约束，很难获得符合地质实际裂缝空间分布. 针对库车露头剖面裂缝较为发育的情况下，利用三维激光扫描技术获得点云数据建立丰富的裂缝知识库，进而为随机模拟方法建立三维空间裂缝模型提供合理的约束.

5.1　库车河剖面裂缝建模

根据三维激光扫描数据所提取的裂缝统计分析认为，裂缝走向主要集中于280°~350°，呈现双峰态的特征其中峰值主要集中于320°~350°（占74%），290°~310°（占21%），其余的占5%. 可以简化为双峰正态分布：WN _a（θ，δ）=0.74WN（338°，10°）+0.21WN（300°，8°），表现为两个裂缝簇. 根据野外实际观测每个裂缝簇内单条裂缝之间存在一定程度的角度偏转，这是裂缝自身发育不稳定所决定的. 本次模拟利用Fisher分布来描述研究区裂缝簇内裂缝之间的不稳定性，由于库车河剖面裂缝簇内单条裂缝之间相对较为集中，根据野外实测数据的统计分析，其中的K值选择138能很好展示这一特性（图16）. 最后在激光扫描获得的裂缝位置、裂缝大小、裂缝密度等信息的约束下，运用随机建模的方法对重构的结构剖面进行了三维空间随机模拟（图17），结果显示：裂缝三维空间展布与野外实际观测到的裂缝信息一直，吻合程度可高达85%以上，数据处理范围是人工测量的25倍，所消费的时间仅仅是人工处理的1/10，大大提高了工作效率（表1）. 三维随机模型不仅表现出了左边裂缝密度低，右边裂缝密度高特征. 同时裂缝簇的延伸长度在每个区域都展现的极为直观，精细将裂缝之间的交切、单条裂缝的弯曲程度刻画的淋漓尽致. 为了更为直观分析裂缝对流体渗流的影响，应用Oda方法对每条裂缝的渗流能力进行了计算分析（Oda，1987；苏小鹏，2020），然后叠加计算在裂缝簇上，获得各个裂缝簇综合渗流效应. 结果显示，裂缝密集程度高的裂缝簇综合渗流效应高，对两侧围岩的影响范围要宽，右侧裂缝簇导致的渗流带宽度可达3.1 m，延伸长度平均为2.6 m；中间区域所展现裂缝簇渗流带宽度1.2 m，延伸长度平均为1.4 m；左侧裂缝簇渗流带宽度约0.9 m，延伸长度平均仅仅为0.33 m左右. 由此可见裂缝簇内裂缝密集程度越高，对流体渗流影响就越大，流体在裂缝簇渗流扩展的宽度，影响的范围就越广（图18）. 我们可以预测经历相同构造作用的地下巴西盖所发育裂缝规律以及对流体渗流状态的影响也应该相同，这为后期指导地下储层裂缝分布规律及流体影响有较大的推动作用.

5.2　建模结果与实际检验

为了检验建模的有效性，我们对库车河剖面局部范围且能进行人工实测裂缝的区域进行了对比. 结果显示，实测裂缝43条，模拟裂缝43条，其中裂缝方位吻合程度高达91%，局部由于岩性的变化，导致裂缝方位小范围的跳跃，致使其存在误差（图19）. 可以说该方法的可靠程度还是比较高的.

6 结论

（1）三维激光扫描技术对露头裂缝信息提取是对传统手工实测剖面获得裂缝信息的方法进一步深化. 该方法的应用能够将获得的裂缝信息量提高几十倍，较好的丰富了传统裂缝地质知识库. 同时裂缝海量数据处理方法的不断提升，很好的满足不同剖面、不同岩性、不同构造部位点云数据拼接、去噪、切割等处理的需要，为合理有效建立数字化三维剖面模型和定量获取裂缝参数提供了保障，极大提高了裂缝属性解析精准度，为地下储层裂缝发育规律分析提供直观可靠的参考证据.

（2）对深层致密砂岩气藏来说，裂缝的发育是改善其物性条件的关键因素. 裂缝属性的定量表征是深层致密砂岩气藏研究的核心，三维空间裂缝随机模拟一直是定量表征裂缝空间属性的热点. 但裂缝随机模拟的难点在于计算分析需要大量的裂缝样本数据来支撑分布函数参数的获取，同时还需要明确裂缝发育区域来约束空间裂缝生成. 这都是测井资料、地震数据体、构造计算很难精准定量描述清楚的，而三维激光扫描所获得的点云数据则很好满足了三维裂缝空间随机模拟的数据需要，同时裂缝随机模拟计算则在点云数据的支持下能够完美真实的展现裂缝结构变异性和裂缝属性非均质性.

（3）库车河剖面野外裂缝数据较为丰富，应用三维激光扫描技术采集裂缝数据，并对该数据集进行精确处理，与人工实测裂缝数据进行对比，计算分析了两者之间的误差，以及计算面积和数据处理的时间：基于三维激光扫描技术的裂缝随机模拟方法在数据处理能力远远高于人工实测计算裂缝属性数据，是人工的25倍以上；而时间上的花费仅仅是人工的1/10. 同时其获得的裂缝数据与人工获得裂缝数据误差基本控制在10%左右. 因此基于三维激光扫描技术的裂缝随机模拟方法具有数据处理能力强、效率高、误差低的特点，具有较高的实用工程应用价值.

参考文献

[1] Baecher, G.B.,1983.Statistical Analysis of Rock Mass Fracturing. Journal of the International Association for Mathematical Geology,15(2):329-348.

[2] Bao, Y.C., Liu, Q.H., Du, X.F., et al.,2021.Division of Glutenite Lithofacies Based on the Trielement of Gravel-Matrix-Fracture. Earth Science, 46(6):2157-2171 (in Chinese with English abstract).

[3] Bellian, J.A., Kerans, C., Jennette, D.C., et al.,2005.Digital Outcrop Models: Applications of Terrestrial Scanning Lidar Technology in Stratigraphic Modeling. Journal of Sedimentary Research,75(2):166-176.

[4] Cao, T., Xiao, A.C., Wu, L., et al.,2017. Automatic Fracture Detection Based on Terrestrial Laser Scanning Data: A New Method and Case Study. Computers & Geosciences,106:209-216.

[5] Casini, G., Hunt, D.W., Monsen, E., et al., 2016. Fracture Characterization and Modeling from Vortual Outcrops. AAPG Bulletin,100(1):41-61.

[6] Dang, F.N., Han, W.T., Zheng, Y.N.,2007. 3D Numerical Simulation of Failure Process of Concrete. Chinese Journal of Computational Mechanics, 24(6):829-833 (in Chinese with English abstract).

[7] Deng, S.G., Li, Z.Q.,2009. Simulation of Array Laterolog Response of Fracture in Fractured Reservoir. Earth Science, 34(5):841-847 (in Chinese with English abstract).

[8] Ding, W.L., Yin, S., Wang, X.H., et al.,2015. Assessment Method and Characterization of Tight Sandstone Gas Reservoir Fractures. Earth Science Frontiers,22(4):173-187 (in Chinese with English abstract).

[9] Dong, S.Q., Zeng, L.B., Xu, C.S., et al., 2018. Some Progress in Reservoir Fracture Stochastic Modeling Researc. OGP, 53(3):625-642 (in Chinese with English abstract).

[10] Dowd, P.A., Xu, C., Mardia, K.V., et al.,2007.A Comparison of Methods for the Stochastic Simulation of Rock Fractures. Mathematical Geology,39(7):697-714.

[11] Feng, C., Zhang, L.,2019.Application of Delaunay Triangulation Method in Noise Map.Geomatics & Spatial Information Technology,42(10):80-82 (in Chinese with English abstract).

[12] Gong, L.,Gao, M.Z., Zeng, L.B., et al.,2017. Controlling Factors on Fracture Development in the Tight Sandstone Reservoirs: A Case Study of Jurassic-Neogene in the Kuqa Foreland Basin. Natural Gas Geoscience,28(2):199-208 (in Chinese with English abstract).

[13] He, D.F., Zhou, X.Y., Yang, H.J., et al.,2009. Geological Structure and its Controls on Giant Oil and Gas Fields inKuqa Depression, Tarim Basin. Geotectonica et Metallogenia,2009,33(1):19-32(in Chinese with English abstract).

[14] Huang, Y.Y., Wang, G.W., Song, L.T., et al.,2022. Fracture Logging Identification and Effectiveness Analysis of Shale Reservoir of the Permian Fengcheng Formation in Mahu Sag, Junggar Basin. Journal of Palaeogeography, 24(3):540-554 (in Chinese with English abstract).

[15] Ivanova, V.M., Sousa, R., Murrihy, B., et al.,2014.Mathematical Algorithm Development and Parametric Studies with the Geofrac Three-Dimensional Stochastic Model of Natural Rock Fracture Systems. Computers & Geosciences,67:100-109.

[16] Jia,C.Z.,1999. Structural Characteristics and Oil and Gas Accumulation Rules in Tarim Basin.Xinjiang Petroleum Geology,20(3):177-183 (in Chinese with English abstract).

[17] Jia, C.Z., Zou, C.N., Li, J.Z., et a.,2012. Assessment Criteria, Main Types, Basic Features and Resource Prospects of the Tight Oil in China. Acta Petrolei Sinica, 33(3):343-350 (in Chinese with English abstract).

[18] Karim, F.M., Durlofsky, L.J., Aziz, K.,2004. An Efficient Discrete Fracture Model Applicable for General-Purpose Reservoir Simulators. Society of Petroleum Engineers Journal,9(2):227-236.

[19] Kemeny, J., Post, R.,2003.Estimating Three-Dimensional Rock Discontinuity Orientation from Digital Images of Fracture traces. Computers & Geosciences,29(1):65-77.

[20] Lato, M.J., Diederichs, M.S., Hutchinson, D.J., 2010.Bias Correction for View-Limited Lidar Scanning of Rock Outcrops for Structural Characterization. Rock Mechanics and Rock Engineering,43(5):615-628.

[21] Lehmer, E.,1938. On Congruences Involving Bernoulli Numbers and the Quotients of Fermat and Wilson. Ann. Math,39:350-360.

[22] Li, Y., Hou, G.T., Hari, K.R., et al.,2018. The Model of Fracture Development in the Faulted Folds: The Role of Folding and Faulting. Marine and Petroleum Geology,89:243-251.

[23] Li, Y., Qi, J.F., Shi, J., et al.,2017. Characteristics of Mesozoic Basin in Kuqa Depression, Tarim Basin. Geotectonica et Metallogenia,41(5):829-842 (in Chinese with English abstract).

[24] Lin, L.B., Yu, Y., Nan, H.L., et al.,2021.ReservoirTightening Process and Its Coupling Relationship with Hydrocarbon Accumulation in the Fourth Member of Upper Triassic Xujiahe Formation in the Western Sichuan Depression, Sichuan Basin. Oil & Geology,42(4):816-828 (in Chinese with English abstract).

[25] Lu, H.F., Jia, D., Chen, C.M., et al.,1999. Nature and Timing of the Kuqa Cenozoic Structures. Earth Science Frontiers,6(4):215-221 (in Chinese with English abstract).

[26] Lv, W.Y., Zeng, L.B., Liu, J., et al.,2016.Fracture Research Progress in Low Permeability Tighe Reservoirs. Geological Science and Technology Information,35(4):74-83 (in Chinese with English abstract).

[27] Marr, D., Hildreth, E.,1980.Theory of Edge Detection. Proceedings of the Royal Society of London. London, UK.

[28] Oda, M., Hatsuyama, Y., Ohnishi, Y.,1987.Numerical Experiments on Permeability Tensor and Its Application ti Jointed Granite at Stripa Mine. Journal of Geophysical Research,92:8037-8048.

[29] Ren, Y.P.,2019. Terrestrial Lidar Based Fracturing and Folding Mechanism Analyzing and Fracture Network Modeling(Dissertation). Zhejiang University, Hangzhou(in Chinese with English abstract).

[30] Shi, C.Q., Li, Y., Yuan, W.F., et al.,2021. Characteristics on Reservoir Architecture and Quality of Tight Sandstone Reservoirs: Taking Jurassic Ahe Formation in Dibei Area of Kuqa Foreland Basin As an Example. Journal of China University of Mining & Technology,50(5):877-892 (in Chinese with English abstract).

[31] Su, X.,2020.Flow Behavior in Granite Fractures and Mechanism of Hydraulic Fracture Propagation(Dissertation). Chongqing University, Chongqing(in Chinese with English abstract).

[32] Tang, Y., Mei, L.F., Tang, W.J., et al.,2010. Property Analysis and Stochastic Modelling of Fractured Reservoirs. Journal of Southwest Petroleum University(Science & Technology Edition),32(4): 56-66 (in Chinese with English abstract).

[33] Tian, Y.P., Liu, X., Li, X., et al.,2011. Finite Element Method of 3-D Numerical Simulation on Tectonic Stress Field. Earth Science Frontiers,36(2):375-380 (in Chinese with English abstract).

[34] Wang, K., Zhang, R.H., Zeng, Q.R., et al.,2022. Characteristics and Formation Mechanism of Lower Cretaceous Deep and Ultra-Deep Reservoir in Bozi-Dabei Area, Kuqa Depression. Journal of China University of Mining & Technology,51(2):257-274 (in Chinese with English abstract).

[35] Wang, Z., Xue, H.Y., Sun, D.,2020. Automatic Stitching Method and Precision Analysis of 3D Point Cloud Based on Visual Tracking Technology. Bulletin of Surveying and Mapping,2020(S1): 164-167(in Chinese with English abstract).

[36] Xu, J.S., Wang, J.X., Chen, X.Q., et al.,2022. Effects of Poisson Ratio on In-Situ Stress Field Near the Jiali Fault Along the Sichuan-Tibet Railway. Earth Science Frontiers,47(3):818-830 (in Chinese with English abstract).

[37] Yang, H.J., Zhang, R.H., Yang, X.Z., et al.,2018. Characteristics and Reservoir Improvement Effect of Structural Fracture in Ultra-Deep Tight Sandstone Reservoir: A Case Study of Keshen Gasfield, Kuqa Depression, Tarim Basin. Natural Gas Geoscience,29(7):942-950 (in Chinese with English abstract).

[38] Zeng, Q.L., Zhang, R.H., Lu, W.Z., et al.,2017. Fracture Development Characteristics and Controlling Factors Based on 3D Laser Scanning Technology: An Outcrop Case Study of Suohan Village, Kuqa Foreland Area, Tarim Basin. Natural Gas Geoscience, 28(3):397-409 (in Chinese with English abstract).

[39] Zhu, X.M., Pan, R., Zhu, S.,F., et al.,2018. Research Progress and Core Issue in Tight Reservoir Exploration. Earth Science Frontiers,25(2):141-146 (in Chinese with English abstract).

[40] 鲍怡晨,刘强虎,杜晓峰,等,2021.基于砾石-基质-裂缝三元素的砂砾岩岩相划分.地球科学,46(6):2157-2171.

[41] 党发宁,韩文涛,郑娅娜,等,2007.混凝土破裂过程的三维数值模型.计算力学学报,24(6):829-833.

[42] 邓少贵,李智强,2009.裂缝性储层裂缝的阵列侧向测井响应数值模拟.地球科学,34(5):841-847.

[43] 丁文龙,尹帅,王兴华,等,2015.致密砂岩气储层裂缝评价方法与表征.地学前缘, 22(4):173-187.

[44] 董少群,曾联波,Xu C.S.,等,2018.储层裂缝随机建模方法演技进展.石油地球物理勘探,53(3):625-642.

[45] 冯驰,张黎明,2019. Delaunay三角剖分法在噪声地图中的应用.测绘与空间地理信息,42(10):80-82.

[46] 巩磊,高铭泽,曾联波,等,2017.影响致密砂岩储层裂缝分布的主控因素分布——以库车前陆盆地侏罗系-新近系为例.天然气地球科学,28(2):199-208.

[47] 何登发,周新源,杨海军,等,2009.库车坳陷的地质结构及其对大油气田的控制作用.大地构造与成矿学,33(1):19-32.

[48] 黄玉越,王贵文,宋连腾,等,2022.准噶尔盆地玛湖凹陷二叠系风成组页岩储集层裂缝测井识别与有效性分析.古地理学报,24(3):540-554.

[49] 贾承造,1999.塔里木盆地构造特征与油气聚集规律.新疆石油地质, 20(3):177-183.

[50] 贾承造,邹才能,李建忠,等,2012.中国致密油评价标准、主要类型、基本特征及资源前景.石油学报,33(3):343-350.

[51] 李勇,漆家褔,师俊,等,2017.塔里木盆地库车坳陷中生代盆地性状及成因分析.大地构造与成矿学, 41(5):829-842.

[52] 林良彪,余瑜,南红丽,等,2021.四川盆地川西坳陷上三叠统须家河组四段储层致密化过程及其与油气成藏的耦合关系.石油与天然气地质,42(4):816-828.

[53] 卢华复,贾东,陈楚铭,等,2017.库车新生代构造性质和变形时间.地学前缘,6(4):215-221.

[54] 吕文雅,曾联波,刘静,等,2016.致密低渗透储层裂缝研究进展.地质科技情报.35(4):74-83.

[55] 任宇鹏,2019.基于激光雷达的褶皱裂缝力学系统及裂缝分布模型构建（博士毕业论文）.杭州：浙江大学.

[56] 史超群,李勇,袁文芳,等,2021.致密砂岩储层构型特征及评价：以库车前陆盆地迪北地区侏罗系阿合组为例.中国矿业大学学报, 50(5):877-892.

[57] 苏小鹏,2020.干热岩裂缝渗流及储层水力压裂缝网扩展机理研究（博士毕业论文）.重庆：重庆大学.

[58] 唐永,梅廉夫,唐文军,等,2010.裂缝性储层属性分析与随机模拟.西南石油大学学报（自然科学版）, 32(4):56-66.

[59] 田宜平,刘雄,李星,等,2011.构造应力场三维数值模拟的有限单元法.地球科学,36(2):375-380.

[60] 王珂,张荣虎,曾庆鲁,等,2022.库车坳陷博孜-大北地区下白垩统深层-超深层储层特征及成因机制.中国矿业大学学报, 51(2):257-274.

[61] 王智,薛慧艳,孙迪,等,2020.基于视觉追踪技术的三维点云自动拼接方法及精度分析.测绘通报,(S1):164-167.

[62] 许俊闪,王建新,陈兴强,等,2022.泊松比对川藏铁路嘉黎断裂附近地应力场的影响.地球科学,47(3):818-830.

[63] 杨海军,张荣虎,杨宪彰,等,2018.超深层致密砂岩构造裂缝特征及其对储层的改造作用：以塔里木盆地库车坳陷克深气田白垩系为例.天然气地球科学, 29(7):942-950.

[64] 曾庆鲁,张荣虎,卢文忠,等,2017.基于三维激光扫描技术的裂缝发育规律和控制因素研究——以塔里木盆地库车前陆区索罕村露头剖面为例.天然气地球科学.28(3):397-409.

[65] 朱筱敏,潘荣,朱世发,等,2018.致密储层研究进展和热点问题分析.地学前缘,25(2):141-146.

基金资助

国家十三五重大专项(2017ZX05008-001;2017ZX05008-005)

湖北省教育厅科学研究计划项目(T201905)