第13章：IVUS 图像的 Gamma 混合模型（A Gamma Mixture Model for IVUS Imaging）

1. Introduction

本章开篇即将动脉壁中的胆固醇与血细胞团块定义为人所熟知的"易损斑块"（vulnerable plaque），并指出其容易引发心肌梗死、脑卒中等心血管事件[1,2]；此外狭窄本身也是斑块衍生的问题之一。在所有常用于检测冠状动脉粥样硬化斑块的成像技术中，冠状动脉造影长期被视为标准手段，但其在精确测量狭窄程度以及刻画斑块性质方面能力有限——具体而言，多数急性冠脉综合征患者在造影下仅表现为轻度或轻度的冠脉管腔阻塞[2]。鉴于这些局限，检测狭窄区域以及识别不同类型斑块的重要性便得以凸显；作为替代方案，血管内超声（intravascular ultrasonography, IVUS）已被证明能够清晰显示动脉壁的内部形态，因此成为评估形态学病变严重程度的便捷方法[2]。

IVUS 数据采集流程与冠状动脉血管成形术所用标准导管类似——将导管送入动脉并推进至待研究节段；导管前端的压电换能器以旋转方式发射声脉冲并收集回波信号 A-line，每个方向对应一组沿深度排列的反射回波，由此得到 360° 横截面的极坐标表示；接着对图像进行插值并按几何关系转换为笛卡尔坐标（图 13.1）。不同类型斑块在声学响应上的特征是定性的：脂质斑块（lipidic plaque）呈现低回声（echolucent）；纤维斑块（fibrous plaque）回声居中；钙化斑块（calcified plaque）呈强回声，并且由于组织内多重反射常伴有声影[2–5]。然而这种定性刻画虽便于直观解读 IVUS 图像，深入理解回声形态学与病理学演化仍需投入大量研究[6]——定量刻画有助于开发或精化斑块检测、风险预测乃至辅助治疗方法的方法。

获取超声回响应的定量刻画需要分析采集过程中超声成像的物理机制。医学超声图像的形成以脉冲包（pulse packet）沿波束向量轴传播开始；脉冲随介质属性而变形，被不同散射深度的散射体散射并产生时延。结果是接收到的反向散射（backscattered）信号被一种特征性的颗粒状噪声——斑点（speckle）——所污染，其强度取决于分辨率单元内的散射体数量及其尺寸[7,8]。这种乘性噪声（其方差依赖于底层信号）同样见于激光[9]和合成孔径雷达（Synthetic-Aperture Radar, SAR）[10]等采用相干辐射的成像模态。

斑点主要依赖于组织的微观结构，因此其统计量可被用作组织学描述符[11]；其统计特性强烈依赖于散射密度（即每个分辨率单元内散射体的数量与强度分布），由此产生的散斑噪声可分为四类：完全发展的散斑（散射体数量大且无确定性分量，可用 Rayleigh 分布建模[12,13]）；完全分辨的散斑（散射体数量大且存在确定性分量，可用 Rice 分布建模[14]）；部分发展的散斑（散射体数量小且无确定性分量，可用 K 分布建模[15,16]）；以及部分分辨的散斑（散射体数量小且存在确定性分量，可用 K-Homodyne 分布建模[17]）。其中完全发展的散斑是最常见的散斑形成模型——它假设组织或区域由大量散射体组成，其行为如回声扩散器，这些散射体源于尺度与超声波长相当或更小的结构（例如组织实质中存在微观尺度声阻抗变化）。针对这一模型有多种近似，其中最常用的单参数概率模型是用 Rayleigh 概率密度函数（probability density function, PDF）描述包络数据中的像素强度[7,13,18]。除此之外还有其它被提出的刻画散斑的分布，其中最引人注目的是 Nakagami 分布[19]——它是一个双参数分布，可视为 Rayleigh 分布的推广；[20]中即基于 Nakagami 分布提出了刻画反向散射回声的模型。

2.1 Overview

本节聚焦于完全发展的散斑。在完全发展的散斑区域，采集到的信号可以用 Rayleigh 分布建模[13]。然而，为了生成最终的笛卡尔图像，这些原本服从 Rayleigh 分布的数据必须经过插值处理，因此最终图像将不再服从 Rayleigh 分布。本节的目的即在于建立一个能够刻画插值过程后最终分布的概率模型[7,26]。作者指出，若不做对数压缩（这在本节始终作为假设），则文献[18,25]中所给出的结论仍然成立，并且作者所建立的插值概率模型可以方便地推广到其他分布。一些基于最大似然与最大后验的超声图像滤波与分割技术依赖于对图像中不同区域的精确统计建模[25]；这种模型通常根据声学物理分析与超声探头可得信息推导[7,26]，但由于采集过程中的全部信息（例如增益、对比度调整）通常无法获得，因此必须作出若干假设。

2.2 Interpolation Model

超声图像由若干"线"或向量按顺序排列构成[26]——这些向量经过包络检测后形成图像中的线，每条线代表来自不同深度的散射波的时间记录。成像过程以脉冲包沿波束向量轴的传播开始，脉冲形状随介质特性变化；传播的脉冲被位于不同散射深度的物体散射并产生时延；反射回波被换能器接收，并在假定声速恒定的前提下由传播时间估算散射体的深度；这些被截取的波在换能器表面以适当的加权与时延进行积分以实现聚焦与波束形成。最终包络记录的幅度通常会进行对数压缩，但这一步骤取决于具体的超声设备。到此步骤，当观察到完全发展的散斑区域时，Rayleigh 概率分布常被采纳[9]；随后所有线被插值以从按几何排布的图像线形成完整的笛卡尔图像[7,26]。

