液态CO₂瞬变致裂技术是一种新型的无炸药爆破方式^[1-2]，通过可重复利用的致裂管充装液态CO₂，加热使液态CO₂瞬间气化，产生强大的膨胀压力，快速释放的高压CO₂气体在岩体内产生拉应力场，从而破碎岩体.该方法在能量利用率上效率更高，且安全性高、震动小、粉尘低，大大减小了对岩体、周边建筑物及环境的损伤破坏^[3].目前，关于液态CO₂瞬变致裂的理论研究和实践多集中在瓦斯型煤炭开采中，对于其在岩石开挖中的应用研究较少，且对其致裂效果评价研究更为稀少，限制了CO₂致裂开挖效益的充分发挥.

在液态CO₂瞬变致裂的科研和生产工作中，准确而又迅速地获取致裂岩石的块度大小及分布，是调整致裂参数和施工工艺，实现优化生产的重要参考依据.测定岩石块度的方法可追溯到19世纪末，澳大利亚学者Rozival首次用线段法测定出块度组成.目前，块度特征信息的获取多采用大块计数法、筛分称重法、爆堆直接测量法、平面照相法等方式，基本满足了工程实际的需要^[4].但同时也存在着一些比较突出的问题，大块计数法忽略小块石，而筛分称重法则只适用于较小块石的测量，无法综合反映整个致裂爆堆的块度分布，且工作量巨大，劳动强度高^[5]；爆堆直接测量法是在爆堆上设置若干条直线，测量直线上各个岩块尺寸占直线的长度比，其主观影响较大^{[4, 6]}.平面照相法只拍摄典型区域，依次推断整个爆堆的块石特征，资料不全，致使工程技术人员始终处于“管中窥豹”的状态^[6-10].因此，准确、快速地获取致裂爆堆的块石特征信息，进而客观、高效地评价致裂效果，对液态CO₂瞬变致裂的技术水平的进一步提升是非常有益的.近年来，随着中国信息技术的迅猛发展，无人机(unmanned aerial vehicle, UAV)技术取得了长足进步^[11-13].这一技术凭借着飞行灵活、操作简便、成本低廉等优点广泛应用于社会各个领域^[13-14].其中，结合图像处理技术可实现对液态CO₂瞬变致裂爆堆信息的快速获取，信息全面，分辨率高，为相关工作提供了崭新的方法和思路.

在获取了致裂爆堆的块度特征后，如何定量表征岩石的块度分布规律是评价液态CO₂瞬变致裂效果的重要步骤.目前，常用的两类爆破块度描述方法为单一指标描述和分布函数描述.单一指标描述最常用的为“不合格大块率”指标和“平均块度”指标^[15].前者反映了对后续工序(装、云、破碎)有重要影响的那部分(通常必须进行二次破碎)的大块岩石所占比例，而后者则反映了平均破碎程度，该方法主观性强，且工作量巨大，单一指标无法全面表征块度特征，准确性较差.分布函数可描述块度的分布规律，但并不能普遍适用.因此，探寻一种定量表征液态CO₂瞬变致裂后岩石块度分布的方法可为后续致裂影响因素的分析以及参数优化等研究奠定基础.为此，本文将着重探讨如何利用无人机+图像处理技术来快速获取施工现场岩体致裂块度的过程与块度表征方法，进而达到为液态CO₂瞬变致裂效果快速评价提供有效技术支撑之目的.

1 试验

本文依托北京市房山区某安置房基坑开挖工程，采用液态CO₂瞬变致裂技术进行现场试验，试验地区发育石灰岩，岩质坚硬，节理裂隙发育，岩石较完整，风化程度为未风化到微风化.试验采用不同型号(95型、108型)、不同数目的致裂管^[16-17](致裂管如图 1所示，设备参数及布孔参数见表 1、2)进行致裂试验.致裂完成后，采用“无人机摄像+MIPS图像处理”的方式来实现致裂岩体特征(方量、块粒度等)的快速获取.通过数据分析，描述块度分布特征，建立块度分布函数，基于岩石破碎程度，评价液态CO₂瞬变致裂效果.

