节理裂隙广泛分布于自然界当中，其存在较大程度地削弱了完整岩石的力学强度，因而对一些大型岩石工程诸如水工大坝、高放核废料处置库、深部储藏、岩质高边坡等的安全性与稳定性有较大的影响^[1]. 当节理表面不够粗糙、使节理闭合的法向压力水平较低时，节理可能会在扰动作用下发生剪切错动，从而导致工程失稳破坏. 目前，一般认为表面形貌、作用的法向应力水平以及岩石的物理力学特性是影响耦合节理剪切力学行为的3个主要因素^[2]. 考虑到岩石物理力学特性可以采用室内试验进行测定，而法向应力亦可以通过各种各样的监测设备方便地获取，如何表征节理面形貌成为关键性问题.

Barton等^[3]最早提出采用10条标准剖面线来表征节理的形貌特征，通过比对方法确定粗糙度 (joint roughness coefficient，JRC)，然而该方法依赖使用者的经验和主观判断. 许多学者提出了其他的形貌参数来定量化描述JRC. 例如，葛云峰等^[4]提出光亮面积百分比BAP来描述节理面的粗糙度，并建立了BAP与JRC的经验关系.Liu等^[5]考虑节理剪切过程中重叠部分的统计特征，提出2个考虑采样点距的JRC经验公式.Zhang等^[6]则首先考虑节理表面迎剪侧微段的作用，将坡度方均根Z₂修正为Z₂′，然后考虑节理高度起伏特征提出综合参数λ来确定JRC的上下限，结果表明该方法可以较为精确地确定JRC. Tatone等^[7]首先将Barton的10条标准剖面线进行数字化处理，然后在不同采样点距下建立JRC和 $\theta _{\max }^ *$/(C+1)_2D的经验关系（ $\theta _{\max }^ *$为最大视倾角，C为拟合系数），该方法同样具有较好的精度. 班力壬等^[8]将剖面线类比于梳齿结构，采用粗糙度参数c表征面向剪切方向微单元体的平均坡度，综合考虑c的分形特征，提出采用2个分形参数来描述JRC. 综上所述，这些研究方法都是节理粗糙度评估方法的有效探索，但是由于这些研究方法都是二维的，不能充分反映真实节理面的空间几何特性，具有较明显的局限性.

针对上述问题，一些学者在节理粗糙度三维表征方法上进行了一些积极尝试. 例如，Belem等^[9]直接将二维粗糙度参数Z₂推广至三维情形Z_2S. 基于三维的Z_2S，唐志成等^[10]建立了新的节理剪切强度准则，并发现其精度优于经典的Barton模型. 类比于剖面线迹线长度比R_p，EL-Soudani^[11]提出了三维粗糙度指标R_s. Chen等^[12]提出采用变异函数综合表征节理粗糙度JRC. 然而，这些三维粗糙度指标的问题在于，不能区分粗糙度在正、反方向上的差异（相差180°），因而难以与节理的实际受剪状态联系起来.

一般来说，当节理上、下盘发生剪切错动后，只有迎着剪切方向的微单元体发生接触，进而产生滑移、啃断，而背着剪切方向的微单元体则逐步分离，不能抵抗这种相对运动^[13]. 因此，迎着剪切侧的微单元体对节理的粗糙度表征起着关键性的作用. 考虑到这一点，Grasselli等^[14-15]分析节理表面有效抵抗倾角与潜在接触部分的关系，提出了几个新的参数描述节理的粗糙度特征. 这些参数能够有效地区分粗糙度在正反剪切方向上的差异，且不是基于某几条剖面线的平均结果^[16-18]，因此较为广泛地应用于节理剪切强度准则中. 例如，唐志成等^[19-20]结合Grasselli形貌参数，提出了能够反映峰值剪胀角发展演化规律的剪切强度模型，该模型在形式上满足摩尔库伦准则，预测精度高于Grasselli剪切强度模型^[14-15]；Yang等^[21]从节理倾斜试验结果出发，通过经验拟合的方式得到了初始剪胀角，根据峰值剪胀角随法向应力的变化规律，提出了新的峰值抗剪强度模型. Singh等^[22]统计了11种较为合理的峰值抗剪强度准则，通过试验研究发现，基于Grasselli形貌参数的唐志成模型^[19-20]和Yang模型^[21]在形式上和预测精度上均要优于其他模型. 由此可见，Grasselli模型有比较广泛的应用基础.

