用于表征三层生物组织热物性和层厚度的热波技术反演算法

《International Journal of Thermal Sciences》：A thermal wave technology inversion algorithm for characterization of thermophysical properties and layer thicknesses in three-layer biological tissue

　　

**摘要要点**

- 分步反演算法将病态问题降维为三个参数的子问题。
- 层界检测无需预设层厚度。
- 伪热导率谱（PTCS）增强了深度依赖响应在模型配置间的条件分离度。
- 三层体模实验提供了分层检测的概念验证。

**符号表**

*拉丁符号*

- $b$：加热器半宽（m）
- $C_p$：比热容（J kg$^{-1}$ K$^{-1}$）
- $C_R$：DE交叉概率
- $C_V$：体积热容（J m$^{-3}$ K$^{-1}$）
- $d_F$：离散Fréchet距离
- $F$：DE变异因子
- $I_0$：电流幅值（A）
- $k$：热导率（W m$^{-1}$ K$^{-1}$）
- $k_{ps}$：伪热导率（W m$^{-1}$ K$^{-1}$）
- $L_n$：第 $n$ 层的厚度（m）
- $m$：傅里叶变换变量（m$^{-1}$）
- $N$：层数
- $N_f$：测量频率数
- $N_{pop}$：DE种群规模
- $p_0$：单位长度的交流加热功率（W m$^{-1}$）
- $q_n$：复波数
- $r$：重复采集次数
- $R$：电阻（Ω）
- $s$：扫描步序索引
- $s_0$：预热窗口尺寸
- $S$：最终扫描步序索引
- $T$：温度（K）
- $u_n = \sqrt{m^2 + q_n^2}$
- $z_\tau$：变异检测阈值

*希腊符号*

- $\alpha$：热扩散系数（m$^2$ s$^{-1}$）
- $\gamma_n = u_n k_n$
- $\delta$：热穿透深度（m）
- $\varepsilon$：相对地板参数
- $\nu_n = u_n L_n$
- $\omega$：角激励频率（rad s$^{-1}$）
- $\rho$：密度（kg m$^{-3}$）
- $\sigma$：标准差
- $\varphi$：单调耦合

*下标与上标*

- $n$：层索引
- ps：伪（pseudo）
- exp：实验
- sim：仿真
- $\vartheta$：参数向量

---

## 1. 引言

热消融技术——包括激光、射频和微波等模态——通过聚焦能量沉积来摧毁靶向组织区域[1]、[2]。这些手术的成功依赖于精确的生物热传递模型，而该模型的预测取决于两类输入：每层组织的热物性参数（热导率 $k$ 和体积热容 $C_V$）以及结构异质区域之间的几何边界[3]。即便这些输入中存在轻微误差，也会导致严重的临床后果。例如，病理性异质性可使预测的消融体积变化高达36%[4]；比热偏差25%会使预测的消融长度偏移27%[5]。由于组织热物性随每位患者的病理状态而异，且无法可靠地从文献值或离体测量中获得[6]、[7]，因此体内外原位个体化表征至关重要。

传统原位方法——包括自加热热敏电阻[8]、双针探针[9]以及共形3ω传感器[10]——只能返回单一的体有效值，无法在数十微米尺度上分辨微观异质性[11]。频域热波方法提供了一种替代途径：改变激励频率会改变分层介质的热响应，因此宽带频率扫描可包含关于其物性与结构的深度依赖信息[12]。

挑战在于反演问题。从测量到的频率响应中提取分层参数高度病态，因为各层的热导率、体积热容和厚度之间存在强烈耦合。智能优化方法已被证明能有效应对这种病态性：模糊推理方案同时重构了参与介质的热边界条件与物理属性[13]；类似方法在激光诱导热治疗中恢复了表面加热功率与内部温度[14]，或恢复了生物组织的温度依赖物性[15]。这些研究表明热物理量可以从表面或边界测量中稳健提取，但它们均在已知几何域上求解——要么为单一区域，要么为边界已知的介质——无法解析未知的分层结构。当介质被明确分层时，可通过预先给定层厚度来确保可处理性。Gonì和Schmidt[12]在微电子薄膜堆叠中验证了这一策略，Qiu等人[16]近期将其扩展至活体组织，利用预设的解剖学厚度提取了表皮和真皮的热导率。

对外部厚度信息的依赖是这些方法的局限性。获取层厚度需要补充成像手段（OCT、MRI或超声），这会引入配准误差，阻碍独立热诊断的实现，尤其在病理组织边界因患者而异的情况下更是如此。因此，仅从热信号中同时估计热物性和层几何结构的反演方法可减少这种依赖。

本文提出了一种分步算法，旨在不预设被测配置中层厚度的情况下估计层参数与边界。多层频域解析解具有固有递归结构：将正演模型从 $N$ 层扩展到 $N+1$ 层恰好引入三个新未知量。仿照这一递归结构，反演分阶段顺序进行，每阶段通过差分演化求解一个三参数子问题，同时锁定所有先前估计的参数值。这种分解将一个大规模病态优化问题（$3N-1$ 个未知量）替换为一系列低维、有界的三参数子问题。基于导数的变换将原始2ω温度响应映射为频率依赖的伪热导率谱（PTCS），用于增强所建模层配置之间的分离度。为了在不需要外部厚度数据的情况下确定各阶段的光谱边界，自适应频率扫描监测测量PTCS与模型PTCS之间的离散Fréchet距离。该距离的突然增加——称为"损失变异"（loss mutation）——被视作当前模型阶次不足以描述响应的指示，从而触发扩展到下一阶段。当扫描达到最低测量频率而未触发变异时，算法自然终止。

该算法在60个三层动脉粥样斑块配置上进行了数值评估（5个纤维帽厚度 × 4个脂质核厚度 × 3种组分），并在玻璃/Mylar/玻璃三层体模上进行了实验验证。

---

## 2. 数学建模与反演算法

### 2.1 物理模型