设 \(\{X_i\}\) 为独立同分布（IID）的 Rayleigh 随机变量，其概率密度为 \(f_X(x) = \frac{x}{\sigma^2} \exp\left(-\frac{x^2}{2\sigma^2}\right),\ x \geq 0\)（式 13.1）。当采用简单的插值方案（例如二维情形下的双线性插值或三维情形下的三线性插值，这是超声设备出于计算效率最可能采用的方案）时，像素的最终插值结果可写为 \(Y = \sum_{i=1}^n w_i X_i\)，其中权重满足 \(\sum_{i=1}^n w_i = 1\)（式 13.2）。插值得到的随机变量 Y 的概率密度没有闭式表达；近年来已提出多种计算方法[27]，本章采用基于求积法的数值方法以求实现简单——其目的是研究该 PDF 的行为以验证文献中的经验方法。

一种简单的计算 \(f_Y(y)\) 的方法是将其视为各独立随机变量加权 PDF 的卷积；通过特征函数可获得闭式表达，因为 Rayleigh 分布的特征函数已知（尽管并不简单[27]）：\(\Phi_X(t) = E[\exp(itX)] = {}_1F_1(1/2; 1/2; -\sigma^2 t^2/2) + i \cdot \sqrt{\pi} \cdot \exp(-\sigma^2 t^2/2) \cdot \operatorname{erf}(\sqrt{\sigma^2/2} \cdot t)\)（式 13.3），其中 \(\operatorname{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} dt\)（式 13.4），且 \({}_1F_1(a;b;x) = 1 + \frac{a}{b}x + \frac{a(a+1)}{b(b+1)} \frac{x^2}{2!} + \cdots\) 为第一类合流超几何函数（式 13.5）。由此 Y 的特征函数为 \(\Phi_Y(t) = \prod_{i=1}^n \Phi_i(t)\)（式 13.6），其中 \(\Phi_i(t)\) 是各 \(w_i X_i\) 的特征函数；权重 \(w_i\) 作用于 \(X_i\)，因此 Y 可视为带有不同 \(\sigma\) 的 Rayleigh 随机变量之和；PDF 通过对 (13.6) 的数值求积得到。

然而这一分布在实际数据的统计估计中并不实用，因为它涉及的参数太多；更简洁的模型更为可取——具体而言，Rayleigh 变量之和用 Gamma 分布来近似。Gamma PDF 只有两个自由参数且尾部行为与 Y 的 PDF 类似；此外，文献中 Gamma PDF 也被用于此类散斑的建模[25]，但未给出理论证明。图 13.2 显示了 Gamma、Gaussian、Rayleigh 之和的特征函数（a）以及在加权和任意权重下逼近误差随项数增加而变化的行为（b、c）；实验表明即使在可以使用中心极限定理的情形下，Gamma 分布的特征函数仍优于正态分布。为验证这一假设，作者按[28]的方式基于采集模型模拟散斑——扫描图像并将数据以矩阵形式记录，再以[13]的散斑形成模型（将组织视为数量众多且尺寸与波长相当的散射体的集合）进行污染；散斑图样通过不假设任何统计分布的随机游走获得以避免引入偏差。图 13.3 显示了无相干回波且散射体数量足以视为完全发展散斑时的重建图像：插值前图像直方图服从 Rayleigh 分布，而插值后图像的直方图则明显由 Gamma 分布精确拟合。

3. Materials

本章所用的刻画数据集是文献[29]中提出的"金标准"——该 IVUS 数据集由 9 条死后冠状动脉组成，并从中选出 50 张含不同性质斑块的图像；之后将动脉切片以通过组织学分析刻画斑块。图像采集按以下流程进行：将动脉自心脏分离并固定于一中等软度的平面上，通过导管以恒压（约 120 mmHg）充入生理盐水以模拟血压；标记远端、近端、左侧和右侧位置的参考标记；探头通过导管引入，并在斑块对应的位置采集射频（RF）数据。RF 数据采用 Galaxy II IVUS 成像系统（Boston Scientific）配合 Atlantis SR Pro 40 MHz 导管（Boston Scientific）实时采集；为采集和存储 RF 数据，成像系统连接一台配有 12 位 Acquiris 采集卡（采样率 200 MHz）的工作站。

每帧 RF 数据组织为 \(N \times M\) 样本的数据矩阵，其中 \(M = 1024\) 为每条 A-line 的样本数，\(N = 256\) 为旋转超声探头的位置数。斑块的组织学验证流程如下：在与预先标记位置对应的部位切割血管，并通过组织学分析确定斑块组成；通过介入专家与病理学家的合作，由参考位置建立组织学检出斑块与对应 IVUS 图像之间的一一对应；为保持组织与 IVUS 图像区域之间的可靠对应，医疗团队手动完成斑块标注任务，并丢弃无法获得可靠对应的图像。最终数据集包含 50 张不同的图像（每张包含一处或多处斑块），全部具有管腔分割；共识别出 69 处斑块并按以下类型进行组织学表征：30 处钙化、14 处脂质、25 处纤维。

