岩土CT图像边界伪影的机制分析及试样尺度估计解读

第25卷第6期
岩石力学与工程学报 Vol.25 No.6
2006年6月 Chinese Journal of Rock Mechanics and Engineering June，2006
岩土CT图像边界伪影的机制分析及试样尺度估计
范留明1，毛黎明2，张志禹2，丁卫华1
(1. 西安理工大学水利水电工程学院，陕西西安 710048；2. 西安理工大学自动化与信息工程学院，陕西西安 710048)
摘要：伪影是不能正确反映试样内部结构的那部分图像，为了减少对图像的误判断，通常将伪影作为一种干扰进行校正。

边界伪影是与CT图像重构算法有关的一种系统误差，但却包含了试样的边界信息。

基于射线成像常用算法，分析了CT边界伪影的特征和产生机制。

研究表明，边界伪影中CT数极大值位置与试样的边界有一定的关系，对边界伪影的研究有助于准确判断试样在图像中的位置，推算试样的几何尺度。

关键词：岩土力学；CT图像；边界伪影；重构算法；试样尺度
中图分类号：TU 43 文献标识码：A 文章编号：1000–6915(2006)06–1198–06 ANALYSIS OF GEOTECHNICAL CT BOUNDARY ARTIFACT
AND EVALUATION OF SAMPLE SCALE
FAN Liuming1，MAO Liming2，ZHANG Zhiyu2，DING Weihua1
(1. College of Water Resources and Hydropower Engineering，Xi′an University of Technology，Xi′an，Shaanxi 710048，China；
2. College of Automation and Information Engineering，Xi′an University of Technology，Xi′an，Shaanxi 710048，China)
Abstract：An artifact is usually regarded as a measure error，but the boundary artifact is an exception，which occurs in the interface between two different media. As a kind of systematic errors caused by reconstruction algorithm of CT images，the boundary artifact contains some information of the sample size. Based on conventional back-projection reconstruction algorithm of CT image，the mechanism of CT boundary artifacts is analyzed. The study shows that the sample boundary is in correspondence with the position of CT boundary artifact maximum and the study on boundary artifact is beneficial to correctly judge the position of sample in CT image and evaluate the geometrical dimension of sample.
Key words：rock and soil mechanics；computerized tomography image；boundary artifact；reconstruction algorithm；sample size
在各个应力阶段变形破裂过程中的内部细节，为细观岩土力学的研究提供了一条新的测试手段。

通过对反映岩土材料内部细观结构变化CT数的观测、分析、统
计及归纳，在细观尺度上探索岩土损伤演化的机制，以期建立岩土损伤演化本构关系[2]。

按照试验步骤，可将岩土CT分为加荷测量、图像显示和CT识别3个阶段，其中在加荷测量方面的研究最为活跃，取得的成果也最多[3
～10]
1 引言
CT(computerized tomography)，即计算机断层摄影(简称层析)技术是一项基于射线成像的无损检测技术，20世纪70年代后期普及于医学诊断，10 a后推广至工业无损检测，80年代末被岩土工程界所重视[1]。

岩土CT的魅力在于，它能够透视岩土试样
收稿日期：2005–01–10；修回日期：2005–05–08
：加载
作者简介：范留明(1968–)，男，博士，1990年毕业于长安大学应用地球物理系，现任副教授，主要从事岩土工程与地质工程方面的教学和研究工作。

E-mail：fanliuming@
第25卷第6期范留明等. 岩土CT图像边界伪影的机制分析及试样尺度估计• 1199 •
方式从单轴发展到三轴，可静态亦可动态；试样可干可湿、可常温可冻结，体现了加载路径和岩土试样的多样化，并促进了试验仪器的改进[11]。

图像显示从灰度图、等值线图[12]、伪彩增强图[13，14]
发展到
由二维数据重建的三维图像
[15～17]。

CT识别实际上
是从图像中提取信息，旨在建立CT数与岩土损伤特性之间的关系，但研究程度相对薄弱，主要采用基于概率理论的统计分析方法，即使对裂纹的量测也多采用上述方法[18
～20]
，基于CT图像重构理论的
识别方法鲜见报道。

CT识别的一个不利因素是伪影问题。

所谓伪影就是不能正确反映试样内部结构的那部分图像，CT试验中有多种伪影，通常可通过CT机校正。

其中边界伪影与CT图像重构算法有关，是一种系统误差。

对边界伪影分析研究有助于准确识别试样在图像中的位置，据此可估算试样的几何尺度。