测量使用表面安装的微型加热器和频率扫描来探测分层介质。一个半宽为 $b$ 的微型线加热器置于多层生物组织表面，此处将其建模为局部分层系统——三层平面结构（纤维帽、脂质核和外膜；图1）。加热器由正弦电流 $I(t) = I_0 \sin(\omega t)$ 驱动，产生的焦耳加热包含稳态（直流）分量和两倍激励频率处的谐波分量：

$$P(t) = I^2(t) \cdot R = I_0^2 R \sin^2(\omega t) = \frac{I_0^2 R}{2} [1 - \cos(2\omega t)] \tag{1}$$

其中 $R$ 为加热器条的电阻。

2ω温度响应为主要可观测量。由于局部温度变化不超过2 K，在恒定物性近似下[6]，忽略一个测量周期内温度依赖的物性变化。由于加热器半宽仅为500 nm，探测体积高度局部化，组织可视为局部平面多层结构。热扩散长度为：

$$\delta = \sqrt{\frac{\alpha}{2\omega}} \tag{2}$$

其中 $\alpha = k/C_V$ 为热扩散系数，仅作为构建频率采样网格的启发式指标（第3.1.1节）。

**图1.** 2ω测量配置示意图。(a) 三维视图：由正弦电流 $I(\omega)$ 驱动的微型线加热器置于局部三层组织（代表动脉粥样斑块：纤维帽、含特征性脂质液滴的脂质核、外膜）表面。(b) $x$–$y$ 平面二维截面几何结构（$x$：切向；$y$：深度方向）。加热器施加谐波表面热通量；产生的两倍激励频率（2ω）温度振荡构成主要测量可观测量。

### 2.2 正演多层解析模型

在通用 $N$ 层模型中（图2），每层 $n$ 由其厚度 $L_n$、热导率 $k_n$ 和热扩散系数 $\alpha_n = k_n/(\rho_n C_{p,n})$ 表征。表面加热器（半宽 $b$，位于 $y = 0$）以角频率 $2\omega$ 施加谐波热通量 $p_0$；所有其他表面边界绝热。在所有内部界面处施加温度和热通量连续条件，忽略界面热阻，最深一层（$n = N$）视为半无限。

在频域中应用分离变量并施加逆傅里叶变换[17]后，2ω表面温度可表示为傅里叶变换变量 $m$（m$^{-1}$）上的积分。解以下列辅助变量形式表达（$i = \sqrt{-1}$）：

$$q_n^2 = \frac{i \cdot 2\omega}{\alpha_n}, \quad u_n = \sqrt{m^2 + q_n^2}, \quad \gamma_n = u_n k_n, \quad \nu_n = u_n L_n \tag{3}$$

**图2.** 通用 $N$ 层解析模型示意图。坐标原点 $o$ 位于域左上角，$x$ 为切向，$y$ 为正向下方向。表面热通量 $p_0$（红色箭头）施加在 $y = 0$ 处的加热器覆盖区域。每层 $n$（$n = 1, 2, 3, \ldots, N$）具有厚度 $L_n$、热导率 $k_n$ 和热扩散系数 $\alpha_n$。最深层 $N$ 视为半无限。

**三层模型。** 对于三层系统（第三层半无限），2ω温度振幅为：

$$T_{2\omega, 3L}(\omega) = \frac{p_0}{\pi k_1} \int_0^{\infty} \frac{\sin^2(bm)}{(bm)^2 (m^2 + q_1^2)} \cdot \frac{\gamma_1 \cosh\nu_1 (\gamma_3 \sinh\nu_2 + \gamma_2 \cosh\nu_2)}{\gamma_2 \sinh\nu_1 (\gamma_3 \sinh\nu_2 + \gamma_2 \cosh\nu_2) + \gamma_1 \cosh\nu_1 (\gamma_3 \cosh\nu_2 + \gamma_2 \sinh\nu_2)} \, dm \tag{4}$$

**二层模型。** 令 $N = 2$ 并将第二层视为半无限（即 $\nu_2 \to \infty$，消除了对 $\gamma_3$ 的依赖），得到：

$$T_{2\omega, 2L}(\omega) = \frac{p_0}{\pi k_1} \int_0^{\infty} \frac{\sin^2(bm)}{(bm)^2 (m^2 + q_1^2)} \cdot \frac{\gamma_1 \cosh\nu_1}{\gamma_1 \cosh\nu_1 + \gamma_2 \sinh\nu_1} \, dm \tag{5}$$

**单层模型。** 对于均匀半无限介质，2ω温度振幅退化为经典斜截距形式[18]：

$$T_{2\omega, 1L}(\omega) \approx -\frac{p_0}{2\pi k_1} (\ln\omega + \ln(b^2/\alpha_1) - \ln 2) - \frac{i \cdot p_0}{4 k_1} \tag{6}$$

加热器半宽 $b$ 通过 $\sin^2(bm)/(bm)^2$ 核函数进入每个正演模型，因此在每次反演中均予包含。

### 2.3 解析解的递归结构

对 $N$ 层系统所有未知量进行同时反演（$3N-1$ 个参数；表1）严重病态。然而，解析解具有递归结构：比较式(4)–(6)可以看出，每次层扩展恰好引入三个新未知量：新增层的热导率 $k_{N+1}$ 和体积热容 $C_{V,N+1}$，以及紧邻新增层之上的层厚度 $L_N$（该层此前为半无限，因此无需厚度参数）。

这一结构直接决定了算法设计：不求解一个大规模问题，反演分 $N$ 个顺序阶段进行，每阶段求解一个三参数子问题，同时锁定所有先前估计的参数值。本文的三层评估验证了该机制的两次调用（1→2层和2→3层）。

**表1.** 未知量随层数累积。每次层扩展恰好增加三个未知量，使反演可作为一系列三参数子问题进行。

| 层数 | 新增未知量 | 未知量总数 |
|------|-----------|-----------|
| 1 | $k_1, C_{V,1}$ | 2 |
| 2 | $k_2, C_{V,2}, L_1$ | 5 |
| 3 | $k_3, C_{V,3}, L_2$ | 8 |
| ? | ? | ? |
| $N$（半无限） | $k_N, C_{V,N}, L_{N-1}$ | $3N - 1$ |

