基于同步辐射装置定量表征煤孔隙结构非均质性和各向异性

孙英峰，赵毅鑫，王欣，彭磊，孙强

（1.中国矿业大学（北京）能源与矿业学院，北京 100083；2.中国矿业大学（北京）应急管理与安全工程学院，北京 100083；3.共伴生能源精准开采北京市重点实验室，北京 100083；4.中国矿业大学（北京）力学与建筑工程学院，北京 100083；5.中国石油勘探开发研究院，北京 100083）

0 引言

煤层气（CBM）已成为世界天然气资源的重要组成部分[1]。全球可采煤层气总量估计为（1.4～8.5）×1014m3[2]。煤是具有复杂孔隙结构的多孔介质[3]，目前定量表征煤孔隙结构的方法主要是压汞法[4-5]和氮气/二氧化碳低压吸附法（LPA）[6]。压汞法由于高压注入可能破坏煤的孔隙结构或引起孔隙结构的变形[4]。LPA测孔的实验结果需要选用合适的分析模型（如Barrett-Joyner-Halenda（BJH）、Dubinin-Astakhov（D-A）、Horvath-Kawazoe（H-K）模型等）进行数据解释[7-11]，但是每个模型都有其适用的孔隙形状和孔径范围，而且煤的孔隙结构具有复杂的几何形貌和宽泛的孔径范围，所以很难找到一个单一的模型来准确分析煤中完整的孔径分布（PSD），因此，有必要探索新的定量表征煤孔隙结构的方法。

2010年以来，利用小角散射（SAXS）和nano-CT表征煤孔隙结构受到越来越多的关注。新一代实验室nano-CT可以实现纳米级的分辨率[12]。然而，实验室nano-CT没有单色光束，也不能实现高X射线相干，这可能影响成像质量，如射束硬化[13]。由于常规SAXS信号强度较弱，测量时间就要拉长，给实验带来很大的困难。测量一个样品往往需要几个小时或几十个小时，即使采用高功率X射线源，也难以满足许多实验工作的需要。而同步辐射光源不但强度高且光斑尺寸好，可以使用很长的准直系统（可达10 m），从而大大提高实验的分辨率和灵敏度，大幅度节约实验时间（一般样品曝光时间仅需几十秒），简化繁琐的数据修正工作。由于同步辐射光的波长连续可调，在实验中可选择适当的波长来消除SAXS中较难解决的多次散射问题[14]。因此，同步辐射SAXS与常规SAXS相比，具有突出的优点。

国际上应用同步辐射 SAXS对煤的孔隙结构开展了一些研究，Okolo G N等[15]采用低温氮吸附、压汞法和同步辐射SAXS 3种技术研究了4种煤样，指出SAXS能探测的孔径范围更广，且包括闭孔和开孔，因此其测定样品的比表面积和孔隙度最大。Radlinski A P等[16]借助 SAXS和 SANS（小角中子散射）测量了 7组原煤和相应粉末状样品，研究认为球团模型充分代表了煤样的微观孔隙结构。赵毅鑫等[17]采用TEM（透射电子显微镜）和同步辐射SAXS测量了6种不同变质程度的煤样，发现孔隙孔径分布与煤阶无关，但随镜质组含量的变化而变化。Luo L等[18]采用同步辐射SAXS研究了不同变质程度的超细煤粉，得出孔隙表面分形维数随煤阶的增加而减小。王博文等[19]采用同步辐射 SAXS研究表明，低挥发分烟煤脱灰分后，其孔隙总数和体积增加，微孔占比增大，孔表面分形维数和孔比表面积增大，但平均孔径和最可几孔径减小。作者团队2018年以来，应用同步辐射nano-CT研究表征煤孔隙结构的实验方法和图像处理方法，探索基于连通孔隙的 CFD（计算流体动力学）数值模拟方法，CFD数值模拟得到的渗透率与其他学者得到的煤的渗透率在同一数量级[20]。