IVUS 图像直接从原始 RF 信号重建，而非使用 IVUS 设备自身产生的图像。重建算法如[29]所述（图 13.4），依次包括以下步骤：（1）时间增益补偿（Time Gain Compensation），其中 \(TGC(r) = 1 - e^{-\beta r}\)，\(\beta = \log 10 \alpha f / 20\)，\(\alpha\) 为生物软组织的衰减系数（\(\alpha \approx 0.8\) dB/MHz cm，\(f = 40\) MHz[29]），\(f\) 为换能器中心频率（MHz），\(r\) 为距导管的径向距离（cm）；（2）Butterworth 带通滤波器，截止频率为 \(f_L = 20\) MHz 与 \(f_u = 60\) MHz；（3）Hilbert 变换包络恢复；（4）图像下采样以通过线性插值获得各向同性分辨率；（5）对数压缩；（6）数字开发过程（Digital Development Process, DDP）——一种用于增益非线性调整与边缘增强的非线性处理。经过此重建流程后，IVUS 显示图像可通过将极坐标 \((\rho, \theta)\) 插值到笛卡尔网格 \((i, j)\) 简单获得，从而得到 \(256 \times 256\) 像素的非压缩图像（图 13.1a 与图 13.1b）。

传统显示的 IVUS 图像由极坐标表示 \((\rho, \theta)\) 通过插值至笛卡尔网格 \((\rho, j)\) 得到；而本章用于分析的是下采样步骤之后获得的非压缩极坐标图像（图 13.4）。重建过程中的这一阶段涉及线性滤波，因此 Rayleigh 或 Nakagami 模型不再适用；该线性滤波即使在最终分辨率为初始分辨率整数倍的情况下也会默认执行，以保证各向同性分辨率。

4. Statistical Analysis of Envelope Data

本节测试 Gamma 分布在描述含有内部预处理（如线性滤波与插值——由下采样阶段引入，见图 13.4）的散斑时的表现，并与 Rayleigh 分布及其推广 Nakagami 分布进行比较（这两者被认为是描述散斑统计的良好候选）。用于性能测试的两种度量分别为 Kullback-Leibler 散度（\(D_{KL}\)）与累积分布函数（CDF）的无穷范数（即 Kolmogorov-Smirnov 统计量 \(D_{KS}\)）。

\(D_{KL}\) 是两个概率分布之间差异的非对称度量：\(D_{KL}(p_n; f_X) = \sum_{i=1}^N p_n(i) \log \frac{p_n(i)}{f_X(i)}\)（式 13.7），其中 \(p_n\) 为经验 PDF 估计，\(f_X\) 为理论分布（Rayleigh、Gamma 或 Nakagami）。PDF 的近似通过所研究像素邻域的直方图估计，所用邻域大小为 \(11 \times 11\)（这是在同质区域上提供 PDF 估计的合理样本数），直方图的箱数为 \(n = 30\)。Rayleigh 与 Gamma 分布的参数由所定义邻域内数据的最大似然估计给出；Nakagami 分布的参数按[24]计算。

\(D_{KS} = \sup |\hat{F}_n(i) - F_X(i)|\)（式 13.8），其中 \(\hat{F}_n\) 为经验 CDF，\(F_X\) 为 Gamma 或 Rayleigh 分布的 CDF。选择 \(D_{KS}\) 作为补充度量是因为它不依赖于 PDF 估计，并且可以在样本数较少时计算；此外由 Glivenko-Cantelli 定理可知，若样本来自 \(F_X\)，则 \(D_{KS}\) 几乎必然收敛到 0[30]。比较 Gamma、Nakagami 与 Rayleigh 分布时，针对 \(D_{KL}\) 与 \(D_{KS}\) 两种度量在各邻域（\(11 \times 11\)）上分别进行 Welch t 检验：对每张图像（数据集共 50 张）计算平均值，并以 50 张图像为样本进行 Welch t 检验（选择该检验是因为不能假设方差相等）。检验以分布对（Gamma vs Rayleigh、Gamma vs Nakagami、Rayleigh vs Nakagami）进行，原假设为两组总体均值相同。图 13.5 显示了 50 张图像样本在两种度量下的箱线图；从中可见 Nakagami 的表现优于 Rayleigh（Nakagami 是 Rayleigh 的推广，这一结果符合预期），但对于 Gamma 分布，两种度量的均值均低于 Rayleigh 与 Nakagami——这表明 Gamma 分布在描述散斑概率行为方面表现更佳。为确认这一结果，Welch t 检验在 \(D_{KL}\) 与 \(D_{KS}\) 上的 p 值如表 13.1 所示（所有值均极小，原假设必须拒绝）；因此三种分布的拟合显著不同，Gamma 是两种度量下拟合最好的。本节由此确认：当存在内部预处理（线性滤波与插值）时，Gamma 分布对散斑概率本质的描述优于 Rayleigh 或 Nakagami——这与[7,18,25]的结果一致。

5. Gamma Mixture Model