### 2.4 通过伪热导率谱进行特征增强

如图3(a)所示，所选单层、二层和三层示例的原始2ω温度谱差异不足0.04 K，这些配置之间视觉分离度有限。这一观察动机了一个基于导数的变换，旨在增强其分离度，将 $T_{2\omega}$ 映射为频率依赖的伪热导率谱（PTCS）：

$$k_{ps}(\omega) = -\frac{p_0}{2\pi} \left[\frac{d\,\mathrm{Re}(T_{2\omega})}{d\ln\omega}\right]^{-1} \tag{7}$$

"伪热导率"这一术语来源于均匀半无限极限：将式(6)代入式(7)可得 $k_{ps} = k$。限定词"伪"将此极限恒等式与多层情形区分开来——在多层情况下，$k_{ps}(\omega)$ 是基于导数的特征变量而非深度分辨的热导率。随着频率降低，热扩散长度增加，更深处的层逐渐对响应做出贡献。由此产生的低频处对单层平台的偏离反映了下层物性的逐步混合，但个别峰、谷或偏离无法唯一归因于特定层或参数。

相同的变换应用于测量响应和每个候选正演模型响应，反演过程在频率范围内比较完整的PTCS曲线。因此，图3(b)展示的是特征分离而非逐点物理解释。这种深度依赖的频率响应为将扫描划分为连续频率范围并逐阶段引入更深层参数提供了物理基础。阶段边界由损失变异准则决定，而非通过将个别PTCS特征分配给特定层来确定。

**图3.** 三种模型配置的原始2ω温度信号与伪热导率谱（PTCS）对比：单层半无限介质（$k_1 = 0.476$ W·m$^{-1}$·K$^{-1}$）、二层系统（纤维帽 $L_1 = 100\,\mu$m 置于纯脂质核上，$k_2 = 0.200$ W·m$^{-1}$·K$^{-1}$）、三层系统（额外添加外膜 $k_3 = 0.490$ W·m$^{-1}$·K$^{-1}$，$L_2 = 500\,\mu$m）。(a) 同相2ω温度 vs. $\ln(\omega/\mathrm{rad\,s}^{-1})$：三条曲线紧密重叠；插图显示差异不足0.04 K，说明原始信号中层配置的视觉分离度有限。(b) PTCS $k_{ps}(\omega)$ vs. $\ln(\omega/\mathrm{rad\,s}^{-1})$：单层曲线（蓝色）为 $k_1$ 处的水平线；二层（绿色）和三层（红色）曲线分离更为明显，尤其在低频处。其峰谷被视为曲线特征，不进行单独的分层物理解释。

### 2.5 基于Fréchet距离的目标函数

反演最小化测量PTCS $P = k_{ps,exp}$ 与模型PTCS $Q = k_{ps,NL}$ 之间的差异。采用离散Fréchet距离 $d_F$[19] 作为距离度量，该度量在约束对应点必须遵循相同顺序的前提下，衡量两条曲线之间的最坏情况偏差：

$$d_F(P, Q) = \min_{\varphi \in \mathcal{M}(p,q)} \max_{(i,j) \in \varphi} \|P_i - Q_j\| \tag{8}$$

其中 $P_i$ 和 $Q_j$ 为两条曲线的第 $i$ 个和第 $j$ 个采样点，$\varphi$ 为单调耦合（保持顺序的采样点配对），$\mathcal{M}(p,q)$ 表示 $p$ 点曲线 $P$ 和 $q$ 点曲线 $Q$ 之间所有此类耦合的集合。

选择该度量的动机源于分步反演。每个阶段估计的参数被锁定并传递至下一阶段，因此当后续阶段曲线被评估时，上游估计误差可能表现为分散的局部不匹配。逐点 $l_2$ 目标函数在每个频率处累积残差贡献，因此对这些前向传递误差引起的局部逐点不匹配敏感。加权最小二乘法则会引入权重，从而带来额外的主观选择。

通过允许单调非逐点配对[20]，Fréchet距离比较的是整体有序的PTCS形状，因此更适合本分步反演。

单调配对可以减少与谱特征的有限局部位移相关的惩罚；在存在继承误差的情况下，这种容差是有意为之的。2.6. 逐步差分进化算法 由于 $L$、$k$ 和 $CV$ 之间存在强耦合，基于梯度的方法容易收敛到局部极小值 [21]。因此，本文采用差分进化（DE）——一种基于种群的全球优化器 [22]——以体积热容 $CV=\rho C_p$ 作为优化变量（即控制方程中唯一出现的密度与比热的乘积）。与 2.3 节中的逐步分解一致，每个阶段最多优化三个未知数（图 4）。DE 采用 DE/rand/1/bin 策略 [23]（$N_{pop}=30$，$F=0.8$，$CR=0.9$，最多 500 代，相对变化小于 $1\times10^{-6}$ 时提前终止）。搜索边界覆盖软组织的完整范围 [6]、[24]：$k\in[0.18,0.60]\,\mathrm{W\cdot m^{-1}\cdot K^{-1}}$，$CV\in[1.8,4.7]\,\mathrm{MJ\cdot m^{-3}\cdot K^{-1}}$，$L_1\in[30,350]\,\mu\mathrm{m}$，$L_2\in[150,1400]\,\mu\mathrm{m}$。该算法基于 Python 3.10 实现，使用 SciPy 1.10 进行差分进化和 PyTorch 2.0 进行傅里叶积分求积；源代码可在数据可用性声明中列出的代码仓库中获取。由于 DE 是一种随机优化器，附录 E 中评估了本算法在不同随机种子下的可重复性。 图 4. 物理引导的逐步差分进化（DE）算法的示意流程图。在每个阶段，恰好将三个未知数追加到正向模型中，同时之前估计的所有参数值保持锁定。固定的损失变异准则（2.6.1 节）在无需几何先验输入的情况下实现自主频率轴划分。当自适应扫描达到最低测量频率时，算法自然终止。2.6.1. 带损失变异的自适应频率扫描 每个阶段使用相同的机制：自适应频率扫描从高频率到低频率逐步进行，每步一个频率点，并在每步评估 Fréchet 损失。该损失对近期趋势的突然偏离表明当前模型阶次对于新纳入的频率不够充分，此时算法推进到下一个反演阶段。 设 $\omega_1>\omega_2>\cdots>\omega_{N_f}$ 表示按降序排列的测量频率，设 $\{x^{(s)}\}_{s=1,2,\ldots}$ 表示一个监测序列，此处实例化为扫描步骤 $s$ 时测量响应与模型响应之间的最小 Fréchet 距离 $d_F^*(s)$。 检测并非对原始损失值进行监测，而是对一阶差分 $\Delta x^{(s)}=x^{(s)}-x^{(s-1)}$ 进行操作，该一阶差分代表损失的变化步长速率。在最近的 $s_0=5$ 个一阶差分上的滚动窗口提供滚动均值 $\hat{\mu}_\Delta$ 和滚动标准差 $\hat{\sigma}_\Delta$：