定量表征煤孔隙结构非均质性和各向异性对于揭示气体在煤层中的赋存和运移规律具有重要意义[21-26]。根据文献研究，对于孔隙结构非均质性和各向异性的定量表征方法已经开展了大量的研究[27-36]，但是鲜有基于同步辐射装置定量表征煤孔隙结构非均质性和各向异性的研究。同步辐射 SAXS和同步辐射 nano-CT都具有纳米级的分辨率，同步辐射 SAXS孔径分辨率可以达到1 nm，同步辐射nano-CT的空间分辨率可以达到30 nm，而且对孔隙几何形貌没有要求，这为定量研究煤的孔隙结构提供了可能性。因此，本文利用同步辐射 SAXS和同步辐射 nano-CT，进行定量表征煤孔隙结构非均质性和各向异性研究，以期探索出新的定量表征煤孔隙结构非均质性和各向异性的方法。

1 样品和实验方法

1.1 煤样

实验中两个不同变质程度煤样分别采自山西省忻州窑煤矿 11＃煤层和河北省唐山煤矿的 9＃煤层。煤的工业分析[37]结果见表1。设备采用 MDMDY-300全自动密度仪，通过氦气置换法来测量两个煤样的真密度，忻州窑煤的真密度为1.350 2 g/cm3，唐山煤的真密度为1.535 5 g/cm3。根据中国石油天然气行业标准SY/T 5124—2012[38]，使用Leitz MPV-3显微光度计对抛光样品进行镜质组反射率测量和显微组分定量分析。从镜质组平均最大反射率来看，忻州窑煤为气煤，唐山煤为肥煤，忻州窑煤的镜质组含量更高，两个煤样的惰质组含量非常接近。根据中国石油天然气行业标准SY/T 5163—2010[39]，应用D/MAX 2500 X射线衍射仪对两组煤样进行了X射线衍射（XRD）分析。XRD分析结果显示，唐山煤的黏土矿物含量更高（见表1）。

表1 煤样成分分析

1.2 同步辐射SAXS实验

由于常规的孔隙结构表征方法难以探测到煤中的闭孔，为了全面（包括闭孔和开孔）表征煤中孔隙结构的非均质性，本文探索利用同步辐射 SAXS定量表征煤中孔隙结构非均质性的方法，同步辐射SAXS实验利用北京同步辐射装置的 1W2A小角散射站开展。

同步辐射SAXS光路示意图如图1所示。角分辨率为0.5×10-3rad，能量分辨率约为1×10-3，样品处光强可达到1×1011cps，比常规X光机的强度高1×103～1×104倍，入射 X射线波长为 0.154 nm，光斑尺寸为1.4 mm×0.2 mm。在同步辐射SAXS实验中，实验样品选用粒径为0.18～0.25 mm（60～80目）的煤颗粒，用3M胶带将样品安装在样品仓上，将样品仓安装至综合平台上，启动同步辐射 SAXS装置，设置曝光时间和曝光次数，记录光电二极管读数，保存散射图像数据。

图1 同步辐射SAXS光路示意图

当X射线照射到试样上，如果试样内部存在纳米尺寸（2～100 nm）的密度不均匀区时，则会在入射X射线的周围2°～5°的小角度范围内出现散射X射线。一个电子在不同方向的散射强度由汤姆逊（Thomson）公式决定[40]：

对于多粒子体系，如果粒子之间的相互距离远远大于粒子本身的尺寸，散射强度则可近似为单个粒子散射强度的简单加和。因此，对于M个不相干涉的粒子体系，其散射强度可认为近似满足Guinier定律[41]：

由（3）—（5）式可得，散射强度随散射角度的增加而减小。通过SAXS实验观测可得到q和I(q)的关系，对其进行数据处理，进而获得散射体相关结构信息，如散射体的回转半径、粒径大小、形状、分形维数、界面层厚度等。

1.3 同步辐射nano-CT实验

为了对比验证同步辐射 SAXS定量表征煤中孔隙结构非均质性方法的有效性，同时考虑到同步辐射nano-CT在探测闭孔和获得孔隙结构三维几何形貌方面具有的优势，本文探索利用同步辐射 nano-CT定量表征煤中孔隙结构非均质性的方法，同步辐射nano-CT实验利用北京同步辐射装置的4W1A-X射线成像实验站开展。

