0引言肺脏是人体呼吸系统的关键器官，其主要功能是将空气中的氧气输送到血液，并将血液中的二氧化碳排出至体外。此过程中，肺血管是血液运输和气体交换的媒介。肺血管的分割是肺部结构可视化和疾病定量分析的重要前提。肺血管重塑是慢性阻塞性肺病(慢阻肺)(chronic obstructive pulmonary disease，COPD) 的重要病理学特征，其主要表现有近端血管增粗和远端血管细化消失(Rahaghi等，2019)，如图 1所示。在临床医学研究中，不同半径的肺血管体积占比常用于定量评估肺血管重塑(Matsuoka等，2010；Estépar等，2013)，上述定量评估通常需要人工介入，且很难覆盖所有感兴趣区域。肺血管的自动分割及定量分析，可为医生提供客观评价数据，辅助医生对慢阻肺疾病进行更精准的诊断与分级，具有临床研究与应用价值。 图1 慢阻肺血管重塑表现 Typical manifestations of chronic obstructive pulmonary vascular remodelingFig 1 (箭头所示远端血管变细或消失) (arrows denoting narrowing and pruning of distal vessels)((a) coronal slices of CT scans from COPD patients; (b) maximal intensity projection images of (a)) 然而，在肺血管邻近组织、复杂树状拓扑结构、肺实质病变和成像噪声等因素的影响下，计算机断层扫描(computed tomography, CT)图像中肺血管的自动分割充满挑战。在过去几十年中，许多肺血管分割方法相继提出(边子健等，2018；Moccia等，2018)，不同肺血管分割方法的性能进行了详细比较和分析(边子健等，2018)，Moccia等人(2018)总结了机器学习方法在血管分割方面的应用。本文将血管分割方法划分为3类：基于增长、基于机器学习和基于导数滤波的管树结构分割方法。1) 基于增长的分割方法。CT图像中的肺血管分割可采用基于原始灰度或管状似然强度的区域增长方法(Bülow等，2005)。从肺血管树根节点处的种子点开始，以阈值下限为增长规则，重复合并或以快速行进法(Sethian，2001)演化的波阵面形式扩散种子点标记，形成肺血管树。此外，在血管增强的基础上，Lo等人(2010)、Zhou等人(2012)和Nadeem等人(2021)提出基于血管似然强度和血管轴向方向约束的局部最优路径增长方法来分割肺血管树。2) 基于机器学习的分割方法。机器学习方法主要分为基于手工特征的分类器方法和深度学习方法。Ochs等人(2007)以Hessian矩阵的特征值为特征输入，利用AdaBoost训练得到肺结构分类器。由于此方法采用人工特征，其对数据鲁棒性有限。Kiros等人(2014)使用字典稀疏编码学习方法，通过正交匹配追踪(orthogonal matching pursuit, OMP)算法获得逻辑回归分类器实现肺血管分割。OMP算法需要事先预设稀疏信号的稀疏性，这可能会导致训练过程困难。Zhao等人(2017)提出一种基于稀疏自动编码特征的方法。此方法基于2维切片进行训练，学习3维肺血管结构信息的能力受限。卷积神经网络(convolution neurual network，CNN)具有强大的特征提取与分类能力，可直接以原始图像为输入自动生成待分割目标。基于CNN的深度学习方法已成功应用于医学图像中各种器官和组织的分割(Shen等，2017)。由于肺血管数据标注难度很大，标注成本高，且无公开的肺血管的标注数据集，因此，针对肺血管分割的网络不多。Nardelli等人(2018)率先使用CNN对肺动脉和肺静脉分类，并且利用图分割方法对分类结果进一步优化，剔除空间结构不一致目标。Qin等人(2021)提出一种多任务分割框架，可同时分割气管、肺动脉和肺静脉，此框架融合肺解剖结构空间信息并引入注意力机制增强模型对管树结构的学习表征能力。这是目前公开报道的最新端到端肺血管分割网络，但是尚未在慢阻肺数据集上进行验证。3) 基于导数滤波器的分割方法。最常用的一类血管检测方法是以二阶图像导数(Hessian矩阵)为内核的管状或线状似然滤波器增强方法(Frangi等，1999)，它们假设血管为圆柱形，截面图像灰度服从高斯分布。通过分析和组合Hessian矩阵的特征向量与特征值，获得血管轴向方向和血管似然强度。为了适应管状半径变化，普遍使用归一化的多尺度滤波策略检测管状目标。这类方法使用高斯线性多尺度空间(Lindeberg，1998)，其本质是各向同性扩散方程的解，不可避免地模糊边缘，导致相邻目标相互扩散与融合。有部分研究工作致力于改善上述问题。非线性扩散(Perona和Malik，1990)和保边的各向异性扩散方法(Manniesing和Niessen，2005)可以有效阻止模糊强边缘，但是，容易从弱边缘处发生泄漏。最优方向通量(optimally oriented flux, OOF)(Law和Chung，2008)滤波器仅依赖边缘梯度信息，最大化地限制邻近目标的影响，但没有考虑管状内部灰度的均匀性。双高斯滤波器(Xiao等，2013)在平滑管状目标内部的同时，通过设置目标前景和背景尺度比例参数适当限制邻近目标的干扰。然而，此方法需预先知晓前景和背景的尺度参数。最近，受扩散张量各向异性测量方法的启发，一种基于Hessian矩阵各向异性测量的血管增强方法被提出(Jerman等，2016)，能有效克服血管分叉处响应低和因图像对比度变化引起响应不一致的问题。本文提出一种将导数滤波增强结果融入各向异性正则的连续最大流分割框架的肺血管自动分割方法。肺血管似然强度和轴向方向通过多尺度Hessian矩阵滤波器获得，强度和方向信息分别以数据保真项和各向异性正则项融入到连续最大流模型中。各向异性正则的应用有利于血管边缘沿轴向延生增长，以适应细长结构的分割。本文工作的主要创新点包括：1) 将肺血管轴向信息以各向异性总变分(anisotropic total variation)的形式融入到连续最大流的分割框架，提高细小血管的分割能力。2) 在低剂量CT数据集下进行了分割性能测试，并在多中心临床数据集上(正常剂量)对分割的血管结果进行量化分析，验证了分割算法在不同剂量数据下的鲁棒性。3) 基于血管分割结果，自动提取了全肺的小血管体积占比参数，并在慢阻肺病人的肺血管重塑的定量化分析中得到验证。1方法1.1整体流程本文提出的肺血管分割系统整体流程如图 2所示，主要包括肺体分割、肺血管增强和肺血管分割3个步骤。肺体分割采用基于U-Net网络结构的深度学习方法(Hofmanninger等，2020)；肺血管增强采用基于多尺度Hessian矩阵各向异性度量的增强方法(Jerman等，2016)；肺血管分割使用全局最优的各向异性连续最大流分割方法(Pezold等，2016)。 图2 肺血管分割方法整体流程图 Overview of the proposed pulmonary vasculature segmentation approachFig 21.2基于U-Net的肺体分割肺体分割是肺结构分析的重要步骤，主要为肺血管的后处理提供感兴趣区域。迄今为止，无论是传统的阈值分割方法，还是如今流行的深度学习方法，能适应各种临床肺病病灶(如肺纤维化、毛玻璃影等)干扰的肺体分割方法很少。Hofmanninger等人(2020)指出对于基于深度学习的肺体分割方法，提高分割性能的主要措施在于提高训练数据集的多样性，而不是设计复杂的网络结构和训练策略，原因是肺体的图像对比度高和形状简单。Hofmanninger等人(2020)采用经典的U-Net分割网络，给每一卷积层加上批归一化(batch normalization，BN)，在231例包含多种病例类型和成像参数的CT影像数据(共计62 224个切片)上完成模型训练，所得模型在其测试集上获得了最优分割性能。本文采用此分割模型完成左右肺体分割，对于个别因肺体边界病变引起局部图像灰度过高导致分割出错的数据，在影像科医生(15年以上从业经验)的指导下手动纠正分割结果。1.3多尺度Hessian矩阵肺血管增强大多数管状结构增强滤波器依赖图像灰度的二阶导数来表征结构特征，其潜在假设为：在暗背景下，血管呈明亮、细长的圆柱状结构。基于此假设，众多滤波器计算图像的Hessian矩阵。在图像中的每一体素表示为函数，Hessian矩阵定义为 1 $\begin{array}{c}\boldsymbol{H}(f)=\left[\begin{array}{c}h_{11} & h_{12} & h_{13} \\h_{21} & h_{22} & h_{23} \\h_{31} & h_{32} & h_{33}\end{array}\right]=\\\left[\begin{array}{ccc}\frac{\partial^{2} f}{\partial x_{1}^{2}} & \frac{\partial^{2} f}{\partial x_{1} x_{2}} & \frac{\partial^{2} f}{\partial x_{1} x_{3}} \\\frac{\partial^{2} f}{\partial x_{2} x_{1}} & \frac{\partial^{2} f}{\partial x_{2}^{2}} & \frac{\partial^{2} f}{\partial x_{2} x_{3}} \\\frac{\partial^{2} f}{\partial x_{3} x_{1}} & \frac{\partial^{2} f}{\partial x_{3} x_{2}} & \frac{\partial^{2} f}{\partial x_{3}^{2}}\end{array}\right]\end{array}$ 式中, $h$表示图像的二阶导数，$x_{1}$、$x_{2}$、$x_{3}$表示CT图像的3个维度方向。为了使$f$满足二阶可导的要求，需将原始图像$\boldsymbol{I}$与高斯核进行卷积。原始图像与一系列不同尺度标准差$σ$的高斯核卷积后可建立起多尺度空间。令$\boldsymbol{e}_{1}, \boldsymbol{e}_{2}, \boldsymbol{e}_{3}$为$\boldsymbol{H}(f)$的特征向量，与之对应的特征值为$λ_{1}, λ_{2}, λ_{3}$，且$|λ_{1}|≤|λ_{2}|≤|λ_{3}|$。血管局部可以看做为暗背景下的明亮圆柱状结构，其特征值满足以下特性 2 $\left|\lambda_{1}\right| \approx 0, \lambda_{2} \approx \lambda_{3} \ll 0$ Jerman等人(2016)提出一种基于Hessian矩阵各向异性度量的血管似然增强方法，此方法在血管分叉处可获得正常响应，且对血管截面形状不敏感。此血管增强函数定义为 3 $V= \begin{cases}0 & \lambda_{2} \leqslant 0 \text { 或 } \lambda_{\rho} \leqslant 0 \\ 1 & \lambda_{2} \geqslant \frac{\lambda_{\rho}}{2}0 \\ \lambda_{2}^{2}\left(\lambda_{\rho}-\lambda_{2}\right)\left(\frac{3}{\lambda_{2}+\lambda_{\rho}}\right)^{3} & \text { 其他 }\end{cases}$ 式中，$λ_{ρ}$是受约束的$λ_{3}$，其作用是控制低对比度区域噪声的敏感性，其数学定义为 4 $\lambda_{\rho}= \begin{cases}\lambda_{3} & \lambda_{3}\tau \max _{x} \lambda_{3}(x) \\ \tau \max _{x} \lambda_{3}(x) & 0\lambda_{3} \leqslant \tau \max _{x} \lambda_{3}(x) \\ 0 & \text { 其他 }\end{cases}$ 式中，$τ∈[0, 1]$。从图 3结果中可以看出，非肺血管结构得到了抑制，而肺血管区域得到了增强，有利于提高血管分割的鲁棒性。 图3 肺血管增强示例 Demonstration of pulmonary vessel enhancementFig 3 ((a) original slices of a CT scan; (b) enhancement results of CT slices in (a)) 1.4各向异性最大流肺血管分割由于直接对增强结果进行阈值处理获得的肺血管分割结果不够精准，采用全局最优的连续最大流分割模型进行分割(Yuan等，2010)，并将血管轴向信息以各向异性总变分形式融入其中。3维图像空间$\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} \subset {\bf{R}}^{3}$的二元分割问题可以看做是一个二分类标记问题$u:\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}→{0, 1}$。连续最大流分割是一种与最小割互为对偶的能量最小化模型(Yuan等，2010)。二元分割的最小割问题定义为 5 $\min \limits_{u \in[0,1]} \int_{\boldsymbol \varOmega}(1-u) C_{\mathrm{s}}+u C_{\mathrm{t}}+|\boldsymbol{\nabla} u| C \mathrm{d} x$ 式中, $C_{\rm{s}}$为源节点容量，$C_{\rm{t}}$为汇节点容量，$C$为非端点容量，分别定义为 6 $C_{\mathrm{s}} =\frac{1}{q} \max \left\{\ln \left(\frac{p_{\mathrm{f}}^{v}+\varepsilon}{p_{\mathrm{b}}^{v}+\varepsilon}\right), 0\right\}$ 7 $C_{\mathrm{t}} =\frac{1}{q} \max \left\{-\ln \left(\frac{p_{\mathrm{f}}^{v}+\varepsilon}{p_{\mathrm{b}}^{v}+\varepsilon}\right), 0\right\}$ 8 $q =\ln \left(\max \left\{\hat{p}{}_{\mathrm{f}}^{v}, \hat{p}{}_{\mathrm{b}}^{v}\right\}+\varepsilon\right)$ 9 $C =w \exp \left(-1 / \zeta^{2}|\boldsymbol{\nabla} \boldsymbol{I}|\right)$ 式中，$p^{v}_{\rm{f}}$和$p^{v}_{\rm{b}}$分别为肺血管前景和背景似然概率，${\hat p}={\rm{max}} \{p(x)\}$，$ε$为极小正值，避免零对数和分母为0, $w$为权重参数，$ζ$控制梯度幅值响应的灵敏度。$\int_{\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}}|\nabla u|{\rm{d}}x$为常见的总变分项，其作用是减小分割区域曲面面积，惩罚前景和背景之间的边界跳跃，从而保持分割区域的光滑性。对于各向异性的线状结构分割来说，它并非最优正则项，原因是总变分是各向同性约束，没有考虑分割函数$u$变化的方向性。例如分割管状目标，利用管状结构的形状先验知识，即沿轴向方向无变化，因此，设计沿轴向变化的惩罚力度大于沿垂直方向变化的正则项。换言之，选择更好的各向异性正则项，其定义为 10 $\begin{gathered}A T V(u, \boldsymbol{A})=\int_{\boldsymbol{\varOmega}}\left(\boldsymbol{\nabla} u^{\mathrm{T}} \boldsymbol{A} \boldsymbol{\nabla} u\right)^{1 / 2} \mathrm{~d} x=\int_{\boldsymbol{\varOmega}}\left|\boldsymbol{S}^{\mathrm{T}} \boldsymbol{\nabla} u\right| \mathrm{d} x \\\boldsymbol{A}=\boldsymbol{S} \boldsymbol{S}^{\mathrm{T}}\end{gathered}$ 式中，$\boldsymbol{A}$为正定矩阵，$\boldsymbol{S}$为$\boldsymbol{A}$的矩阵分解分量。将Hessian增强步骤中获得的特征向量$\boldsymbol{e}_{1}, \boldsymbol{e}_{2}, \boldsymbol{e}_{3}$组成正交基$\boldsymbol{V}=[\boldsymbol{e}_{1};\boldsymbol{e}_{2};\boldsymbol{e}_{3}]$，由此可计算出$\boldsymbol{A}=\mathit{\boldsymbol{V\tilde AV}}^{\rm{T}}，\boldsymbol{S}=\boldsymbol{V\tilde A}^{1/2}$，其中$\boldsymbol{\tilde A}={\rm{diag}}(1, a, a)$，$0a≤1$，当$a=1$时，则退化为各向同性总变分正则项。因而各向异性正则化的最小割问题的数学定义表示为 11 $\min \limits_{u \in[0,1]} \int_{\boldsymbol{\varOmega}}(1-u) C_{s}+u C_{t}+\left|\boldsymbol{S}^{\mathrm{T}} \boldsymbol{\nabla} u\right| C \mathrm{d} x$ 最终其对偶问题在增加拉格朗日乘子项后的表达式为(Pezold等，2016) 12 $\begin{gathered}\max \limits_{p_{\mathrm{s}}, p_{\mathrm{t}}, p} \min \limits_{u \in[0,1]} \int_{\boldsymbol{\varOmega}} p_{\mathrm{s}}+u\left({div}(\boldsymbol{S} p)-p_{\mathrm{s}}+p_{\mathrm{t}}\right)-\\\frac{c}{2}\left({div}(\boldsymbol{S} p)-p_{\mathrm{s}}+p_{\mathrm{t}}\right)^{2} \mathrm{~d} x\end{gathered}$ 采用增广拉格朗日乘子算法求解，具体算法流程如算法1所示。算法1：基于增广拉格朗日乘子的各向异性最大流求解算法1) 设置约束变量${\hat \varepsilon }$，步长$γ，c$；计算$C, C_{s}, C_{t}, \boldsymbol{S}$; 初始化$p^{0}, p^{0}_{\rm{s}}, p^{0}_{\rm{t}}, u^{0}$；2) 从$n=0$开始迭代，直到$\frac{1}{{\left| \mathit{\boldsymbol{ \boldsymbol{\varOmega} }} \right|}}\int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {\left| {{\varepsilon ^{n + 1}}(x)} \right|} {\rm{d}}x \hat \varepsilon $；3) while $\frac{1}{{\left| \mathit{\boldsymbol{ \boldsymbol{\varOmega} }} \right|}}\int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {\left| {{\varepsilon ^{n + 1}}(x)} \right|} {\rm{d}}x \hat \varepsilon $ do4) $\tilde{p}^{n+1}=p^{n}+\gamma \boldsymbol{S}^{\mathrm{T}} \boldsymbol{\nabla}\left(\operatorname{div}\left(\boldsymbol{S}p^{n}\right)-p_{\mathrm{s}}^{n}+p_{\mathrm{t}}^{n}-\frac{u^{n}}{c}\right) $；5) $p^{n+1}=\frac{\tilde{p}^{n+1}}{\left|\tilde{p}^{n+1}\right|} \min \left\{\left|\tilde{p}^{n+1}\right|, C\right\} \text { if } \tilde{p}^{n+1} \neq 0 \text { else } 0 $；6) $p_{\mathrm{s}}^{n+1}=\min \left\{\frac{1-u^{n}}{c}+\operatorname{div}\left(\boldsymbol{S}_{p}{ }^{n+1}\right)+p_{\mathrm{t}}^{n}, C_{\mathrm{s}}\right\} $；7) $p_{\mathrm{t}}^{n+1}=\min \left\{\frac{u^{n}}{c}-\operatorname{div}\left(\boldsymbol{S} p^{n+1}\right)+p_{\mathrm{s}}^{n}, C_{\mathrm{t}}\right\} $；8) $ε^{n+1}=c(div(\boldsymbol{S}p^{n+1})-p^{n+1}_{\rm{s}}+p^{n+1}_{\rm{t}}$；9) $u^{n+1}=u^{n}-ε^{n+1}$；10) end while1.5肺血管体积量化分析为了量化评估肺血管重塑程度，如小血管变细或消失、大血管增粗等情形，需要分析血管管径及其体积分布。为了获得肺血管的截面半径，应用距离变换获得肺血管中心线(Selle等，2002)，并给中心线上的每个体素点赋予血管半径值。血管体积可以通过中心线长度和半径值测算得到。血管半径为$r$的血管体积计算公式为：$CSA_{r}=N_{r}×2×{\rm{π}}×r^{2}$，其中$N_{r}$表示中心线数据中半径为$r$的体素数量。本文主要测量截面面积小于5 mm2或10 mm2的血管体积CSA5或CSA10与肺血管总体积，以及小血管体积与肺血管总体积的比值。2实验与结果2.1实验数据在ArteryVein (Charbonnier等，2016)和VascuSynth(Hamarneh和Jassi，2010)公开数据集上对所提出的分割算法进行有效性验证。ArteryVein是一个公开的低剂量CT扫描数据集(管电流30 mAs，管电压120~140 kV)，包含不同严重程度的肺气肿和间质性肺疾病影像。所有影像数据是在全吸气状态下通过CT扫描获得，其横断面分辨率为512×512像素，空间分辨率为0.59~0.83 mm，切片间距均为0.7 mm。ArteryVein数据集含有55例扫描数据，其中对肺血管进行了全标注的有10例，另外45例只进行了部分标注。本文使用肺血管全标注的影像。VascuSynth是一个各向同性的血管仿真数据集，其分辨率为100×100×100，图像灰度分布范围为[0，255]。此数据集共有10组数据，每组数据由12例随机生成包含1~56个分叉点的血管树组成。对每组仿真数据添加不同标准差的高斯噪声，以验证算法的适应性。为了分析慢阻肺病人的肺血管重塑，从4家医院(山东省立医院、长沙市第一人民医院、昆明军区总医院和青岛大学附属医院)收集了1 284例多中心临床肺部CT扫描数据，并对数据进行了匿名化处理，成像设备包含GE、Philips、SIEMENS、TOSHIBA，图像横断面分辨率为512×512像素，层厚包含5.0 mm、1.5 mm、1.25 mm、1.0 mm和0.625 mm。为了适应分割算法以及满足小血管体积测量的需求，将厚层数据(≥1.5 mm)剔除，余下614例，如表 1所示，其中慢阻肺病例352例，非慢阻肺病例262例。慢阻肺组中包含GOLD分级信息的数据有281例，1~4级分别为16、108、108、49例。本文对上述614例临床数据进行肺血管分割及体积占比测量，分析了慢阻肺组和非慢阻肺组之间，以及不同GOLD分级病例之间的肺血管体积占比情况。 表1 实验中不同中心慢阻肺与非慢阻肺临床影像数据使用情况 医院 COPD NonCOPD 筛选条件 长沙市第一人民医院 137 93 Th≤1.5 mm 昆明军区总医院 206 8 Th≤1.5 mm 山东省立医院 9 100 Th≤1.5 mm 青岛大学附属医院 0 61 Th≤1.5 mm Clinical CT scans of COPD and NonCOPD patients from different centers in the experimentTable 1 Th表示层厚。2.2实验参数及平台本文算法参数主要包含肺血管增强算法和分割算法参数。增强算法中的噪声抑制参数$τ=0.5$；尺度参数依据血管管径范围选择(Jerman等，2016)，其中最小尺度$σ_{\rm{min}}=0.5 {\rm{mm}}$，最大尺度$σ_{\rm{max}}=3.5 {\rm{mm}}$，总共离散为15个尺度级别。对于仿真数据实验，分割算法中各向异性正则矩阵中的$a=0.03$，非端点流量中权重$w=0.7$，敏感度$ζ=0.01$，均在部分标注数据上通过网格搜索获得最优值；对于ArteryVein数据实验，获得的最优参数为：$a=0.03$，$w=0.3$，$ζ=0.04$。实验所用的肺血管增强和分割算法均基于并行框架OpenCL实现，实验主机硬件的主要配置为Intel Xeon CPU 8核，RAM 64 GB，NVIDIA Quadro P5000，配备Windows 7操作系统。2.3分割实验结果选用Dice系数作为肺血管分割性能度量指标，其定义为 13 ${ Dice }=\frac{2 \times T p}{F p+F n+2 \times T p}$ 式中，$Tp$表示分割结果与真实血管均包含的体素，$Fp$表示分割结果含有但真实血管不含的血管体素，$Fn$表示分割结果缺失而真实血管含有的血管体素。此度量指标主要衡量分割结果与真实血管之间的重叠程度。鉴于低剂量CT图像密度值噪声的标准差范围为[8, 29] (江一峰等，2010)，对VascuSynth血管仿真数据分别添加标准差$σ$为5、10、15、20、25、30、35的高斯噪声后再进行血管分割。通过对10组数据进行分割性能测试后，所统计的平均Dice系数走势如图 4中曲线所示，$σ=10$时，$Dice=0.8$；$σ=25$时，$Dice=0.73$；$σ=35$时，$Dice=0.69$。图 4曲线表明尽管随着噪声标准差的增加Dice系数降低，但在CT最大噪声干扰条件下，本文血管分割方法的性能保持在$Dice=0.73$附近。 图4 高斯噪声与分割性能关系曲线 The relationship between Gaussian noises and the vasculature segmentation performanceFig 4除仿真数据外，还对ArteryVein低剂量数据集中含有肺血管全标注的CT图像进行了分割性能测试。通过实验从定性和定量两个方面比较了本文各向异性正则化分割方法与各向同性正则化分割方法的性能差异。各向同性正则分割方法可由本文方法通过设置各向异性参数$a=1$退化得到。图 5展示了两种方法对远端细支血管分割性能差异。图 5(a)为各向同性正则分割结果，从黄色圈标记的区域可以看出，末端细支血管已丢失；而从各向异性正则分割方法的分割结果图 5(b)可以看出，这些细支血管得到了保留。对10例CT图像的分割结果以Dice指标进行了量化度量评估，各向异性正则分割方法与各向同性正则分割方法的平均Dice系数分别为0.79和0.76。本文肺血管分割方法的分割精度接近细小血管手动分割的精度(Haft-Javaherian等，2020；Maher等，2020)。值得说明的是基于深度学习的肺血管分割方法(Qin等，2021)在此数据集上获得了Dice为0.83的分割性能(如表 2)，但是其Dice方差值为0.08。 图5 不同正则化分割方法的肺血管分割效果 Comparison of the isotropic and anisotropic regularization based segmentation methodsFig 5((a) segmentation results of isotropic method; (b) segmentation results of anisotropic method) 表2 本文方法与深度学习方法(Qin等，2021)在ArteryVein数据集上的肺血管分割性能比较 测试数据 本文方法 Qin等人(2021) 第1组(1396) 0.77 0.86 第2组(2064) 0.78 0.87 第3组(2604) 0.80 0.89 第4组(3300) 0.81 0.88 第5组(3356) 0.78 0.63 第6组(3476) 0.78 0.70 第7组(3840) 0.78 0.87 第8组(4332) 0.79 0.84 第9组(4584) 0.83 0.88 第10组(4840) 0.80 0.89 平均值 0.79 0.83 方差 0.02 0.08 Comparison of Dice for the proposed approach and a deep learning method (Qin et al., 2021) on ArteryVein datasetTable 2 测试数据栏括号内容为公开测试数据的编号。2.4肺血管体积定量评估结果为了观察慢阻肺病人的肺血管重塑情况，收集了慢阻肺和非慢阻肺病例的临床CT影像和相关临床信息，并对两个实验组CT影像的肺血管分割结果进行定量分析，分别统计小血管占比情况。本文对截面面积小于5 mm2(CSA5)和截面面积小于10 mm2(CSA10)的血管体积进行统计，并计算小血管体积占血管总体积的比例。统计结果如图 6所示，慢阻肺(COPD)组和非慢阻肺(NonCOPD)组之间的小血管体积占比存在差异。 图6 慢阻肺组与非慢阻肺组小血管体积占比统计 Statistics of the proportion of small pulmonary vasculatures in the COPD and the NonCOPDFig 6表 3中列出了两个实验组中小血管体积占比值(CSA5和CSA10)的平均值和中位值，以上统计结果表明，无论是CSA5还是CSA10，COPD组均显著小于NonCOPD组($p0.001$)，相比于NonCOPD组，COPD组存在小血管变细消失的可能性(Rahaghi等，2019)。 表3 慢阻肺组与非慢阻肺组肺血管定量分析结果 分组 指标 统计结果 COPD CSA5 0.341±0.093 平均值±标准差 0.352 中位值 CSA10 0.589±0.074 平均值±标准差 0.588 中位值 NonCOPD CSA5 0.395±0.089 平均值±标准差 0.407 中位值 CSA10 0.656±0.067 平均值±标准差 0.658 中位值 Quantitative results of pulmonary vasculatures in COPD and NonCOPD groupTable 3另外，从以上慢阻肺数据中筛选出含有GOLD分级信息的影像数据，并对其小血管体积占比进行了统计，结果如表 4所示，除了GOLD1和GOLD2之间差异较小外，GOLD3-4级明显低于GOLD1-2级。从图 7中可以看出，随着慢阻肺严重程度升高，小血管体积占比呈现线性下降趋势。以上结果表明，肺内小血管体积占比与慢阻肺严重程度具有一定相关性。 表4 慢阻肺不同分级肺血管定量分析结果 分组 指标 统计结果 GOLD1 CSA10 0.612±0.051 平均值±标准差 0.599 中位值 GOLD2 CSA10 0.600±0.078 平均值±标准差 0.596 中位值 GOLD3 CSA10 0.565±0.067 平均值±标准差 0.567 中位值 GOLD4 CSA10 0.528±0.053 平均值±标准差 0.536 中位值 Quantitative results of pulmonary vasculatures in different COPD stageTable 4 图7 慢阻肺不同分级组小血管体积占比统计 Statistics of the proportion of small pulmonary vasculatures in different COPD stagesFig 73讨论本文方法将多尺度Hessian矩阵管状结构增强信息融入到各向异性连续最大流分割框架中，实现肺血管的自动分割，分别在血管仿真数据集和低剂量临床数据集上进行了分割性能测试，结果表明本文各向异性正则分割方法的分割精度高于各向同性正则分割方法。此外，还将此分割方法应用到慢阻肺肺血管重塑研究中，并且对肺血管分割结果进行了定量化分析，验证了肺血管重塑规律。肺血管轴向方向信息以各向异性正则项融入到连续最大流分割框架，相比于各向同性分割方法，更利于肺血管轴向增长分割。然而，肺血管与气管壁相邻粘连或交错，在增强与分割环节中都未显性考虑。基于Hessian血管增强方法易受病灶和噪声干扰(如图 8)，对于多中心数据肺血管分割任务，需要采用最优的算法适应参数，使其临床实用性降低。虽然深度学习方法(Qin等，2021)可获得较好的分割性能，但需要标注全血管区域，训练集制作耗时费力。因此，在3维影像标注数据获取困难的条件下，将自监督学习或者弱监督学习与本文方法相结合，增强血管检测的准确性和鲁棒性，是本文未来的一个工作方向。 图8 肺血管错误分割示例 Demonstration of pulmonary vessel segmentation errorsFig 8为了验证慢阻肺中肺血管重塑规律，对收集到的可用临床CT影像进行肺血管分割后，统计了小血管体积占比，回顾性地分析慢阻肺组、非慢阻肺组和不同GOLD分级的慢阻肺组的小血管体积占比情况。分析发现慢阻肺小血管体积占比小于非慢阻肺组，这与肺血管重塑病理改变相吻合，即重塑改变主要包括远端血管变细或消失，近端血管管径增大。然而，当前定量评估是针对整个肺内血管进行统计，所统计的肺血管体积并未区分肺动脉和肺静脉。由于慢阻肺是一类高度异质性疾病，若采用以肺叶或肺段来评估肺血管体积变化，并分别统计肺动脉和肺静脉的体积，或许能发现更有价值的肺血管重塑规律。因此，以肺叶或肺段的肺动脉和肺静脉定量分析是本文的又一工作方向。尽管回顾性统计分析可以得出慢阻肺组比非慢阻肺组的小血管体积占比要小，并且GOLD分级越高，比值越低的结论，但并非所有慢阻肺病例都存在肺血管重塑，且非慢阻肺组也有病例可能存在肺血管重塑。由于肺血管重塑的原因可能是异质多样的，并且在病理学上也并未研究透彻，因此在肺血管体积分析的同时，结合更多临床信息，如BMI(body mass index)、家族病史以及是否吸烟等因素，进行多因素相关性分析，为探寻肺血管重塑的诱发因素或许有帮助。本文实验所用CT影像的层厚跨度大(0.625~1.5 mm)，不同层厚数据对细小血管显影效果存在差异，层厚越小，细小血管成像越清晰完整。本文所用小血管的截面面积小于5 mm2(CSA5)或10 mm2(CSA10)。图 6所展示的实验结果表明，CSA10指标更能凸显COPD和NonCOPD组之间的差异，因此在GOLD分级实验中，仅用CSA10指标进行分析。采用公开的深度学习肺血管分割预训练模型(Qin等，2021)对本文部分临床数据做了肺血管分割和截面量化分析，得到与本文方法一致的慢阻肺肺血管重塑结论。本文对血管定量分析的主要目的是探寻与验证慢阻肺肺血管病变的趋势，而非绝对测量血管半径与体积。因此，实验中暂未分析血管增强步骤对肺血管定量分析的影响，以及距离变换的精度对肺血管体积测量准确度的影响。今后，如果需要对血管管径和体积进行绝对测量，则需要进一步进行实验分析与验证。4结论提出一种各向异性正则的连续最大流分割框架实现肺血管的自动分割，在公开数据和仿真数据集下进行了分割性能测试，在临床数据集上对分割的血管结果进行了定量体积占比分析，验证了慢阻肺肺血管重塑的变化趋势，可为临床研究与应用提供一定技术支持。本文肺血管分割的平均精度逊于深度学习方法，因此，为了进一步增强分割的准确性和鲁棒性，将自监督学习与本文方法相结合是未来的工作方向。此外，除了肺血管管径改变外，其曲率、分叉角度和分支模式等形态学属性也可能发生改变，因此，未来工作将研究树状结构的形态学度量方法，以弥补肺内小血管体积占比单一指标的不足，实现更为完备的肺血管形态学改变的描述。