(9) $\hat{\mu}_\Delta = \frac{1}{s_0-1}\sum_{i=1}^{s_0-1}\Delta x(s-s_0+i)$

(10) $\hat{\sigma}_\Delta = \max\left(\sqrt{\frac{1}{s_0-2}\sum_{i=1}^{s_0-1}(\Delta x(s-s_0+i)-\hat{\mu}_\Delta)^2},\ \varepsilon|\hat{\mu}_\Delta|\right)$

其中 $\varepsilon=0.05$ 是一个相对下限，防止当损失平稳演化时 $\hat{\sigma}_\Delta$ 坍缩为零。由于窗口随每一步向前滚动，基线自动适应渐进趋势；只有相对于近期趋势的突然加速才会触发检测。 损失变异检测：在每个步骤 $s>s_0$，通过计算标准化分数将当前一阶差分与滚动基线进行比较：

(11) $z_\Delta(s)=\frac{\Delta x(s)-\hat{\mu}_\Delta}{\hat{\sigma}_\Delta}$

在第一个 $z_\Delta$ 超过阈值 $z_\tau$ 的步骤处宣布损失变异：

(12) $s^*=\min\{s>s_0\ |\ z_\Delta(s)>z_\tau\},\ z_\tau=3.0$

采用单侧准则（$z_\Delta>z_\tau$）是因为所监测的事件是优化失配值在受影响频率（即受未建模层影响的频率）进入扫描时出现向上偏离；向下偏离不被视为模型阶次转变。$z_\tau=3.0$ 的经验值作为三个标准差阈值，用于检测与近期损失演化的显著偏离。窗口大小 $s_0=5$ 作为平衡局部波动率估计与趋势跟踪的短滚动窗口被选定。前 $s_0$ 步被排除在监测之外（预热阶段），因为扫描开始时狭窄的谱窗口会产生不可靠的损失值，这是由于边缘裁剪伪影所致。 边界截断：当在步骤 $s^*$ 检测到变异时，谱边界不设在变异点本身，而设在滚动窗口的中点，即 $\omega_{s^*-\lfloor s_0/2\rfloor}$。这种选择平衡了两个相互竞争的效果：纳入足够的低频数据以获得足够的谱覆盖，同时排除已经受下一层影响的频率。当前阶段产生的谱区域为 $\omega\in[\omega_{s^*-\lfloor s_0/2\rfloor},\ \omega_{upper}]$，其中 $\omega_{upper}$ 继承自前一阶段。将边界设在窗口中点而非 $s^*$ 处，是对滚动窗口引入的滞后所做的实用校正。该偏移量并非作为统计最优边界估计器推导得出。附录 F 报告了 $s_0$ 和 $z_\tau$ 对三种代表性高对比度配置的局部灵敏度检验。2.6.2. 第 1 阶段：浅表层提取 在选定的高频区域内，$2\omega$ 响应被视为主要对第一层敏感。因此，在每个扫描步骤拟合单层模型，由此得到的 $d_F^*(s)$ 作为监测序列。在步骤 $s^*$ 处的损失变异表明谱窗口开始包含第二层的信息。 一旦 Region I 被划定，热导率 $k_1$ 估计为该区域内 PTCS 平台值的平均值：

(13) $k_1=\frac{1}{s^*-\lfloor s_0/2\rfloor}\sum_{s=1}^{s^*-\lfloor s_0/2\rfloor}k_p^s(\omega_s)$

$CV_1$ 然后通过最小化测量同相 $2\omega$ 温度曲线与模型同相 $2\omega$ 温度曲线在完整 Region I 范围内的 Fréchet 距离获得：

(14) $\hat{CV}_1=\arg\min_{CV_1} d_F[\mathrm{Re}(T_{2\omega,\mathrm{exp}}),\ \mathrm{Re}(T_{2\omega,1L}(k_1,CV_1))]|_{\omega\in\mathrm{Region\,I}}$

随后 $\{k_1, CV_1\}$ 在所有后续阶段中锁定。2.6.3. 第 2 阶段及后续阶段：通过 DE 耦合扫描进行递归层扩展 在第 $N$ 阶段（$N\geq 2$），之前估计的参数值被锁定，留下三个新的未知数：$\theta_N=\{L_{N-1}, k_N, CV_N\}$。在每个扫描步骤 $s$，DE 在扩展的谱窗口 $[\omega_{\min}^{(s)},\ \omega_0^{(N)}]$ 上，在生理搜索边界 $\mathcal{B}$（2.6 节）内优化 $\theta_N$，其中 $\omega_0^{(N)}$ 是继承自前一阶段的上边界：

(15) $\hat{\theta}_N^{(s)}=\arg\min_{\theta_N\in\mathcal{B}} d_F[k_{ps,\mathrm{exp}},\ k_{ps,NL}(\theta_N)|_{\omega\in[\omega_{\min}^{(s)},\omega_0^{(N)}]}]$

其中 $k_{ps,\mathrm{exp}}$ 是测量的 PTCS，$k_{ps,NL}$ 是 $N$ 层正向模型预测的 PTCS。优化后的损失 $d_F^*(s)$ 再次作为监测序列。如果在步骤 $s^*$ 检测到变异，则接受并锁定来自中点步骤 $s_m=s^*-\lfloor s_0/2\rfloor$ 的参数向量：