同步辐射nano-CT的能量范围为5～12 keV，光斑尺寸为10 µm×10 µm，可以达到的空间分辨率优于30 nm，对于该同步辐射nano-CT的详细介绍见文献[42]。对于同步辐射 nano-CT成像，样品颗粒的尺寸要求小于 10 μm，因此，首先粉碎煤样，然后选择尺寸接近10 μm的颗粒，并通过图2c所示的仪器粘在大头针的针尖上。由于后续图像校准需要参考点，通过图2d所示的仪器将直径为 0.5 µm 的球形金颗粒安装在样品上，用作后期图像校准的参考点。最后，将粘有样品的大头针夹在如图2b所示的样品转台上。在同步辐射nano-CT成像过程中，样品台以预设的角度步（旋转速度）不连续地旋转以获得不同角度的投影图像。

图2 同步辐射nano-CT成像用到的仪器照片

2 实验结果和讨论

2.1 煤孔隙结构非均质性定量表征

2.1.1 同步辐射SAXS定量表征煤孔隙结构非均质性

同步辐射 SAXS实验所获得的两种不同变质程度煤样的SAXS散射图像如图3所示。

图3 煤的散射图像

为了方便数据分析，通过FIT2D软件将图3所示的SAXS散射图像转为对应的散射曲线，如图4所示。其中，横坐标|q|为倒空间的散射矢量的模，纵坐标为散射强度。

图4 煤的散射曲线

分形物体的X光小角散射强度在分形区可用幂律表示：

α为[0，4]，取曲线线性区域的斜率为d，则α=-d。当0＜α＜3时，散射体为质量分形或孔分形，其分形维数D2=α；当3＜α＜4时，散射体为表面分形，其分形维数D1=6-α。因此，若曲线存在线性区域，则表明存在分形现象，进而可通过α大小判断分形为表面分形还是孔分形[43]。

图5为两个煤样的|q|和散射强度I对应的双对数曲线。通过在不同|q|值范围下，对两个煤样的双对数曲线作切线的方法，可得相应|q|值范围段煤样孔隙的分形维数。

在低|q|值区域（-3～-1 nm-1）曲线的分形特征为表面分形（D1），由此可分析孔隙表面的光滑度和平整度，其分形维数描述孔隙表面光滑度和平整度，其值越大，则说明孔隙表面越不规则。高|q|值区域（-1～1 nm-1）曲线的分形特征表现为孔分形（D2），其分形维数可分析孔隙空间分布的非均质性，其值越大，说明孔隙空间分布的非均质性越强。由图5的表面分形（D1）可以看出，忻州窑煤孔隙表面更不规则，由图5的孔分形（D2）可以看出，忻州窑煤孔隙空间分布的非均质性更强。

2.1.2 同步辐射nano-CT定量表征煤孔隙结构非均质性

同步辐射 nano-CT获得煤的成像以后，图像处理流程主要包括图像校准、三维重构、兴趣区域（ROI）选取、噪声滤除和图像分割，如图6所示。图像处理具体方法见文献[20]。

图像分割以后，对孔隙结构的连通性进行了分析。连通孔隙指ROI立方体从一个面到对面连续连通的孔隙，两个煤样中的连通孔隙如图7a和 7b所示。通过分析计算，忻州窑煤连通孔隙的体积是 1.71 µm3，占总孔隙体积2.47 µm3的69%，而唐山煤连通孔隙的体积是2.70 µm3，占总孔隙体积3.24 µm3的83%，这说明唐山煤有更多的连通孔隙。

图5 煤样分形求解对应的双对数曲线

图6 同步辐射nano-CT图像处理

图7 连通孔隙的几何形貌

煤样ROI中不同子块孔隙度之间的差异可以表征孔隙结构非均质性，本文中煤样ROI尺寸均为200像素×200像素×200像素，像素尺寸为 0.014 59 μm×0.014 59 μm× 0.014 59 μm。煤样 ROI中不同子块孔隙度之间的差异可以通过计算各个子块孔隙度的相对标准偏差来定量表征：

