引言

生物组织、多孔岩体、土壤、颗粒物质、胶体等复杂介质(complex media)，是介于理想固体和液体之间的复杂物质状态，其力学和物理性质复杂，通常表现出记忆性和路径依赖性，难以用一般的经典力学模型来描述.此外，复杂介质中的能量耗散往往呈现出幂律特征，即其试验数据的拟合公式具有幂函数形式^[1-2].

大量的室内和室外实验发现，声波在复杂介质中传播时，其衰减系数呈现出频率的幂律依赖关系^[3-7]

$ \alpha =\alpha _0 \omega ^\eta $

(1)

其中，$\omega $为角频率，$\alpha $为衰减系数，$\alpha_{0}$和$\eta$是实验数据拟合得到的参数, $\eta$介于0和2之间.

Szabo等^[8]给出了一组不同材料的频率与衰减系数之间关系的实验数据，如图 1所示.图 1将衰减系数和频率之间关系刻画在双对数坐标上，直线的斜率即为参数$\eta$.对于理想固体和液体，$\eta $为0或者2.例如，水中超声波耗散的$\eta $值为2，图 1所示的YIG (钇铟石榴石)作为一种理想固体，其剪切波和纵波的$\eta $值都为2.对于复杂介质，其力学性质介于理想固体和液体之间，所以参数$\eta $在0和2之间变化.例如，肝脏纵波在1$\sim$100 MHz内$\eta =1.3$，蒸馏水和石油的$\eta$值都在1和2之间，海底沉积物对应的$\eta$为1.

气体中能对振动起到阻尼作用的只有黏滞性和热传导.经典的声波控制方程包括阻尼波动方程和热黏性波动方程两类.阻尼波动方程的形式如下

$\dfrac{1}{c_0^2 }\dfrac{\partial ^2P}{\partial t^2} + \dfrac{2\alpha _0 }{c_0 }\dfrac{\partial P}{\partial t} = \nabla ^2P $

(2)

其中，$\alpha_{0}$为衰减系数，$c_{0}$为波传播速度，$P$为声场变量的振幅，例如，声速或者声压.该模型可以由经典的Maxwell黏弹性本构关系推导而来，其描述了与频率无关的声衰减.

声波在黏性流体中传播时，其控制方程的形式如下^[9-10]

$\dfrac{1}{c_0^2 }\dfrac{\partial ^2P}{\partial t^2} - \dfrac{\mu }{c_0^2 }\dfrac{\partial }{\partial t}\left( {\nabla ^2P} \right) = \nabla ^2P $

(3)

其中，$\mu $为材料参数.该波动方程可由经典的Kelvin黏弹性本构关系推导得到.如果考虑低频近似，即$\left( {\mu \omega / c_0^2 } \right)^2 \ll 1$，则式(3)可以退化为一种近似热黏性波方程^[10]

$\dfrac{1}{c_0^2 }\dfrac{\partial ^2P}{\partial t^2} - \dfrac{\mu }{c_0^2 }\dfrac{\partial ^3P}{\partial t^3} = \nabla ^2P $

(4)

该模型描述了频率二次方依赖的声衰减，即$\alpha \propto \omega ^2$.

从以上公式可以发现，经典的阻尼波方程和近似热黏性波方程分别描述了与频率无关和频率二次方依赖的声衰减，但无法描述频率非整数次依赖的衰减.为了研究声波衰减的频率任意次幂律依赖现象，许多学者提出了相应的模型.例如，基于频域分析的声波模型^[11-18]，多松弛模型^[19-20]等.然而，这些模型由于需要消耗大量的计算成本^[21-23]而未被广泛接受.

分数阶导数算子已被广泛研究，其是一种能描述具有历史记忆性和空间依赖性力学行为的有力建模工具^[24-25].采用分数阶导数建模，以修正传统的达西定律、傅里叶热传导、黏弹性本构关系，可以很好地解释传统模型无法解释的``反常''现象.此外，分数阶导数的阶数有明确的物理意义.例如，扩散方程中的分数阶阶数可以区分正常扩散，超扩散和次扩散^[26]；黏弹性本构方程中的分数阶阶数可以定性描述黏性和弹性的比重^[27-28].类似地，分数阶导数模型也可以被用于描述任意阶频率幂律依赖衰减的现象.目前，已有许多学者在这方面做了大量的工作，提出了相应的分数阶导数声波模型.本文旨在对这些已有的模型作一个总结，从时间和空间分数阶两个方面概述了不同的分数阶导数声波模型，并讨论了相关模型之间的联系和区别，阐述了频率依赖指数变化范围的统计解释，最后对尚未涉及的声波传播耗散机制的建模领域提出了研究设想，期望能起到抛砖引玉的效果.

1 在软物质中传播的分数阶导数声波模型 1.1 时间分数阶导数声波模型

时间分数阶导数声波模型大致可分为唯象模型和本构模型两类.唯象模型主要是基于声衰减的频率幂律依赖现象，采用分数阶微积分理论修正已有的声波方程而获得的；本构模型的构建主要是从分数阶黏弹性本构方程推导得到.

1.1.1 唯象模型