(16) $\hat{\theta}_N\leftarrow\hat{\theta}_N^{(s_m)}$

然后算法进入第 $N+1$ 阶段，参数向量为 $\theta_{N+1}$。如果在扫描到达最低测量频率之前没有发生变异（自然终止），则接受最终扫描步骤的参数向量：

(17) $\hat{\theta}_N\leftarrow\hat{\theta}_N^{(S)}$

其中 $S$ 表示最后一个扫描步骤；不再推断更深层。 在本次三层评估中，第 2 阶段检测到变异并进入第 3 阶段。第 3 阶段到达最低频率 $\omega_{N_f}$ 而未触发变异，因此算法在三层处自然终止。3. 数值与体模评估 3.1. 数值模型 3.1.1. 有限元模拟 在 COMSOL Multiphysics 6.1（COMSOL Inc.，瑞典）中构建了时域有限元模型，以在交变表面加热条件下为三层动脉粥样硬化斑块（纤维帽、脂质核、外膜）生成合成 $2\omega$ 数据。利用双边对称性，仅对横截面的右半部分建模（图 5）。第 1 层和第 2 层表示为矩形，第 3 层表示为半径 2 mm（超过最大渗透深度的 5 倍）的四分之一圆，以近似半无限条件 [25]。 加热器（半宽 500 nm）放置在表面。每层中的热传导遵循具有恒定各向同性属性的瞬态扩散方程，由表面条带 $x\in[0,b]$ 上的焦耳热源驱动；所有其余边界为绝热或对称，内部界面处温度和热通量连续。为解析深层结构，激励按热渗透深度均匀采样而非按频率采样，在 0.08 Hz 至 1249 Hz 范围内给出近似对数频率间隔，集中在低频区域以增大分辨率。控制方程、热源定义、网格和深度均匀采样详见附录 C。 图 5. 用于数值评估的二维 COMSOL Multiphysics FEM 域示意图。仅对轴向横截面的右半部分建模，利用双边对称性。第 1 层（纤维帽，厚度 $L_1$）和第 2 层（脂质核，厚度 $L_2$）表示为矩形子域。第 3 层（外膜）建模为半径 $r=2\,\mathrm{mm}$ 的四分之一圆，以复现半无限基体条件。金条加热器（红色条带，半宽 $b=5\times10^{-4}\,\mathrm{m}$）放置在域表面的左上角。插图：加热器区域的放大视图。3.1.2. 参数空间 通过组合五种纤维帽厚度（$L_1=50-250\,\mu\mathrm{m}$）、四种脂质核厚度（$L_2=250-1000\,\mu\mathrm{m}$）和三种脂质核成分（0%、50% 和 75% 纤维），构建了 60 种配置的参数空间，覆盖临床上报道的范围 [26]、[27]、[28]、[29]。脂质核的热导率从 0.200（纯脂质）到 0.407 W·m?1·K?1（75% 纤维），产生的第 1/2 层对比度从 2.4× 到 1.2×。完整的热物性参数和配置范围见附录 C 中的表 C.7。3.1.3. 数据处理 在每个激励频率下，模拟两个完整的热周期（$N_{cyc}=2$），每个周期 160 个采样点。$2\omega$ 温度振幅通过离散傅里叶变换从采样波形中提取，对应于相应谐波分量的频率。每个频率处的同相分量 $T_{2\omega,\mathrm{Re}}(\omega_n)$ 产生用于通过公式 (7) 计算 PTCS 的数据集。 进入公式 (7) 的单位长度交流焦耳加热振幅为：

(18) $p_0=\frac{1}{2}\frac{I^2\rho_{\mathrm{Au}}h_{\mathrm{Au}}}{S_{\mathrm{Au}}^2}$ 3.2. 实验装置 样品组装：三层体模由浮法钠钙玻璃片（第 1 层，0.4 mm）、透明双向拉伸 Mylar 聚酯薄膜（第 2 层，0.3 mm）和厚浮法钠钙玻璃块（第 3 层，3 mm）组装而成，后者作为半无限背衬。每层的横向延伸范围（60×15 mm）均大于扫描中最大的热扩散长度。实验前在均匀静载荷下将各层压合在一起，以最小化界面空隙；外层玻璃层具有相同的热物性，而中间的 Mylar 薄膜具有较低的热导率，从而在每个内界面处产生热导率对比。加热与检测：溅射在顶部玻璃表面的铟锡氧化物（ITO）条带（宽 300 μm，长 30 mm）作为焦耳加热源（图 6），由波形发生器（33600A，Keysight Technologies，美国）以 11 V 峰峰值电压驱动。$2\omega$ 表面温度响应通过长波红外热像仪（A655sc，FLIR Systems，美国；640×480 像素，噪声等效温差 <30 mK）记录。制造商规定的测量精度为 ±2 °C 或 ±2% 读数（取较大者）[30]。为在长波波段提供正确的发射率，成像表面涂覆一层均匀的薄黑色丙烯酸层。激励在 1000 至 10 mHz 之间按对数间隔逐步进行 20 个频率，这是热渗透深度从第一层玻璃内部经过 Mylar 薄膜进入背衬的频段。信号提取与重复：在每个频率下，数据采集仅在表面达到准稳态（5 s 窗口内直流温度漂移低于 5%）后才开始。温度在横跨加热器条带中间三分之一的感兴趣区域上进行空间平均，以减少端部效应的影响，同相 $2\omega$ 振幅通过参考驱动信号的离散傅里叶变换提取（3.1.3 节）。完整的频率扫描在三个独立组装和排列的实验中重复进行，因此实验间的展幅反映了测量噪声与运行间装配变异的综合效果。在每个实验的每个频率下，三次重复振幅测量取平均后用于反演。对于实验内这三次重复，样本标准差为 $s_f$，则平均振幅的标准误为 $s_f/\sqrt{3}$。三个独立组装实验之间的运行间标准差单独保留；它包含残余测量变异和装配变异，因此不解读为纯仪器噪声。钠钙玻璃的参考物性由热常数分析仪（TPS3500，Hot Disk AB，瑞典）独立测得（$k_{\mathrm{glass}}=1.137\,\mathrm{W\cdot m^{-1}\cdot K^{-1}}$，$CV_{\mathrm{glass}}=1.001\,\mathrm{MJ\cdot K^{-1}\cdot m^{-3}}$），Mylar 薄膜的参考物性取自文献 [31]（$k_{\mathrm{Mylar}}=0.155\,\mathrm{W\cdot m^{-1}\cdot K^{-1}}$，$CV_{\mathrm{Mylar}}=1.629\,\mathrm{MJ\cdot K^{-1}\cdot m^{-3}}$）。 代表性的红外帧（图 7a）显示了加热器条带中间三分之一的热场（ROI，黑色矩形）。ROI 内的空间平均温度（图 7b）在缓慢变化的直流偏移上呈现出稳定的 $2\omega$ 振荡，与规定的准稳态准则一致；用作反演输入的同相 $2\omega$ 振幅通过离散傅里叶变换从该序列中提取。 图 6. 三层体模的实验装置。堆叠的玻璃/Mylar/玻璃样品由交流电源驱动的 ITO 条带加热，全场表面温度由长波红外热像仪记录并传输至计算机进行 $2\omega$ 振幅提取。 图 7. (a) 和 (b)。

