0 前言

       在组学和大数据时代，整合影像模态相同但仪器平台、技术参数、检查时间等不同中心和站点的影像数据将成为常态^{[1, 2, 3]}。联合多站点影像数据成功的关键取决于跨站点图像的可比性^[4]，而MRI影像数据受仪器设备、扫描序列等成像环境的影响。不同MRI扫描仪的特性，例如制造商、场强、非线性梯度场和纵向漂移等的不同会增加脑结构、扩散和功能影像数据的偏差和变异性，数据中通常包含非生物混淆变量^{[4, 5]}，从而影响跨站点研究下游分析的一致性和可重复性^{[5, 6, 7, 8, 9]}。因此，影像组学特征面临可重复性和稳定性较差的问题，这将阻碍影像组学模型的探索以及成功转化为更好的诊断、分类或监测工具^[4]。在神经影像学中使用术语“扫描仪效应”或“站点效应”来指代这类非生物变异^[4]。与CT和PET图像分别以kBq/mL和Hounsfield Unit（HU）标化单位表示其强度数值不同，MR图像无标准的强度灰度值，缺乏组织特异性的绝对强度数值含义。因此，在比较基于强度的特征值时，需要患者之间MR图像强度的标准化及校准。既往关于乳腺体模^[7]、脑肿瘤^[9]、前列腺癌^[10]等研究中，已提出几种协调CT、PET-CT影像组学特征的方法，例如混合白色条纹（hybrid white stripe, hWS）、直方图匹配（histogram matching, HM），ComBat法也已在影像组学研究中得到验证，但上述方法只能协调部分特征，且尚未在健康人群脑MRI数据校准中进行前瞻性的验证及比较。

       本研究中，首先，利用青少年健康受试者的脑MRI结构影像学数据即三维T1加权成像（three-dimensional T1 weighted image, 3D-T1WI）、扩散张量成像（diffusion tensor imaging, DTI）验证N4偏置场校正方法。其次，在此基础上分析HM标准化技术对N4偏置场校正后的MR图像强度分布的影响。第三，利用ComBat协调进一步减少影像组学特征的“扫描仪效应”。本研究通过将MR图像N4偏置场校正、HM标准化与ComBat协调相结合，成功减少了扫描仪和站点效应，同时保留了与生物学相关的变异性，从而确保跨扫描仪汇集MRI影像组学研究的更高敏感度和特异性。

1 材料与方法

1.1 研究对象

       本研究于2021年9月至2021年12月间通过广告从当地初级和高级中学招募青少年健康志愿者。由1位精神科医师对参与者进行临床访谈及心理量表测试。纳入标准：（1）年龄13～18岁，性别不限；（2）右利手；（3）流调用抑郁自评量表（Center for Epidemiological Survey, Depression Scale, CES-D）^{[11, 12]}评分小于16分。排除标准：（1）既往或当前患有符合《精神疾病诊断与统计手册》轴I障碍^[13]疾病的发病史及存在的任何精神障碍；（2）有严重颅脑损伤史或神经系统疾病史（如脑部创伤、脑肿瘤、脑出血、脑部炎症及癫痫等）；（3）患有重大或慢性躯体疾病；（4）有精神活性物质滥用或依赖史；（5）MRI检查禁忌证（幽闭恐惧症、体内铁磁性植入物、妊娠等）。

       受试者在访谈及量表评估后3周内接受3次不同设备脑MRI扫描。本研究遵守《赫尔辛基宣言》，并获得昆明医科大学第三附属医院伦理委员会批准（批准文号：KYCS202107）。在充分了解MRI扫描目的和过程后，获得所有受试者监护人的同意并签署知情同意书。

1.2 脑MRI结构影像采集、预处理