基于频率幂律依赖现象，Szabo^[11]从热黏性波模型和阻尼波模型的频域形式出发，提出了声波模型的一般化频域模型，然后将频域声波模型转化为时域上的卷积积分声波方程

$\dfrac{1}{c_0^2 }\dfrac{\partial^2 P}{\partial t^2} + \dfrac{2\alpha _0 }{c_0 } {\rm L}_\eta \left( P \right) = \nabla ^2P $

(5)

其中，耗散项${\rm L}_{\eta}(P)$为一个卷积算子.

由于Szabo提出的声波模型含有超奇异性，数值上难以处理，Chen等^[29]引入分数阶微积分理论对其进行修正，推导出了一种修正的Szabo模型，其控制方程为

$\dfrac{1}{c_0^2 }\dfrac{\partial^2 P}{\partial t^2} + \dfrac{2\alpha _0 }{c_0 \cos \left( {\eta \pi / 2} \right)}\dfrac{\partial ^{\eta + 1}P}{\partial t^{\eta + 1}} = \nabla ^2P $

(6)

其中，$\eta$为分数阶导数的阶数，其范围在0和2之间.不难发现，修正的Szabo模型介于阻尼波模型($\eta=0$)和近似黏性波模型($\eta = 2$)之间，即这两种整数阶导数的波动方程是该分数阶声波方程的特例.紧随其后，Liebler等^[30]提出了一种求解修正Szabo方程的有效方法，数值结果与实验测量数据能很好地吻合. Chen等^[31]和Zhang等^[32]分别将此类声波控制方程用于二维和三维的医学超声成像，以识别肿瘤的位置和大小.如图 2所示，肿瘤的位置和大小与声压曲线凸起的位置和大小相对应.

Szabo模型和修正的Szabo模型都满足频率幂律依赖的衰减$\alpha = \alpha _0 \omega^\eta$.然而，这种关系只是在频率比较小的时候才能满足.为了使控制方程对于任意频率都能满足频率幂律依赖关系，Kelly等^[33]提出了分数阶导数幂律声波模型

$\begin{gathered} \frac{1}{{c_0^2}}\frac{{{\partial ^2}P}}{{\partial {t^2}}} + \frac{{2{\alpha _0}}}{{{c_0}\cos \left( {\pi \eta /2} \right)}}\frac{{{\partial ^{\eta + 1}}P}}{{\partial {t^{\eta + 1}}}} + \\ \frac{{\alpha _0^2}}{{{{\cos }^2}\left( {\frac{{\pi \eta }}{2}} \right)}}\frac{{{\partial ^{2\eta }}P}}{{\partial {t^{2\eta }}}} = {\nabla ^2}Pn \\ \end{gathered} $

(7)

并成功采用格林函数法求解了该方程.从图 3可以看出，幂律声波在脂肪组织($\eta=1.5$)中的衰减明显比在肝脏($\eta=1.139$)中大. Meerschaert等^[34]从随机行走的角度给出了分数阶导数项的物理解释和方程(7)的统计解.

不同于Szabo的研究方式，Wismer^[35]从频率的幂律依赖现象出发，直接对经典的黏性波动方程进行修正，得到了一种修正的声波控制方程

$\dfrac{1}{c_0^2 }\dfrac{\partial^2 P}{\partial t^2} - \tau^{\eta - 1} \dfrac{\partial^{\eta - 1}}{\partial t^{\eta - 1}}\left( {\nabla ^2P} \right) = \nabla ^2P $

(8)

其中，分数阶阶数$<\eta-1<1$，$\tau^{\eta- 1}$为松弛时间常数.该声波方程的衰减系数$\alpha \propto \omega ^\eta $.

Ochmann等^[36]根据幂律衰减现象推导得到了一类广义的Burgers方程

$ \dfrac{1}{c_0^2 }\dfrac{\partial ^2P}{\partial t^2} + \left( {\eta + 1} \right)\nabla \left( {\nabla ^2P} \right) + \dfrac{2D}{c_0 }\dfrac{\partial ^{\eta + 1}P}{\partial t^{\eta + 1}} = \nabla ^2P $

(9)

其中，$D$为常数，$P$为质点移动速度，该模型可与以下一维非线性波动方程统一起来^[37]

$ \dfrac{\partial ^2P}{\partial t^2} + c_0^2 \left( {\gamma + 1} \right) \dfrac{\partial ^3P}{\partial x^3} + 2Dc_0 \dfrac{\partial ^{\gamma + 1}P}{\partial t^{\gamma + 1}} = c_0^2 \dfrac{\partial ^2P}{\partial x^2} $

(10)

该非线性模型中，$P$为质点振幅，$\gamma $为比热容比.

1.1.2 本构模型

与理想固体和流体不同，软物质的力学行为具有时间和空间上的长程相关性，其通常可以被当作黏弹性材料研究.分数阶微积分是一种能描述黏弹性材料力学行为的强有力工具.目前，这一观点已得到广泛认可^{[1, 25, 38-43]}.从分数阶黏弹性本构出发，研究黏弹性材料中的声波耗散机制，是不同于唯象建模的另一种研究思路.