FLIR A655sc红外相机记录的代表性测量数据。（a）全场表面温度分布。黑色矩形框定感兴趣区域（ROI），即ITO加热条中央三分之一区域，在此区域上计算空间平均温度?Δ???。（b）ROI内空间平均温度序列?Δ???(??)，在直流偏移上叠加稳定的2??振荡。准稳态的定义为：在5 s观察窗口内，直流分量漂移低于5%。

**4. 结果与讨论**

**4.1 整体重建精度**

全部60个案例均无需先验厚度信息即可处理。八个参数从无噪声合成响应中顺序估计，逐例结果见附录A。第一层参数???和????,?从高频PTCS平台区估计，观测平均误差分别为0.07%和1.98%。本节其余部分聚焦于构成该方法核心评估的六个参数：{???, ???, ????,?, ???, ???, ????,?}。

累积误差分布（图8）显示出明显的阶段依赖模式：对于阶段2参数（???, ???, ????,?），分别有95.0%、83.3%和48.3%的案例误差低于10%阈值。对于阶段3参数（???, ???, ????,?），这些比例分别降至61.7%、20.0%和25.0%，与热波敏感性随深度增加而衰减的预期一致。总体而言，66.7%的参数-案例组合误差低于10%。

下载：下载高分辨率图片（696KB）下载：下载完整尺寸图片
图8. 所有60个数值测试案例中参数重建误差的全局经验累积分布函数（ECDF）。红色虚线垂直线标记10%误差阈值；绿色点线标记每个参数的第90百分位数。参数???和????,?（第一行）聚集在零附近，表明在此无噪声测试中表浅属性的条件敏感性较强。阶段2参数（???, ???, ????,?）分布紧凑，而阶段3参数（???, ???, ????,?）误差分布更宽，反映了深层光谱敏感性下降的趋势。

**4.2 递归提取的逐层分析**

两个提取阶段表现出不同的精度特征，误差从阶段2到阶段3递增，因为热波探测到介质更深处。

**阶段2：第一层厚度和第二层参数（1层到2层扩展）**

算法产生的{???, ???, ????,?}估计值，平均误差分别为2.26±2.48%、6.91±5.85%和10.28±10.04%（第90百分位数：4.60%、10.70%、14.93%）。

**阶段3：第二层厚度和第三层参数（2层到3层扩展）**

算法产生的{???, ???, ????,?}估计值，平均误差分别为8.92±7.25%、16.61±7.32%和17.94±10.81%（第90百分位数：17.94%、22.45%、28.27%）。

**生理搜索边界作为隐式正则化的作用**

软组织搜索范围（??≤0.60 W?m?1?K?1）充当基于物理的正则化：即使光谱信息不足以唯一确定，也将最深层估计约束在物理上可接受的范围内。在60个案例中的27个案例中，估计的???收敛到该上界；我们称之为边界约束案例。在其余33个案例中，优化器收敛到内点解（内点收敛案例）。在全部60个测试数值配置中均检测到第三层，随后估计???和????,?。阶段2误差在两个子集之间无显著差异（Welch ??检验：???、???和????,?的??值分别为0.530、0.096和0.383；表2）。在当前样本中，两个子集之间阶段2误差无统计学显著差异。

在内点收敛子集中，???的平均误差为11.84%。边界约束子集的???精度较低（11.89% vs. 6.49%，??=0.005），这反映了???与???之间的耦合：当最深层的光谱信息有限时，???的不确定性会传播到???。

下载：下载高分辨率图片（397KB）下载：下载完整尺寸图片
图9. 按反演阶段分组的所有60个数值测试案例中六个核心参数重建误差的小提琴图和条带图。蓝色阴影：阶段2参数{???, ???, ????,?}（1层到2层扩展）；红色阴影：阶段3参数{???, ???, ????,?}（2层到3层扩展）。各个案例以抖动条带点显示。箱体标记显示四分位距（IQR）、第10至90百分位须线、中位数（水平横线）和均值（菱形）。红色虚线标记10%误差阈值。阶段3分布系统性宽于阶段2，与热波衰减随深度增加的趋势一致。

下载：下载高分辨率图片（188KB）下载：下载完整尺寸图片
图10. 在全部60个自洽合成数据配置上，使用离散Fréchet度量评估的PTCS和原始信号目标函数对三个深层参数的敏感性。一个参数受到固定相对偏移的扰动，其他参数保持在真值；更大的目标值表示对该单次扰动的更陡峭条件惩罚，而非联合可辨识性的保证。箱图显示中位数、四分位距和第5至95百分位须线（省略异常值）。（a）深层厚度???，（b）深层电导率???，（c）深层体积热容????,?。在每个偏移量上，PTCS目标（蓝色）比原始信号目标（红色）高约一个数量级。