1.2.1 MRI影像采集

       在昆明医科大学第三附属医院使用3台不同场强及头线圈的磁共振仪对每位受试者分别进行3次MRI扫描，1号扫描仪（MRI1）为Philip 1.5 T（Achieva，Philips Medical Systems，Netherlands），配备15通道相控阵线圈；2号扫描仪（MRI2）为Philip 3.0 T（Elition，Philips Medical Systems，Netherlands），配备32通道相控阵线圈；3号扫描仪（MRI3）为GE 3.0 T（Pioneer，GE Healthcare，USA），配备21通道相控阵线圈。检查前准备包括：（1）受试者取仰卧位，泡沫垫固定头部以减少成像过程中头部的移动；（2）佩戴橡胶耳塞减少噪音影响；（3）嘱受试者身体保持不动并闭眼，头脑保持清醒、不可以睡眠。依次进行3D-T1WI、DTI、液体衰减反转恢复序列（fluid-attenuated inversion recovery, FLAIR）扫描，参数见表1及表2。扫描FLAIR序列以排除脱髓鞘、肿瘤、血管畸形等病变。序列完成后，即刻由神经影像医生查看图像质量，对有运动或其他伪影的序列进行重新扫描。

点击查看表格

表1  3D-T1WI扫描序列及参数
Tab. 1  Sequences and parameters of 3D-T1WI protocol

点击查看表格

表2  DTI扫描序列及参数
Tab. 2  Sequences and parameters of DTI protocol

1.2.2 3D-T1WI预处理

       使用基于统计参数映射软件的CAT 12^[14]工具包完成。主要步骤：（1）从原空间标准化到蒙特利尔标准脑空间坐标系统（montreal neurological institute, MNI）^{[15, 16]}；（2）利用N4 Bias-Field-Correction算法^[17]进行偏置场校正；（3）将空间标准化后的图像分割成灰质、白质和脑脊液，并根据自动解剖标记图谱（automated anatomical labeling, AAL）将受试者图像分割为116个脑区；（4）把空间归一化过程中产生的变形场作用到分割结果上，消除空间标准化过程中因受试者大脑个体差异而产生的体积变化。由1名具有7年神经影像学经验的放射科主治医师使用ITK-SNAP 4.0.0 beta软件^[18]检查纠正分割错误的脑组织区域。

1.2.3 DTI预处理

       使用FSL即FMRIB's Software Library软件^[19]对DTI数据进行涡流及头部运动校正、脑分割和张量模型拟合等预处理。来自所有参与者的扩散张量数据通过DTI-TK软件^[20]使用基于迭代张量的配准算法，结合局部纤维方向以驱动单个纤维束对齐创建本研究特定的模板。通过ANTs可变形配准将约翰霍普金斯大学ICBM-DTI-81白质图谱的各向异性分数（fractional anisotropy, FA）等图配准到研究特定模板。通过对方向编码彩色分数各向异性图像和各个空间中的包裹图集进行目视检查来进行质量控制。

1.3 偏置场校正

       在图像预处理过程中，使用N4 Bias-Field-Correction图像过滤器对原始图像的偏置场不均匀性进行估计，并配置不同的参数以优化校正效果。N4图像过滤器分辨率级别的数目和每个级别的迭代次数默认设置为四个级别和五十次迭代^{[7, 21]}。通过三个、四个、五个及六个分辨率级别，50或100次迭代，结合默认掩模或全掩模的参数交叉组合来检测不同设置水平时N4算法的校正效果。

       为了评估N4算法校正相似组织类型内强度不均匀性的能力，在每位受试者的3D-T1WI脑白质区使用LIFEx软件^[22]勾画8个（双侧半卵圆中心2个、双侧额叶深部白质2个、双侧枕叶深部白质2个、双侧颞叶深部白质2个）含900个体素的球形感兴趣体积（volume of interest, VOI），依次为VOI1、VOI2、VOI3、VOI4、VOI5、VOI6、VOI7、VOI8，同时避免任何灰质区域，测量其强度，通过比较不同校正的球形小区域的平均强度的变异系数来评估N4算法性能。

1.4 MRI直方图归一化

       MR图像在偏置场校正后，进一步采用分段线性HM归一化将像素值0～10000标准化到（-1，1）之间。HM法包括两个阶段：首先，HM学习标准直方图的界标，然后将图像直方图的界标非线性映射到标准直方图的界标以对齐强度分布。HM独立应用于两个序列，代码改编自REINHOLD等^[23]和SHAH等^[24]定义的十分位地标和标准比例尺。通过定性比较强度直方图对齐来评估归一化对校正受试者个体间和设备间变异性的影响。