2 CT机工作原理
基于透射射线理论的CT图像重构技术已经成熟，并得到广泛的应用，其基本思想为：被测物体放置在射线源与探测器之间，射线源所发出的射线穿透被测物体后必然引起射线速度、强度、频率等物理量数值的变化，这些数据将会被探测器所检测到。

在每一个方向上，都会有一组射线穿透被测物体，被测物体包含在这组射线所组成的几何区域中，所测数据集称为此方向上的CT投影(见图1)。

通过转动或平移改变射线源(或探测器)位置(见图2)，则可以得到不同方向的CT投影，据此可重构CT图像。

在CT装置中，射线源可以是(超)声波、电磁波、微波、核磁共振(NMR)，X射线以及其他粒子流，其中X射线应用最为广泛。

X射线可以穿透非金属，不同波长的X射线有不同的穿透能力，而不同物质对同一波长的X射线的吸收能力也不相同，从射线源发出的X射线穿透
图1 X射线几何投影图
Fig.1 Geometrical projection by X-ray
图2 CT机工作原理图 Fig.2 CT working principle
物体后的衰减服从比尔(Beer)定理：
I=I0exp[−(μ1l1+μ2l2+L+μnln)] (1)
式中：μi，li分别为第i层物质的衰减系数和厚度。

式(1)反映了初始强度为I0的射线源穿过n层物质衰减后的强度I。

对于被检物体的某一截面，实际上是一定厚度薄片体，CT检测的目的就是要计算截面内的衰减系数μ (x，y)。

为了显著反映物质衰减系数的差别，通常用水的衰减系数
μw
做参照定义CT数(用H表示)，
据此表示表征物质对X射线的吸收特性，H与衰减系数μ存在如下关系：
H=μ−μw
μ μ×1 000 (2)
w
3 CT图像重构算法
根据投影数据重建(重构、再生)对象的过程称为CT图像重建(重构)，是一类射线成像方法。

与其他射线成像方法不同的是，X射线的折射率→1，在
CT图像重建，可以忽略其折射特性，这也是X射线CT成像算法优于其他CT算法的主要原因。

比较成熟的CT图像重建算法[21]有联立方程组法、反投影法、傅里叶变换法、卷积法以及逐步逼近法等等，其中卷积法在CT机上得到最广泛的应用。

3.1 卷积法
卷积法又称卷积滤波法或卷积反投影法，重建
表达式为
f(x，y)=∫
2π+∞
∫
g(ξ，θ)h(xcosθ+
ysinθ−ξ)dξdθ (3)
其中，
• 1200 • 岩石力学与工程学报 2006年
h(ρ)=∫ +∞
−∞
|r|e2πr(xcosθ+ysinθ)jdr (4)
为了将式(4)写成卷积形式，令
ρ=xcosθ+ysinθ (5) g(ρ，θ)h(ρ)=∫ +∞
−∞g(ξ，θ)h(xcosθ+
ysinθ−ξ)dξ (6)
由式(3)～(6)，可得到二维CT图像重建公式卷积形式的积分表达式：
f(x，y)=∫
2πg(ρ，θ)h(ρ)
dθ (7)
式中：g为投影函数，对于X射线CT，有
g=1
Lln(I0/I) (8)
式中：L为射线穿过岩土试样的距离，这里忽略了
X射线在空气中的衰减。

h(ρ)称为ρ
滤波器，它是造成CT图像边界效应
的主要原因，下面着重分析ρ
滤波器及其产生计算
误差的机制。

3.2 ρ
滤波器的特性
利用卷积法重建CT图像，可以分为两步进行，首先对投影数据进行卷积运算，然后按照不同方向
θ角做求和运算。

在实际数值计算中，需要进行离
散处理，而且被测对象长度是有限的，假如采样间隔为d，则频率范围为(－0.5～0.5) d－1
，于是ρ
滤波
器为
h(ρ)=∫
12d −1|r|e2πrρjdr (9)
2d
对式(9)进行积分，并整理后，得
h(ρ)=1⎡sin(πρ/d)sin2(πρ/2d)⎤
2d2⎢−⎣πρ/d2(ρ/2d)2 (10)
⎦
其频率响应函数为
H(r)=|r| (−0.5d−1≤r≤0.5d−1) (11)
ρ
图3 ρ 滤波器函数Fig.3 ρ filter function
图4 ρ 滤波器的频率响应Fig.4 Frequency characteristics of ρ filter
由于滤波器在截止频率|0.5 d－
1|处突然间断，将
会在该点附近产生较大的振荡，从而引起较大的噪音。