此外，节理形貌表征结果与采样点距密切相关. 唐志成等^[18]研究评估指标 $ \theta _{\max }^ * /(C + 1) $与采样点距的关系，发现随着采样点距的增加， $ \theta _{\max }^ * /(C + 1) $不断减小. 葛云峰等^[23]研究大规模结构面的采样点距效应，得到了类似的结论. 目前，关于小尺度、低采样点距的节理面采样效应的研究还较少.

本研究对天然红砂岩粗糙节理开展了形貌扫描试验，在Grasselli模型^[14-15]基础上，重点研究节理表面有效倾角与剪切方向上潜在接触部分的关系，提出新的三维粗糙度指标. 将改进模型与其他模型进行详细的比较，分析其合理性，并讨论新指标的各向异性特征. 最后，结合改进的形貌描述方法，分析形貌的采样点距效应.

采集3组形貌各异的天然红砂岩节理用于形貌扫描. 节理尺寸均为100 mm × 100 mm × 50 mm（长×宽×高），记为J1、J2、J3. 扫描试验在武汉大学岩土与结构工程安全湖北省重点实验室的三维工业扫描仪X3上进行（见图1(a)、(b)）. 该扫描系统精度约为0.025 mm，单次扫描区域的面积约为280 mm×250 mm. 该仪器发射的激光对人眼安全，测量物件大小范围约为0.1~6.0 m，使用时手持，携带较方便.

系统与软件Geomagic Studio兼容（见图1(b)、(c)），可以进行数据的初步处理. 在扫描时，扫描头距离待测物体为30~40 cm，通过扫描背景布中的标记棱镜识别物体位置. 在每个扫描视野中，一般应保证不少于3个标记点以确保扫描结果的有效性. 在扫描时点距设定为0.2 mm. 对于每个节理，选择面积大小为90 mm×90 mm的区域用于形貌分析（见图1(d)、(e)）. 3个节理面的形貌如图2所示.

根据Tatone等^[24]提出的方法将获得的点云变为系列微单元体，如图3所示. 图中，S为剪切方向，v为微单元体的外法线，v₀为剪切平面的外法线，v₁为v在剪切平面上的投影方向，β为S和v₁间的夹角，θ、θ^*分别为微单元体的倾角和有效倾角. 由于有效倾角θ^*与剪切方向密切相关，θ^*的统计特征可以较好地解释节理粗糙度的方向相关性. 对于单个微单元体，θ可以通过立体几何分析得到，而β与选定研究的剪切方向有关，表达式如下：

(1) $ \cos\;\beta = \frac{{{\boldsymbol{S}} \cdot {{\boldsymbol{v}}_1}}}{{\left| {\boldsymbol{S}} \right| \times \left| {{{\boldsymbol{v}}_1}} \right|}} \text{，} $

(2) $ \cos \;\theta = \frac{{{\boldsymbol{v}} \cdot {{\boldsymbol{v}}_0}}}{{\left| {\boldsymbol{v}} \right| \times \left| {{{\boldsymbol{v}}_0}} \right|}} . $

通过θ和β可以得到单个微单元体θ^*的计算表达式：

(3) $ \tan \;{\theta ^ * } = - \tan \;\theta \cos \;\beta . $

节理表面数理构建主要过程如下. 1) 将三维激光扫描获取的点云进行初步处理. 利用逆向工程软件Geomagic Studio删除点云中的杂质，并对节理面进行裁剪以选择相应的分析区域. 在删除杂质点时，文件以.asc格式进行处理. 在裁剪分析区域时，须利用Geomagic Studio自带的划分网格功能将文件以.stl格式进行处理. 此外，在利用逆向工程软件处理表面时，往往会对孔洞进行填充. 在本研究中须避免这样的操作，因为这样做会导致粗糙表面形貌失真. 在处理后须将结果以.asc格式进行存储，以便进行下一步操作. 2) 将上述点云以矩阵形式导入Matlab作为初始数据，然后采用Delaunay算法将点变成面. Delaunay剖分算法实际上建立了点云和面的关系，这种关系信息存储在一个检索矩阵中. 利用检索矩阵，可以找到每个三角形. 在本研究中，点云是三维的，因此须将点云投射到平面上，以便进行二维Delaunay网格划分. 随后，建立二维Delaunay网格和三维网格的对应关系. 3) 三角形单元几何信息的计算和剪切过程中潜在接触部分的检索. 例如，利用海伦-秦九韶公式计算每个三角形单元的面积，根据立体几何关系确定每个单元的真倾角和有效倾角. 表面所有三角形单元是利用Delaunay剖分得到的检索矩阵进行遍历的. 最后，给定不同的门槛值θ^*，然后寻找所有三角形单元中有效倾角大于θ^*的单元，即为潜在接触部分.