Caputo^[44]早在1967年即从分数阶Kelvin黏弹性本构方程推导出分数阶Kelvin波动方程，其形式如下

$\dfrac{1}{c_0^2 }\dfrac{\partial^2 P}{\partial t^2} - \tau _\sigma ^\eta \dfrac{\partial ^\eta }{\partial t^\eta }\left( {\nabla ^2P} \right) = \nabla ^2P $

(11)

其中，$\eta$为分数阶阶数，且$0<\eta <1$，$\tau _\sigma ^\eta $为表征材料物理性质的一个常数.其衰减系数和频率之间的关系可见图 4(a)^[45-48].需要注意的是，图中的$\alpha $即为本文的分数阶阶数的表示符号$\eta$；图中的$\alpha_{k}$即为本文的衰减函数$\alpha $.

从形式上看，分数阶Kelvin模型与Wismer提出的模型(式(8))可以记作一个统一的形式，其分数阶阶数在0和1之间(Wismer模型的分数阶阶数为$\eta-1$, $1 <\eta<2$).经典的热黏性波动方程(式(3))，是该统一形式的一个特例($\eta=1$). Caputo等^[49]采用有限差分法分析了该模型在医学超声上的应用.最近，Zhao等^[50]采用格林函数法比较了修正Szabo模型(式(6))，幂律声波模型(式(7))和分数阶Kelvin波动模型(式(11))，并用数值方法证明了在输入条件为脉冲函数时，这3种模型几乎是等价的，如图 5所示.需要注意的是，图中的Szabo模型即本文的修正Szabo模型.

类似地, Konjik等^[51]、Näsholm等^[52]都推导了分数阶Zener波动方程

$ \dfrac{1}{c_0^2 }\dfrac{\partial^2 P }{\partial t^2} - \tau _\sigma ^\eta \dfrac{\partial ^\eta }{\partial t^\eta }\left( {\nabla ^2P} \right) + \dfrac{\tau _\varsigma ^\beta }{c_0^2 }\dfrac{\partial ^{\beta + 2}P}{\partial t^{\beta + 2}} = \nabla ^2P $

(12)

其中，$\eta$和$\beta $为两个不同的分数阶阶数，且$0 <\eta$，$\beta<1$，$\tau _\sigma ^\eta $和$\tau _\varsigma ^\beta $为表征材料物理性质的常数. 图 4(b)反映了其衰减系数和频率之间的关系^[45-48]. Holm等^[52]分析了分数阶Kelvin模型和Zener模型的衰减系数的逼近关系.此外，多松弛模型、幂律衰减和本构型时间分数阶波动方程是相互等价的^{[46, 53]}.

Buckingham^[54-55]发现在海底沉积物中传播的纵波和横波同样具有频率的幂律依赖现象.因此，他采用一种具有应变硬化性质的修正Maxwell模型来描述颗粒之间的剪切机制.该学者将时间依赖的脉冲形状函数应用到Navier-Stokes方程中，得到波动方程如下

$\left. \begin{gathered} {\nabla ^2}\psi = \frac{1}{{c_0^2}}\frac{{{\partial ^2}\psi }}{{\partial {t^2}}} - \hfill \\ \left( {\frac{{{\lambda _{\text{p}}}}}{{{\rho _0}c_0^2}}} \right)\frac{\partial }{{\partial t}}{\nabla ^2}\left( {{h_{\text{p}}}\left( t \right) * \psi } \right) - \hfill \\ \frac{4}{3}\frac{{{\gamma _{\text{s}}}}}{{{\rho _0}c_0^2}}\frac{\partial }{{\partial t}}{\nabla ^2}\left( {{h_{\text{s}}}\left( t \right) * \psi } \right) - \hfill \\ \frac{{{\eta _{\text{s}}}}}{{{\rho _0}}}{\nabla ^2}\left( {{h_{\text{s}}}\left( t \right) * \boldsymbol{A}} \right) - \frac{{\partial \boldsymbol{A}}}{{\partial t}} = 0 \hfill \\ \end{gathered} \right\}$

(13)

其中，速度向量$v = \nabla \psi + \nabla \times {\boldsymbol {A}}$，$*$表示卷积，$\psi$是速度势函数，${\boldsymbol {A}}$是向量场函数. $h\left( t \right) = \dfrac{\theta }{\xi _0 }\left( {1 + \dfrac{\theta }{\xi _0 }t} \right)^{ - E / \theta }$，$E$为弹性模量，$\xi_{0}$为初始黏滞系数，$\theta $为黏滞系数随时间的变化率，下标s和p分别代表剪切波和纵波，$\rho_{0}$为材料的体积密度，$\lambda_{\rm p}$和$\gamma_{\rm s}$分别为纵波和剪切波对应的应力松弛系数. Pandey等^[56]将上述方程中的纵波和横波参数统一为一个固定参数(不区分s和p)，并以逆幂律函数$h\left( t \right) = \tau ^{ - 1}\left( {\dfrac{t}{\tau }} \right)^{ - \eta }$逼近$h\left( t \right) = \tau ^{ - 1}\left( {1 + \dfrac{t}{\tau }} \right)^{ - \eta }$，得到含有分数阶导数项的纵波方程