图3所示为ρ滤波器[21]，图4所示为其频率响应，由图可见，在±0.5d
－1
处存在一截断。

图5所示
是ρ滤波器与矩形窗卷积结果，矩形窗长度(40 d)为两倍滤波器长度(20 d)。

由图可见，在矩形窗两端产生对称的振荡波形。

中间出现约一个滤波器长的平直段((－10～10) d)，将两边振荡波形分开。

这一现象表明，影响范围与滤波器长度有关。

只有在超过滤波器长度之外的区域，才能够完全摆脱滤波器长度的影响，滤波结果是可靠的。

尽管后来提出许多改进方法，虽然在一定程度上改善了ρ滤波器在边缘附近的振荡特性，但是却无法完全消除由此引起的误差。

15
10
5
μ0
－5 －10 －15
－40 d
－20 d
020 d 40 d
x
图5 ρ 滤波器产生的边缘效应
Fig.5 Boundary effect of
ρ filter
4 CT数值模型算例
为了进一步分析CT图像边界伪影的产生机制，选择均匀圆柱截面试样进行CT 重建。

重建函数为
X射线衰减系数μ (x)，投影函数依式(8)，显然投影函数为一矩形窗函数，且在各个投影方向上一致，即与投影方向无关。

考虑到圆的对称性，只需要对
x轴(y = 0)进行重建计算即可。

离散后重建公式为
第25卷第6期范留明等. 岩土CT图像边界伪影的机制分析及试样尺度估计• 1201 •
Nx
2
f(x)=
∑h(kd) (12)
k=−
N2k
5 CT测试实例
图7(a)～(c)是某一断面在单轴荷载分别为18.85，20.59
，31.03 MPa作用时的CT图像，试验设备为SIEMENS-PLUS
螺旋CT扫描仪器，CT图像规格为512×512。

样品为直径50 mm，高100 mm的砂岩试样，峰值强度为31.93 MPa，图中色标为医用CT数(水的CT值为0)。

图8(a)～(c)是该砂岩试样沿水平直径方向的CT数分布图，分别对应于图7(a)～(c)中白线指示的位置，位于图像255行，相当于圆截面的直径。

图7中，水平方向为试样的位置(像素单位)，垂直方向为CT数Hμ，原点位于图像中心处。

在图8(a)～(c)中，CT数分布总体上为一矩形函数，高值约为1 600～1 700，为砂岩的CT数；低值约为－1 000～－900，接近于空气的CT数。

在两种介质的分界处，也出现与图6类似的振荡现象，只是振荡幅度相对减小，主要表现为极大值
与高CT数以及极小值与低CT数之间差别较小，提高了成像精度，但极值点仍然清晰可辨。

CT数极大值均出现在左侧的－172和右侧165处(像素单位)，与图7中高灰阶白色的试样边界相对应。

极大值位置并不严格对称，表明试样圆心与图像中心并不一致，有微小偏移。

由于图8(a)～(c)中CT数极大值的位置都是一致的，可以推断该点是CT图像的特征点能够作为测量试样的参考点。

以0.15 mm /pixel计(据实测资料)，试样直径约为50.55 mm，与实际尺度非常吻合。

μ(x)=∑f[xcos(kΔθ)]Δθ (13)
其中，
⎧≤x w−Nd⎪N |x|
Nx=⎨
⎪⎩(xw−2x)/d |x|＞xw−Nd
式中：d为距离，Δθ为角度采样间隔，N为滤波器长度，xw为窗口宽度。

取d = 0.1，N = 20，xw = 4Nd，
μ= 1.0重建结果如图6所示，图中横轴表示成像
范围，纵轴表示衰减系数μ。

由图6可知，重建图
像与图5极其相似，在矩形窗两端产生对称的振荡现象，极大值位于－40 d和40 d处，两者距离约为ρ滤波器长度的4倍，与矩形窗长度一致，相当于圆形试样的直径大小。

－40 d和40 d外侧的平直段μ值为0，与真空中X射线情形一致；其内侧平直段μ值接近于1，反映了均质试样的衰减系数。