整个现场试验工作流程为：试孔(爆孔)钻探→致裂气体充装→致裂器安装(→监测系统布设→无人机遥感)→起爆网络连接→启爆(→无人机致裂效果获取)→爆后检查与效果评价→提管.试验过程中应特别注意钻孔后的填塞作业，填塞材料一般采用钻屑、干燥细石粉，充分填实，防止“飞管”.

2 基于图像处理的致裂岩石特征信息获取

致裂岩石块度的研究是评价液态CO₂瞬变致裂效果的关键因素，岩石的破碎程度直接影响着后续工序，精准的获得岩石块度的信息对优化致裂参数、提高生产效率、致裂质量管理及控制成本具有重要意义.常用的块石特征获取方法主要有筛分法和平面照相法.筛分法是目前采用的最精确方法，但筛分法对作业干扰较大，劳动力投入多，即使是处理小型的爆堆也需要高额的费用，同时，筛分法针对小块石具有高准确率，但对于尺寸较大的岩石无法进行试验，准确性低.因此，筛分法通常仅用于校正数据或稀有矿物致裂开采，在其他工程中应用较少.平面照相法具有操作迅速、费用低廉等优势，但只能拍摄典型区域，依次类推整个爆堆得块石特征，资料不全，若为了获得整体块度分布，则需处理大量照片，无疑降低了效率，同时其自身的技术误差也影响了评定结果的准确性.本文采用“无人机摄像+图像处理”的方式进行块度信息的获取.小型无人机操作灵活、视觉可控，结合自行研发的图像处理软件能够更准确更全面的获取致裂块度信息.

2.1 无人机摄像

无人机摄像系统主要由飞行平台、传感器稳定云台、飞行控制系统和相机等几部分构成.其中，无人机搭载的相机一般为固定镜头的非量测型的普通数码相机，为焦距值固定的定焦镜头^[11].在进行拍摄前，首先需根据致裂范围、精度要求等合理规划无人机拍摄路线，由此保证影响资料的可靠性.航线规划主要涉及航向重叠率、旁向重叠率、航高、航速、影响分辨率、测区范围等参数的设置，如图 2所示^[13-14].

拍摄过程中，无人机从拍摄区域的一角起飞，到达合适飞行高度后，沿着既定航带飞行，相邻两幅照片的拍摄范围需保持一定的航向重叠.随后，依次完成不同航带的飞行，两航带之间照片拍摄范围保持一定的旁向重叠.其中保证一定的航向重叠率与旁向重叠率尤为重要，由此避免无人机飞行平台姿态不稳时出现的拍摄遗漏，保证照片的顺利拼接，且降低照片畸变对获取数据的准确性的影响，无人机的路径示意图如图 3所示^[14].

本次实验使用的小型无人机型号为“大疆精灵Phantom 4 Pro”，技术性能完全可以满足本文致裂实验要求，其主要技术参数如下：1英寸2 000万像素摄像传感器、可拍摄4 K/60 fps视频，并以14张/s的速度拍摄静态照片，5向环境感知，机械快门，双编码方式，30 min续航能力，7 km飞行距离.由于无人机自带摄像机与专业航空摄影测量摄像机有一定差别，故拍摄过程中航向和旁向重叠率分别不低于80%和60%为宜，其中无人机高度3 m，拍摄时采用自动拍照模式，飞行速度0.3 m/s，自动拍照间隔为0.3 m/s.

2.2 图像处理技术

数字图像处理技术是针对图像实现高效的复原、校正、增强、统计分析、分类和识别等的加工技术.在许多情形中，图像通过肉眼观察存在模糊以及不可见的特征，而依靠适合的图像预处理，则能够将原本模糊的图像演变得更加清晰^{[6, 18-22]}.同时在数字图像处理的基础上，能够针对肉眼无法观察的信息进行分类统计.依靠计算机识别能够实现高效的检索、匹配以及识别操作^[23].这些优势非常适合运用在液态CO₂瞬变致裂块石特征、结构要素的分析、统计上.本次试验采用自主研发的数字图像处理软件系统(MIPS)进行图像处理，其基本技术路线如图 4所示^[24].