如图8所示，当煤样 ROI在（X，Y，Z）方向上分别被等分成s段，煤样ROI被分成s3个子块，各个子块孔隙度的相对标准偏差H随分割段数s的变化而变化。所以有必要提出一个与子块数量无关的参数来表征煤的三维孔隙结构的非均质性。虽然曲线随着分割段数s的增加而上升，但是曲线变得越来越平缓。为了获得一个常数来表征煤的三维孔隙结构的非均质性，用下面的函数来拟合图8中的数据，函数H（s）的极值是一个常数：

图8 子块孔隙度相对标准偏差与分割子块数量的关系

HV值可以表征煤样三维孔隙结构的非均质性，HV值越大，孔隙结构的非均质性越强。从（9）式可以发现该方法得到的非均质性值消除了分割子块数量的影响。

忻州窑煤：A=0.210 47，B=0.065 62，R2=0.999 41，HV=3.21；唐山煤：A=0.175 47，B=0.064 63，R2=0.994 02，HV=2.71。

如图8所示，两个煤样的非均质性值可以很好地用（8）式拟合。计算结果表明，忻州窑煤的HV值为3.21，而唐山煤的HV值为2.71；可见，忻州窑煤比唐山煤孔隙结构表现出更高的非均质性，也可以从图8中两条曲线的对比中得到相同的结论。

关于两个煤样孔隙结构的非均质性，忻州窑煤孔隙结构的非均质性要比唐山煤强，同步辐射 SAXS和同步辐射 nano-CT得到的结果一致，这也验证了这两种表征煤孔隙结构非均质性方法的有效性。虽然孔隙结构非均质性分级标准的建立需要对大量煤样进行统计，但是借助同步辐射SAXS和同步辐射nano-CT分别获得的孔分形（D2）和HV值为定量描述孔隙结构的非均质性提供了手段，这两个指标可以实现对多个样品的孔隙结构非均质性进行排序，也为孔隙结构非均质性分级标准（强-中-弱）的建立提供了可能。

根据文献研究，煤孔隙结构的非均质性随镜质组含量的增加而增加，随着惰质组含量的增加而降低[44]，与灰分含量呈现线性正相关关系[45]。忻州窑煤的镜质组含量高于唐山煤，两个煤样的惰质组和灰分含量接近，忻州窑煤孔隙结构的非均质性强于唐山煤，这与前人的相关研究结果一致。

2.2 煤孔隙结构各向异性定量表征

考虑到孔隙结构各向异性和渗透能力各向异性之间的对应关系，孔隙结构各向异性的定量表征可以通过运用LBM（格子Boltzmann方法）方法计算渗透率张量来实现，孔隙结构的各向异性通过渗透率张量的特征值和特征向量来表征。渗透率张量的特征值反映了由相应特征向量确定的 3个方向上的渗透率。通过特征值的对比分析可以找到孔隙结构渗透性最好的方向，渗透率张量特征值中最大值对应的特征向量确定的方向上渗透性最好，而且特征值中的最大值也代表了孔隙结构可以达到的最佳的渗透率，为对比不同样品的渗透性提供了参考。

本文中用到的LBM方法采用单松弛时间Bhatnagar-Gross-Krook（BGK）模型[46]，本模型通过求解一个离散的 Boltzmann方程来模拟格子节点中迁移和碰撞的流体颗粒分布的变化，然后宏观流体流动用虚拟流体颗粒的整体运动来近似。在方向ei上的流体颗粒分布函数fi在每一个时间步通过以下的方程来更新：

式中feq是截面平衡分布，定义为：

τ是与流体动力黏度υ相关的参数，在 Chapman-Enskog方程中表示为：

ρ和v是宏观的密度和速度，它们由颗粒分布的变化决定：

根据达西定律，渗透率张量可以由下式计算[47]：

通过3个独立的模拟，在3个流动方向（X，Y，Z）的两个相对面上施加压力梯度并在剩余 4个面上施加无流动边界条件，计算出渗透率张量的3个对角值。

在渗透率张量是对称的和正定的假设下，达西定律可以扩展为：

如（16）式所示，如果得到v1、v2、v3，渗透率张量的非对角值可以由渗透率张量已经得到的对角值计算。

为了获得v1、v2、v3，进行了另一组模拟，在这组模拟中，压力梯度同时施加在3个正交方向上。图9展示了4组LBM模拟的流速流线图，可以发现流线分布与压力梯度的方向一致。