151050
－5－10－15－80
－40 0
μ
40 80 x
图
6 CT
数值模型图
Fig.6 Numerical model image of CT
512
(a)
(b) (c)
图7 不同荷载阶段砂岩试件在同一扫描截面处的CT图像
Fig.7 CT images of sandstone sample in the same scanning section with different stresses
• 1202 • 岩石力学与工程学报 2006年
2 000
2 000
Hμ
2 000
1 000Hμ
Hμ
1 000
1 000
－1 000
－300 d
0 x
300 d
－1 000
－300 d
0x
300 d
－1 000
－300 d
0x
300 d
(a) (b) (c)
图8 砂岩试件沿水平直径方向的CT数分布图
Fig.8 Distribution of CT data of sandstone sample in horizontal diameter direction
6 结语
CT边界伪影是一种与图像重构算法有关的误差，对图像分析产生不利的影响，通常作为一种误差需要校正。

但是，它又包含了试样的边界信息，CT数极大值位置与试样的边界有一定的关系，据此可估算试样的几何尺度。

当然，由于各类CT机对CT图像重建算法不尽相同，加之各种误差干扰，CT数极大值位置并不一定与试样边界严格对应，但是两者之间必然存在一定的相互关系，通过必要的误差校正，仍然可以判断试样在图像中的位置，推算试样的几何尺度。