2.2.1 比例尺标定

拍摄前，事先放置了一个地面比例标志(50 cm*30 cm*30 cm)，因此在计算机标定过程中，按地面标志进行换算，即可获取图像比例.

2.2.2 图像处理

首先进行图像预处理，主要实现方式为图像滤波.MIPS软件主要通过均值滤波来完成图像滤波，均值滤波是一种线性空间滤波，算法简单，计算速度快，本次实验采用的滤波窗口为3*3.图像处理为图像的二值化.用0~255之间的256个灰度级别来表示预处理后图像各个像素的明暗程度，采用为OSTU算法进行阈值选取，即按图像的灰度特性将图像分为两类，当类间方差与类内方差的分离度最大时为最佳阈值.得到二值化图像后，针对干扰因素较大的部分进行人工处理，图像切割及剔除，提高准确性.处理后的图像如图 5所示(以108多孔试验为例).

2.2.3 信息获取

对处理后的数字图像进行分析统计，通过比例尺，由像素点信息转化成实际信息，从而获得致裂块石的特征信息.

2.3 方法对比

常规的块石特征获取方法主要为筛分法和平面照相法，由于现场实验中人员、时间的限制，未能进行现场筛分称重，无法与无人机拍摄的照片进行对比分析.本文选取一个小区域采用常规的平面照相法进行块石特征获取(以108多孔为例)，从而进行对比，结果如图 6所示，二者大致相同，但平面摄像结果中小块石比例较多，“无人机+MIPS”技术的结果更符合实际.因此，基于“无人机摄像+MIPS图像处理”技术的块石特征快速获取方法对于工程实践具有重要意义.

3 块度分析结果与模型建立 3.1 块度分布特征

由于利用上述图像处理技术仅仅得到的是致裂岩块投影面积的信息数据，当各粒级内岩块数目足够多时，以岩石碎块的平面投影尺寸(如面积或周长)作为检测致裂块度的几何特征量，可反映实际的致裂效果.岩块的平面尺寸与体积之间具有一定的数量关系，目前较为常用的为公式^[25]为:

$ \begin{array}{l} \;\;a = \frac{{\frac{P}{{\rm{ \mathsf{ π} }}} + \sqrt {\frac{{{P^2}}}{{{{\rm{ \mathsf{ π} }}^{\rm{2}}}}} - \frac{{4S}}{{\rm{ \mathsf{ π} }}}} }}{2}, \\ \;\;b = \frac{{\frac{P}{{\rm{ \mathsf{ π} }}} - \sqrt {\frac{{{P^2}}}{{{{\rm{ \mathsf{ π} }}^{\rm{2}}}}} - \frac{{4S}}{{\rm{ \mathsf{ π} }}}} }}{2}, \\ D = 1.16b\sqrt {1.35\frac{b}{a}} , \\ \;\;\;\;\;\;\;V = D \times S. \end{array} $

(1)

式中：S为岩块面积；P为岩块周长；D为岩块的筛网尺寸；a为最佳椭圆最大半径；b为最佳椭圆的最小半径；V为岩块体积.

采用式(1)进行体积换算.首先，对换算后的体积进行简单分析，图 7表示液态CO₂瞬变致裂不同方案岩石块石总体积对比图.主要由于致裂管中CO₂储量不同，所释放能量变化较大，故多孔致裂的致裂块度总体积明显大于单孔致裂体积，且108型致裂试验致裂方量均大于95型致裂方量.图 8为4次试验结束后其块度组成，对比两次多孔试验，108型多孔试验块度分布较为均匀，大块石含量较95型较小，同时前者致裂块石总体积也多于后者.综上所示，108型多孔致裂试验致裂效果优于95型多孔致裂试验.两次单孔试验结果对比中，95型单孔大块石含量较小，致裂效果较好，原因考虑为108型单孔致裂试验的试验点岩石风化较为严重，岩石较为松散，对力的传导具有削减作用，故致裂效果变差.