$\begin{gathered} {\nabla ^2}\psi = \frac{1}{{c_0^2}}\frac{{{\partial ^2}\psi }}{{\partial {t^2}}} - \\ \Gamma \left( {1 - \eta } \right)\left( {\frac{{{\lambda _{\text{p}}}}}{{{\rho _0}c_0^2}} + \frac{4}{3}\frac{{{\gamma _{\text{s}}}}}{{{\rho _0}c_0^2}}} \right){\tau ^{\eta - 1}}\frac{{{\partial ^\eta }}}{{\partial {t^\eta }}}{\nabla ^2}\psi \\ \end{gathered} $

(14)

和横波方程

$\dfrac{\partial ^{2 - \eta }{\boldsymbol {A}}}{\partial t^{2 - \eta }} = \Gamma \left( {1 - \eta } \right)\dfrac{\gamma _{\rm s} }{\rho _0 }\tau ^{\eta - 1}\nabla ^2 {\boldsymbol {A}} $

(15)

通过比较式(11)与式(14)，可以发现，这两个波动方程具有相似的形式，而式(15)是一类分数阶扩散-波方程.如此，Pandey等^[56]将分数阶Kelvin波动方程与海底沉积物中的声波传播方程紧密地联系起来，并赋予分数阶阶数以物理意义$\eta=E/\theta $.然而，Pandey等^[56]在推导分数阶纵波方程时，将横波和纵波的松弛时间设定为同一参数，这一做法不可取.因此，本文作者认为采用推广的分数阶Kelvin模型来描述该问题比较合适，即两个Soctt-Blair元件与一个弹簧元件并联^[25].基于此推广的分数阶Kelvin模型得到的波动方程为

$ \dfrac{1}{c_0^2 }\dfrac{\partial^2 P }{\partial t^2} - \tau _{\rm p}^{\eta _1 } \dfrac{\partial ^{\eta _1 }}{\partial t^{\eta _1 }}\left( {\nabla ^2P} \right) -\tau _s^{\eta _2 } \dfrac{\partial ^{\eta _2 }}{\partial t^{\eta _2 }}\left( {\nabla ^2P} \right) = \nabla ^2P $

(16)

其中，$\eta_{1}$和$\eta_{2}$分别为描述纵波和剪切波松弛系数的分数阶阶数，$\tau _{\rm s}^{\eta _2 } $和$\tau _{\rm p}^{\eta _1 } $为两个材料常数.如此，可避免在推导式(14)时采用剪切波与纵波的$h(t)$相同的假设.

Prieur等^[57]采用分数阶Kelvin本构模型，与分数阶热传导模型相结合，得到的波动方程如下

$\begin{gathered} {\nabla ^2}P - \frac{1}{{c_0^2}}\frac{{{\partial ^2}P}}{{\partial {t^2}}} + {L_{\text{v}}}\frac{{{\partial ^{\eta - 1}}}}{{\partial {t^{\eta - 1}}}}{\nabla ^2}P - \frac{{{L_{\text{t}}}}}{{c_0^2}}\frac{{{\partial ^{\eta + 1}}P}}{{\partial {t^{\eta + 1}}}} = \hfill \\ \;\;\;\;\;\;\; - \frac{\beta }{{{\rho _0}c_0^4}}\frac{{{\partial ^2}P}}{{\partial {t^2}}} \hfill \\ \end{gathered} $

(17)

式中，$\eta$为分数阶阶数，其他参数可参见文献[57]，等式左边第三项和第四项为衰减项，等式右边为非线性项.低频条件下，这个方程的衰减系数为$\alpha \propto \omega ^\eta $，即该方程的声衰减也是频率幂律依赖的.另外，此方程可以看作是Westervelt方程、Burgers方程和Khokhlov-Czabolotskaya-Ckuznetsov方程的一个推广.

1.2 空间分数阶声波方程

声波传播时，其衰减过程不仅与传播的时间相关，也与传播的距离和位置相关.受启发于热黏性波动方程的衰减项$\dfrac{\partial }{\partial t}\left( {\nabla ^2P} \right)$和近似热黏性波动方程的衰减项$\dfrac{\partial ^3P}{\partial t^3}$之间的低频近似的等价关系，Chen等^[58]提出了一种空间分数阶拉普拉斯声波模型

$\dfrac{1}{c_0^2 }\dfrac{\partial^2 P}{\partial t^2} + \dfrac{2\alpha _0 }{c_0^{1 - \eta } }\dfrac{\partial }{\partial t}\left( { - \nabla^2 } \right)^{\eta / 2}P = \nabla ^2P $

(18)