根据上述方法遍历节理表面可以获取每个微元体的几何信息. 在给定剪切方向后，背离剪切侧的微单元体（θ^*<0°）并没有参与到抵抗节理剪切错动的过程中，因此在确定节理粗糙度指标时应不予考虑. 在Grasselli^[14-15]方法中，先设定一个门槛值θ^*，那么可以检索到所有有效倾角大于门槛值θ^*的微单元体，这些微单元体即是在剪切过程中的潜在接触部分. 不同的门槛值对应着不同的潜在接触部分.

对于如图4(a)所示的剪切方向，当门槛值分别为0°、5°、10°、15°、20°时，节理J1对应的潜在接触部分（黑色部分）分别如图4(b)~(f)所示. 门槛值越大，潜在接触部分越少.

根据Grasselli^[14-15]的研究，节理的潜在接触部分用潜在接触面积比A_θ*来描述. 如图5所示为不同有效倾角门槛值θ^*对应的A_θ*. 可以看到，A_θ*为有效倾角门槛值θ^*的函数，Grasselli^[14-15]指出两者关系可以表示为

(4) $ {A_{{\theta ^ * }}} = {A_0} {\left( {\frac{{\theta _{\max }^ * - {\theta ^ * }}}{{\theta _{\max }^ * }}} \right)^C} . $

式中：A₀为最大可能接触面积比； $\theta _{\max }^ * $为最大有效倾角；C为粗糙度参数，通过拟合得到. Tatone等^[24]分析发现，采用非线性拟合的效果较好，建议将式（4）的积分视为粗糙度指标：

(5) $ \int^{\theta _{\max }^ * } _0 {A_0}{\left( {\frac{{\theta _{\max }^ * - {\theta ^ * }}}{{\theta _{\max }^ * }}} \right)^C}{\text{d}}{\theta ^ * }{\text{ = }}\frac{{{A_0}\theta _{\max }^ * }}{{C{\text{ + }}1}} . $

Tatone等^[24]还发现A₀近似为常数，因此进一步考虑将 $\theta^*_{\max} $/(C+1)作为评估参数.

然而上述方法^{[14-15, 24]}尚存一些缺点. 一些学者针对存在的问题进行了一些改进工作^[2,25-28].

1） $\theta^*_{{\rm{max}}}$不适合用于描述形貌^[2]. 这是因为即使 $\theta^*_{{\rm{max}}}$产生计算误差，参数C的存在使得 $\theta^*_{{\rm{max}}}$/(C+1)几乎不受影响^[2]. 此外，班力壬等^[25-27]也强调，在计算过程中若不剔除较大的离散倾角， $\theta^*_{\max} $的确定也会受到影响. 鉴于此，Tian等^[2]将式（4）中的 $\theta^*_{\max} $用90°来代替，并提出了新的关系：