1.5 MRI影像组学特征提取

       在MR直方图归一化之后，同一序列获得3个MRI影像集，包括原始影像、N4偏置场校正影像、HM归一化的N4校正影像。将分割后的脑灰质、白质数据导入LIFEx软件^[22]，按照IBSI指南^[25]和RQS标准^[26]，在提取特征之前，使用最近邻插值将所有6×3个（3台设备×2个序列×3次转换）MRI影像集进行重采样，单体素分辨率为1 mm×1 mm×1 mm。对于影像组学特征计算，采用绝对离散化方法，对于每个序列，在每个步骤分别计算区域内的最小和最大强度以确定强度范围，并将区域内强度的平均标准差定义为固定的bin大小。最后提取各个脑灰质和白质区域的形状特征（5个）、一阶直方图特征（12个）、二阶纹理特征（46个）共三大类63个影像特征。

1.6 MRI影像组学特征ComBat协调

       由于即使在N4偏置场校正、HM直方图归一化后，同一受试者的3台MRI设备所采集的图像的影像组学特征值也可能不同，因此对于每个序列，直方图归一化后提取的影像组学特征使用ComBat方法在3台设备之间进行协调。JOHNSON等^[4]在基因表达分析的背景下，引入以下ComBat位点和比例（L/S）调整模型，见公式（1）。

       在上述函数中，αv是体素v处的整体特征值，如皮层厚度、体积或FA值等，X是感兴趣的协变量（例如性别、年龄）的设计矩阵，βv是对应于X的回归系数的体素特定向量。假设误差项εijv服从均值为零且方差为σv2的正态分布。γiv和δiv分别表示站点 i对体素v的加性和乘性站点效应。ComBat使用经验贝叶斯框架来改进参数估计γiv和δiv的方差。超参数γi，τi2，λi，θi是根据JOHNSON等^[4]所述的数据经验估计的，场地效应参数γiv*和 δiv*是使用条件后验方法计算的。因此，最终的影像特征值ComBat协调模型见公式（2）。

       ComBat分别应用于T1和DT1测量的特征值。

1.7 统计学分析

       影像特征值使用R软件进行统计分析，P＜0.05认为差异有统计学意义。8个白质VOI脑区影像数据，通过N4偏置场校正、HM归一化及ComBat协调，最终产生2个序列（T1、DTI）×4次转换（原始、N4偏置场校正、HM归一化、ComBat协调）共8个数据集，每个VOI数据集包含63+16共计79个影像组学特征，使用Kruskal-Wallis检验评估每个影像特征统计分布的设备间差异。通过定性比较ComBat协调方法对校正线圈间变异性的影响。

       为了方便计数和直观显示，根据Kruskal-Wallis检验的P值结果将影像特征分为以下五个P值范围：P＜10^-5、10^-5≤P＜10^-3、10^-3≤P＜0.01、0.01≤P＜0.05和P≥0.05，并列举通过N4偏置场校正、HM归一化、ComBat协调后每个范围内的特征数量。

       为了检验ComBat协调是否会降低影像特征的判别能力，进一步评估ComBat对区别T1WI序列中两种脑组织（灰质和白质）的影响，分别在每位受试者的T1WI脑白质区（右侧半卵圆中心、左侧颞叶深部白质）、灰质区（右侧尾状核头、左侧背侧丘脑）使用LIFEx软件勾画4个含600个体素的3D球形VOI区域，依次为VOIc1、VOIc2、VOIc3、VOIc4，使用Wilcoxon检验比较上述脑灰质与白质ComBat协调前后影像组学特征。

2 结果

       本研究结果依次展示脑MRI多尺度结构影像组学特征在不同场强、不同扫描参数设置时MRI设备之间的扫描仪效应及处理方法。首先，以T1WI图像上测量的皮层厚度和DTI图像上获得的FA值为例展示并说明MRI设备的“扫描仪效应”；然后，评估N4偏置场校正、HM归一化及ComBat协调方法在纠正扫描仪效应方面的表现；最后，评估ComBat协调前后脑灰白质影像组学特征的生物判别能力。

2.1 人口学资料和心理量表测试

       最终纳入受试者23例，一般人口学资料和心理量表测试结果如表3所示。男性和女性受试者的年龄（t=1.090，P=0.316）、受教育年限（t=-0.638，P=0.574）、CES-D评分（t=-0.670，P=0.510）差异均无统计学意义（P＞0.05）。