其中，$\eta$为分数阶阶数，$\left( { - \nabla^2 } \right)^{\eta / 2}$为分数阶拉普拉斯算子.特别当$\eta=0$时，该方程退化为阻尼波模型(式(2)); 当$\eta=2$时，该模型即为热黏性波方程(式(3)).张晓棣^[59]用数值方法研究了空间分数阶Laplace声波模型与幂律声波模型、修正Szabo声波模型和分数阶Kelvin声波模型之间的区别，如图 6所示.当衰减指数($\eta=1.7$)比较大时，4种模型的波形几乎相同，即此时4种模型的衰减、速度频散几乎相同，且与声源中心频率的大小无关(如图 6(a)和图 6(b))；当衰减指数($\eta=1.3$)比较小时，幂律声波模型、修正Szabo声波模型和分数阶Kelvin声波模型之间的差异随着声源中心频率的增大而愈发明显.另外，由于空间分数阶拉普拉斯声波模型不能描述声波的频散，因此其传播速度明显慢于其他模型(如图 6(c)和图 6(d)).需要注意的是，图中的time fractional wave equation即本文的分数阶Kelvin模型.

空间分数阶拉普拉斯声波方程满足频率幂律依赖关系$\alpha = \alpha _0 \omega^\eta $.但是，如上文所述，该方程描述的波速是无频散的，所以Treeby等^[60]提出了一种既满足因果律(Kramers-$\!$-Kronig关系)又能刻画速度频散的色散关系$k = \dfrac{\omega }{c_0 } - \dfrac{\alpha _0 \left( { -{\rm i}} \right)^{\eta + 1}\omega ^\eta }{\cos \left( {\pi \eta / 2} \right)}$，其中$k$为复波数，进而推导出了修正的空间分数阶导数耗散声波模型

$\begin{gathered} \frac{1}{{c_0^2}}\frac{{{\partial ^2}P}}{{\partial {t^2}}} + \frac{{2{\alpha _0}}}{{c_0^{1 - \eta }}}\frac{\partial }{{\partial t}}{\left( { - {\nabla ^2}} \right)^{\eta /2}}P + \hfill \\ \;\;\;\;\;\;2{\alpha _0}c_0^\eta \tan \frac{{\pi \eta }}{2}{\left( { - {\nabla ^2}} \right)^{(\eta + 1)/2}}P = {\nabla ^2}P \end{gathered} $

(19)

传统的非局部弹性理论采用权函数表征物质点之间力学性质的相互影响关系.随着距离的增大，权函数逐渐减小并趋于0.此外，权函数与材料内部特征长度相关，特征长度越小，则权函数影响范围越小，反之亦然.常见的权函数有格林函数、高斯函数等，形式上表现为一类指数函数.

分数阶导数是一个非局部算子，在时间上表征历史依赖性，在空间上刻画长程相关性.分数阶黏弹性模型是最常见的时间分数阶模型之一，大量的实验也证明了这类模型的有效性.目前，已有研究学者将经典的非局部弹性理论与分数阶黏弹性模型相结合，提出了非局部的分数阶黏弹性理论.

若不考虑时间域上的异同，仅讨论不同位置的物质点对某一点力学性质的影响，理论上可以采用空间分数阶导数来进行非局部弹性理论的建模.然而，经典的分数阶导数的权函数是一个幂函数，且与特征长度无关.

庞国飞^[61]在幂函数中引入特征长度对其进行修正，并与Riesz势的概念有机统一起来，提出了空间依赖的分数阶非局部算子.该分数阶非局部算子的定义为

