功能性梯度皮肤组织中瑞利波传播的全球敏感性分析：一种分数阶三相滞后热粘弹性模型

《Materials Today Communications》：Global Sensitivity Analysis of Rayleigh Wave Propagation in Functionally Graded Skin Tissue: A Fractional Three-Phase Lag Thermo-Viscoelastic Model

　　马阿兹·阿里·汗 | 卡迪贾·M·阿布阿尔纳贾 | 阿德南·贾汉吉尔 | 埃马德·E·马哈茂德 | 乌斯曼·里亚兹 | 穆拉特·亚亚拉克
巴基斯坦白沙瓦库尔塔巴科学与信息技术大学物理与数值科学系

**摘要**
本文提出了一个完整的分数阶三相滞后热粘弹性模型，用于研究功能梯度的人体皮肤组织中的瑞利波传播。所提出的框架结合了非局部弹性、分数阶记忆效应以及深度依赖的材料梯度——这是将多物理场理论与生物波建模相结合的一个关键缺失环节。此外，本文使用基于方差的Sobol方法进行了全局敏感性分析（GSA），以量化主要参数（即弹性非局部性（?1）、热非局部性（?2）、分数阶数（α）、梯度参数α*、相位滞后τq、τT、τv以及静水应力（P）对相位速度、衰减、穿透深度和比热损失的影响。所有相位速度均以参考速度c0=10m/s（皮肤典型值）进行归一化，以获得无量纲结果；相应的物理值范围为100?1000m/s，与弹性成像测量结果一致。研究结果表明，波的特性主要受热非局部性（?2）控制（特别是在热损失方面），穿透深度受分数阶数（α）影响，而相位速度和衰减则是所有参数复杂协同作用的结果。GSA为模型简化及参数优先级提供了有力的支持，有助于优化临床应用中的诊断弹性成像和热疗技术。

**1. 引言**
弹性动力学作为连续介质力学的一个分支，研究弹性材料在动态载荷下的复杂行为，特别是波的传播和材料的时变响应。该学科对于阐明机械波（包括体波和表面波）在固体介质中的传播和相互作用至关重要。弹性动力学理论的发展已有150多年历史，其长波长领域的理论公式已通过实验得到验证。如今，在工程和科学的许多领域中，动态弹性理论都是不可或缺的，应用于机械振动、冲击载荷、地球物理学和生物物理学等领域[1]。
表面波（如瑞利波、洛夫波和斯通利波）的传播是弹性动力学领域的重要研究课题，特别是在生物组织（如人类皮肤）中的应用[2]。这些波在材料界面或表层内部表现出独特的传播特性。例如，瑞利波涉及表面的椭圆颗粒运动，对表面性质具有高度敏感性，因此成为无创诊断方法（如超声弹性成像）的宝贵工具。洛夫波是水平偏振的剪切波，存在于两种固体之间，其传播特性取决于这两种介质的性质。

在热弹性理论的发展过程中，对其进行了不断的改进。虽然经典弹性理论适用于等温条件，但它忽略了热效应——这在处理生物组织时是一个严重的缺陷。Biot[3]基于傅里叶的热传导定律和运动方程发展了热弹性耦合理论。然而，该模型仍假设热波具有无限速度。为了解决这一问题，Cattaneo[4]和Vernotte[5]提出了双曲热传导模型，确保了热波的有限传播速度。进一步的推广包括具有一个松弛时间的Lord–Shulman（LS）模型[6]以及具有两个松弛时间的Green–Lindsay（GL）模型[7]。Tzou[8]、[9]进一步发展了双相位滞后（DPL）模型，随后又发展了三相位滞后（TPL）模型，这些模型考虑了微观结构效应的相互作用。同时，Green和Naghdi GN提出了基于熵平衡和热位移的三热弹性模型GN-I、GN-II、GN-III[10]、[11]、[12]。Roy Choudhuri[13]将这些概念扩展为完整的三相位滞后理论，考虑了由于热流、温度梯度和热位移引起的滞后。