致裂块体块径级配曲线，如图 9所示，该曲线呈下凹型，大块石体积分数较高，其中曲线中出现多个水平段，表示缺少该块径区间的块石，即块石粗的较粗，细的较细，级配不良.

为了更好地表示致裂块体的均匀程度和级配的优劣，选取d₁₀(小于此块径的块石含量占总含量的10%，也称有效块径)、d₃₀(小于此种粒径的块石含量占总含量的30%)、d₆₀(小于此粒径的块石含量占总含量的60%，也称控制块径)3种特征块径进行表示，见表 3.

采用d₅₀表示致裂块石平均块度，d₅₀越大，致裂块石越大.不均匀系数C_u越大，致裂块石越不均匀，即大块石和小块石尺寸相差越悬殊.曲率系数C_c表征块径级配累积级配曲线的斜率是否连续.同时，当块堆中所缺少的中间块径大于连续级配曲线的d₃₀时，曲率系数变小，小于d₃₀时，曲率系数变大.综上所述，结合致裂块石级配曲线，多孔致裂试验中，108型不均匀系数C_u较小，且C_c较大(水平段出现在大于d₃₀处)，故其级配较好，综合平均块度及块石体积，其致裂效果最好.单孔致裂试验中，95型单孔试验块石级配不良，但其平均块度较小，主要是由于总体积与致裂点岩石性质不同导致.其具体致裂效果比较需增加试验进行进一步验证.

3.2 块度分布分形模型

致裂岩块在几何形状和块度分布两方面都具有统计自相似性，从分形角度描述岩体破碎过程即是：岩体最初破碎成有限个形状相似的块体，其中的一部分在荷载作用下各自破碎成与原块体形状相似的次一级块体，这些次一级块体中的一部分又被破碎成形状相似的更次一级块体，此过程的每一次重复都产生更小的块体^[26-28].由此，本文从分形的角度分析液态CO₂瞬变致裂破碎岩石的分布.引入分形维数D这一反映岩石块度分布的特征值，建立致裂块度分形分布模型.

令岩样为零级碎块，各级碎块均以破坏概率P破碎为n个子碎块，且相邻级碎块尺度的相似比为1/r，则对i级碎块总数为^[26]

$ {N_i} = (1 - P)n\left[ {1 + nP + {{(nP)}^2} + \cdots + (nP){i^{ - 1}}} \right], $

由分形的基本定义可知：

$ \frac{{{N_{i + 1}}}}{{{N_i}}} = {\left( {\frac{1}{r}} \right)^D}. $

因为i> >1，且假定nP>1，则可得到：

$ D = \log (np)/\log \left( {\frac{1}{r}} \right). $

设致裂过程中岩体的破碎模型总体积不变，每次构成时，将原有块体以破坏概率P、相似比r分裂成下一级块体，重复进行直到无穷，如图 10所示，生成一系列大小不同的具有相似形状的碎块.

该过程中次一级岩块的量为^[29]

$ N = np = {\left( {\frac{1}{r}} \right)^3}P, $

由此，

$ D = \frac{{\log N}}{{\log (1/r)}} = 3 - \frac{{\log P}}{{\log r}}. $

(2)

设岩体最初由M个线尺寸为x的最大块体组成，各块体的体积V=C_vx³(C_v为体积形状系数)，则所有源块体经过k次分形构造后产生的第k次块体的线尺寸为^[27]

$ {x_k} = {r^k}x(k = 0, 1, 2, \cdots ), $

(3)

其个数为

$ {N_k} = {\left[ {{{\left( {\frac{1}{r}} \right)}^3}P} \right]^k}(1 - P)M, $

(4)

由式(3)、(4)得第k级块体的体积v_k为

$ {v_k} = {C_v}x_k^3{N_k} = {C_v}{x^3}{P^k}(1 - P)M. $

线尺寸小于或等于x_i的块体总体积V_i为