${{\text{I}}^\gamma }u\left( x \right) = \left\{ \begin{gathered} \int {\frac{A}{{\Gamma \left( l \right)}}{{\left| {x - y} \right|}^{ - n + \gamma \left( l \right)}}u\left( y \right){\text{d}}y} \hfill \\ \;\;\;\;\;\;\;\left| {x - y} \right| \leqslant l{\mkern 1mu}, \gamma < n \hfill \\ 0, \;\;\;\;\;\left| {x - y} \right| > l \end{gathered} \right.$

其中，$A$为归一化常数, $l$为材料内部特征尺寸，$\gamma =n-1+kl$，$n=1$，2，3代表维数，常数$k$满足$kl <1$，且$\gamma <n$以保证函数的衰减性.所提出的权函数满足衰减性，非负性，归一性等.数值实验表明，采用该权函数的本构模型能很好地解释应力集中问题.基于此定义，该作者提出了空间非局部Maxwell和Kelvin黏弹性本构模型，并导出了相应的非局部Kelvin黏弹性波方程

$\dfrac{1}{c_0^2 }\dfrac{\partial ^2P}{\partial t^2} - \hat {\tau } {\rm I}^\gamma \dfrac{\partial }{\partial t}\left( {\nabla ^2P} \right) = \nabla ^2P $

(20)

和非局部修正Maxwell黏弹性波方程

$ \dfrac{1}{c_0^2 }\dfrac{\partial ^2P}{\partial t^2} + \dfrac{\rho }{\eta } {\rm I}^{ - \gamma } \dfrac{\partial P}{\partial t} + \dfrac{1}{c_0^2 } {\rm I}^{ - \gamma _1 } \dfrac{\partial ^2P}{\partial t^2} = \nabla ^2P $

(21)

这两种方程对应的衰减系数分别为$\alpha \propto \omega ^{2 - \gamma }$和$\alpha \propto \omega ^\gamma $，$0 < \gamma <1$，满足频率的幂律依赖关系.

1.3 时间和空间分数阶波动方程的联系

Treeby等^[62]利用低频近似的色散关系$k=\omega /c_{0}$将时间分数阶和空间分数阶波动方程有机结合起来.该学者将分数阶Kelvin波动方程(式(11))作傅里叶变换后，代入色散关系，再作傅里叶逆变换，可推导得到修正的空间分数阶波动方程(式(19)).这类分数阶波动方程可通过$k$波分裂法进行求解^[63].

Carcione等^[64]根据频率依赖的体积模量, 得到频域上的应力应变关系，并推导得到了一类扩散-波方程

$\dfrac{\partial ^{2 - 2\eta }P}{\partial t^{2 - 2\eta }} = c_0^2 \omega _0^{ -2\eta } \nabla ^2P $

(22)

随后，Carcione^[65-66]采用傅里叶伪谱法对其进行了数值分析.类似地，Zhu等^[67-68]基于黏弹性本构得到形如式(22)的扩散-波方程，再利用上述色散关系，得到了一类与品质因子$Q$相关的空间分数阶导数波动方程

$ \dfrac{1}{c_0^2 }\dfrac{\partial^2 P }{\partial t^2} - \tau \dfrac{\partial }{\partial t}\left( { - \nabla ^2} \right)^{\eta + 1 / 2}P - \gamma \left( { - \nabla ^2} \right)^{\eta + 1}P = 0 $

(23)

其中，$\eta={\rm arctan} (1/Q)/\pi $，$\eta$和$\gamma $均为与频率相关的数.将式(23)与物理方程和动量守恒方程相结合即可得到相应的波动方程.该波动方程可以看作是以下两个方程的组合

$\dfrac{1}{c_0^2 }\dfrac{\partial ^2P}{\partial t^2} = \gamma \left( { - \nabla ^2} \right)^{\eta + 1}P $

(24)

$\dfrac{1}{c_0^2 }\dfrac{\partial ^2P}{\partial t^2} = \nabla ^2P + \tau \dfrac{\partial }{\partial t}\left( { - \nabla ^2} \right)^{\eta + 1 / 2}P $

(25)

式(24)主要控制声波的频散，而式(25)与空间分数阶拉普拉斯声波模型类似(式(18))，主要控制声波的耗散，等式右边第2项为其耗散项. 图 7给出了4种波动控制方程下波传播130 ms后相位和波幅示意图.与(a)区相比，(b)区的波幅明显降低；(c)区的波幅未有明显变化，但相位滞后了；(d)区的波幅降低并且相位滞后了.该图验证了式(23)在描述波频散和耗散方面的有效性.

将上述的耗散声波方程统一记作如下形式

$\dfrac{1}{c_0^2 }\dfrac{\partial^2 P}{\partial t^2} - L\left( P \right) = \nabla ^2P $

(26)

其中，$L(P)$为耗散项，则这些方程耗散项的对比可见表 1.

2 在多孔介质中传播的分数阶导数声波模型

多孔介质(如土壤、砂岩、泡沫等)是多相物质共存的一种组合体，以固相为骨架，骨架孔隙中充满液体、气体或者液气混合物.多孔介质有多种分类方法，如果按固体骨架的性质分类，可以将其分为刚性固体骨架多孔介质和骨架可变形的多孔介质.多孔介质往往具有非均匀性和各向异性的特点，声波在其中传播时，传播过程和耗散机理的复杂性促使越来越多的研究人员对其进行深入研究.特别当介质中存在断面、裂隙等界面层时，波传播的机理更加复杂^[69].

声波在各向异性的介质中传播时，其衰减依赖于空间各个位置.基于此考虑，Meerschaert等^[70]提出了一种可以描述各向异性的分数阶Stokes波动方程

$\dfrac{1}{c_0^2 }\dfrac{\partial ^2P}{\partial t^2} = \Delta ^{\overrightarrow r / 2}P + \tau \dfrac{\partial }{\partial t}\Delta ^{\overrightarrow r / 2}P $

(27)

其中，空间依赖项的定义为

$\Delta ^{\overrightarrow r / 2}P\left( {\overrightarrow x } \right) = \sum\limits_{j = 1}^d {\dfrac{\partial ^{r_j }}{\partial x_j^{r_j } }} P\left( {x_1, x_2, \cdots, x_d } \right) $

具体可参见文献[71-72].如果此方程中的拉普拉斯算子退化为一般的拉普拉斯算子, 并用于描述各向同性材料, 则该方程与热黏性波动方程(式(3))具有相同的形式.

各向同性多孔介质中声波传播可以由两个复变量来描述，即，特征阻抗和传播系数.实验发现，特征阻抗和传播系数的实部和虚部都与频率呈幂律相关^[73]，如图 8. Chen等^[74]和胡帅^[75]将特征阻抗的实部和虚部的幂律指数设定为同一值, 并对实验数据进行拟合，提出了刻画在均匀多孔介质中传播的声波模型

$\dfrac{1}{c_0^2 }\dfrac{\partial ^2P}{\partial t^2} + \dfrac{\eta }{c_0^2 }\dfrac{\partial ^{2-\eta }P}{\partial t^{2 - \eta }} = \nabla ^2P $

(28)

其中，$P$为质点速度.特别地，当$\eta=1$时，该方程退化为经典的阻尼波方程.同时，该声波方程满足Kramers-Kronig方程，其衰减系数为$\alpha \propto \omega ^{1 - \eta }$.

声波在刚性固体骨架中传播时，研究学者发现需要在欧拉方程和质量守恒方程中分别引入频率依赖的动曲率和动压缩率来描述声波在充满空气的多孔介质中的传播^[76-77]. Fellah等^[78]将动曲率和动压缩率分别进行低频和高频近似，推导出了低频和高频对应的波动方程

$\left.\!\!\begin{array}{l} \dfrac{\partial ^2P}{\partial x^2} - a\dfrac{\partial ^2P}{\partial t^2} - d\dfrac{\partial P}{\partial t} = 0 \, , \ \ \hbox{low frequency} \cr \dfrac{\partial ^2P}{\partial x^2} - A\dfrac{\partial ^2P}{\partial t^2} - C\dfrac{\partial P}{\partial t} -\cr \qquad B\int_{ - \infty }^t \dfrac{\partial ^2P / \partial t'^2}{\sqrt {t - t'} } {\rm{d}}t' = 0\, , \ \ \hbox{high frequency} \end{array}\!\!\right\} $