本节提出 Gamma 混合模型（Gamma Mixture Model, GMM）并描述其参数与系数的计算方法。使用 GMM 的动机在于：斑块的回声形态可能源自不同回声组分的贡献，而各组分服从不同的分布；在 Gamma 散斑假设下，GMM 自然地成为这一描述框架。设 \(X = \{x_i\}\)（\(1 \leq i \leq N\)）为超声图像某一区域内的一组样本（像素强度）。由于下采样阶段削弱了相邻像素之间的相关性，这些样本可视为 IID 随机变量。GMM 假设这些变量来自 \(J\) 个分布的混合：\(p(x_i | \theta) = \sum_{j=1}^J \pi_j f_X(x_i | \theta_j)\)（式 13.9），其中 \(\theta = (\theta_1, \ldots, \theta_J; \pi_1, \ldots, \pi_J)\) 为 GMM 的参数向量，\(\theta_j\) 为第 \(j\) 个 PDF 的参数（本章中 Gamma 分布的参数为 \(\alpha_j\) 和 \(\beta_j\)），其定义为 \(f_X(x | \alpha, \beta) = \frac{x^{\alpha-1} e^{-(x/\beta)}}{\beta^\alpha \Gamma(\alpha)}\)（式 13.10），其中 \(\Gamma(x) = \int_0^\infty t^{x-1} e^{-t} dt\)（\(x > 0\)）为 Euler Gamma 函数；约束条件 \(\sum_{j=1}^J \pi_j = 1\) 必须成立以保证 \(p(x_i | \theta)\) 为定义良好的概率分布。

IID 样本的联合分布为 \(p(X | \theta) = \prod_{i=1}^N p(x_i | \theta)\)（式 13.11）。当引入隐藏离散随机变量 \(Z = \{Z_i\}\)（取值于 \(\{1, \ldots, J\}\)，含义为 \(Z_i = j\) 表示样本 \(x_i\) 属于第 \(j\) 类分布）时，采用期望最大化（Expectation-Maximization, EM）算法以最大化对数似然函数。设 \(\theta^{(n)}\) 为第 \(n\) 次迭代的参数估计，期望步计算对数似然函数 \(L(\theta; X, Z)\) 关于隐藏变量的期望：\(Q(\theta | \theta^{(n)}; X) = E_{Z | \theta^{(n)}, X}[L(\theta | X, Z)]\)（式 13.12）；极大化步通过最大化 \(Q(\theta | \theta^{(n)}; X)\) 得到新估计 \(\theta^{(n+1)}\)；迭代至停止条件 \(\|\theta^{(n+1)} - \theta^{(n)}\| < TOL\)（预设阈值）满足为止。

一般分布的 EM 算法并非新方法，可参考文献[23,31]；对于 GMM，Webb 在[32]中首先给出了推导，[24]给出了另一种类似推导；为完整起见本章给出自己的推导。对数似然关于隐藏随机变量的期望（在已知数据 \(x\) 与前一次估计 \(\theta^{(n)}\) 时）为：\(Q(\theta | \theta^{(n)}; x) = \sum_{i=1}^N \sum_{j=1}^J p(Z_i = j | x_i, \theta^{(n)}) \log p(x_i | \theta_j) + \log \pi_j\)（式 13.13）。

其中 \(p(Z_i = j | \theta)\) 是 \(x_i\) 属于第 \(j\) 类的概率，记为 \(\pi_j\)。另一方面，\(p(Z_i = j | x_i, \theta^{(n)})\) 可由贝叶斯定理导出：\(p(Z_i = j | x_i, \theta^{(n)}) = \frac{p(x_i | \theta_j^{(n)}) p(Z_i = j | \theta^{(n)})}{f_X(x_i | \theta^{(n)})}\)（式 13.14），其中与式 (13.9) 一致 \(\sum_{j=1}^J p(x_i | \theta_j^{(n)}) p(Z_i = j | \theta^{(n)}) = \sum_{j=1}^J \pi_j f_X(x_i | \theta_j^{(n)})\)（式 13.15）。由于 (13.13) 由两项组成，可分别对每一项独立最大化。对依赖于 \(\pi_j\) 的项可借 Lagrange 乘子法（在约束 \(\sum_{j=1}^J \pi_j = 1\) 下）以直接的方式优化，Lagrange 函数为 \(\ell(\lambda, \theta) = \sum_{i=1}^N \sum_{j=1}^J \gamma_{i,j} \log \pi_j + \lambda \left(\sum_{j=1}^J \pi_j - 1\right)\)（式 13.16），其中 \(\gamma_{i,j} = p(Z_i = j | x_i, \theta^{(n)})\)。对每个 \(\pi_j\) 求导并令之为 0，可得 \(\sum_{i=1}^N \gamma_{i,j} = \lambda \pi_j\)（式 13.17）。对 \(j\) 求和得 \(\lambda = N\)，从而最大化 Lagrange 函数（以及对数似然中的相应项）的 \(\hat{\pi}_j\) 值为 \(\hat{\pi}_j = \frac{1}{N} \sum_{i=1}^N \gamma_{i,j} = \frac{1}{N} \sum_{i=1}^N p(Z_i = j | \theta^{(n)})\)（式 13.18）。

接下来可对依赖于 \(\theta_j = (\alpha_j, \beta_j)\) 的 (13.13) 中的项分别关于 \(\beta_j\) 和 \(\alpha_j\) 求极大。首先关于 \(\beta_j\) 求偏导并令之为 0：\(\frac{\partial}{\partial \beta_j} \left[ \sum_{i=1}^N \gamma_{i,j} \log p(x_i | \theta_j) \right] = 0\)（式 13.19），其中 \(p(x_i | \theta_j)\) 的对数似然为 \(\log p(x_i | \theta_j) = (\alpha_j - 1) \log x_i - \frac{x_i}{\beta_j} - \alpha_j \log \beta_j - \log \Gamma(\alpha_j)\)（式 13.20），由此可得 \(\sum_{i=1}^N \gamma_{i,j} \left( \frac{x_i}{\beta_j^2} - \frac{\alpha_j}{\beta_j} \right) = 0\)（式 13.21），从而 \(\beta_j\) 的估计为 \(\hat{\beta}_j = \frac{1}{\alpha_j} \cdot \frac{\sum_{i=1}^N \gamma_{i,j} x_i}{\sum_{i=1}^N \gamma_{i,j}}\)（式 13.22）。