为了分析三维孔隙结构的各向异性，计算了两个煤样ROI的渗透率张量的特征值和特征向量：

忻州窑煤的渗透率张量计算结果是：

图9 渗透率张量计算得到的流速流线图

特征值是k1=-4.10×10-6μm2，k2=0.42×10-6μm2，k3=5.33×10-6μm2，相应特征向量是：

唐山煤的渗透率张量计算结果是：

采用同样的方法，计算出渗透率张量的特征值是k1=-5.30×10-6μm2，k2=0.54×10-6μm2，k3=6.69×10-6μm2，相应的特征向量是：

可以看出唐山煤的渗透率张量的特征值比忻州窑煤对应的特征值大，这和唐山煤有更多的连通孔隙的相符合。

孔隙结构各向异性分析的目的是找到渗透性最好的方向，但是值得注意的是，两组煤样的渗透率张量特征值都有负值，表明在相应特征向量所确定的方向上存在U形连通孔隙。渗透率张量特征值中最大的值对应的特征向量所确定的方向上渗透性最好，在该方向上连通孔隙的截面积也最大，渗透率张量特征值中第 2大的值对应的特征向量所确定的方向上渗透性第2好，在该方向上连通孔隙的截面积也第 2大，在图7a和图7b中，相应特征向量所确定的方向被标出，发现渗透率的各项异性和孔隙结构几何相貌的各向异性之间有明显的相关性。特征值比较大的方向，连通孔隙的横截面积也比较大。

2.3 表征煤孔隙结构不同方法的对比分析

为了体现同步辐射SAXS和同步辐射nano-CT方法的特点和优势，本文将同步辐射 nano-CT和同步辐射 SAXS测得的孔径分布与低温氮气吸附解吸和核磁共振冻融测得的结果进行了对比（见图10）。低温氮气吸附解吸和核磁共振冻融实验和数据处理方法在作者之前发表的论文中做了详细介绍[48]。从图10a可以发现，孔径小于100 nm时，低温氮气吸附解吸和核磁共振冻融测得的孔体积比同步辐射 nano-CT要高，这可以由同步辐射 nano-CT分辨率的局限性解释，一些孔径小于同步辐射 nano-CT分辨率的孔隙不能被同步辐射 nano-CT探测到，但是可以被低温氮气吸附解吸和核磁共振冻融探测到，而当孔径大于100 nm时，同步辐射 nano-CT测得的孔体积比其他两种方法要大，这是由于低温氮气吸附解吸和核磁共振冻融探测不到闭孔，而同步辐射nano-CT可以探测到闭孔所导致的。同步辐射 nano-CT在探测闭孔方面具有优势，但是由于分辨率的限制，在探测孔径较小的孔隙方面存在不足。

图10 不同方法测得的煤孔径分布对比

从图10b可以发现，同步辐射SAXS和低温氮气吸附解吸的测试结果具有很好的一致性，在0～40 nm和40～80 nm孔径范围内分别具有两个峰。在20～25 nm孔径范围内，同步辐射SAXS测得的孔隙体积要大于低温氮气吸附解吸，这可以由同步辐射 SAXS可以测得闭孔，而低温氮气吸附解吸测不到闭孔来解释。在其他范围内同步辐射 SAXS测得的孔体积小于低温氮气吸附解吸，出现这种现象的原因还有待进一步研究。在笔者开展的应用同步辐射 SAXS和核磁共振冻融测试页岩的孔隙结构中，也出现在部分区域同步辐射 SAXS测得的孔体积大于核磁共振冻融，而在其他区域小于核磁共振冻融的情况[41]。虽然同步辐射SAXS测得的孔径范围不能覆盖全部的孔径范围，但是如图11所示，根据核磁共振冻融测试的结果，孔径小于80 nm的孔体积占孔径500 nm以下的孔体积的74%；所以，同步辐射 SAXS所得的孔径范围对于分析煤的孔隙结构还是具有一定意义的。

