Seismic Facies Analysis Based on Spectral Clustering with Waveform Characteristic Vector
-
摘要:
本文提出一种基于地震沉积学原理沿层提取地震波形特征向量，并以谱聚类（spectral clustering）分析进行地震相划分的方法。谱聚类能够处理非线性的数据结构和高维数据的聚类问题，但其相似度矩阵的构建和谱分解的计算较为复杂，需要较高的计算资源和时间成本。为提高谱聚类算法的效率和可扩展性，本文提出将Mini-batch K-means算法与谱聚类算法结合起来的MKSC算法，在提高谱聚类算法精度的同时大大降低谱聚类空间的复杂度。经过对数值模拟、地球物理模型数据和实际地震资料的处理分析，证明该方法在沉积相划分、沉积相特征识别方面的效果明显，是一种具有良好应用前景的新型沉积特征分析工具。
Abstract:Based on the principle of seismic sedimentology, the feature vectors of seismic waveforms are extracted along stratum slices, and spectral clustering analysis is introduced to classify seismic facies. Spectral clustering is an unsupervised machine learning algorithm. Its essence is to simplify the expression of high-dimensional seismic data in the form of feature vectors, which belongs to the process of dimensionality reduction. Considering the traces with specific time windows in the seismic work area as nodes of the graph and the similarity between traces as the weight of the edges, a graph model can be constructed. Spectral clustering must determine the best segmentation method to complete the segmentation of the graph, so that different types of sedimentary characteristics can be distinguished. Physical model and actual data processing and analysis demonstrate that this method is capable of dividing sedimentary facies characteristics and is a new kind of facies analysis tool for reservoir classification, which has good application prospects.
-
冠状动脉疾病的病死率居高不下，在全球范围内造成了巨大的经济负担和人文负担[1]。随着多层螺旋CT的发展，冠状动脉CT血管造影（coronary artery CT angiography，CCTA）已经成为临床无创评估及诊断冠状动脉疾病的首选检查方法[2-3]。然而CCTA检查过高的电离辐射也受到人们更多的关注[4]。
目前有多种降低CCTA辐射剂量的方法，如前瞻性心电门控扫描、自动管电压和管电流技术、应用迭代重建算法等[2, 5-8]。应用β受体阻滞剂降低患者心率也可间接降低辐射剂量[4]。控制心率对图像质量至关重要，对低心率患者使用较窄的扫描时间窗也可降低患者的辐射剂量[9]。
本研究旨在探究不同心率CCTA检查的最佳采集时相，优化不同心率患者CCTA的心电图扫描时间窗，比较不同心率CCTA的全心动周期（时间窗采集范围为0%~100%）与缩窄扫描时间窗后图像质量和辐射剂量的差异。同时部分患者以DSA结果作为金标准，比较两组患者CCTA检查对冠状动脉及其分支狭窄程度的诊断准确性。
1. 材料与方法
1.1 一般资料
本研究经我院伦理委员会批准，全部患者签署检查知情同意书。前瞻性收集我院自2023年1月至11月行CCTA的患者
1000 例，其中男513例，女487例，年龄58岁（49，66），心率范围37~117 bpm，所有患者均未接受舌下含服硝酸甘油。根据患者CCTA扫描时的心率将患者分为A1和B1亚组（BMP＜51 bpm）、A2和B2亚组（51 bpm≤BMP≤55 bpm）、A3和B3亚组（56 bpm≤BMP≤60 bpm）、A4和B4亚组（61 bpm≤BMP≤65 bpm）、A5和B5亚组（66 bpm≤BMP≤70 bpm）、A6和B6亚组（71 bpm≤BMP≤75 bpm）、A7和B7亚组（76 bpm≤BMP≤80 bpm）、A8和B8亚组（81 bpm≤BMP≤85 bpm）、A9和B9亚组（BMP＞85 bpm）。
排除标准：①碘对比剂过敏；②不能配合屏气；③妊娠；④心脏起搏器植入；⑤心、肾功能不全。
1.2 扫描方案
所有患者在uCT 968 CT扫描仪（联影医疗科技股份有限公司）上进行扫描。扫描参数如下：
管电压：自动管电压调制技术；管电流：ECG mA调制，150~180 mA，管电压及管电流具体数值根据患者体型进行调节；机架转速0.25 s/r；SFOV：20 cm×20 cm；探测器准直宽度：120~160 mm，根据患者心脏大小决定。
两组患者均采用前瞻式心电门控，扫描范围从气管隆嵴下至心底下0.5 cm处。A组患者行全心动周期（时间窗采集范围为0%~100%）进行图像采集。扫描结束后，以5%为间隔重建0%~95%间期的20组图像，根据右冠状动脉、左冠状动脉前降支和回旋支的图像质量确定最佳的一组图像，再以1%的间隔重建以上图像的前后各4组图像，选择以上9组图像中质量最佳的一套图像，作为各亚组患者图像的最佳采集时相，以最佳采集时相的95%置信区间作为优化扫描时间窗，即最佳重建时相（平均值±标准差）×2[10]。B组各亚组患者采用A组优化后的扫描时间窗行CCTA检查。
所有患者使用高压注射器（Ulrichmissouri-XD2001）将碘对比剂（碘佛醇，320 mgI/mL）以4.0 mL/s的速率注射，注射总量为0.8 mL/kg，注射完成后以4.0 mL/s速率注射30 mL生理盐水。
扫描采用自动扫描触发技术，将6 mm2感兴趣区（region of interest，ROI）放置于四腔心层面的降主动脉，ROI内CT值达到132 HU后延迟6 s启动扫描。
原始图像的层厚及层间距均为0.5 mm，图像重建矩阵：512×512，使用C-SOFT-AA滤波函数和混合迭代重建技术（KARL 3 D）进行图像重建，迭代权重为8。
1.3 图像分析
客观评价。测量升主动脉中心，右冠状动脉、左前降支及左回旋支的近段以及同一层面的心包脂肪组织CT值和标准差（standard deviation，SD），测量时避开病变和伪影。计算所有图像的信噪比（signal-to-noise ratio，SNR）、对比度噪声比（contrast-to-noise ratio，CNR），计算公式：
$$ {\mathrm{SNR}}=血管平均{\mathrm{CT}}值/血管平均{\mathrm{SD}}值， $$ (1) $$ {\mathrm{CNR}}=\frac{(血管平均{\mathrm{CT}}值-心包脂肪平均{\mathrm{CT}}值)}{心包脂肪{\mathrm{SD}}值}。$$ (2) 记录所有患者的容积CT剂量指数（volume CT dose index，CTDIvol）和剂量长度乘积（dose length product，DLP），并计算有效剂量（effective dose，ED），
$$ {\mathrm{ED}}={\mathrm{DLP}} \times K，$$ (3) 其中，成人心脏系数K=0.014 mSv/mGy·cm[11]。
冠脉血管狭窄程度的计算公式为：
$$ \frac{(冠脉血管狭窄段直径-狭窄处的直径)}{狭窄血管直径} \times 100 \% 。$$ (4) 狭窄程度＜50%为轻度狭窄，50%~70%为中度狭窄，70%~100%为重度狭窄，100%为完全闭塞，冠脉狭窄≥50%认为是阳性。以DSA结果为金标准，评估CCTA结果对冠状动脉节段狭窄的诊断准确性，并计算CCTA的敏感性和特异性。
$$ 敏感性=\frac{真阳性数}{(真阳性数+假阴性数)} \times 100 \% ，$$ (5) $$ 特异性=\frac{真阴性数}{真阴性数+假阳性数)} \times 100 \%。$$ (6) 主观评价。按照SCCT18段分段标准，评价左前降支、左回旋支和右冠状动脉的图像质量。图像由两位具备10年以上诊断经验的放射医师采用双盲法进行评价：1分，血管伪影严重，无法诊断；2分，血管部分节段伪影明显，影响诊断；3分，血管部分节段轻度伪影，不影响诊断；4分，血管无伪影。最终取两名医师评分值的平均值。主观评分值≥3分认为满足临床诊断要求[12]。
1.4 DSA检查
术前常规向患者简要介绍DSA过程，进行术前沟通，减轻患者紧张情绪。
患者取平卧位，连接心电监护，常规消毒右侧上肢，铺无菌巾，肝素盐水冲洗所有介入器材，取右侧桡动脉途径。1%利多卡因局麻。采用Seldinger穿刺法穿刺成功后置入6F动脉防漏鞘管，注入普通肝素
2500 U、硝酸甘油200 μg。透视下送入导丝，应用5F造影导管行冠状动脉造影。1.5 统计学分析
使用SPSS 26.0软件进行统计学分析。
计数资料采用频数和百分比（n，%）表示，组间比较采用卡方检验。采用Shapiro-Wilk检验分析连续变量数据的正态性，不符合正态分布的连续变量以M(Q1, Q3)表示，并采用Mann-Whitney检验分析差异；符合正态分布且方差齐性的变量以
$ ({\bar{x }}\pm s) $ 表示，并采用单因素方差分析比较组间差异。用 Kappa检验评估观察者间的一致性。k＜0.40为一致性较差，0.4≤k＜0.75为一致性一般，k≥0.75为一致性较好。P＜0.05为差异有统计学意义。
2. 结果
2.1 一般情况分析
两组患者均顺利完成CCTA检查。
结果显示：两组患者的年龄、性别、身高和体重均无统计学差异。
2.2 A组中各亚组CCTA最佳采集时相及扫描时间窗
A组中各亚组CCTA的最佳采集时相及扫描时间窗见表1和图1。
表 1 A组各亚组患者CCTA的最佳采集时相及扫描时间窗Table 1. Optimal reconstruction phase and acquisition time window of CCTA for subgroups in Group A心率/bpm 最佳采集时相 最佳扫描时间窗 ＜51 73±6 61~85 51~55 76±4 68~84 56~60 76±3 70~82 40±3 34~46 61~65 76±3 70~82 40±3 34~46 66~70 76±3 70~82 42±3 36~48 71~75 77±6 65~89 44±3 38~50 76~80 46±5 36~56 76±4 68~84 81~85 46±4 38~54 ＞85 48±5 38~58 ![]()
图 1 4名不同心率患者使用优化后扫描时间窗得到的CCTA心电图采集时相图及图像Figure 1. Four cases of different heart rates using the optimized acquisition window2.3 客观评价
两组患者图像的客观评价值见表2，均无统计学差异。
表 2 A组和B组CCTA图像客观评价及辐射剂量表Table 2. Objective evaluation of CCTA images and radiation doses for groups A and B组别 统计检验 A组 B组 F/Z P 降主动脉 CT值a 457.48±60.59 466.80±61.77 3.457 0.063 标准差a 18.42±2.79 17.95±2.24 5.182 0.230 右冠状动脉 CT值b 439.50（391.60，489.50） 442.85（403.00，492.43） −0.821 0.412 标准差a 16.24±2.25 15.58±1.79 15.912 0.593 左冠状动脉前降支 CT值 433.38±60.85 426.56±60.04 0.644 0.423 标准差b 15.75（13.13，17.40） 15.65（13.43，17.30） −1.205 0.228 左冠状动脉回旋支 CT值b 437.05（401.17，477.87） 426.45（383.25，472.03） −0.433 0.665 标准差a 14.88±3.06 15.57±2.61 3.064 0.082 SNRb 28.88（25.81，3.163） 27.40（24.40，30.79） −1.869 0.062 CNRb 33.66（30.35，36.78） 31.96（29.48，35.34） −1.640 0.101 CTDIvolb/mGy 30.13（27.56，32.67） 17.89（16.24，19.36） 11.952 ＜0.05 DLPb/mGy·cm 414.18（380.58，456.61） 247.52（217.20，270.18） 11.750 ＜0.05 EDb/mSv 5.80（5.33，6.39） 3.47（3.04，3.78） 11.750 ＜0.05 注：a为$(\bar x\pm s )$，b为$M(Q_1,\;Q_3) $。 两组患者的CTDIvol中位数、DLP中位数、ED中位数见表2。统计学结果显示：B组CTDIvol中位数和DLP中位数、ED中位数较A组分别下降40.62%、40.24%和40.17%。
各亚组患者的ED中位数见表3。统计学结果显示：A组和B组各亚组的ED中位数差异有统计学意义。
表 3 A组和B组各亚组患者的辐射剂量Table 3. Radiation doses for subgroups of groups A and B心率/bpm ED中位数/mSv Z P 人数 A组 B组 A组 B组 ＜51 7.36（7.01，7.61） 3.21（2.84，4.01） −5.945 ＜0.05 31 60 51~55 6.57（6.38，6.87） 3.99（3.22，4.17） −7.924 ＜0.05 73 70 56~60 6.13（5.92，6.40） 3.79（3.61，3.88） −10.160 ＜0.05 136 103 61~65 5.68（5.49，5.87） 3.54（3.49，3.61） −8.916 ＜0.05 86 92 66~70 5.46（5.26，5.64） 3.33（3.18，3.46） −7.440 ＜0.05 76 57 71~75 5.18（5.09，5.46） 3.30（3.08，3.37） −6.420 ＜0.05 43 50 76~80 4.92（4.53，5.07） 3.11（2.77，3.19） −4.899 ＜0.05 24 30 81~85 4.64（4.55，4.81） 2.07（2.01，3.05） −3.762 ＜0.05 17 16 ＞85 4.44（4.43，4.62） 1.98（1.98，1.99） −3.606 ＜0.05 14 22 本研究将32例患者的546节段冠脉进行诊断，其中包含A组患者279节段，B组患者267节段。在A组中CCTA检查的敏感性为83.64%（46/55）特异性为94.64%（212/224），B组中CCTA检查的敏感性为84.91%（45/53）特异性为92.99%（199/214），两组的敏感性和特异性无明显差异。CCTA结果和DSA进行比较，结果表明：两组患者CCTA检查对冠状动脉分支狭窄的敏感性与特异性无明显差异。
2.4 主观评价
两位放射医师对两组患者CCTA图像的主观评分值见表4。
表 4 A组和B组CCTA图像主观评价表Table 4. Subjective evaluation of CCTA images for groups A and B分组 A B $\chi^2 $ P 右冠状动脉 3分 247（49.4） 236（47.2） 0.126 0.722 4分 253（50.6） 264（52.8） 左冠状动脉前降支 3分 354（70.8） 326（65.2） 0.827 0.363 4分 146（29.2） 174（34.8） 左冠状动脉回旋支 3分 349（69.8） 315（63.0） 1.100 0.294 4分 151（30.2） 185（37.0） 统计学分析显示：两组患者CCTA图像的主观评分值无统计学差异。两名放射医师的评分一致性较好，k值为0.894。
3. 讨论
患者心率高低是CCTA检查成功的关键因素之一[13]。前瞻性心电门控技术只在心脏运动相对较缓慢的时间进行数据采集，在心电图上的扫描时间窗较普通单心动周期模式缩小，在扫描时间以外的时间减少曝光剂量甚至不曝光，所以是降低CCTA检查辐射剂量的有效手段之一[13]。
有研究证明：低心率（BMP＜61 bpm）患者CCTA图像的最佳采集时相位于舒张中末期，高心率（BMP＞75 bpm）患者CCTA图像的最佳采集时相位于收缩末期；中心率（61 bpm≤BMP≤75 bpm）患者最佳采集时相位于收缩末期或舒张中末期[9-14]。但本研究结果与上述结论有部分不符。
本研究中A3亚组（56 bpm≤BMP≤60 bpm）中有12.50%的患者最佳采集时相位于收缩期；A7亚组（76 bpm≤BMP≤80 bpm）中有29.41%的患者最佳采集时相位于舒张期。究其原因可能为A3亚组患者心率较缓慢，收缩期时间较长，所以收缩期也存在心脏运动较缓慢的时相。A7亚组患者在舒张期仍存在最佳时相，在此心率下仍可进行成像。
本研究还发现80 bmp为CCTA图像从舒张期重建到收缩期重建的临界心率。另外，随着心率增快，最佳重建时相逐渐向后偏移。
为保证CCTA检查的成功率和图像质量，进行CCTA检查时往往采用全心动周期或较宽的时间窗进行扫描[13]。这导致患者的辐射剂量大大增加。
本研究结果显示：A组和B组的图像质量无统计学差异。但B组的ED中位数较A组下降42.56%，这是因为B组扫描时间窗较A组缩窄，患者接受X线曝光时间缩短。B1亚组、B2亚组（BMP＜55 bpm）、B8亚组和B9亚组（BMP＞80 bpm）的辐射剂量降低尤为明显，较A1亚组、A2亚组、A8亚组和A9亚组分别下降56.39%、39.27%、55.39%和55.41%。这是因为B1、B2、B8以及B9仅有一个扫描时间窗，而B组中其他各组均有两个扫描时间窗，且在两个扫描时间窗之间仍有低剂量X线进行曝光。
本研究对部分患者CCTA结果和DSA进行比较，两组患者CCTA检查对冠状动脉分支狭窄的敏感性与特异性无明显差异，此结果表明缩窄采集时间窗不会降低对冠状动脉狭窄的诊断效能。
本研究的不足：①各亚组间样本量不相同，高心率患者样本量较少，后续应增加高心率的样本数量；②本研究属于单中心研究；③未对所有患者的CCTA结果和DSA结果进行比对验证。
综上所述，不同心率患者具有不同的最佳采集时相。低心率患者（BMP＜55 bpm）可将扫描窗缩窄至舒张中末期；中心率患者（55 bpm≤BMP≤80 bpm）可将扫描窗缩窄至收缩末期或舒张中末期；高心率患者（BMP＞80 bpm）可将扫描窗缩窄至收缩末期。通过优化CCTA扫描时间窗可在保证图像质量前提下大大降低患者的辐射剂量。
-
![]()
图 1 无向图及其邻接矩阵特征向量
Figure 1. Feature vectors of the undirected graph and its adjacency matrix
![]()
图 2 相同数据集不同聚类方式分类结果
Figure 2. Classification results of the same dataset by different clustering methods
![]()
图 3 物理模型的采集尺寸及参数示意图
Figure 3. Schematic of the acquisition size and parameters of the physical model
![]()
图 6 物理模型偏移剖面及基于地震沉积学原理解释的层位
Figure 6. Migration section of the physical model and horizon interpretation based on the seismic sedimentology principle
![]()
图 8 物理模型第2至第5层砂、泥岩空间展布及形态
Figure 8. Spatial distribution and morphology of sand and mudstone in the second to fifth layers of the physical model
![]()
图 9 为各层的地层切片均方根振幅地震相
Figure 9. Slices of root mean square amplitude seismic phase of stratums
![]()
图 10 各层MKSC 聚类地震相分布
Figure 10. MKSC clustering seismic facies distribution of stratums
![]()
图 11 目标区基于地震沉积学原理解释的层位
Figure 11. Horizon interpretation based on seismic sedimentological principles in the target area
![]()
图 14 分流河道中心部位不同水动力条件下的相变示意图
Figure 14. Schematic diagram of phase transition in the central part of distributary channel under different hydrodynamic conditions
![]()
图 16 H3c段MKSC聚类属性平面图和沉积相平面图
Figure 16. MKSC clustering attribute plan and sedimentary facies plan of H3c
![]()
图 17 H3b段MKSC聚类属性平面图和沉积相平面图
Figure 17. MKSC clustering attribute plan and sedimentary facies plan of H3b
表 1 不同聚类分析算法计算效率比较
Table 1 Comparison of the computational efficiency of different clustering analysis algorithms
算法类型 Mini-batch K-means Spectral clustering MKSC 计算时间/s 0.17 0.21 0.17+0.09 
下载: 导出CSV
-
[1] ZENG H, BACKUS M M, BARROW T K, et al. Stratal slicing, Part I: Realistic 3-D seismic model[J]. Geophysics, 1998, 63(2): 502−513. doi: 10.1190/1.1444351
[2] ZENG H, HENRY C S, RIOLA P J. Strata slicing; Part II, Real 3-D seismic data[J]. Geophysics, 1998, 63(2): 514−521. doi: 10.1190/1.1444352
[3] ZENG H, HENTZ F T, WOOD J L. Stratal slicing of Miocene-Pliocene sediments in vermilion block 50-Tiger Shoal Area, offshore Louisiana[J]. The Leading Edge, 2012, 20(4): 408−418.
[4] ZENG H, AMBROSE A W, VILLALTA E. Seismic sedimentology and regional depositional systems in Mioceno Norte, Lake Maracaibo, Venezuela[J]. The Leading Edge, 2001, 20(11): 260−269.
[5] BHATTACHARYA S, CARR R T, PAL M. Comparison of supervised and unsupervised approaches for mudstone lithofacies classification: Case studies from the Bakken and Mahantango-Marcellus Shale, USA[J]. Journal of Natural Gas Science and Engineering, 2016, 33: 1119−1133. doi: 10.1016/j.jngse.2016.04.055
[6] 刘爱群, 陈殿远, 任科英. 分频与波形聚类分析技术在莺歌海盆地中深层气田区的应用[J]. 地球物理学进展, 2013, 28(1): 338-344. LIU A Q, CHEN D Y, REN K Y, Frequency decomposition and waveform cluster analysis techniques Yinggehai Basin gas field in the deep area of application[J]. Progress in Geophysics, 2013, 28(1): 0338-0344. (in Chinese).
[7] 白博, 舒梦珵, 康洪全, 等. 基于伪阻抗体波形聚类的贝壳灰岩储层预测方法[J]. 物探化探计算技术, 2015, 37(6): 724-727. BAI B, SHU M C, KANG H Q, et al. Coquina reservoir prediction method of seismic pseudo-impedance waveform clustering[J]. Computing Techniques for Geophysical and Geochemical Exploration, 2015, 37(6): 724-727. (in Chinese).
[8] 李辉, 罗波, 何雄涛, 等. 基于地震波形聚类储集砂体边界识别与预测[J]. 工程地球物理学报, 2017, 14(5): 573-577. LI H, LUO B, HE X T, et al. Boundary identification and prediction of sand body based on seismic waveform[J]. Chinese Journal of Engineering Geophysics. 2017, 14(5): 573-577. (in Chinese).
[9] 徐海, 都小芳, 高君, 等. 基于波形聚类的沉积微相定量解释技术研究−以中东地区 X油田为例[J]. 石油物探, 2018,57(5): 744−755. doi: 10.3969/j.issn.1000-1441.2018.05.014 XU H, DU X F, GAO J, et al. Quantitative interpretation of sedimentary microfaacies based on waveform clustering: A case study of X oilfield, Middle East[J]. Geophysical Prospecting for Petroleum, 2018, 57(5): 744−755. (in Chinese). doi: 10.3969/j.issn.1000-1441.2018.05.014
[10] 刘仕友, 宋炜, 应明雄, 等. 基于波形特征向量的凝聚层次聚类地震相分析[J]. 物探与化探, 2020,44(2): 339−349. DOI: 10.11720/wtyht.2020.1153. LIU S Y, SONG W, YING M X, et al. Agglomerative hierarchical clustering seismic facies analysis based on waveform eigenvector[J]. Geophysical and Geochemical Exploration, 2020, 44(2): 339−349. DOI: 10.11720/wtyht.2020.1153. (in Chinese).
[11] 杨随心, 耿修瑞, 杨炜暾, 等. 一种基于谱聚类算法的高光谱遥感图像分类方法[J]. 中国科学院大学学报, 2019,36(2): 267−274. doi: 10.7523/j.issn.2095-6134.2019.02.015 YANG S X, GENG X R, YANG W T, et al. A method of hyperspectral remote sensing image classification based on spectral clustering[J]. Journal of University of Chinese Academy of Sciences, 2019, 36(2): 267−274. (in Chinese). doi: 10.7523/j.issn.2095-6134.2019.02.015
[12] NG A Y, JORDAN M I, WEISS Y. On spectral clustering: Analysis and an algorithm[J]. Advances in Neural Information Processing Systems, 2001, 14: 849−856.
计量
- 文章访问数: 206
- HTML全文浏览量: 103
- PDF下载量: 50