表2. 内点收敛案例（??=33）与边界约束案例（???达到生理上界，??=27）之间阶段2和阶段3重建误差的比较。均值±SD（%）。

| 子集 | 阶段2 | | | 阶段3 | | |
| | ??? | ??? | ????,? | ??? | ??? | ????,? |
| 内点（??=33） | 2.1 ± 2.7 | 5.7 ± 4.9 | 9.3 ± 11.4 | 6.5 ± 5.5 | 11.8 ± 6.8 | 17.1 ± 11.8 |
| 边界约束（??=27） | 2.5 ± 2.2 | 8.3 ± 6.6 | 11.5 ± 7.9 | 11.9 ± 8.0 | 22.5 ± 0.0 | 19.0 ± 9.3 |
| 全部（??=60） | 2.3 ± 2.5 | 6.9 ± 5.8 | 10.3 ± 10.0 | 8.9 ± 7.3 | 16.6 ± 7.3 | 17.9 ± 10.8 |
| ??值 | 0.530 | 0.096 | 0.383 | 0.005* | <0.001* | 0.488 |

*底部行报告双侧Welch ??检验??值；显著结果（??<0.05/6≈0.0083，Bonferroni校正）。

**与内点收敛的观测关联**

在测试网格中，与第一层厚度???的关联最强：???=50 μm的92%案例收敛到内点，因为较薄的第一层消耗更少的穿透预算，为更深层留下更多光谱范围。层间对比度???/???起次要作用，内点收敛率从???/???=0.41时的65%降至0.69–0.83时的50%。

各参数误差分布（图9）显示这种阶段依赖模式，阶段3箱体相对于阶段2向更高误差偏移。

**4.3 条件PTCS特征空间敏感性**

上述阶段3误差增长有一个可能的来源：信息到达优化器的PTCS特征空间。我们通过比较PTCS目标与基于原始温度信号ln[Re(?????)]构建的参考目标来分析这一点，后者保留了完整的未微分幅值；固定度量以隔离特征空间的影响。使用自洽合成数据（在真值处解析模型，因此不涉及前向模型差异），每个深层参数（???、???或????,?）受到固定相对偏移（?20%、+20%、+50%）的扰动，其他参数保持在真值。所得目标函数上升衡量条件单次敏感性；它本身不能建立三个参数同时变化时的联合可辨识性。

在全部60个案例中，PTCS目标施加的惩罚显著高于原始信号目标（图10）。对相同扰动的中位Fréchet响应大约大一个数量级，???的因子为9.6–18.4，???为34.9–43.7，????,?为27.8–28.8，两个分布在真值两侧的每个偏移量上均保持良好分离。该差异在????,?中依然存在，这是本测试中误差最大的参数（表2）。

在此条件比较下，PTCS中的微分相对于原始信号特征空间增加了所测试深层参数扰动的惩罚。在此特征变换后剩余的联合参数耦合在附录D中单独评估。

还进行了局部条件CRLB和分块条件Fisher信息分析，以表征测量信息的深度依赖性和阶段内参数耦合。其公式和完整结果见附录D。这些分析是局部测量级诊断，而非噪声数据下非线性反演的端到端验证。

**4.4 组织成分的影响**

为评估脂质核成分是否影响重建精度，对三种纤维含量组（0%、50%、75%）在匹配几何配置下进行了配对Welch ??检验。经Bonferroni校正18次比较（??_Bonf=0.05/18≈0.0028）后，仅有一项比较达到显著性：???在0%与75%之间（??<0.001；表3，图11）。另外两项???比较（0% vs. 50%，??=0.035；50% vs. 75%，??=0.036）在??=0.05下名义显著，但未通过校正。校正后其他参数未检测到统计学显著的成分差异（??>0.10）。

???敏感性反映了高纤维含量时界面对比度的减弱：当纤维分数接近75%时，Layer 1/2电导率比???/???=0.407/0.476=0.85，使界面在PTCS中变得不清晰，模糊了损耗突变边界。

表3. 六个核心参数配对成分比较的双侧独立样本Welch ??检验??值。

| 参数 | 0% vs. 50% | 0% vs. 75% | 50% vs. 75% |
|------|-----------|-----------|------------|
| ??? | 0.035 | <0.001* | 0.036 |
| ??? | 0.105 | 0.485 | 0.637 |
| ????,? | 0.777 | 0.660 | 0.951 |
| ??? | 0.816 | 0.774 | 0.941 |
| ??? | 0.413 | 0.456 | 0.966 |
| ????,? | 0.692 | 0.787 | 0.350 |

*显著性在Bonferroni校正水平??_Bonf=0.05/18≈0.0028（6参数×3比较）下评估；标记为显著结果。

下载：下载高分辨率图片（186KB）下载：下载完整尺寸图片
图11. 按脂质核纤维含量分组的六个核心参数平均重建误差（±1 SD）（蓝色：0%，红色：50%，绿色：75%）。水平虚线标记10%误差阈值。

**4.5 几何配置的影响**

在此有限测试设计中，未检测到重建误差与规定层厚度之间的统计学显著关联。将每个参数的误差与???和???相关后，经Bonferroni校正（图12）无显著关联，在???–???网格上的平均误差热力图（图13）在20个测试几何配置中也未显示明显的视觉系统性盲区。这些空结果不能建立测试网格之外与几何无关的结论。

下载：下载高分辨率图片（628KB）下载：下载完整尺寸图片
图12. 六个核心参数（???, ???, ????,?, ???, ???, ????,?）重建误差相对于规定层厚度???（左列）和???（右列）的散点图。所有60个案例（三种脂质核成分合并）以均匀蓝色符号显示。实线显示普通最小二乘回归；阴影带表示95%置信区间。显著性在Bonferroni校正水平??_Bonf=0.05/12≈0.0042下评估；无面板达到该阈值；在测试样本中未找到统计学可检测的线性关联。