$ {V_i} = \mathop {\lim }\limits_{j \to \infty } \sum\limits_{k = i}^j {{v_i}} = {C_v}{x^3}{P^i}M. $

所以，线尺寸小于或等于x_i的块体体积占总体积的比例为

$ {y_i} = \frac{{{V_i}}}{{{V_k}}} = {P^i}. $

(5)

由式(3)可得i=log(x_i/x)/log r，结合式(2)、(5)，则有

$ {y_i} = {\left( {\frac{{{x_i}}}{x}} \right)^{3 - D}}. $

(6)

由此，即为分形维数与致裂块度各个级别尺寸筛下累积量之间的关系，也是致裂岩石块度分布的分形模型，其中x为岩块的最大尺寸.式(6)是从分形构造推导中获得的块度分布，该式与爆破中常用的G-G-S经验分布函数在形式上完全相同.分形模型的优点是不用考虑原岩的赋存状态、地质构造条件、致裂参数等，只要从形成的破碎岩石中求出分形维数D，即可了解其块度分布，且分维的概念更明确，适用性更加广泛.

根据上述理论分析，表 4、图 11即为该4组液态CO₂瞬变致裂试验的筛下体积分数及分形结果值.

由表 4结果可看出，尽管4次试验参数不尽相同，但致裂后的岩石块度仍具有良好的分形结构.分别对比两次单孔试验和两次多孔试验，不难发现，致裂后小块度占比越高，分形维数越大；反之，越小.其中108单孔实验分形维数为0.183，而95单孔则为0.970，存在异常，考虑原因为95单孔实验炮孔位置黏土含量较高，为土石堆积层，强风化层，岩石性质极差，导致分形维数增大，后续增加现场实验进一步验证上述结论.分形维数D是评价岩石液态CO₂瞬变致裂块度分布的定量指标.

以95型单孔致裂试验为例，绘制分形模型分布曲线与筛分后实际块度分布曲线，如图 12所示.

块度分形模型误差计算公式为

$ e = \sqrt {\frac{{\sum\limits_{i = 1}^n {{{\left[ {\left( {{y_{i实际}}\;\; - {y_{(i - 1)实际}}\;\;} \right) - \left( {{y_{i分形}\;\;} - {y_{(i - 1)分形}}\;\;} \right)} \right]}^2}} }}{n}} , $

式中y_i实际、y_i分形分别为实际块度与分形模型中，小于第i级岩块尺寸的体积累积体积分数，其余同理.利用上述公式4次试验误差计算结果见表 5.

根据误差分析，分形模型与实际块度分布平均误差在4%，误差较小，分型模型精度较高，能够满足工程实际应用的要求

虽然4次试验岩石块度均具有良好的分形结果，但在不同的致裂参数条件下，致裂后岩石块度的分形维数不同，所以，将分形维数D引入液态CO₂瞬变致裂块度分布函数中的意义在于可将致裂前的函数参数、致裂布孔参数与致裂后的块度分布指数等联系起来，建立岩石块度分维值与致裂参数的关系，为控制液态CO₂致裂块度分布、致裂前块度估算及预报提供了一个新的思路.

4 结论

1) 通过无人机与MIPS图像处理技术，分别对块度的投影面积、周长等参数进行统计分析，实现了对液态CO₂瞬变致裂后块石宏观结构信息的快速获取.

2) 根据分形理论，建立液态CO₂瞬变致裂块度分形模型.运用分形维数刻画块度的分布情况，小块度占比越高，分形维数越大，反之越小.因此，分形维数D是评价岩石液态CO₂瞬变致裂块度分布的定量指标.

3) 本文将分形维数D引入液态CO₂瞬变致裂块度分布评价中，为致裂块度分布的优化、致裂前块度估算及预报提供了一个新的思路.但影响其块度分布的原因众多，包括致裂能量、布孔参数、岩石破碎程度等，后续研究中可采用分形维数评价岩石块度分布，建立分形维数与上述影响因素的定量关系，为后续设计合理致裂方案提供理论依据.