(29)

其中，$a$，$d$，$A$，$B$和$C$为与材料相关的常数.容易发现，低频声波的控制方程与阻尼波控制方程(式(2))类似，其衰减是与频率无关的.高频声波的控制方程，可以看作一个分数阶阶数为3/2的分数阶波动方程.随后，这组声波方程被用于研究多孔介质的散射^[79]，测定多孔介质的曲率、黏度、热特征长度、孔隙度等传播参数^[80-81].此外，该模型还被成功应用于人体骨骼超声成像^[82]和双层多孔介质^[83]等方面.

考虑声波在可变形的固体骨架中传播时，动态渗透率是频率依赖的^[76].将其代入广义的达西定律中，可以得到一组修正的Biot-JKD模型

$\left.\!\!\begin{array}{l} \sigma _{ij, j} + F_i^b = \rho \ddot {u}_i + \rho _f \ddot {w}_i \cr - p_{, i} + F_i^f = \rho _f \ddot {u}_i + m\ddot {w}_i +\cr \qquad \dfrac{k_0 }{\eta }a_{\rm s}^{ - 1 / 2} \left( { {\rm D} + a_{\rm s} } \right)^{1 / 2}\dot {w}_i \end{array}\!\!\right\} $

(30)

其中，$(D+a_{\rm s})^{1/2}$为移位的分数阶导数(shifted fractional derivative)，式中各参数可参见文献[84-85].该模型已被成功应用于一维、二维各向同性和横观各向同性多孔介质中^[86-88].

3 衰减指数$\eta$分布在[0，2]之间的统计解释

大量的实验证明，声波的能量耗散具有幂律特征，由式(1)可知，其幂律衰减指数$\eta$在[0, 2]范围内变化.从图 1可以看出，幂律衰减指数$\eta$只与材料本身相关，与频率无关.事实上，电介质、振动阻尼、热黏性流体的耗散行为也都存在幂律依赖现象. Chen^[89]首次从统计学的角度阐明了声衰减指数$\eta$在[0, 2]范围内变化的原因，并将其与介质的分形结构联系在一起.

若将式(18)右侧的拉普拉斯项去掉，再对时间$t$作一次积分，即可将空间分数阶拉普拉斯波动方程近似为一个广义的扩散方程^{[11, 89]}.为不失一般性，记空间分数阶扩散方程的形式如下

$\left. \begin{array}{l} \frac{{\partial P}}{{\partial t}} + \kappa {\left( { - \frac{{{\partial ^2}}}{{\partial {x^2}}}} \right)^{\eta /2}}P = 0\\ P\left( {x, 0} \right) = \delta \left( x \right), - \infty < x < \infty \end{array} \right\} $

(31)

其中，$\kappa $为扩散系数，$\delta (x)$为狄拉克函数.该方程的解为

$P\left( {x, t} \right) = \dfrac{1}{t^{1 / \eta }}w_\eta \left( {\dfrac{x}{t^{1 / \eta }}} \right) $

(32)

式中

${w_\eta }\left( \xi \right) = \frac{1}{{2\pi }}\int_{ - \infty }^\infty {{{\rm{e}}^{ - {\rm{i}}q\xi }}} {W_\eta }\left( k \right){\rm{d}}k $

其中，$\xi = x / t^\eta$，$W_\eta \left( k \right) = \exp \left( { - \kappa k^\eta } \right)$是稳态Lévy分布概率密度函数进行傅里叶变换后的形式.特别地，当$\eta=1$时，稳态Lévy分布为柯西分布；当$\eta=2$时，该分布退化为正态分布.所以，这两种分布是稳态Lévy分布的特例.特殊情况$\eta=2$时，式(31)描述的是正常扩散，该式的解对应的统计分布是正态分布.相应介质多为晶体等各向同性均质材料(例如，金属等).从随机行走的角度来看，该方程描述了粒子的布朗运动，即粒子的等待时间满足泊松过程，跳跃步长则是正态分布的.一般情况下，式(31)描述的是反常扩散过程，该式的解服从指标为$\eta$的稳态Lévy分布.相应介质为多相非均质材料(例如，生物体、岩土、颗粒材料等).从随机行走的角度来看，该方程描述的粒子等待时间仍然满足泊松过程，而跳跃步长则服从Lévy分布.