**2. 数学公式化**
2.1. 具有非局部弹性的本构关系
皮肤组织被建模为一个均匀的、各向同性的、具有非局部效应的热粘弹性半空间。根据Eringen的非局部弹性理论[27]，热粘弹性材料的本构关系可以表示为：
$$
(1)
(1-\epsilon_1^2\nabla^2)\sigma_{ij}=-P(\delta_{ij}+\omega_{ij}+\{2\mu\bar{e}_{ij}+(\lambda^\bar{e}\nabla r-\gamma^\bar{T})\delta_{ij}
$$
其中：
• $\sigma_{ij}$ 是柯西应力张量
• $e_{ij}=\frac{1}{2}(u_{ij}+u_{ji})$ 是应变张量
• $u_{i}$ 是位移分量
• $T$ 是相对于参考温度$T_0$的温度变化
• $\delta_{ij}$ 是克罗内克delta函数
• $\epsilon_1$ 是考虑长程相互作用的弹性非局部参数
• $P$ 是静水应力
• $\omega_{ij}=\frac{1}{2}(u_{ji}-u_{i}$ 是旋转张量，表示介质的局部刚性体旋转

我们使用的粘弹性效应表达式为：
$$
(2)
^* (\lambda^\bar{\mu}\nabla\gamma^\bar{k}+\lambda^\bar{k}^*)R+i\omega(\lambda^\bar{\mu}\nabla k^*)I
$$

2.2. 运动方程
忽略体力作用时的运动方程为：
$$
(3)
\sigma_{ij,j}=\rho u_{i潜能}
$$
其中 $\rho$ 是组织密度，$u_{i}$ 是位移分量，点表示时间微分。

2.3. 广义生物热方程
考虑灌注和代谢效应的生物组织能量守恒方程为：
$$
(4)
\overline{\rho_t}\frac{c_t}{c}\frac{\partial T}{\partial t}+\gamma^\bar{T}_0\nabla\cdot\frac{\partial u}{\partial t}=-\nabla\cdot q+\omega_b\rho_b c_b T_b-T+Q
$$
其中：
• $\rho_t, c_t$ 是组织密度和比热
• $\gamma^\bar{=}\gamma+\omega\gamma^*$ 是复杂的热弹性耦合常数
• $q$ 是热流矢量
• $\omega_b, \rho_b, c_b$ 是血液灌注率、密度和比热
• $T_b$ 是动脉血液温度
• $Q$ 是代谢产生的热量

2.4. 具有分数阶导数的三相位滞后热传导
热流矢量定义为：
$$
(5)
q=-(\k\nabla T+k^*\nabla v)
$$
其中 $k^*$ 是热导率。根据Green–Naghdi II型热弹性理论，热导率$k^*$表示为：
$$
k^*=c_t\lambda+\frac{2\mu}{4}
$$
$ct$ 是热波速度，$\lambda$ 和 $\mu$ 是拉梅常数。为了同时考虑局部和非局部热效应，热传导方程可写为：
$$
(6)
\overline{\nabla^2}\frac{1-\epsilon_2^2}{\nabla^2}q=-(\k^*\nabla T+k^*\nabla v)
$$
其中 $\epsilon_2$ 是热非局部参数。

为了纳入记忆效应，我们使用Caputo分数阶导数（$0<\alpha\leq1$）：
$$
(7)
D_{\alpha}c_f(t)=1\Gamma(1-\alpha)\int_0^t f'(\tau)\tau^{-\alpha}d\tau
$$
带有Caputo导数的分数阶TPL热传导方程和生物热方程为：
$$
(7)
\tau_q^{*}(1-\epsilon_2^2\nabla^2)(\rho_t c_t\frac{\partial^2 \theta}{\partial t^2}+\gamma^\bar{T}_0\nabla\cdot\frac{\partial^2 u}{\partial t^2}+\omega_b\rho_b c_b\frac{\partial \theta}{\partial t)}=k^*\left[(\tau_k+\tau_k\tau_T+\tau_v\right)\nabla^2 \theta'+\nabla^2 \theta]
$$
其中：
• $\tau_q$ 是热流相位滞后
• $\tau_T$ 是温度梯度相位滞后
• $\tau_v$ 是热位移梯度相位滞后
• $k, k^*$ 是热导率系数

3. 问题解决方法
3.1. 问题陈述与假设
引入笛卡尔坐标系（$x, y, z$），其中皮肤表面位于$z=0$，下层组织占据半无限区域$z\geq0$。假设瑞利波沿正x方向在表面上传播。

在功能梯度介质中，材料属性不是均匀的，而是随位置平滑变化。这种空间变化在研究表面波时尤为重要，因为沿传播方向的轻微不均匀性就会显著影响应力分布和波的特性。与大多数假设深度依赖梯度的研究不同，我们考虑了沿传播方向x的横向梯度。这种选择在伤口愈合、肿瘤生长或妊娠纹等情况下具有生理学意义，因为这些情况下皮肤的硬度会显著变化。如果忽略这种梯度，得到的波模型可能无法准确捕捉介质的实际机械响应[38]。

条件（8）将（1）、（3）、（7）修改为：
$$
(9)
(1-\epsilon_1^2\nabla^2)\sigma_{ij}=-P(\delta_{ij}+\left[2\mu\bar{0}e_{ij}+(\lambda^\bar{0}\nabla r-\gamma^\bar{0}T\right)\delta_{ij}\times g(x)
$$
$$
(10)
\sigma_{ij,j}=\left(1-\epsilon_1^2\nabla^2\right)\times\left[\left(-P(\delta_{ij}+\omega_{ij}\right)+\left(2\mu\bar{0}e_{ij}+(\lambda^\bar{0}\nabla r-\gamma^\bar{0}T\right)\delta_{ij}\right]\times g(x)
$$
$$
(11)
\tau_q^{*}(1-\epsilon_2^2\nabla^2)(\rho_t c_t\frac{\partial^2 \theta}{\partial t^2}+g(x)\gamma^\bar{T}_0\nabla\cdot\frac{\partial^2 u}{\partial t^2}+\omega_b\rho_b c_b\frac{\partial \theta}{\partial t})=g(x)k_0\left[(\tau_k+\tau_k\tau_T+\tau_v\right)\nabla^2 \theta'+\nabla^2 \theta]
$$
为了以简单可行的方式考虑这种效应，假设材料参数沿x方向通过指数梯度函数变化。因此，功能梯度表示为$g(x)=e^\alpha x$，其中$\alpha^*$是控制空间变化率的不均匀性参数。这种选择保持了分析的简洁性，同时有效考虑了材料梯度对瑞利波传播的影响。测试表明，幂律梯度的敏感性排名在质量上几乎相同（相位速度的差异<5%）。

由于平面应变假设，所有机械和热场变量被认为与y坐标无关。因此，位移矢量限制在$x-z$平面内，表示为：
$$
(12)
u=(u,0,w)(x,z,t)
$$
使用方程（12）代入方程（9-11），得到：
$$
(13)
\left(-P+(\lambda^\bar{0}+\mu^\bar{0})\times e^\alpha x\right)\frac{\partial^2 u}{\partial x^2}+\left((\lambda^\bar{0}+\mu^\bar{0})\alpha \times e^\alpha x\right)\frac{\partial u}{\partial x}+\left(-P^2+(\lambda^\bar{0}+\mu^\bar{0})e^\alpha x\right)\frac{\partial^2 w}{\partial x\partial z}+\alpha^\alpha \lambda^\bar{0}e^\alpha x\frac{\partial w}{\partial z}\left[-P^2+2\mu^\bar{0}e^\alpha x\right]\frac{\partial^2 u}{\partial z^2}-\rho\epsilon_{12}\frac{\partial^4 u}{\partial x^2}-\rho\epsilon_{12}\frac{\partial^4 u}{\partial z^2}\left[-P^2+2\mu^\bar{0}e^\alpha x\right]\frac{\partial^2 u}{\partial z^2}+\alpha^\alpha \lambda^\bar{0}e^\alpha x\frac{\partial w}{\partial z}+\alpha^\alpha 2\mu^\bar{0}e^\alpha x\frac{\partial u}{\partial z}\left[-P^2+(\lambda^\bar{0}+2\mu^\bar{0}e^\alpha x\right]\frac{\partial^2 w}{\partial z^2}-\gamma^\bar{0}e^\alpha x\frac{\partial \theta}{\partial z}\right]=\rho\frac{\partial^2 u}{\partial t^2}-\rho\epsilon_{12}\frac{\partial^4 u}{\partial x^2}-\rho\epsilon_{12}\frac{\partial^4 w}{\partial z^2}\frac{\partial t^2}\tau_q^{*}\left[\rho_t c_t\frac{\partial^2 \theta}{\partial t^2}+\gamma^\bar{T}_0\left(\frac{\partial^3 u}{\partial t^2}+\frac{\partial^3 w}{\partial t^2}\frac{\partial x+}\frac{\partial^3 w}{\partial t^2}\right)\times e^\alpha x+\omega_b\rho_b c_b\frac{\partial \theta}{\partial t}\right]-\tau_q^{*}\epsilon_{22}\left[\rho_t c_t\frac{\partial^4 \theta}{\partial t^2}+\alpha^\bar{\gamma^\bar{T}_0\left(\frac{\partial^5 u}{\partial t^2}+\frac{\partial^5 w}{\partial t^2}\frac{\partial x^2}\right)\times e^\alpha x+\gamma^\bar{T}_0\left(\frac{\partial^5 u}{\partial t^2}+\frac{\partial^5 w}{\partial t^2}\frac{\partial x^2}\right)e^\alpha x+\omega_b\rho_b c_b\frac{\partial \theta}{\partial t}\right]-\tau_q^{*}\epsilon_{22}\left[\rho_t c_t\frac{\partial^4 \theta}{\partial t^2}+\alpha^\bar{\gamma^\bar{T}_0\left(\frac{\partial^5 u}{\partial t^2}+\frac{\partial^5 w}{\partial t^2}\frac{\partial x^2}\right)e^\alpha x+\omega_b\rho_b c_b\frac{\partial \theta}{\partial t}\right]
$$

3.2. 问题求解
假设形式为[33]的时谐解：
$$
(16)
\Xi=\Xi^*e^{i\xi(x+mz-\ct)}
$$
其中：
• $\Xi=(u,w,\theta)$ 是复波动数（$\xi=\xi_R+i\xi_I$），$c$ 是相位速度，$m$ 控制深度衰减。衰减定义为：沿传播方向x的空间衰减为$\alpha_{att}=\xi_I$（nepers/m）。图中绘制的是无量纲衰减因子$AF=\xi_I/\xi_R$。穿透深度为$\delta=\frac{1}{|\xi_I|}$。

3.3. 以势能为基准的控制方程
将方程（16）代入方程（13–15）并消除应力，得到三个耦合方程：
$$
(17)
\begin{align}
(a_1+a_2m^2)u^{*}+a_3m^w^{*}-a_4\theta^{*}=0 \\
(a_5m^2)u^{*}+a_6+a_7m^2)w^{*}-a_8m\theta^{*}=0 \\
(a_9+a_{10}m^2)u^{*}+a_{9}m+a_{10}m^3)w^{*}+a_{11}+a_{12}m^2)\theta^{*}=0
\end{align}
$$
为了完整性，系数$a_1\cdots a_{12}$在附录I中给出。

3.4. 特征方程
对于非平凡解，系数矩阵的行列式必须为零：
$$
(20)
A_m^6+B_m^4+C_m^2+D=0
$$
特征方程中的系数$A\cdots D$在附录I中给出。

3.5. 通解形式
方程（20）得到六个根$m_j$（$j=1,2,3,4,5,6$），其中三个根满足$R(m_j)\geq0$，对应于随深度衰减的表面波。通解为：
$$
(21)
\Xi=(1\delta_j l_j)A_j e^{i\xi(x-c_t)
$$
其中$A_j$是根据边界条件确定的常数：
$$
(22)
d_j=-a_{14}m^3-a_{13}m^{16}m^2+a_{15} \\
(23)
l_j=a_{17}m^6+a_{18}m^4+a_{19}m^2+a_{20}a_{21}m^4+a_{22}m^2+a_{23}a_{13}=a_1a_8-a_4a_5; \quad a_{14}=a_8a_2; \quad a_{15}=-a_4a_6; \quad a_{16}=a_8a_3-a_4a_7; \quad a_{17}=a_{14}a_{10}; \quad a_{18}=-a_{16}a_{10}+a_{13}a_{10}+a_{14}a_9; \quad a_{19}=-a_{15}a_{10### 地程求解与应用边界条件后得到一个齐次系统：
$$
\begin{cases}
\sum_{j=1}^{3} k_{ij} A_j = 0, \\
K_{1j} = (-P - (\lambda_{\bar{0} + 2\mu_{\bar{0}} + P)(i\xi)m_j d_j e^{\alpha^{\star}x} + \lambda_{\bar{0}}i\xi e^{\alpha^{\star}x} - \gamma_{\bar{0}}l_j e^{\alpha^{\star}x) = 0, \\
K_{2j} = (-P^2 + e^{\alpha^{\star}x} - \mu_{\bar{0}})(m_j + d_j)(i\xi) = 0, \\
K_{3j} = (-i\xi m_j + h)l_j = 0.
\end{cases}
$$
对于非零的 $ A_j $，行列式必须为零：
$$
\det(K_{ij}) = 0
$$
方程 (30) 就是该方程，其解给出了色散关系 $ c = c(\xi) $ 和其他波动特性。

### 全局敏感性分析 (GSA)
全局敏感性分析用于量化不确定输入参数对皮肤组织中热粘弹性波响应的影响。与仅检查标称点附近扰动的局部敏感性方法不同，GSA 探索了全部允许的参数空间，因此更适合非线性和耦合的多物理系统，在这些系统中参数相互作用起着重要作用 [39]。
在本研究中，采用了基于方差的 Sobol 框架将模型输出的总方差分解为来自各个输入参数及其高阶相互作用的贡献。这种方法能够全面评估主要效应和相互作用效应。

### 4.1. 问题表述
设模型输出表示为
$$
Y = f(X_1, X_2, \ldots, X_n)
$$
其中 $ X_i, i = 1, 2, \ldots, n $ 表示假设在单位超立方体 $ [0, 1]^n $ 上均匀分布的独立输入参数。在可积性的假设下，函数 $ f(X) \in L^2([0, 1]^n) $，可以使用 Sobol–Hoeffding 展开将其分解为：
$$
f(X) = f_0 + \sum_{i=1}^n f_i(X_i) + \sum_{i 0.02 $ 且 $ p < 0.01 $，则认为相互作用显著。

### 表 1. 仿真中使用的皮肤组织属性
| 参数 | 符号 | 值 | 单位 |
|----------------|------------------|-----------------|-------------------|
| 组织密度 ρ_t | $ \rho_t $ | 1190 | Kg/m^3 |
| 组织密度 ct | $ \rho_c $ | 3770 | J/kg |
| 热导率 k* | $ k^* $ | 8.063 \times 10^{11} | W/mK |
| 分数阶参数 α | $ \alpha $ | 0 ≤ α ≤ 1 | |
| 剪切模量 μ | $ \mu $ | 4.456 \times 10^6 | Pa |
| Lame 常数 λ | $ \lambda $ | 9.15 \times 10^7 | Pa |
| 血液比热 cb | $ c_b $ | 3650 | J/kg |
| 血液密度 ρ_b | $ \rho_b $ | 10^6 | kg/m^3 |
| 组织导电率 k | $ k $ | 0.246 | W/mK |
| 动脉温度 Tb | $ T_b $ | 310 | K |

### 5. 数值结果和讨论
#### 5.1. 材料参数
- 弹性非局部性：$ \epsilon_1 \in [0.1, 0.7] $
- 热非局部性：$ \epsilon_2 \in [0.1, 0.7] $
- 分数阶：$ \alpha \in [0, 0.1] $
- 相位滞后：$ \tau_q, \tau_T, \tau_v \in [0.1, 1.0] $ s
- 频率范围：$ \omega \in [1, 100] $ MHz（临床弹性成像范围）
- 相位速度归一化：所有相速度均通过 $ c_0 = 10 \, m/s $（放松状态下的人体皮肤中的典型剪切波速度）进行归一化。物理相速度恢复为 $ c = c^\* \times c_0 $。例如，$ c^\* = 40.5 $ 对应于 450 m/s。未经归一化的原始计算值为 100–1000 m/s，与文献一致。
- 非局部参数的物理基础：弹性非局部参数 $ \epsilon_1 = e_0a $，其中 $ a $ 是胶原纤维间距（10–100 μm），$ e_0 \approx 0.39 $；因此 $ \epsilon_1 $ 在 0.1–0.7 范围内对应于物理上的 0.01–0.07 mm。热非局部参数 $ \epsilon_2 $ 表示热载体（声子或间质液）的平均自由路径，也在 10–100 μm 的数量级。这些值与皮肤微观结构一致。

使用 MATLAB 软件，我们图形化地展示了不同分数参数 $ \alpha $、非局部性参数 $ \epsilon_1, \epsilon_2 $ 对相速度、比热损失、穿透深度和衰减系数的变化。

#### 图 2. 随非局部弹性参数 $ \epsilon_1 $ 和热参数 $ \epsilon_2 $ 变化的相速度
相速度色散曲线展示了在具有分数三阶段滞后热弹性的功能梯度皮肤组织中瑞利波的频率依赖传播。随着角频率 $ \omega $ 的增加，相速度的一致减小证实了这种复杂介质的强色散特性，这是由于材料梯度、粘弹性阻尼和非局部热记忆的共同作用。参数研究表明，增加分数阶数（$ \epsilon_1 $ 或 $ \epsilon_2 $）会显著降低整个频率范围内的波速。从物理上看，这意味着更强的分数记忆效应（与热和机械滞后相关）增强了能量耗散机制，从而减缓了波的传播。实际上，对于高频超声弹性成像等诊断技术来说，这意味着一个基本的权衡：更高的频率提供了更好的分辨率，但由于这些内在的组织特性，会导致更大的速度减小和衰减。因此，准确的体内表征需要模型精确捕捉这些分数阻尼和梯度参数，以便将波测量结果与组织健康状况相关联。

#### 图 3. 随非局部弹性参数 $ \epsilon_1 $ 和热参数 $ \epsilon_2 $ 变化的衰减
衰减谱描述了弹性非局部性和热非局部性参数在控制功能梯度皮肤组织中瑞利波耗散中的主要作用。随着角频率 $ \omega $ 的增加而单调增加，表明高频波周期更广泛地与组织的非局部微观结构相互作用，从而导致更大的能量损失。任意一个非局部参数的增加都会显著增加衰减，表明弹性和非热部门内的空间相互作用增强了波能量的扩散耗散。这些参数的存在直接通过引入超出经典局部理论的空间记忆效应来影响波与组织的耦合。实际上，这为高频诊断波的穿透深度设定了一个严格的限制，因为更高的分辨率被强烈的频率依赖性衰减所抵消。因此，为了准确地进行体内组织表征，需要精确量化 $ \epsilon_1 $ 和 $ \epsilon_2 $。

#### 图 4. 随非局部弹性参数 $ \epsilon_1 $ 和热参数 $ \epsilon_2 $ 变化的穿透深度
该图表示了在具有分数三阶段滞后模型下，由弹性非局部参数 $ \epsilon_1 $ 和热非局部参数 $ \epsilon_2 $ 决定的瑞利波在功能梯度皮肤组织中的穿透深度损失。从图中可以看出，当初始深度较浅时，穿透深度的下降非常陡峭，表明由于组织的粘弹性阻尼和非局部相互作用，波能量紧密地限制在表面附近。从曲线中可以推断，增加 $ \epsilon_1 $ 或 $ \epsilon_2 $ 会进一步减少有效穿透深度，表明由于更高的非局部性导致空间耗散增加。从物理上解释，可以确定分数非局部算子与明显的深度依赖性衰减相关，这限制了基于表面的诊断波的探测深度。对于实际弹性成像，更高的频率或具有更强非局部特性的组织会大大降低成像深度。因此，选择正确的波频率和准确估计 $ \epsilon_1 $ 和 $ \epsilon_2 $ 对于平衡分辨率和穿透深度至关重要。

#### 图 5. 随非局部弹性参数 $ \epsilon_1 $ 和热参数 $ \epsilon_2 $ 变化的声波损耗
从图中可以看出，由于瑞利波传播的比热损失与角频率 $ \omega $ 成线性关系，与非局部参数 $ \epsilon_1 $ 和 $ \epsilon_2 $ 成非线性关系。随着频率的增加，热量损失的增加验证了快速的周期性变形导致能量更有效地转换为其热形式，因为粘弹性和热机械耦合更为有效。这两个非局部参数的较大值导致更高的耗散，从物理上讲意味着在弹性和非热层面有更大的空间相互作用负责每个波周期的更高能量产生。这种关系在热基疗法中很重要，因为可控的波频率可以引起局部加热，在诊断成像中需要最小化热损失以防止组织损伤。因此，正确建模这些参数对于提高基于波的生物医学应用的安全性和有效性至关重要。从物理角度来看，这意味着遗传阻尼和结构异质性都是限制基于表面基波的波能量传播深度的因素。在弹性成像应用方面，这导致了一个重要的基本权衡：通过提高频率来增加弹性成像的空间分辨率会导致成像深度的直线下降。因此，诊断协议的优化需要平衡频率选择以及组织固有的衰减特性和记忆效应。下载：下载高分辨率图像（144KB）下载：下载全尺寸图像图9. 分数阶α和渐变参数α*的表面声波（SHL）。比热损失谱定义了依赖于分数阶α和渐变参数α*的瑞利波的热机械能量转换效率。随着频率的增加，热损失的增强表明了组织对粘弹性和热记忆的不可逆工作的加速。分数阶α的增加强调了这种耗散，因为提升的时间非局域性延长了热松弛过程，从而增加了每个周期的累积热量产生。渐变参数α*的变化改变了材料属性在空间中的分布，从而影响了机械能量随深度的差异吸收及其转化为热量的过程。从物理上讲，这意味着热量产生不仅是一个频率依赖的属性，还与组织的结构——即梯度和遗传响应——内在相关。在实际的治疗应用中，如可控热疗或热消融，这种关系允许通过调整波频率并根据组织分级来定位热剂量。另一方面，在诊断成像领域，减少这种加热对于安全性至关重要，例如，这需要创建预测通过这两个参数的热量损失的模型。图10解释了静水压力（P）对功能梯度皮肤组织中瑞利波传播相位速度的强烈依赖性。随着（P）值的增加，色散曲线的上升位移支持了这样一个假设：施加的压缩应力使介质变硬，导致所有频率范围内的波传播速度增加。这种行为可以通过静水压力对组织粘弹性矩阵的压缩作用来解释，这种压缩作用减少了微观孔隙率，同时提高了有效弹性模量。从物理角度来看，这些结果突出了静态应力配置与动态波特性之间的关键相互作用，特别是在分数阶三维滞后模型的背景下。在临床弹性成像的背景下，结果表明波速的诊断评估并不是组织本身的固有属性，而是直接受到探针外部压力或内在预应力（例如来自进展中的病变）的影响。因此，为了准确地在体内表征组织，需要纳入局部应力场，以便正确解释速度数据，并将明显的硬化归因于病理变化或机械加载。下载：下载高分辨率图像（121KB）下载：下载全尺寸图像图10. 带有静水压力P的PV。下载：下载高分辨率图像（117KB）下载：下载全尺寸图像图11. 带有静水压力P的AF。衰减系数谱显示，随着静水压力（P）的增加，波能量耗散显著减少，从而证明了静态应力状态与动态耗散机制之间的直接耦合。这种衰减减少可以通过压缩应力对粘弹性组织基质的影响来解释，这种压缩作用减少了负责将波能量转化为热量的微观摩擦和内部散射点。在分数阶三相滞后范式内，这种压缩随后降低了材料每个周期的遗传阻尼能力。从物理上讲，这意味着更大的、更具弹性的、耗散更少的介质是一个预应力介质。在实际弹性成像中，探针压力可以通过最小化衰减来增加信号的穿透深度及其清晰度，但这会改变被测量组织的内在属性。因此，诊断协议必须包括标准化或其他校正措施，以区分真实的病理衰减和由机械加载引起的人为修改/减少（即病理变化）。下载：下载高分辨率图像（120KB）下载：下载全尺寸图像图12. 带有静水压力P的PD。穿透深度谱显示，当静水压力（P）进入组织时，瑞利波的穿透深度显著增加。这种效应是由于压缩应力增加了粘弹性基质的有效硬度，从而减少了材料的固有阻尼能力，并降低了波能量随深度的减少率。从物理角度来看，施加的压力减少了微观空隙和内部摩擦，使得传播路径更具弹性且耗散更少。在实际弹性成像标准中，这一事实提出了一个重要方面：诊断探针通常施加的接触压力可能会人为增加穿透深度，从而干扰组织的真实衰减特性。因此，使用标准化的测量协议或应力补偿模型在区分皮肤内在结构属性和测量过程中外部影响下的机械调制效果方面发挥着不可或缺的作用。下载：下载高分辨率图像（145KB）下载：下载全尺寸图像图13. 带有静水压力P的SHL。图表显示，静水压力P的增加系统性地减少了瑞利波传播过程中的比热损失，从而表明压缩应力是一种降低粘弹性滞后和内部摩擦的方式，这些因素负责将波能量转化为热量。从机制上讲，施加的应力用于压缩组织的多孔、耗散性微观结构，从而有利于组织的更弹性响应和较少地耗散能量。实际上，这一观察表明，在弹性成像过程中探针接触压力会极大地影响组织感受到的热负荷，这对治疗性超声来说是一个重要的安全问题。因此，在基于波的诊断和治疗中，必须明确考虑局部应力状态，以便准确预测热剂量，防止在未受应力组织区域产生局部热量。下载：下载高分辨率图像（80KB）下载：下载全尺寸图像图14. 涉及模型参数对输出QoI的影响。5.2. 全局敏感性分析结果全局敏感性分析定量评估了所有模型参数对功能梯度皮肤组织中瑞利波传播的影响。结果确定热相位滞后τT和分数阶α是最重要的单个因素，它们控制着波的耗散和记忆依赖行为。热图忽略了非局部参数?1、?2和渐变参数α*之间的显著非线性相互作用，这表明它们对波特性的综合效应是协同的，而不是叠加的。相比之下，如静水压力（P）以及松弛时间τq和τv等参数的一阶敏感性指数相对较小。这种层次结构为模型简化和实验设计提供了重要路线图，最重要的结论是：正确表征热记忆和分数阶动力学是模型预测准确性的前提条件。下载：下载高分辨率图像（283KB）下载：下载全尺寸图像图15. 重要参数的排序。5.2.1. 参数排序和因子固定全局敏感性分析显示，参数排序和固定显示了不同的、依赖于输出的层次结构，这有助于简化模型和实验设计。相位速度对弹性非局部参数?1（24%）和粘弹性相位滞后τv（21%）最为敏感，这揭示了空间记忆和粘性松弛效应组合的重要性。此外，分数阶α（11%）和热非局部性?2（13%）也具有相当的敏感性。在速度的初始近似中，可以将敏感性较低的参数赋予其标准值。衰减因子受到几乎所有参数的相对均匀贡献（每个参数12%-13%），反映了复杂的协同机制，没有主导贡献。因此，不建议进行因子固定，而需要同时校准所有参数以实现准确的衰减建模。穿透深度主要由分数阶α（72%）决定，使其成为最重要的参数。热滞后τT（17%）是次要因素，而所有其他因素的影响都小于5%。为了简化专门用于深度预测的模型，可以固定这些参数。比热损失主要由热非局部参数?2（48%）控制，因此将其确定为热机械能量转换的主要控制因素。渐变参数α*（15%）起次要作用，而所有其他参数（如相位滞后τq、τv、τT和弹性非局部性?1）的个体影响可以忽略不计（<11%）。在预测热损失时，这些参数可以固定在名义值上。5.3. 额外验证5.3.1. 对梯度函数的敏感性我们比较了指数梯度eα*x和幂律1|βx)n在相同刚度变化范围内的表现。所有频率下的相位速度相差<5%，Sobol指数变化<1%（在蒙特卡洛误差范围内）。因此，指数形式是稳健的。5.3.2. 可变相位滞后的影响尽管τq、τT、τv可能取决于水合作用或深度，但全局敏感性分析表明τT对穿透深度有中等影响（17%），而τq、τv的敏感性较低（<5%）。因此，恒定滞后假设对于排名目的来说是可接受的。5.3.3. 各向异性预应力初步的双轴应力（Px,Pz）测试表明，敏感性排序保持不变；P的总效应指数从<5%增加到约7%，仍然低于主导参数（α、?2、?1）。6. 与文献和实验数据的比较和验证为了评估所提出的分数阶三相滞后热粘弹性模型的预测能力，我们将我们的归一化相位速度结果与现有的实验测量和文献中的理论预测进行了比较。本研究中的所有相位速度都通过参考速度c0=10 m/s进行了归一化，得到无量纲值cˉ=c/c0。我们模拟得到的相应物理相位速度范围为100–1000 m/s，具体取决于非局部参数（?1、?2）、分数阶α、渐变参数α*和静水压力P的值。6.1. 与实验弹性成像数据的比较Nenadic等人[26]使用机械振动器和超声跟踪技术报告了人类前臂皮肤的表面波速为8–35 m/s。Han等人[25]通过光学相干弹性成像测量了猪皮肤的15–45 m/s。Grinspan等人[21]在骨骼肌中获得了10–40 m/s的速度，其机械性质与皮肤相似。我们的物理速度（100–1000 m/s）高于这些典型值，因为我们的模型包括了几种在简单弹性成像解释中经常被忽视的硬化机制：•静水预应力P（高达20 kPa）压缩组织，增加了其有效弹性模量。•分数阶α<0.1引入了强烈的记忆效应，提高了高频动态模量。•非局部弹性（?1、?2）增加了空间相互作用，特别是在较高频率（MHz范围内）时可以提高波速。•横向梯度（α*）改变了阻抗剖面，也影响了测量到的速度。当我们将P设为0，α设为0.1（弱记忆），?1=?2=0.1（小非局部性），α*=0（无梯度）时，我们的模型恢复了20–50 m/s的相位速度，这与实验范围非常吻合。因此，我们的框架可以通过调整这些生理相关的参数来再现报告的皮肤波速的完整谱。更高的值（100–1000 m/s）对应于老化皮肤、纤维化组织或高频治疗超声等条件，其中硬化和记忆效应变得占主导。6.2. 与理论模型和极限情况的比较在经典极限（?1=?2=0，α=1，α*=0，τq=τT=τv=0，且没有粘弹性）下，我们的方程（30）简化为各向同性弹性半空间的众所周知的瑞利波色散关系。在这种情况下，无量纲相位速度cˉ=c/c0接近经典值cR/cS=0.9194（其中cS=μ/ρ是剪切波速度）。在我们的归一化假设下，c0=10 m/s，并使用名义剪切模量μ=4.456×10^6 Pa和密度ρ=1190 kg/m3，cS≈61.2 m/s，因此经典的瑞利波速度约为56.3 m/s，对应的cˉ≈5.6。当关闭所有非局部效应和分数效应时，我们的数值求解器能够将这个值再现到0.5%的范围内，验证了我们实现方法的正确性。我们的结果也与Ezzat等人[14]的分数热弹性模型以及Biswas和Mukhopadhyay[22]及Tiwari和Kumar[23]的非局部TPL研究结果相符。在没有功能梯度和侧向不均匀性的情况下，我们的色散曲线与[22]中报道的非局部热弹性半空间的色散曲线相匹配。侧向梯度的加入（α*>0）使色散曲线向上移动（相位速度增加），并引入了额外的频率依赖性，这与Kumar等人[43]和Khan等人[38]的分级波分析结果一致。

6.3. 衰减和热损失预测的验证
在MHz频率下直接实验验证人体皮肤的衰减和比热损失更具挑战性，因为难以在体内隔离热耗散和粘弹性耗散。然而，我们计算出的衰减因子（AF = I(ξ)/R(ξ)）在1–100 MHz频段内的范围为0.05到0.5。这些值与软生物组织中测量的损耗正切值一致：例如，Madsen等人[41]报告肝脏和肾脏的衰减系数为0.5–1.0 dB/cm/MHz，这对应于10 MHz时的AF值约为0.1–0.2。我们的比热损失值（10–500 W/m3在10 MHz）位于低强度治疗超声的报告范围内[42]。热损失对?^2的强依赖性（48%敏感性）是一个新的预测，需要使用具有预设热非局部长度的受控 phantom 材料进行实验测试。

7. 临床意义
本研究的结果对临床诊断和治疗应用有几项直接的影响，这些影响源自全局敏感性分析（第5.2节）和参数研究（图2、图3、图4、图5、图6、图7、图8、图9、图10、图11、图12、图13）。

7.1. 弹性成像解释和参数优先级
传统的弹性成像主要根据剪切模量来解释波速。我们的全局敏感性分析（GSA）显示，相位速度对弹性非局部参数?1（24%）和粘弹性相位延迟τv（21%）最为敏感，而不仅仅是局部剪切模量。因此，临床系统应该设计为从多频率色散测量中提取?1和τv。?1的增加可能表明胶原纤维间距增大（例如，在水肿或老化的皮肤中），而τv的增加则表明粘弹性松弛作用增强。这些参数可以作为纤维化、硬皮病或淋巴水肿早期诊断的新生物标志物。

7.2. 穿透深度和频率选择
穿透深度主要由分数阶α（72%敏感性）决定，热延迟τT（17%）也有次要贡献。这意味着组织的记忆依赖性粘弹性行为是限制瑞利波探测深度的主要因素。对于无创成像，选择操作频率需要权衡：更高的频率可以提高空间分辨率，但会降低穿透深度。我们的结果显示，如果α较低（记忆性强，如年轻健康的皮肤），穿透深度已经很浅，因此进一步提高频率几乎没有额外好处。相反，如果α较高（记忆性弱，例如在光损伤的皮肤中），更高的频率仍可能提供足够的深度。因此，临床医生应首先估计α（例如，通过低频蠕变测试或从衰减曲线的斜率来确定），然后再选择弹性成像频率。

7.3. 热疗安全性和加热控制
比热损失主要由热非局部参数?^2控制（48%），其次是梯度参数α*（15%）。对于治疗性超声（加热、消融），热扩散的空间范围（非局部长度）比经典热导率更为重要。这意味着治疗计划应结合术前成像或患者皮肤类型的?^2估计（例如，较厚的真皮层会增加?^2）。为了避免意外的热损伤，应用频率和强度应调整，以使热损失保持在安全阈值以下（例如，对于表层治疗，<100 W/m3）。我们的模型提供了一个定量地图：对于给定的?^2和α*，可以预测热损失变得过高的频率。

7.4. 探针压力标准化
静水应力P对相位速度、穿透深度和比热损失的敏感性较低（<5%），对衰减的影响也仅约为7%。因此，在弹性成像过程中探针接触压力的小变化不太可能显著改变测量得到的波速或深度，只要压力保持在10 kPa以下。这简化了临床协议，减少了严格压力控制的需要，这在实践中往往很困难。然而，在比较不同操作者或设备的测量结果时，仍然建议记录施加的压力并在必要时使用我们的模型进行校正。

7.5. 基于模型的诊断和逆问题
我们的集成框架可以反向应用：给定多频率的相位速度、衰减和穿透深度测量数据（例如来自临床弹性成像系统），可以通过解决逆问题（例如使用贝叶斯优化或神经网络）来估计关键参数?1、α、?^2、α*。这将使得无需活检就能无创地表征皮肤微观结构（胶原间距、记忆效应、热非局部长度和侧向刚度梯度）。这样的工具可以作为软件模块集成到商用超声或光学相干弹性成像设备中。

7.6. 个性化热疗计划
热损失对?^2的强依赖性表明，应使用患者的特定?^2值来制定热疗计划。例如，?^2较高的患者（例如，皮肤厚且散射强烈的患者）在同一声学强度下会经历更大的局部加热；因此应降低功率。相反，?^2较低的患者（皮肤薄且散射较少的患者）可能需要更高的功率来达到相同的热剂量。我们的模型为这种个性化调整提供了合理的依据。

8. 结论
本研究通过开发一个集成的分数三相延迟热粘弹性模型，成功解决了已识别的研究空白，该模型适用于功能梯度皮肤组织中的瑞利波传播。该模型同时考虑了非局部弹性、分数记忆效应、空间不均匀性（侧向梯度）和热机械耦合，这是之前未曾报道过的组合。数值模拟揭示了?1、?^2、α、α*在调节波色散、衰减、穿透和热生成方面的不同作用。全局敏感性分析首次提供了参数影响的定量排名：分数阶α主导穿透深度（72%），热非局部性?^2主导比热损失（48%），弹性非局部性?^1主导相位速度（24%）。衰减结果是由复杂的协同作用产生的，没有单一的主导参数。至关重要的是，静水应力P对大多数输出的影响较低（<5%），表明探针压力效应次于组织的内在记忆和非局部特性。这些发现为基于波的组织表征提供了新的理论基础，并为弹性成像和热疗提供了实际指导。未来的工作将包括深度依赖的相位延迟、各向异性预应力以及在病理组织中的实验验证。

未引用的参考文献[28]CRedi
作者贡献声明：
Murat Yaylac?：数据整理、调查、资源管理、验证、撰写——审稿与编辑。
Usman Riaz：监督、方法论、调查、形式分析、数据整理、概念化、撰写——审稿与编辑。
Maaz Ali Khan：撰写——初稿、可视化、验证、软件开发、项目管理、形式分析、数据整理、撰写——审稿与编辑。
Adnan Jahangir：可视化、验证、监督、软件开发、资源管理、方法论、调查、撰写——审稿与编辑。
Khadijah M. Abualnaja：概念化、形式分析、方法论、验证、撰写——审稿与编辑。
Emad E. Mahmoud：形式分析、方法论、软件开发、可视化、撰写——审稿与编辑。