点击查看表格

表3  人口学资料和心理量表测试
Tab. 3  Demographic data and psychological test of subjects

2.2 影像特征扫描仪效应表征及校正

2.2.1 扫描仪效应表征

       图1展示了每例受试者所有脑区皮层厚度测量值，箱线图（图1A）显示MRI1（Philip 1.5 T）皮质厚度值整体向下移动，较MRI2（Philip 3.0 T）和MRI3（GE 3.0 T）测量值小，而且其变异性也增加。图1B显示的3个箱线图总结了每台设备中位皮质厚度的分布，ANOVA显示3台设备测量的皮质厚度中值不同（P＜1.4×10^-5）。采用Tukey-HSD（honestly significant difference test）检验进一步两两比较发现MRI1皮质厚度测量中值较另外2台测量值小，MRI2与MRI3测量值之间差异无统计学意义。采用Shapiro-Wilk检验评估皮层厚度测量值是否符合正态分布，P值均大于0.05（MRI1：P=0.731；MRI2：P=0.230；MRI3：P=0.182）。由于测量值的大小也可能受到不同扫描仪的影响，进一步Bartlett球形检验估计扫描仪之间的方差是否均匀（同质化）。为了明确并避免年龄和性别带来的效应，回归发现年龄和性别对测量值无影响（P=1.807）。

       图2展示了每例受试者所有48个白质区域FA值，箱线图（图2A）显示MRI1（Philip 1.5 T）FA值较MRI2（Philip 3.0 T）和MRI3（GE 3.0 T）测量值整体稍下移，而MRI3（GE 3.0 T）变异性较高。ANOVA和Tukey-HSD分析显示MRI1的FA值较另外2台值小，MRI2与MRI3之间差异无统计学意义（P＞0.05）（图2B）。Shapiro-Wilk检验P值均大于0.05（MRI1：P=0.794；MRI2：P=0.811；MRI3：P= 0.182）。Bartlett球形检验及年龄和性别回归结果与皮层厚度测量类似（P＞0.05）。

点击查看大图

图1  3台MRI扫描仪皮层厚度测量。1A：按设备分组的皮质厚度箱线图，每个箱线图代表1例受试者90个皮层脑区皮质厚度测量值的分布；1B：各设备中位皮质厚度箱线图；MRI1平均皮质厚度中位数较小，变异性较高。MRI1为Philip 1.5 T；MRI2为Philip 3.0 T；MRI3为GE 3.0 T。
Fig. 1  Cortical thickness for three MRI scanners. 1A: Boxplots of the cortical thickness sorted by scanner. Each boxplot represents the distribution of the 90 cortical regions for one subject; 1B: Boxplots of the median cortical thickness, grouped by scanner. The mean cortical thickness of MRI1 is small and has higher variability. MRI1: Philip 1.5 T; MRI2: Philip 3.0 T; MRI3: GE 3.0 T.

点击查看大图

图2  3台MRI扫描仪FA值测量。2A：按设备分组的FA值箱线图，每个箱线图代表1例受试者48个白质区域FA值的分布；2B：各设备平均FA值箱线图，MRI1平均FA值较小，而MRI3平均FA值变异性较高。FA为各向异性分数；MRI1为Philip 1.5 T；MRI2为Philip 3.0 T；MRI3为GE 3.0 T。
Fig. 2  FA values for three MRI scanners. 2A: Boxplots of the FA values sorted by scanner. Each boxplot represents the distribution of the 48 regions of white matter for one subject; 2B: Boxplots of the FA values, grouped by scanner. The mean FA value of MRI1 is small, while the mean FA value of MRI3 has higher variability. FA: fractional anisotropy; MRI1: Philip 1.5 T; MRI2: Philip 3.0 T; MRI3: GE 3.0 T.

2.2.2 扫描仪效应校正

2.2.2.1 偏置场校正

       从定性的角度分析，与原始图像相比，N4算法的默认参数设置（4个分辨率级别，每个级别50次迭代和使用由阈值定义的默认掩模）显示对于脑MR图像的偏置场校正不是最佳的。在50次迭代后，校正效果趋于稳定，增加迭代次数并不能改善校正效果。当使用5级（50次迭代）、全掩模进行N4偏置场校正时，脑组织区域的强度变异性最低。不同的校正参数设置时，N4算法对各设备8个白质VOI区域内强度均值的变异系数的影响如图3所示。