Saichev等^[90]指出，为了保证概率密度函数为正，Lévy分布的指数必须在(0, 2]之间.也就是说，在$\eta$大于2时，材料在统计上来说是不稳定的；当$\eta$为0时，此时的耗散是与频率无关的.如此，即从统计上解释了$\eta$的变化范围.

将式(1)中的衰减函数写成如下形式

$\eta = \dfrac{\ln \left( {\alpha / \alpha _0 } \right)}{\ln \omega } $

(33)

式(33)揭示了频率幂律依赖衰减的自相似性. Mandelbrot^[91]和Sato^[92]就Lévy概率密度函数的自相似性，指出了Lévy稳态分布与分形之间的内在联系.基于此，$\eta$也可以被认为是与分形维数相关的(前文提到衰减指数$\eta$为Lévy分布的指数).所以$\eta$只与材料的微观和细观结构有关.例如，不同人体组织的$\eta$表征了组成该生物材料的大分子的统计几何性质.

目前关于分数阶模型与材料微结构分形特征之间关系的研究成果多是定性的，定量结论很少.因而，分数阶模型一般多是唯象模型，缺少坚实的物理基础.主要的验证手段是实验数据的拟合.统计热力学可能是建立两者之间联系的一个有效方法.分数阶导数模型的统计力学背景比较清楚，研究也较深入；但目前材料分形结构与其热力学行为关系的研究也还不成熟.

另一方面，函数分形维与分数阶导数的关系的研究较多.函数的分形维数和分数阶微积分的阶数之间呈线性关系. Mandelbrot等^[93]的研究表明分数阶布朗运动函数的分形维数与Hurst指数之间的关系为$D_{\rm f} = 2 - H$，其中$D_{\rm f}$为分形维数，$H$为分数阶布朗运动的Hurst指数，$H - \dfrac{1}{2}$为分数阶阶数. Berry^[94]在其研究中指出分形维数和衰减指数如下关系，$D_{\rm f} = \dfrac{5 - \alpha }{2}$，而分数阶阶数$\eta$与衰减指数$\alpha $的关系为$\eta = \alpha / 2$^[95-96].由此可发现，分形维数可由分数阶阶数线性表示. Su及其合作者^[97-99]也证明了，对于Weierstrass函数而言，其分形维数与分数阶阶数之间的关系也是线性的.所以，分数阶阶数与分形维数之间是呈线性相关的，而他们之间普适性的定量关系，目前尚未有报道.

此外，Henry等^[100]指出实际问题中自相似性只表现在一个有限的尺度范围内，这与图 1所观察到的实验现象也是吻合的.

4 结论和展望

复杂介质，如土壤、生物组织等，是固、液、气体的多相混合物，其物理和力学性质介于理想固体和液体之间.声波在此类复杂介质中传播时，其幂律衰减的指数介于0和2之间.经典的声波控制方程只能描述与频率无关或频率平方依赖的声黏性阻尼衰减现象，而分数阶导数模型能够刻画声波的频率任意次幂律依赖衰减.分数阶阶数的物理意义也很明确，即为双对数坐标下衰减系数与频率所呈直线的斜率.

目前的主要研究思路是针对声衰减的幂律依赖现象展开研究，分析所得的波动方程是否满足因果律(Kramers-Kronig关系)，最后采用数值算法进行仿真分析.本文概述了现有的各种分数阶导数波动模型，分析和比较了这些模型之间的异同点.可以发现，这些模型主要分为唯象模型和本构模型两大类.最后，本文介绍了幂律衰减指数介于0和2之间的统计解释，并将频率依赖黏性衰减与介质的分形结构联系起来.

尽管分数阶导数波动方程在刻画声波的黏性耗散研究方面有了较深入的进展，然而引起声波衰减的因素还有其他情况.例如，反射、散射等.本文讨论的声波模型主要考虑的是声黏性吸收衰减.考虑到声波在传播过程中对空间位置的依赖性和历史的记忆性，如何利用分数阶导数建立散射、反射等引起的衰减声波控制方程是未来的一个研究课题.

目前，研究声波在多孔介质中传播的问题时，主要依赖于Johnson-Koplik-Dashen (JKD)模型.该模型中包含较多的参数，考虑到分数阶导数建模的优势在于控制方程的参数少、形式简单，如何对传统的JKD模型进行优化或者建立等效模型也是一个值得研究的课题.

分数阶导数算子全局相关的性质，会导致在数值求解分数阶导数波动方程时，计算成本和存储量很高.分形导数作为一个局部算子，目前已被成功用于模拟扩散和阻尼振动过程.发展相应的分形导数波动方程可以大大降低仿真成本.