将 (13.22) 代入 (13.13) 并对 \(\alpha_j\) 求导：\(\frac{\partial}{\partial \alpha_j} \left[ \sum_{i=1}^N \gamma_{i,j} \log p \left( x_i \Big| \alpha_j; \frac{\sum_{k=1}^N \gamma_{k,j} x_k}{\sum_{k=1}^N \gamma_{k,j}} \right) \right] = 0\)（式 13.23），重新整理后得到：\(\sum_{i=1}^N \gamma_{i,j} \log x_i - \sum_{i=1}^N \gamma_{i,j} \log \frac{\sum_{k=1}^N \gamma_{k,j} x_k}{\sum_{k=1}^N \gamma_{k,j}} + \sum_{i=1}^N \gamma_{i,j} \log \alpha_j - \sum_{i=1}^N \gamma_{i,j} \psi(\alpha_j) = 0\)（式 13.24），其中 \(\psi(x) = \Gamma'(x)/\Gamma(x)\) 为 Digamma 函数。最终得到 \(\log \alpha_j - \psi(\alpha_j) = \log \frac{\sum_{i=1}^N \gamma_{i,j} x_i}{\sum_{i=1}^N \gamma_{i,j}} - \frac{\sum_{i=1}^N \gamma_{i,j} \log x_i}{\sum_{i=1}^N \gamma_{i,j}}\)（式 13.25）；该式没有闭式解，但可通过数值方法获得（因为 \(f(x) = \log x - \psi(x)\) 行为良好）。由使对数似然极大化的 \(\hat{\alpha}_j\) 值出发，\(\hat{\beta}_j\) 可由 (13.22) 直接得到。

6. Results

本节对 GMM、RMM（Rayleigh 混合模型）与 NMM（Nakagami 混合模型）在组织类别拟合方面的表现进行比较——分别考察它们在不同组分数下对不同组织类型的拟合能力，两种度量 \(D_{KL}\) 与 \(D_{KS}\) 均被采用。\(D_{KL}\) 通过估计经验 PDF（300 箱）来计算；由于这一度量不仅强依赖于样本数，还依赖于箱的选择，因此在不同箱选择下可能产生偏差。相比之下 \(D_{KS}\) 仅依赖于样本，对拟合性能的描述更为鲁棒。所有混合模型均使用三个组分进行拟合的一个示例见图 13.6，结果汇总于表 13.2。

结果表明 GMM 在几乎所有情形下对两种度量（\(D_{KS}\) 与 \(D_{KL}\)）都给出更好的拟合；但也存在值得评论的例外——具体而言，对于管腔（Lumen）四到七个组分的拟合，RMM 在 \(D_{KL}\) 度量上给出更好的结果，而 GMM 在 \(D_{KS}\) 度量上表现更佳。这一现象源于从样本估计 PDF 的困难：管腔类别的值域动态范围过大，PDF 在较低值附近更为集中，从而形成一个难以通过分箱近似精确描述的窄峰；该峰附近 PDF 估计值的波动会降低 \(D_{KS}\) 的精度；而 \(D_{KL}\) 仅使用管腔样本本身，因而能给出更优的误差度量。在少数情形下 NMM 与 RMM 的 \(D_{KS}\) 略低于 GMM，但 RMM 与 NMM 相对于 GMM 的 \(D_{KS}\) 差异几乎可以忽略（最坏情形下约为 0.0005），因此可以认为在这些情形下 GMM 与其他两者的拟合相当。

这些结果可用于选择 GMM 的组分数：随着组分数增加，度量以更慢的速度下降；图 13.7 显示了三种组织类别（钙化、脂质、纤维）以及管腔下 GMM 与 RMM 两种度量的下降速率——下降速率由 \(r_{\text{rate}}(n) = D(n-1) - D(n)\) 计算，其中 \(D(n)\) 为 \(D_{KL}\) 或 \(D_{KS}\) 度量，\(n\) 为混合模型组分数。全局最大值在 \(n = 3\) 处取得，因此可以将三组分视为描述每种斑块与管腔的最优组分数。

7.1 Formulation

本节展示 GMM 所提供的概率图（由贝叶斯定理得到）的应用潜力。GMM 概率图的形式为 \(p(Z_i = j | x_i, \theta^{(n)}) = \frac{p(x_i | \theta_j^{(n)}) p(Z_i = j | \theta^{(n)})}{f_X(x_i | \theta^{(n)})}\)（式 13.26）；这些概率图可用于确定组织类别之间的最可能轮廓，这一信息对开发分割或滤波算法很有价值。作为概率图潜力的一个示例，本节提出一种新颖而快速的超声图像滤波方法——它考虑了原图像的任何线性滤波或插值。

设 \(X_i\) 为 IID 随机变量，服从参数为 \(\sigma\) 的 Rayleigh 分布。经过插值/滤波阶段，这些随机变量按下式变换：\(Y = \sum_{i=1}^4 w_i X_i = \sigma \sum_{i=1}^4 w_i R_i\)，其中 \(\sum_{i=1}^4 w_i = 1\)（式 13.27），其中 \(R_i\) 为 IID 归一化 Rayleigh 随机变量（\(\sigma_R = 1\)）；权重 \(w_i\) 在二维插值情形下由两个 IID 均匀随机变量 \(U_1, U_2 \in [0, 1]\) 按以下方式构造：\(w_1 = U_1 U_2\)，\(w_2 = U_1 (1 - U_2)\)，\(w_3 = (1 - U_1) U_2\)，\(w_4 = (1 - U_1)(1 - U_2)\)。注意新随机变量 Y 不再服从 Rayleigh 分布，而是如第 2 节所示呈 Gamma 类分布；借助 GMM 可以方便地获得描述随机变量概率行为的分布，从而能够有效地计算与各组织类别对应的随机变量的局部条件矩。

线性最小均方误差估计（Linear Minimum Mean Square Error Estimation, LMMSE）滤波器可按下式定义[33]：\(\hat{\sigma} = E\{\sigma\} + C_{\sigma Y} C_{YY}^{-1} (Y - E\{Y\})\)（式 13.28），其中 \(C_{YY}^{-1}\) 是 Y 的协方差，\(C_{\sigma Y}\) 是互协方差。两者以及 \(E\{\sigma\}\) 均可通过对 (13.27) 直接计算得到：\(E\{\sigma\} = \sqrt{\frac{2}{\pi}} \cdot E\{Y\}\)（式 13.29）；\(C_{YY} = E\{Y^2\} - E\{Y\}^2\)（式 13.30）；\(C_{\sigma Y} = E\{(Y - E\{Y\})(\sigma - E\{\sigma\})\} = \sqrt{\frac{2}{\pi}} (E\{\sigma^2\} - E\{\sigma\}^2)\)（式 13.31）；由此 \(E\{Y^2\} = E\{\sigma^2\} \left[ E\{R_i^2\} + E\{R_i\}^2 \right] = E\{\sigma^2\} \cdot \frac{4 - \pi/2}{\pi/2}\)（式 13.32），而 \(\sigma\) 的方差为 \(\text{var}\{\sigma\} = E\{Y^2\} \cdot \frac{\pi/4}{8 - \pi/2} - E\{Y\}^2\)（式 13.33）。

由于 \(\text{var}\{\sigma\}\) 可通过对 \(Y\) 的局部矩以与[34,35]相似的方式计算得到，但该等式不直接依赖于 \(\text{var}\{Y\}\)，因此基于局部矩对 \(\text{var}\{\sigma\}\) 的估计可能产生负值——这是一个依赖于图像（\(Y\) 的取值）的不良特性。作为对该问题的解决方案，作者利用概率图在修正的 LMMSE 等式中获得一个合适的 \(C\)：\(\hat{\sigma} = E\{\sigma\} + C \cdot (Y - E\{Y\})\)（式 13.34）——这是一个线性估计。现在定义一个新的常数 \(C\)，使其具备与 LMMSE 中 \(C_{\sigma Y} C_{YY}^{-1}\) 类似的理想性质：该系数应当描述协方差与互协方差之间的统计差异；当系数接近 0 时，归一化随机变量 Y 的统计行为接近 \(\sigma\)；相反当系数接近 \(1\)（定义域上界），估计退化为 \(\hat{\sigma} = Y\)（\(\sigma\) 的无偏估计，因为 \(E\{\hat{\sigma}\} = E\{Y\} = \sigma\)）。为定义 \(C\)，作者提出取各组织类概率图的梯度——这样在统计行为相同的区域系数将接近 0；反之在异质区域系数将接近 1，估计应保留图像细节。形式上所提系数定义为 \(C_i = \frac{\sum_{j=1}^J |\nabla^\sigma p(Z_i = j | y_i; \theta)|}{\max_j |\nabla^\sigma p(Z_i = j | y_i; \theta)|}\)（式 13.35），其中 \(Z_i\) 为表示 \(Y_i\) 属于第 \(j\) 类的隐藏随机变量；算子 \(\nabla^\sigma\) 表示对与高斯核（标准差为 \(\sigma\)）卷积后的图像所计算的梯度算子；下标 \(i\) 对应图像的每个像素，因此该系数的定义对每个像素独立给出。式 (13.35) 的定义域为 \([0, 1]\)，因为在概率图上用单位步长有限差分实现的离散梯度总是有界于 \([0, 1]\)。

为获得对 \(\sigma\) 的更佳估计，概率图还可用于将 \(\sigma\) 计算为各组织类别局部邻域局部均值的贡献——即每个组织类别条件下 \(\hat{\sigma}(m, n)\) 的局部均值。计算条件局部矩的公式为 \(\hat{Y}(m, n) = \sum_{j=1}^J P(Z = j) \hat{Y}(m, n) | Z = j\)（式 13.36）。当以这种方式计算时，仅属于各组织类别的像素参与各类别局部均值的估计，从而得到对局部矩的更准确估计。这一结果可以视为计算局部矩的一种非各向同性方式——因为当计算局部矩时各像素并非被同等对待。局部条件矩的主要优势在于其使方法具备鲁棒性：换言之，计算局部矩时更大的窗口可能包含不同组织（不同统计量），从而引入对矩计算的偏差；而引入条件概率允许丢弃或削弱属于其他组织的像素的影响。

将属于各组织类别的概率引入局部矩计算后，\(\hat{Y}(m, n) | Z = j\) 可按下式计算：\(\hat{Y}(m, n) | Z = j = \frac{\sum_{(m', n') \in (m, n)} Y(m', n') p(Z(m', n') = j | Y; \theta)}{\sum_{(m', n') \in (m, n)} p(Z(m', n') = j | Y; \theta)}\)（式 13.37），其中 \(p(Z(m', n') = j | Y; \theta) = \frac{\pi_j p(Y(m', n') | \theta_j)}{\sum_{j' = 1}^J \pi_{j'} p(Y(m', n') | \theta_{j'})}\)（式 13.38）。权重 \(P(Z(m, n) = j)\) 是像素 \((m, n)\) 邻域属于第 \(j\) 类组织的概率，由 \(P(Z(m, n) = j) = \frac{\prod_{(m', n') \in (m, n)} p(Z(m', n') = j | Y; \theta)}{\sum_{j' = 1}^J \prod_{(m', n') \in (m, n)} p(Z(m', n') = j | Y; \theta)}\)（式 13.39）给出。

作为示例，对实际 IVUS 图像的下采样阶段进行所提方法的滤波（图 13.8a 显示原图，图 13.8b 显示滤波后图像——两幅图像仅在对数压缩下用于显示）；所用滤波参数为邻域半径 7、\(\sigma = 1\)、GMM 组分数 12——选择 12 个组分是为了描述所有可能的斑块（每类斑块各 3 组，管腔 3 组）。

7.2 Robustness of the Method

在计算局部矩时引入条件矩的优势在于矩估计的鲁棒性——它允许使用任意大小的邻域而无须担心通常各向同性局部均值所带来的模糊效应（实际上各向同性局部均值相当于与常值进行卷积，从而表现为低通滤波）。为检验两个参数变化的影响，图 13.10 显示了由图 13.8a 滤波得到的一些结果，并附上残差噪声模态的相对误差。邻域窗口的尺寸 \(D\) 在 \(\{3, 113\}\) 范围内变化，细节参数 \(\sigma\) 在 \(\{1, 9\}\) 范围内变化；结果显示对于所有参数组合，所提滤波方法均能得到增强后的图像，并且避免了由更大邻域带来的模糊效应。在小邻域情形下，局部强度的变化得以保留——这是滤波最严格的设置，适用于希望良好增强图像但同时保留散斑纹理的情形。随着邻域尺寸 \(D\) 增大，局部矩的估计使用更多来自同类（同组织类别、相同概率分布）的样本进行计算，从而得到更准确的估计；矩的变化更为平滑但始终保持组织之间的边缘。即使 \(D = 113\)（图像尺寸为 \(256 \times 256\)）这种极大尺寸，滤波后的图像也不存在因矩计算不准确而导致的模糊效应——这是由按组织类别条件计算矩的鲁棒性所保证的。

8. Conclusions

本章分析了 IVUS 图像的采集过程以及更好拟合散斑统计性质的概率分布。作者强调线性滤波或插值过程对概率模型的影响（这些过程在采集阶段常被施加）；结果表明完全形成的散斑分布被变换为 Gamma 类分布；该分布在 IVUS 图像情形下与最广泛使用的两个完全形成散斑模型（Nakagami 与 Rayleigh）进行了比较——结果显示 Gamma 分布在 IVUS 图像情形下拟合更佳。一个 GMM 被提出来拟合由不同回声组分贡献而成的回声形态；该 GMM 与 NMM 和 RMM 进行了比较，结果证明 GMM 提供了优于 NMM 或 RMM 的拟合。作为 GMM 在 IVUS 图像上的应用，本章提出了一种细节保持滤波方法——它利用了 GMM 所提供的概率信息：具体而言，GMM 给出每个像素属于各组织类别的概率，这一信息为计算局部矩提供了一种有用且鲁棒的方式（即使在极大邻域下也是如此）；此外概率图还可用于检测组织类别之间最可能的边缘。局部矩的鲁棒计算与边缘概率的组合，可用于定义对完全形成散斑滤波底层参数的估计——这种估计在估计的同时保持细节。该滤波方法在未经对数压缩变换的实际 IVUS 图像上进行了测试；通过将 \(\sigma\) 估计直接除以原图像得到的残差噪声，其模态接近理想值 \(\sqrt{\pi/2}/2\)（参数为 1 的 Rayleigh 分布的均值）——这一结果表明所提滤波方法的一致性。

本章个人批注

本章是定量 IVUS 系列中专门讨论统计建模的一章，与此前以成像临床问题为切入的章节（Ch10–Ch12 围绕 CT，Ch9 关于 MRI 易损斑块）形成显著反差——本章的"内容"几乎完全是概率模型与 EM 推导，作者把 IVUS 的"定量"明确等价于"对散斑统计的精确刻画"。这一定位决定了本章对前置数学知识（Rayleigh/Rice/K/Nakagami/Gamma 分布、特征函数、EM 算法、LMMSE）有较多依赖，对临床医生的可读性也相对受限；但它为后续使用 GMM 进行斑块组织分类、图像滤波乃至分割任务提供了统一的方法论基础。Ch14–Ch17 转入超声颈动脉斑块与回声分析，可以预期作者将从统计建模返回到临床信号——但本章建立的"Gamma 优于 Rayleigh/Nakagami"的结论很可能会在后续章节作为既有事实被引用。

作者选择 Gamma 而非 Nakagami 的核心论据是 §2.2 中的插值论证：原始 Rayleigh 分布在经过双线性/三线性插值后不再保持 Rayleigh 分布，其特征函数的尾部形态更接近 Gamma。这是一条很优雅的论证——它把"哪个 PDF 更好"的统计拟合问题转化为"采集流水线中哪一步引入了线性组合"的信号处理问题。但需要注意：这一论证的成立依赖于"插值等价于对 Rayleigh 变量进行加权和"的假设（式 13.27），而这一等价在加性白高斯噪声主导的极限下由中心极限定理保证；在 IVUS 的低散射密度（小数散射体）情形下，原始信号服从 K 分布而非 Rayleigh，插值后是否仍可由 Gamma 准确描述是开放问题——这也是作者在 §1 中给出 K 分布作为另一类散斑模型的原因。从这一视角看，Gamma 在本数据集中的胜出并不能完全推广到全部 IVUS 成像条件。

在 §5 的 EM 推导中作者给出了 \(\hat{\alpha}_j\) 的隐式方程 (13.25)：\(\log \alpha_j - \psi(\alpha_j) = \log(\bar{x}_j) - \overline{\log x}_j\)，其中 \(\bar{x}_j\) 为加权样本均值、\(\overline{\log x}_j\) 为加权样本对数均值。这一等式的右侧在每步 E 步后可直接计算；左侧作为 \(\alpha_j\) 的单调函数（当 \(\alpha_j > 0\) 时）可以通过 Newton 法或简单二分法求解。作者将 \(\hat{\beta}_j\) 的闭式表达 (13.22) 与 \(\hat{\alpha}_j\) 的数值解并列呈现，这是一种很典型的"两参数 Gamma M-step"模式——可以与 SciPy 的 gamma.fit、Stata 的 gammafit 或 R 的 MASS::fitdistr 在内部实现的算法对照。值得注意的是作者使用全局 12 个组分（每类斑块各 3 + 管腔 3）对一张图像建模——这一选择来自 §6 的组分数曲线（图 13.7）外推以及 §7 滤波的实际需要，与 §3 数据集中"每张图像含 1–多处斑块"的组织学描述并不直接对应；这意味着 GMM 在此处的"组分"并不严格对应"组织类别"，而是带有冗余表达能力的统计基——这一现象与 ICA、字典学习中的过完备分解在精神上相通。

§7.2 关于"邻域尺寸高达 113 时仍无模糊效应"的论证条件矩鲁棒性，作者的论据本质上是：当邻域跨越组织边界时，条件概率会按 \(p(Z=j | Y; \theta)\) 自动降低边界像素的权重——这等价于一种隐式的、各向异性的、依赖数据驱动的窗口函数。这一观察在图像去噪文献中并非全新（与非局部均值、BM3D 等方法的精神相通），但在 IVUS 散斑语境下的明确陈述值得注意：作者把"细节保持"从一个工程目标转化为概率图梯度幅值的形式化定义（式 13.35），这一思路在医学超声图像的 despeckling 文献中较少见。

最后注意一处排印层面的小问题：§5 中 M-step 关于 \(\alpha_j\) 的推导作者写了 \(\log p(x_i | \theta_j)\) 形式（式 13.20），但紧接着在推导 (13.24) 时把 (13.22) 的 \(\hat{\beta}_j\) 直接代入了对数似然，这一做法在 EM 标准教材中通常需要解释"为何 \(\beta_j\) 的条件极大值就是全局极大值"——作者省略了这一论证，直接给出了最终 (13.25)。这在工程导向的论文中是常见做法（只要数值实验能复现即可），但在严格的统计推断视角下属于跳步。

与上下章的衔接（一段话）

本章是 Part IV Ultrasound Atherosclerosis Imaging 的开篇章节，标志着全书的关注点从基于电离辐射的 CT（Ch10–Ch12）与基于非电离辐射的 MRI（Ch3–Ch9）切换到基于机械波的 IVUS/超声（Ch13–Ch20）。在前一章（Ch12）中，作者把"定量 CT"界定为几何参数（PL、MLD、%DS、%AS、%AV、RI）与成分分类（CP/HD-NCP/LD-NCP）两套方法，并明确承认几何测量有"宽符合性限度"、成分分类受 HU 重叠困扰；Ch13 则把"定量 IVUS"重新定义为"散斑统计的精确建模"——这种定义本身就体现了成像物理的差异：CT 的对比度由组织密度差决定、可逐点读出 HU，而 IVUS 的对比度由散射体密度决定、必须以概率分布描述像素强度。在本卷书中，Ch13 与 Ch7（"Carotid Plaque Stress Analysis"）共同代表了"方法学深度"型章节——它们的共同特点是使用大量数学/统计工具来回答一个临床问题；从后续衔接看，Ch14–Ch17 将从本章的统计建模返回到更临床化的"超声颈动脉斑块 → 卒中预测 / 狭窄定量 / 成分分类"语境，但本章所建立的 Gamma 模型很可能以"工具箱"的方式在 Ch16（超声组织学/生化成分）、Ch17（自动 IMT CALEX）、Ch19（超声斑块分割）等章节中再次被引用——读者可以在这些章节留意作者是否重新引入 GMM/NMM/RMM 作为对照方法。