点击查看大图

图3  N4不同校正参数设置时各设备8个白质VOI区域强度均值变异系数。对于T1WI和DTI，MRI1的图像区域强度均值变异系数较大；与原始图像相比，N4算法的默认参数显示对于脑MR图像的偏置场校正并不是最佳的设置；当使用5级分辨率、全掩模进行N4偏置场校正时，脑组织的区域强度变异性最低。VOI为感兴趣体积；T1WI为T1加权成像；DTI为扩散张量成像。
Fig. 3  The mean coefficient of variation over eight VOI of the white matter across different N4 corrections settings. For T1WI and DTI, the mean coefficient of variation of region intensity for MRI1 is large; Compared with the original images, the default parameter of N4 algorithm for the bias field correction of brain MR is not the best setting; When N4 bias field correction was performed using 5-level (50 iterations) full mask, the intensity variation coefficient of brain tissue among the three scanners was the lowest. VOI: volume of interest; T1WI: T1 weighted image; DTI: diffusion tensor imaging.

2.2.2.2 MRI结构图像HM归一化及影像组学特征ComBat协调

       图像后处理流程中各环节不同的后处理方法对于T1WI和DTI两个序列的效果表现相似。在校正之前的原始图像中，强度直方图的分布范围及峰值均未对齐，尤其是MRI1（Philip 1.5 T），强度分布范围明显大于来自MRI2和MRI3设备的强度。N4校正锐化了强度峰值，但各设备的强度峰型并未对齐。HM归一化与N4校正相结合，进一步优化对齐了各强度峰，使峰值对齐在相同值附近。

       为了显示图像处理流程中各环节对影像组学特征值的影响，图4显示了从3个设备的T1WI图像上提取的灰度级长矩阵纹理特征-短行程低灰度强调（GLRLM_SRLGE）特征分别在3个处理阶段的统计分布。图4表明HM归一化有助于重新对齐各设备特征的分布（图4C），但需要ComBat进一步协调以对齐所有三个设备的分布（图4D）。

       表4展示了Kruskal-Wallis检验P值在特定范围内的影像组学特征数量。例如，在T1WI原始图像中，63个特征中有56个特征在3台设备间的差异具有统计学意义，其中有43个特征P＜10^-5，5个特征的P值范围在10^-5≤P＜10^-3，6个特征的P值范围在10^-3≤P＜0.01，2个特征的P值范围在0.01≤P＜0.05，仅有7个特征的P值大于0.05。在T1WI和DTI序列中观察到相同的趋势：当从N4偏置场校正数据、HM归一化N4校正数据、ComBat协调后数据计算时，具有差异的特征数量逐渐减少。ComBat协调对于大幅减少3台设备具有差异的特征数量至关重要。

点击查看大图

图4  各设备T1WI图像提取的GLRLM_SRLGE纹理特征强度的统计分布。4A：原始T1WI图像中强度分布；4B：N4（全掩模、5级）偏置场校正后强度分布；4C：HM归一化N4校正的强度分布；4D：HM归一化N4校正并进一步由ComBat协调的强度分布。T1WI为T1加权成像；HM为直方图匹配；GLRLM_SRLGE为灰度级长矩阵纹理特征-短行程低灰度强调。
Fig. 4  Statistical distributions across scanners of the GLRLM_SRLGE feature extracted from T1WI. 4A: Intensity distribution in raw T1WI; 4B: Intensity distribution in N4 corrected (full mask, 5 levels) T1WI; 4C: Intensity distribution in HM normalized N4 corrected T1WI; 4D: Intensity distribution in HM normalized-N4 corrected T1WI and harmonized by ComBat. T1WI: T1 weighted image; HM: histogram matching; GLRLM_SRLGE: Gray-Level Run Length Matrix_Short Run Low Gray Level Emphasis.

点击查看表格

表4  5个不同P值范围内的MRI影像组学特征数量
Tab. 4  Number of radiomic features in the 5 different ranges of P-values