(6) $ {A_{{\theta ^ * }}} = {A_0} {\left( {\frac{{{{90}^ \circ } - {\theta ^ * }}}{{{{90}^ \circ }}}} \right)^{C'}}. $

式中：C′为参数C的修正. 然而，将C′代入节理剪切强度准则中，剪胀角量纲往往不符合要求.

2）班力壬等^[8]注意到，应当以有效倾角门槛值等于0°时的潜在接触部分面积之和，而非节理表面所有微元面积总和^[28]作为面向剪切侧的节理单元总量. Liu等^[28]基于这一考虑建立一个精度更高的公式如下：

(7) $ A_{{\theta ^ * }}^ + = {{{A_{{\theta ^ * }}}}}/{{{A_0}}} = \exp\; \left( { - {{\left( {{{{\theta ^ * }}}/{{{{\bar \theta }^ * }}}} \right)}^n}} \right). $

式中： ${\bar \theta }^ * $、n均为待定参数，Liu等^[28]将 $\bar \theta^* $/n视作形貌评估参数； $A_{{\theta ^ * }}^ +$为初始接触面积比，当θ^*=0°时， $A_{{\theta ^ * }}^ + =1$，公式仅考虑面向剪切方向的微元的作用. 根据Liu等^[28]的这一做法，式（4）须改写为

(8) $ A_{{\theta ^ * }}^ + = \frac{{A{}_{{\theta ^ * }}}}{{{A_0}}} = {\left( {\frac{{\theta _{{\text{max}}}^ * - {\theta ^ * }}}{{\theta _{\max }^ * }}} \right)^C}. $

式（6）须改写为

(9) $ A_{{\theta ^ * }}^ + = \frac{{{A_{{\theta ^ * }}}}}{{{A_0}}} = {\left( {\frac{{{{90}^ \circ } - {\theta ^ * }}}{{{{90}^ \circ }}}} \right)^{C'}}. $

由以上分析可以看出，改进工作关键主要包括2点：1）改进模型不应考虑 $\theta^*_{{\rm{max}}}$，且应具有较好的精度；2）改进后的模型应采用重新定义的初始接触面积比描述微单元体的接触特性；评估指标的量纲应为角度.

鉴于此，本研究提出了改进模型来描述初始接触面积比和有效倾角门槛值的关系：

(10) $ A_{{\theta ^ * }}^ + = \exp\;\left( {1 - \exp \; {{{\left( {{{{\theta ^ * }}}/{{\theta _{\rm{r}}^ * }}} \right)}^k}} } \right) . $

式中： $\theta^*_{\rm{r}}$、k均为拟合参数. 类似于Liu等^[28]将参数 $\bar {\theta^*} $、n组合为 $\bar {\theta^*} $/n的做法，本研究将 $\theta^*_{\rm{r}}$/k视为新的粗糙度指标. 这是因为根据式（10）， $\theta^*_{\rm{r}}$越大、k越小，则 $ A_{{\theta ^ * }}^ + $越大，式（10）下方面积就越大. 根据文献[24]，曲线下方面积越大意味着更多较陡倾角被考虑在内，则节理越粗糙. $\theta^*_{\rm{r}}$/k越大，节理越粗糙，两者正相关，这种组合可以视为粗糙度表征. 此外， $\theta^*_{\rm{r}}$的量纲为角度，而k无量纲，因此 $\theta^*_{\rm{r}}$ /k满足量纲分析的要求.

结合试验获取的形貌数据，自行编制了相关的Matlab程序来分析节理表面的几何特征，计算过程如图6所示. 对于每个节理，先规定某一剪切方向为正向（例如规定图7中向左为正方向，则图8中向右为反方向），然后分别分析其正向和反向剪切的粗糙度特征，结果如图7、8所示.可以看出，相较于式（7）~（9），式（10）的精度总是最佳（相关系数R²总大于0.999）. 当有效倾角门槛值约为0°~10°时，4个公式与试验结果基本吻合，差异并不明显；当有效倾角门槛值约为10°~25°时，改进模型与试验结果较吻合，而式（9）明显比试验结果偏大，式（7）、（8）则有不同程度的偏离. 式（9）虽然没有考虑 $\theta^*_{{\rm{max}}}$但是与实测结果相比有较大的差距，表明Tian等^[2]的改进仍存在一些不足. 此外，式（7）有时候比式（8）精度高（见图7(b)、8(b)），有时候比式（8）精度还要低（见图7(a)、(c)，图8(a)、(c)），表明这种改进亦不够充分. 相较而言，本研究新模型的精度较好， $\theta^*_{\rm{r}}$/k可能成为一种新的粗糙度表征方式.

上、下盘节理沿不同方向相对运动往往会产生大小不同的阻力，这是由节理形貌的各向异性所导致的. 因此，一个合理的粗糙度指标必须能够反映粗糙度在各个方向上的差异. 对于每个红砂岩节理，首先规定某一方向为0°剪切方向，然后以10°为间隔，逆时针选择36个方位作为剪切方向（见图9），通过编程计算每个分析方向的形貌参数（包括新指标 $\theta^*_{\rm{r}}$/k、Liu等^[28]提出的 $ {\bar \theta ^ * }/n $、Tatone等^[24]提出的 $\theta^*_{{\rm{max}}}$/(C+1)、Tian等^[2]提出的C′），如图10所示. 可以看出，沿不同剪切方向的 $\theta^*_{\rm{r}}$/k构成一个椭圆，根据文献[24]的分析，新的评估指标 $\theta^*_{\rm{r}}$/k具有各向异性特征. 具体来看，新的粗糙度指标具有以下2个特点：1）新粗糙度指标 $\theta^*_{\rm{r}}$/k随着剪切方向的变化连续变化，没有出现突变或者锯齿状^[12]，这种特性和Tian等^[2]、Liu等^[28]、Grasselli^[14-15]的研究结果基本上是一致的，同时也符合经验判断. 2）相较于JRC^[3]、Z₂^[10]等指标，新指标可以有效地体现粗糙度在相反方向上的差异. 因此，新指标 $\theta^*_{\rm{r}}$/k适合用于表征节理形貌特征.

不同采样点距d下节理J1的形态如图11所示. 可以看到，随着采样点距增大，节理面逐渐变“模糊”，主要起伏特征未发生改变（起伏和凹陷较大的区域，一阶粗糙度），但是细部信息（二阶粗糙度）逐渐被忽略，且有变平滑的趋势.

为了研究新评估指标的采样特征，对每个节理每隔60°选择一个分析方向（见图12），计算不同采样点距下的 $\theta^*_{\rm{r}}$/k. 结合文献中常见的取样距离，从0.2~3.0 mm选择共计12种点距，计算结果如图12所示. 可以看出，总体上，随着采样点距的增大，各个方向上的形貌均逐渐变平缓，粗糙度越来越小，和图11所反映的情况基本一致. 这种情况主要是因为，随着采样点距的增大，形貌中波动较大的起伏度部分(可以称之为一阶起伏度)被保留，而一些细部的不平整度结构(可以称之为二阶起伏度)则被忽略；而粗糙度是起伏度和不平整度的叠加^[29]，因而节理总体的粗糙度降低.

此外，在某些情况下出现了一些异常采样点距特征. 如图12所示，当采样点距为1.0 ~3.0 mm时，随着采样点距小幅度增大，粗糙度评估指标偶尔出现小幅度增大. 这主要是因为，在少数情形下，由于采样间距的特殊性，虽然细部不平整度被忽略，但是起伏度却增大了，总体粗糙度可能略微增大. 如图13所示为一种特殊情况. 当采样点距为0.1 mm时，情形1粗糙度tan⁻¹(Z₂)=22.5°；当采样点距为0.2 mm时，情形1变为情形2，节理粗糙度tan⁻¹ (Z₂)=22.9°，略有增加.

因此，在实际工程中，建议在描述节理粗糙度时采用统一的采样点距，使得在各种条件下的研究结果具有可比性. 此外，考虑到异常采样点距情形主要发生在1.0 ~ 3.0 mm，建议采样点距取在1.0 mm以内，以避免这种异常点距效应的出现.

（1）本研究改进模型与试验结果更加贴近，形式也更为合理，新指标 $\theta^*_{\rm{r}}$/k可以为粗糙度指标 $\theta^*_{{\rm{max}}}$/(C+1)的改进提供一种新的思考.

（2）新指标是三维指标，能够充分描述粗糙度的各向异性. 该指标在不同的剪切方向上连续变化，且能反映正反方向上的敏感性. 采用新指标评估节理面的粗糙度是可行的.

（3）对于小尺度节理，粗糙度随采样点距增大总体上呈减小趋势，这主要是因为采样点距的增大会导致一些二阶粗糙度特征被忽略. 在少数特殊情形下，采样点距增大会导致粗糙度表征指标暂时增大，这主要是因为，在某些情况下，采样点距增大虽然会忽略二阶粗糙度，但是同时也放大了一阶粗糙度特征.

（4）改进模型在粗糙度表征上更为精细，仍须开展更多的扫描试验，以进一步验证新模型对于起伏度差异较大的节理的适用性. 新模型在工程实践中的应用还有待进一步检验.

致　　谢

扫描试验得到了唐少辉博士的指导和马建武硕士的协助，在此表示谢意！