总之，与常规的孔隙结构表征方法相比，虽然同步辐射 SAXS可探测的孔径范围较小，同步辐射nano-CT难以探测到一些孔径小于其分辨率的孔隙，但二者的优势在于可以探测到闭孔，且同步辐射nano-CT可以获得孔隙结构的三维几何形貌。另外，两种方法的结合使用可以互相弥补各自在探测孔径较大和较小的孔隙时存在的劣势。

图11 核磁共振冻融测得的累计孔体积

3 结论

针对定量表征煤孔隙结构非均质性和各向异性的研究存在的不足，探索了利用同步辐射 SAXS和同步辐射 nano-CT定量表征煤孔隙结构非均质性和各向异性的方法。在利用同步辐射 SAXS定量表征煤孔隙结构非均质性中，由孔分形D2定量表征两个煤样孔隙结构的非均质性，忻州窑煤的孔分形维数为 2.74，唐山煤的孔分形维数为 1.69。在利用同步辐射 nano-CT定量表征煤孔隙结构非均质性中，把煤样ROI分成一定数量的子块，通过计算各个子块孔隙度的相对标准偏差的极值来定量表征孔隙结构的非均质性，忻州窑煤的孔隙结构非均质性值为 3.21，唐山煤的孔隙结构非均质性值为 2.71。两种手段获得的煤孔隙结构的非均质性具有一致性，即忻州窑煤孔隙结构的非均质性强于唐山煤。在利用同步辐射 nano-CT定量表征煤孔隙结构各向异性中，考虑到孔隙结构各向异性和渗透能力各向异性之间的对应关系，孔隙结构各向异性的定量表征通过运用 LBM 方法计算孔隙结构的渗透率张量来实现，利用渗透率张量的特征值和特征向量来定量表征孔隙结构的各向异性，孔隙结构的形貌验证了本文提出的定量表征煤孔隙结构各向异性方法的有效性。

致谢：对中国科学院高能物理研究所李志宏和黄万霞在同步辐射SAXS和同步辐射nano-CT实验中给予的帮助、北京市计算中心李娇娇在数值计算中提供的帮助表示衷心感谢。

符号注释：

a——试样到接收器的距离，m；A——拟合参数，无因次；B——拟合参数，无因次；c——光速，299 792 458 m/s；cs——声速，346 m/s；d——斜率，无因次；D——孔径，nm；D1——表面分形维数，无因次；D2——质量或孔分形维数，无因次；e——电子电荷，1.6×10-19C；ei——方向向量；E1，E2，E3——特征向量，无因次；fi——流体颗粒分布函数；feq——截面平衡分布函数；H——各个子块孔隙度的相对标准偏差，无因次；HV——煤样三维孔隙结构的非均质性值，无因次；ie——常量，7.9×10-20cm2；Ie——单个电子对光子的散射强度；Ie（θ）——单个电子在不同方向对光子的散射强度；I（q）——散射强度；I0——入射X射线光通量，phs/（mm2·s）；k1，k2，k3——渗透率张量的特征值，μm2；kij——渗透率张量，μm2；m——电子质量，9.109 56×10-31kg；M——不相干涉的粒子体系数量，个；n——一个粒子中的电子数目，个；N——子块数量，个；p——流体压力，Pa；q——散射矢量，nm-1；r——颗粒半径，nm；R——粒子半径，nm；Rg——粒子的回转半径，nm；s——煤样ROI在（X，Y，Z）方向上分别被等分成的段数，无因次；t——时间，s；Δt——时间增量，s；2θ——散射角；v——宏观速度，m/s；vi(x)——速度分布函数；V——SAXS辐照的体积，mm3；Vp——孔隙体积，cm3/g；VΩ——积分区域的体积，m3；wi——权系数，无因次；x——流体颗粒的位置，m；α——常量，无因次；γ0——相关函数；λ——波长，nm；μv——流体的运动黏度，m2/s；ρ——宏观密度，kg3/m3；ρm——基质电子密度，eA-3；ρp——微孔电子密度，eA-3；τ——松弛时间，s；υ——流体动力黏度，Pa·s；φs——样品的孔隙度，%；φi——第i个子块的孔隙度，%，φ——煤样整体的孔隙度，%；Ω——积分区域；▽x——x方向梯度算子，Pa/m。