下载：下载高分辨率图片（657KB）下载：下载完整尺寸图片
图13. 六个核心参数在???–???几何空间中的平均重建误差热力图（阶段2：???, ???, ????,?，第一行；阶段3：???, ???, ????,?，第二行）。每个单元格报告跨三种脂质核成分的平均误差。色标为参数特定（绿色：低误差；红色：高误差）。在测试参数网格（???=50–250 μm，???=250–1000 μm）中未识别到明显的视觉系统性几何盲区。

**4.6 三层模型评估**

将分步程序应用于第3.2节描述的三层玻璃/聚酯薄膜/玻璃模型。在三次实验中，在20个从1000降至10 mHz的激发频率下测量了同相2??幅值。测量幅值及其与三层解析预测的偏差见附录B。该程序过渡到三层模型并在每次实验中返回第二边界???的估计值。

表4报告了三次实验的平均估计参数以及运行间标准偏差，后者衡量每个参数在重复组装和测量下的可重复性。完整相对误差以平均误差±运行间标准偏差报告，而非仅以平均值报告。在FLIR A655sc的指定测量精度下，表4中的实验反演误差总体可接受。然而，????,?的运行间标准偏差为43.9%，远大于其他参数。

表4. 三层模型（浮法玻璃/聚酯薄膜/浮法玻璃）的估计参数，报告为三次独立实验的均值±运行间标准偏差，与参考值比较。

| 阶段 | 参数 | 参考值 | 平均估计值 | 平均相对误差/SD（%） |
|------|------|--------|-----------|---------------------|
| 1 | ??? (W?m?1?K?1) | 1.137 | 1.139 | +0.2 ± 2.4 |
| | ????,? (MJ?m?3?K?1) | 1.001 | 0.969 | ?3.2 ± 4.1 |
| 2 | ??? (μm) | 400 | 353 | ?11.8 ± 5.2 |
| | ??? (W?m?1?K?1) | 0.155 | 0.157 | +1.1 ± 11.9 |
| | ????,? (MJ?m?3?K?1) | 1.629 | 1.721 | +5.6 ± 18.4 |
| 3 | ??? (μm) | 300 | 314 | +4.7 ± 12.6 |
| | ??? (W?m?1?K?1) | 1.137 | 1.177 | +3.5 ± 11.4 |
| | ????,? (MJ?m?3?K?1) | 1.001 | 1.070 | +6.9 ± 43.9 |

**4.7 计算成本与局限性**

完整的三层反演（阶段1–3）在标准工作站（Intel Core i7-12700，2.10 GHz）上运行约90 s。分阶段分解大致为阶段1约1 s，阶段2约35 s，阶段3约54 s。由于每次扫描步骤中??_pop=30个DE候选者独立评估，预计8线程并行执行可将运行时间减少至约15 s。实测串行运行时间和估计并行运行时间在原则上可能适用于离线或预处理分析；本研究未评估临床适用性和患者特异性性能。

当前工作存在以下局限性：

(i) 无法实时体内操作。当前采集和串行反演流程速度不够快，无法在治疗过程中持续追踪热物性参数的变化。

所报告的运行时间仅适用于离线或预处理分析。

（二）缺乏活体验证。在本研究中，构建并集成一个足够小且具有可控组织接触条件的表面加热器-传感器组件超出了研究范围。因此，验证仅限于合成组织模型和分层仿体。

（三）概念验证仿体对比度。玻璃- Mylar-玻璃样品提供了明显的层间对比度，用于测试逐步过程。其第1层/第2层电导率比约为1.137/0.155=7.3，远低于早期玻璃-XPS双层样品的约38倍比值，但仍高于合成组织构型的1.2-2.4范围。因此，该实验仅支持所测试构型下的可行性，而非通用临床性能。

（四）忽略界面热接触电阻。正演模型在每个内部界面处强制执行完美热接触。组装后的仿体中有限的界面电阻或不完美耦合的界面可能因此会在估计的层参数中引入偏差。

（五）深度相关信息与噪声敏感性。数值评估使用了无噪声的合成数据。附录D中的条件Cramér-Rao和块条件FIM分析提供的是局部测量信息诊断，而非端到端的含噪声反演验证。尽管三层实验误差在所给定的相机精度下总体可接受，但CV,3的运行间标准差为43.9%，表明该参数的可重复性降低。

**5. 结论**

本研究使用合成生物组织模型和分层仿体对逐步热波反演进行了评估。在所测试的三层构型中，单次表面频率扫描即可估计层属性和边界，无需预先指定层厚度。多层正演模型本质上具有递归性，每次层扩展恰好增加三个未知量；这将完整的3N?1参数反演分解为一系列低维、有界的三参数子问题。在所测试的微扰比较中，基于导数的伪热导率谱（PTCS）产生了比近乎重叠的原始温度曲线更大的条件目标变化，且损失突变准则检测到了所测试三层构型中的界面，引导了从一个阶段到下一阶段的顺序扩展。

在60种三层动脉粥样硬化斑块构型上的无噪声数值评估表明，在施加的生理约束下，第1、2、3阶段的平均误差分别低于2%、10.5%和18%。三层玻璃/Mylar/玻璃仿体实验在所给定的红外相机测量精度下产生了总体可接受的误差；然而，CV,3的运行间标准差为43.9%。

在所测试的构型范围内，结果表明该方法有望降低对外部提供的层几何信息在边界检测和受限属性估计方面的依赖。在考虑使用导管兼容的薄膜3ω传感器进行血管内实施之前，需要进一步评估受控测量不确定性、额外的独立可观测量、经验证的先验信息以及具有生物代表性的仿体。

**CRediT作者贡献声明**

Zixun Zhan：原始草稿撰写、可视化、验证、软件、方法论、调查研究、形式分析、数据整理。Shiqing Zhao：审阅与编辑、验证、指导、方法论、调查研究、资金获取、概念构思。Aili Zhang：审阅与编辑、指导、项目管理、资金获取、概念构思。

**科学写作中生成式人工智能的声明**

在本文准备过程中，作者使用了Claude（Anthropic）来提高英文散文的清晰度和可读性。使用此工具后，作者对内容进行了审查和编辑，并对发表文章的内容承担全部责任。

Zixun Zhan | Shiqing Zhao | Aili Zhang