2.3 脑MRI结构影像组学特征鉴别脑灰白质

       对于两种脑组织类型灰质（VOIc2、VOIc3）和白质（VOIc1、VOIc4），对每个序列及每个特征进行Wilcoxon秩和检验。协调前，由T1WI及DTI图像上提取的79个特征中分别有29个和5个特征在两种脑组织类型之间存在差异（P＜0.05），而在ComBat协调后，分别有52个和12个特征在两种脑组织类型之间存在差异（P＜0.05）。图5以一阶直方图平均绝对偏差（mean absolute deviation of histogram, HISTO_MAD）特征为例，显示了ComBat对该特征值平均强度的影响，降低了扫描仪间的影响。

点击查看大图

图5  ComBat协调前后T1WI序列脑灰质和白质一阶直方图HISTO_MAD特征的平均强度比较。5A：ComBat协调前HISTO_MAD特征的平均强度；5B：ComBat协调后HISTO_MAD特征的平均强度；ComBat协调降低了脑灰白质HISTO_MAD特征扫描仪效应并保留其生物性差别。HISTO_MAD为一阶直方图平均绝对偏差；VOI为感兴趣体积；VOIc1为右侧半卵圆中心；VOIc2为右侧尾状核头；VOIc3为左侧背侧丘脑；VOIc4为左侧颞叶深部白质。
Fig. 5  Mean intensity of the first order statistics HISTO_MAD feature in gray and white matter before and after ComBat harmonization across scanners. 5A: The mean intensity of HISTO_MAD feature before ComBat; 5B: The mean intensity of HISTO_MAD feature after ComBat. ComBat harmonization on the mean intensity reduces scanner effect and retains its biological differences. HISTO_MAD: Histogram_ Mean Absolute Deviation; VOI: volume of interest; VOIc1: Right centrum ovale; VOIc2: Right caudate head; VOIc3: Left Dorsal thalamus; VOIc4: White matter of the left temporal lobe.

3 讨论

       本研究通过分析3台扫描仪获得的脑T1WI和DTI数据，发现“扫描仪效应”对MRI数据的统计结果确实存在影响，同时还发现既往为乳腺体模^[7]或脑肿瘤^{[9, 10]}MRI开发的校准方法并不完全适用于无实质性占位的脑MRI结构图像。本研究证实并推荐适合于脑MRI结构影像的一般校准过程需包括三个步骤，N4偏置场校正可以减少相似区域中的局部/区域不均匀性（图像内变化），HM强度归一化可以减少采集间的变化，以及进一步ComBat协调以使各扫描仪之间的结果具有可比性。每一步骤都能够有效解决不同来源的可变性，三者是必要的和互补的；且在校准非生物因素的同时，保留了脑灰白质之间的生物差异性。本流程的提出为多扫描仪采集的影像组学研究的高效标准化提供了帮助。

3.1 脑MRI结构影像N4偏置场校正

       本研究通过比较不同参数设置时N4偏置场估计，得出一组适用于无器质性病变的脑MRI结构图像偏置场校正的最佳参数。根据本研究结果，建议使用包含全掩模（与阈值定义的默认掩模不同）并在五个分辨率级别而不是四个上执行优化，每个级别进行50次迭代。使用为脑肿瘤^{[9, 10]}MRI优化的默认参数低估了偏置场的变化，导致强度不均匀。默认掩模和五个分辨率级别的缺点是在中心部分低估了偏置场，在校正后的图像中会产生超信号效应，增加了背景信号强度异质性，因此增加了平均强度的变异系数。通过定性分析发现，所提出的脑MRI结构影像最佳N4参数设置会最大程度降低不均匀区域的平均强度变异系数。在N4校正图像上，组织内部的强度直方图显示其峰值在最大结构的平均值附近显著锐化。需要注意的是估计的偏置场不仅取决于MR扫描仪，还取决于线圈参数、序列类型以及靶部位在线圈内的定位，所以即使采用相同的线圈，偏置场也会跨序列出现明显的波动^[4]。