参考文献(References)：
[1] 杨更社，谢定义，张长庆，等. 岩石损伤特性的CT识别[J]. 岩石
力学与工程学报，1996，15(1)：48–54.(Yang Gengshe，Xie Dingyi，Zhang Changqing，et al. CT analysis of mechanical characteristics of damage propagation of rock[J]. Chinese Journal of Rock Mechanics and Engineering，1996，15(1)：48–54.(in Chinese))
[2] 葛修润，任建喜，蒲毅彬，等. 岩石材料宏细观试验研究[M]. 北
京：科学出版社，2004.(Ge Xiurun，Ren Jianxi，Pu Yibin，et al. Macro and Meso-test Study of Rock Material[M]. Beijing：Science Press，2004.(in Chinese))
[3] 任建喜，葛修润. 单轴压缩岩石损伤演化细观机制及其本构模型研
究[J]. 岩石力学与工程学报，2001，20(4)：425–431.(Ren Jianxi，Ge Xiurun. Study of rock meso-damage evolution law and its constitutive model under uniaxial compression loading[J]. Chinese Journal of Rock Mechanics and Engineering，2001，20(4)：425–
[4] 杨更社，谢定义，张长庆，等. 单轴受力状态下岩石损伤扩展本构
关系的CT识别研究[J]. 岩石力学与工程学报，1998，17(3)：279–285.(Yang Gengshe，Xie Dingyi，Zhang Changqing，et al. Study of the constitutive relation of damage propagation of rock under uniaxial compression by CT[J]. Chinese Journal of Rock Mechanics and Engineering，1998，17(3)：279–285.(in Chinese))
[5] 马巍，吴紫旺，蒲毅彬. 冻土单轴蠕变过程中结构变化的CT动
态监测[J]. 冰川冻土，1996，18(4)：306–311.(Ma Wei，Wu Ziwang，Pu Yibin. Monitoring the change of structures in frozen soil in uniaxial creep process by CT[J]. Journal of Glaciology and Geocryology，1996，18(4)：306–311.(in Chinese))
[6] 马巍，吴紫旺，蒲毅彬，等. 冻土三轴蠕变过程中结构变化的
CT动态监测[J]. 冰川冻土，1997，19(1)：52–57.(Ma Wei，Wu Ziwang，Pu Yibin，et al. Monitoring the change of structures in frozen soil in triaxial creep process by CT[J]. Journal of Glaciology and Geocryology，1997，19(1)：52–57.(in Chinese)) [7] 李晓军，张登良. 路基填土单轴受压细观结构CT监测分析[J]. 岩
土工程学报，2000，22(2)：205–209.(Li Xiaojun，Zhang Dengliang. Monitoring change of structure of road foundation soil in uniaxial compression test with CT[J]. Chinese Journal of Rock Mechanics and Engineering，2000，22(2)：205–209.(in Chinese))
[8] 蒲毅彬，朱元林. CT用于与冻结、岩及冰的无损动态试验研究[J].
自然科学进展，1998，8(2)：251–253.(Pu Yibin，Zhu Yuanlin. A real-in-time CT undamaged testing study of freezing，rock and ice[J]. Progress in Natural Science，1998，8(2)：251–253.(in Chinese)) [9] 葛修润，任建喜，蒲毅彬，等. 煤岩三轴细观损伤演化规律的CT
动态试验[J]. 岩石力学与工程学报，1999，18(5)：497–502.(Ge Xiurun，Ren Jianxi，Pu Yibin，et al. A real-in-time CT triaxial testing study of meso-damage evolution law of coal[J]. Chinese Journal of Rock Mechanics and Engineering，1999，18(5)：497–502.(in Chinese))
431.(in Chinese))
第25卷第6期范留明等. 岩土CT图像边界伪影的机制分析及试样尺度估计• 1203 •
[10] 葛修润，任建喜，蒲毅彬，等. 岩石疲劳损伤扩展规律CT细观分
析初探[J]. 岩土工程学报，2001，23(2)：191–195.(Ge Xiurun，Ren Jianxi，Pu Yibin，et al. Primary study of CT real-time testing of fatigue meso-damage propagation law of rock[J]. Chinese Journal of Rock Mechanics and Engineering，2001，23(2)：191–195.(in Chinese)) [11] 葛修润，任建喜，蒲毅彬，等. 岩石细观损伤演化规律的CT实时
试验研究[J]. 中国科学(E辑)，2000，30(2)：104–111.(Ge Xiurun，Ren Jianxi，Pu Yibin，et al. Study on the law of minute damnification spread of rock by CT real time test[J]. Science In China(Series E)，2000，30(2)：104–111.(in Chinese))
[12] 陈蕴生，李宁，李爱国，等. 非贯通节理介质细观损伤演化CT
分析[J]. 岩石力学与工程学报，2000，19(6)：702–706.(Chen Yunsheng，Li Ning，Li Aiguo，et al. Analysis of meso-damage progress of non-interpenetrated jointed media by using CT[J]. Chinese Journal of Rock Mechanics and Engineering，2000，19(6)：702–706.(in Chinese))
[13] 蒲毅彬，陈万业，廖全荣. 陇东黄土湿陷过程的CT结构变化研
究[J]. 岩土工程学报，2000，22(1)：49–54.(Pu Yibin，Chen Wanye，Liao Quanrong. Research on CT structure changing for damping process of loess in Longdong[J]. Chinese Journal of Rock Mechanics and Engineering，2000，22(1)：
49–54.(in Chinese))
[14] 范留明，李宁，丁卫华. 数字图像彩色增强方法在岩土CT图像
分析中的应用[J]. 岩石力学与工程学报，2004，23(13)：2 257–2 261.(Fan Liuming，Li Ling，Ding Weihua. Application of pseudo color enhancement method to geotechnical CT image analysis[J]. Chinese Journal of Rock Mechanics and Engineering，2004，23(13)：2 257–2 261.(in Chinese))
[15] Kawakata H，Cho a，Kiyama T，et al. Three-dimensional observations
of faulting process in Westerly Granite under uniaxial conditions by CT X-ray Scan[J]. Tectonophysics，1999，313：293–305. [16] Ueta K，Tani K，Kato T. Computerized X-ray topography analysis of
three-dimensional fault geometries in basement-induced wrench faulting[J]. Engineering Geology，2000，56：197–210.
[17] Bay B K，Smith T S，Fyhrie D P，et al. Digital volume correlation：
Three-dimensional strain mapping using X-ray topography[J]. Experimental Mechanics，1999，39：217–226.
[18] 任建喜，葛修润. 单轴压缩岩石损伤演化细观机制及其本构模型研
究[J]. 岩石力学与工程学报，2001，20(4)：425–431.(Ren Jianxi，Ge Xiurun. Study of rock meso-damage evolution law and its constitutive model under uniaxial
compression loading[J]. Chinese Journal of Rock Mechanics and Engineering，2001，20(4)：425–431.(in Chinese))
[19] 丁卫华，仵彦卿，蒲毅彬，等. 基于X射线CT的内部裂纹宽度测
量[J]. 岩石力学与工程学报，2003，22(9)：1 421–1 425.(Ding Weihua，Wu Yanqing，Pu Yibin，et al. Measurement of crack width in rock interior based on X-Ray CT[J]. Chinese Journal of Rock Mechanics and Engineering，2003，22(9)：1 421–1 425.(in Chinese)) [20] 崔中兴，仵彦卿，吕惠民，等. 砂岩裂宽变化与CT数变化的相对
关系[J]. 岩石力学与工程学报，2004，23(20)：3 484–3 488.(Cui Zhongxing，Wu Yanqing，Lu Huimin，et al. Correlation between crack width and variety of CT values in sandstone sample[J]. Chinese Journal of Rock Mechanics and Engineering，2004，23(20)：3 484–3 488.(in Chinese))
[21] 赵荣椿. 数字图像处理导论[M]. 西安：西北工业大学出版社，
1995.(Zhao Rongchun. Introduction to Digital Image Processing[M]. Xi′an：Northwestern Polytechnical University Press，1995.(in Chinese))。