3.2 脑MRI结构影像HM归一化

       为了减少受试者间的差异，需要在偏置场校正后对MRI影像进行归一化^[27]。本研究中跨采集和线圈的强度直方图的共同对齐突出了强度归一化的影响，结果也证实直方图匹配的归一化方法的可行性及高效性。既往研究证实混合条纹法确实压缩了一系列值内的所有强度，但未能对齐组织特异性峰^[28]。另有研究结果表明^[29]，在前列腺影像组学特征归一化中，分段线性直方图在重新调整强度分布方面具有显著而可靠的效果。然而，直方图匹配取决于选择用于提取标准直方图的图像集。在临床环境中，重要的是使用来自各种扫描仪和不同生物变异性的患者影像来识别稳健的指标，同时在新患者的影像中也能够稳定表现^[24]。

3.3 脑MRI结构影像组学特征ComBat协调

       本研究结果表明，结合直方图匹配归一化的N4校正并不能完全消除扫描仪效应。虽然上述校准环节都不同程度减少了3台MRI设备之间有差异的特征数量，但仍然难以协调所有影像组学特征。这与既往在脑肿瘤^[9]和前列腺癌^{[10, 29]}影像组学研究中观察到的趋势一致。因此，本研究除了对MRI扫描仪内的成像数据进行协调，另外关注的重点是跨扫描仪的影像数据的校准。事实上，即使对于在同一扫描仪上获得的影像，参与者之间由技术参数和物理环境引起的非生物变化也存在，这需要一种不依赖于站点或扫描仪变量的协调技术。

       ComBat技术^[4]可以通过估计来自控制区域的不需要变化的潜在因素，且可以扩展到站点内协调，其优势是能够扩大规模以进行大型神经影像学研究^{[9, 30]}。事实上，ComBat算法与成像特征的数量呈线性关系，因此适用于在数量通常可以达到数百万级别的体素水平执行图像分析。在分类结果的情况下，例如疾病状态，可以仅使用来自健康人群的参与者来估计仪器或站点效应参数。这种方法对于将临床研究与高异质性疾病相结合特别有用，例如对于自闭症和创伤性脑损伤的诊断和鉴别^{[7, 22, 30]}。但是，对于小样本量研究，区分疾病异质效应和部位效应可能比较困难。因此，可通过使用相对更稳定的健康人群来标准化数据以提高性能。

       既往研究对每个站点或中心分别使用缩放或Z分数对特征进行归一化^[31]，但与ComBat不同的是，这些方法无法对可能影响特征的可能协变量进行建模。SAINT等^[7]使用2个扫描仪和3个乳房专用线圈采集2个乳腺体模的T1WI和T2WI图像，校正了图像偏置场的不均匀性，并使用Z分数或直方图匹配进一步归一化，然后提取影像组学特征并通过ComBat方法在6个线圈之间进行协调。通过利用模拟正常乳腺组织及致密病变的凝胶中计算的两组影像特征值比较，对整个校准流程进行定性和定量评估。研究结果显示校准流程明显降低了采集内和采集间的差异，ComBat协调将6个线圈之间有差异的影像特征的百分比从偏置场校正和归一化后的87%降低到3%，同时提高了影像组学特征对于2种病变分类的性能。该流程的实现为减少乳腺MRI影像数据的变异性，并进行稳健的多扫描仪影像学研究奠定了基础，但作者仍强调需要根据患者检查部位及具体数据模态进行评估，不能盲目直接采用。虽然ComBat在降低扫描仪效应方面具有重要作用，但也需要N4校正和HM归一化来减少无法通过ComBat仿射变换准确补偿的图像内和采集间变化^{[7, 10]}。因此，组合校正对于校正不同来源的变异性至关重要。

3.4 本研究的局限性

       首先，本研究脑MRI结构影像数据采集是在同一机构分别使用3台MRI扫描仪进行的，故需要进一步验证该校准框架扩展到不同站点及中心影像数据的效果。其次，本试验中仅纳入23名健康志愿者，样本量小，且仅采集了结构影像数据，未验证其在脑功能磁共振出血等与时序有关的影像数据校准中的价值。

4 结论

       总之，本研究结果表明，汇总多台设备和多个序列参数设置所获取的脑MRI结构影像组学研究是可行的，且有必要对采集的数据进行校准。利用N4偏置场校正、HM归一化与ComBat协调的联合校准方法，能够有效减少脑MRI结构图像的“扫描仪效应”，同时保留与生物学相关的变异性，从而确保跨扫描仪、多中心MRI影像组学研究的可重复性和准确性。