王振龙 王学林* 胡于进 师 洪 程华茂

1(华中科技大学机械科学与工程学院，武汉 430074)

2(华中科技大学同济医学院附属协和医院耳鼻咽喉科，武汉 430022)

引言

人耳的结构十分复杂而且精细，从空气传导过程来看，外界的声音信号传递到中枢听觉系统需要经过外耳道、中耳、内耳等器官。为了解声音在人耳内传递的机理以及人耳的结构与其功能的关系，已经提出了很多中耳和耳蜗传声分析模型。有限元数值方法由于能较准确地模拟复杂的几何形状，较容易用统一的算法处理声音与结构的耦合，因而在耳声传递分析中得到了很大的发展［1－13］。中耳声传导分析模型主要模拟声音由外耳道、鼓膜、听骨链至镫骨足板传递机制［1－5］。在这类分析的有限元模型中，耳蜗对中耳的作用常以集中质量和弹簧阻尼单元代替［6］，当前普遍采用的耳蜗模型，主要计算压力自卵圆窗到耳蜗顶端的传递和基底膜沿长度方向的运动响应［6－10］。耳蜗分析模型的输入均理想化为蹬骨足底板的活塞状运动，较少考虑中耳与耳蜗的耦合。中耳有限元模型将耳蜗对中耳的作用简化为与频率无关的弹簧和阻尼单元［1－3，5］，这与实验测量的耳蜗输入阻抗不完全相符［16－19］。在文献［4］中，其模型包含内耳部分，但没有讨论耳蜗模型。Gan等建立了中耳与耳蜗耦合的有限元模型［12］，耳蜗采用简化的直腔模型，但其耳蜗体积为220 μL，该体积与内耳中液体空间(包括前庭和耳蜗等，总容积为 204 μL ～ 228 μL)相当［14］，远大于人耳蜗中前庭阶与鼓阶的容积之和。为了使计算的中耳压力增益合理，文献中还将流体的体积模量取为220 GPa［12］，这个值较正常水的体模量高出2个数量级。此外，在高频范围上，计算的鼓膜和蹬骨底板的位移过小。在建立该模型时还发现，如采用真实的耳蜗前庭阶与鼓阶容积，该模型的蹬骨底板位移计算结果在高频时会变得更小。这表明在建立耳声传递计算模型的时候，中耳和耳蜗模型的匹配还有待深入讨论，建立精确的中耳有限元模型就要考虑耳蜗对中耳作用的频率相关性。同时，合适地选择耳蜗的流体计算模型，也是改进计算精度的重要研究内容。

本研究利用高解析CT扫描影像数据，构建外耳道和中耳的几何模型和有限元模型，同时建立了真实大小的简化耳蜗模型;通过外耳道、中耳和被动耳蜗的集成，完成了从外耳道到耳蜗的传声模拟，并与相关文献中中耳转换功能的实验数据进行了比较。

1 有限元模型

建立的有限元模型包括外耳道、中耳听骨链，中耳韧带肌腱和简化耳蜗。

1.1 中耳模型

利用华中科技大学同济医学院附属协和医院采集的无听力损伤病史的颞骨螺旋CT扫描影像数据集(Siemens Sensation 16)，通过图像处理软件和CAD软件，获得中耳几何模型，并通过有限元前处理软件Hypermesh进行网格划分，获得中耳的三维有限元模型。主要建模过程参照文献［1－2］，只是获得影像数据的途径不同。

中耳有限元模型的材料参数，包括鼓膜、砧骨、镫骨，锤骨、砧骨-锤骨关节、锤骨柄的密度和杨氏模量，参考文献［2］和［13］的数据。另外，经过计算分析过程的反复验证，最终将砧骨-镫骨关节杨氏模量取为5.5 MPa。中耳各个部分的泊松比都取0.3，瑞利(Raleigh)阻尼系数 α 和 β分别假设为 0 s－1和0.75 ×10－4s［12］。中耳结构的韧带、肌腱以及鼓膜和镫骨环形韧带都假定为弹性材料，材料参数列于表1。图1为建立的集成外耳道、中耳和耳蜗的有限元模型。

表1 有限元模型材料参数Tab.1 Material properties of ligaments and tendons in FE model

图1 中耳与耳蜗有限元模型Fig.1 The FE model of human ear with external ear canal，middle ear components and simplified cochlea

1.2 耳蜗模型

采用耳蜗力学分析中常用的非螺旋结构的耳蜗模型(见图2，其中鼓阶进行了半透明处理)，将耳蜗简化处理为非螺旋状的、充满液体的双腔导管结构，主要包括前庭阶、鼓阶、基底膜和蜗孔。基底膜长度取32 mm，基底膜在宽度和厚度上的尺寸都是线性变化的，顶端和底端的宽度分别为0.1 mm和0.5 mm，顶端和底端的厚度分别为 7.5和2.5 μm［12］;前庭阶和鼓阶的截面积自蜗底向蜗顶线性减小，容积分别为48 mm3和47 mm3。根据文献［14］的测量数据，人耳蜗前庭阶与鼓阶的总容积约为 75.9 μL［14］。根据文献［23］，人耳前庭阶的截面积在底端最大处约为0.7 cm2，在距底端约3 mm处阶梯状减小，此后前庭阶的截面积在0.01～0.018 cm2间变化。采用阶梯的前庭阶截面，其截面积与文献［23］的尺寸基本相同，耳蜗模型体积与人耳蜗体积基本一致。本研究建立的耳蜗模型没有考虑柯蒂氏器的作用，即没有考虑耳蜗的主动机制，为一被动模型。

图2 简化耳蜗模型的结构Fig.2 The structure of simplified cochlea

耳蜗模型中的结构材料属性包括卵圆窗、圆窗、基底膜的密度和弹性模量，上述结构的密度均取为1 200 kg/m3，泊松比均为0.3，阻尼系数 β均为0.75×10－4。卵圆窗(或者镫骨环)的杨氏模量为4 MPa，圆窗膜的杨氏模量为0.2 MPa。基底膜的弹性模量沿着基底膜长度方向变化［15］，在本研究中，基底膜的弹性模量由底部的25 MPa线性地减小到10 mm处的15 MPa，然后按指数函数规律减小到顶端的1.8 MPa。基底膜的阻尼系数β是按照双线性变化的，首先由底端的 1.0×10－6线性地增加到5.0×10－6，然后再线性地增长到基底膜顶端的5.0×10－5。

1.3 结构与声传递耦合模型

在耳声传递计算模型中，有外耳道内声场与鼓膜、鼓膜及听小骨与中耳腔内声场的耦合，也有耳蜗前庭阶和鼓阶内液体分别与两个圆窗、基底模间的耦合。假设流体是无黏性的，并且满足下面的控制方程，即

在本研究中，声音在空气中的传播速度以及空气的密度分别取343 m/s和1.21 kg/m3;耳蜗中的流体也采用式(1)的描述［6］，淋巴液中声音的传播速度取为1 400 m/s，淋巴液密度取1 000 kg/m3。

计算采用有限元分析软件 ANSYS 10.0，鼓膜、听小骨、肌腱和韧带、关节软组织以及耳蜗中的基底膜采用实体单元，流体结构包括外耳道和耳腔内的空气、耳蜗内的外淋巴液都使用 Fluid30单元模拟。通过定义的流固耦合面，可以将声压作为力边界条件耦合到结构分析中，其中包括:声音单元与鼓膜、听小骨以及韧带表面的耦合，鼓膜两侧分别与中耳腔和外耳道内的流体耦合，卵圆窗和圆窗的两侧分别与耳道和耳蜗内流体的耦合，基底膜与两侧的前庭阶和鼓阶内的淋巴液的耦合。

在人耳听力阈0.1～10 kHz范围内，在外耳道距离鼓膜2 mm处施加90 dB SPL的压力时，分析中耳在声音传播过程中的增益以及耳蜗内基底膜的波动传播情况。通过分析计算，可以获得鼓膜、镫骨足板和基底膜沿长度方向的位移，前庭阶以及鼓阶内的压力值，将有限元分析结果与相同条件下的实验数据比较，可以验证模型的合理性，为模型的进一步利用奠定基础。

3 计算结果与验证

3.1 鼓膜振动位移及声压增益

图3(a)所示为模型计算鼓膜脐部(TM)和镫骨足板(FP)的位移响应，分别对应空心三角形标记的曲线和空心矩形标记的曲线，图3(b)表示的是对应的相位曲线。

图(3)中虚线表示的是Gan等对人耳的颞骨进行实验测量的平均位移曲线［2］;模型的分析结果与Gan的实验数据相比具有较好的一致性，鼓膜和镫骨足板的有限元分析位移曲线在小于4 kHz的时候，与实验的平均位移曲线有很好的符合。

图4为模型计算获得的中耳压力增益(无标记粗实线)，压力增益定义为前庭阶靠近卵圆窗处的压力值与外耳道内鼓膜附近的压力的比值。图4(a)表示中耳压力增益的幅值，图4(b)表示中耳压力增益的相位。在1 000 Hz时，计算的中耳压力增益为21 dB、2 000 Hz时，中耳压力增益为23 dB。

在图4中，带圆点细实线表示Aibara等的实验结果的平均值［16］，两条黑色细实线分别表示上限和下限。图4中也列出了 Puria、Hüttenbrink与 Hudde和 Nakajima 的实验结果曲线［17—19］。

图3 模型的耳膜和镫骨足板处的位移和相位角。(a)幅值 ;(b)相位Fig.3 Displacement magnitude and phase around ear drum and stapes footplate of the model.(a)amplitude;(b)phase

图4 中耳的压力增益。(a)幅值;(b)相位角Fig.4 The middle ear pressure gain.(a)amplitude;(b)phase

3.2 基底膜的运动

由外耳道内90 dB SPL的声压引起的基底膜由蜗底至蜗顶的位移曲线如图5所示。

图5 不同频率上基底膜在长度方向的位移。(a)位移幅值;(b)位移的标准化形式Fig.5 Distribution of displacement magnitude of the BM from the base to apex at frequencies of 500–10 kHz.(a)displacement amplitude;(b)standardized displacement style

图5(a)包含了0.5～10 kHz范围内的7个不同声音频率下基底膜位移(dBM)分布曲线。图5(b)为用镫骨足板位移规范化处理的基底膜沿长度方向振动位移(dB)。频率在0.5～10 kHz变化时，对应的基底膜位移的峰值位置依次由靠近顶端的位置变化到靠近底端的位置。如图5(a)所示，7个频率处的峰值在1.38～129 nm之间变化;相应地，图5(b)中被动基底膜的最大位移相对于足板位移放大约为25 dB。

为了说明基底膜沿长度方向的频率选择特性的模拟结果，将沿基底膜的位移(dBM)相对于镫骨足板(dFP)的位移进行规范化，用dBM/dFP表示，并绘于图6。图6(a)为不同频率下基底膜位移相对于镫骨足板位移的比值，图6(b)表示与图6(a)对应的相位角的曲线。从图6(a)中可明显观察到每条曲线都有各自的峰值，而且是有限的并具有单一峰值，基底膜长度方向振动峰值的位置随着频率的减小由蜗底部逐渐移动到顶部。随着频率的增加，曲线的变化越来越陡，而且峰值区域越窄。

图6 基底膜位移与镫骨足板位移的比值。(a)幅值;(b)相位Fig.6 The BM displacement normalized with respect to the footplate displacement.(a)magnitude;(b)phase

从图6中也可以建立有限元模型计算的基底膜位置-频率关系图(place-frequency map)，相应的结果绘于图7。实线表示本模型计算所得曲线，虚线为 Greenwood等的模型计算所得曲线［22］，带标记的曲线为 Békésy［20］实验测量获得的曲线［22］。

图7 基底膜上不同频率的响应峰值位置Fig.7 The place of maximum response on the BM vs.frequency

图8表示耳蜗前庭阶内靠近基底膜处由耳道内90 dB SPL的声压引起的压力分布，压力曲线描述了0.5～10 kHz频率范围内8个频率点沿基底膜长度方向的压力大小。由耳蜗底部到顶部，前庭阶内的压力逐渐减小。

在图8(a)中，当横坐标轴大于1.6 cm或者在前庭阶靠近顶部的一半时，压力变化逐渐减小并趋于稳定，而且随着频率的增加达到稳定越快。图8(b)所示为前庭阶压力的相位角变化，曲线显示，随着频率的增加，压力的相位角变化也随之变大。当频率在200～600 Hz范围内时，压力的相位角都是正值;当频率大于等于800 Hz时，压力的相位角都是由低频的正角度变化到高频的负角度。

4 讨论

本研究建立的中耳与耳蜗集成的有限元模型，一个重要的考虑就是希望合理地反映了耳蜗与中耳的相互作用。从图3可看出，建立的模型得到了与实验结果较为一致的中耳响应曲线。与文献［12］的结果相比，镫骨足板位移的高频计算结果有了很大的改进。

图8 不同频率上前庭阶在基底膜长度方向的压力。(a)幅值;(b)相位Fig.8 Distribution of the sound pressure in scala vestibuli along the BM length at frequencies of 400 Hz～10 kHz.(a)amplitude;(b)phase

耳蜗的输入阻抗是中耳计算中不易准确决定的参数，而前庭阶的截面积是耳蜗输入阻抗最重要的影响因素［23］。本研究采用人耳蜗尺寸的模型，得到的中耳压力增益与相关文献的实验结果基本一致。在图4中，计算结果与Aibrara和Nakajima的实验数据曲线相比较，在0.4～10 kHz的频率范围内很好符合，因此可以较准确地体现中耳的声压增益。我们注意到，Nakajima等的实验采用了最新的高精度光纤压力传感器［19］，测量中对标本的破坏较少，因而有较高的可靠性。在图4中，也画出了采用文献［12］的模型，取正常液体体模量的计算结果，其低频范围的计算结果偏低。本模型在低频计算上有了较大的改进。对中耳压力增益模拟，只有在中耳和耳蜗集成模型下才能计算得到。目前普遍采用的中耳声音传递计算有限元模型将耳蜗简化为集中质量和弹簧阻尼单元，不可能分析以中耳压力增益为指标的中耳转换功能，因此本研究建立的模型使对中耳转换功能的模拟更加深入。同时，部分频率的中耳压力增益和相位的计算结果与实验结果有一定的差别，在其他类型的模型上也有相似的问题，按照文献［23］对耳蜗输入阻抗的模型分析，相位计算结果的误差可能是由于采用单调减少的前庭阶截面模型所造成的。

从图6(a)中，可明显观察到每频率对应的基底膜位移曲线都有各自的峰值，基底膜长度方向振动峰值的位置随着频率的减小由蜗底部逐渐移动到顶部。随着频率的增加，曲线的变化越来越陡，而且峰值区域越窄。这种变化与von Békésy的实验测量［20］是一致的，与其他模型分析推断［21－22］以及对动物的耳蜗实验观察［24］相类似。

耳蜗基底膜位置－频率关系图定义了沿耳蜗长度方向基底膜速度峰值的位置，它是输入频率的函数，通过对比分析，可以看出模型结果与实验结果及其他经验模型结果一致。基底膜的刚度对基底膜位置-频率关系和基底膜位移曲线的形态有重要的影响。在耳蜗振动解析计算模型中，有很多基底膜刚度模型，但这些刚度模型在有限元数值计算中无法直接应用。在有限元数值计算中，必须确定基底膜弹性模量沿耳蜗长度方向的变化，在缺乏实验数据的情况下，基底膜弹性模量的确定非常困难，只有通过基底膜位置-频率的关系和基底膜位移曲线的形态近似估计，因此这是一个非常耗时的计算。由于采用了不同类型的流体单元和不同密度的网格，文献［12］采用的基底膜弹性模量用于本研究的模型，得不到合理的基底膜位移曲线的形态，也得不出正确的基底膜位置-频率关系图。通过模型的计算分析，基底膜刚度、阻尼等的变化会影响响应峰值的位置以及响应幅值的大小，但刚度对基底膜位置-频率关系的影响较为显著。

5 结论

本研究建立了包括外耳道、中耳、中耳腔和耳蜗结构的有限元集成模型，运用声音-结构耦合动力学分析，模拟中耳声音传递功能，计算鼓膜和镫骨足板对耳道内压力的位移响应、中耳的声压增益。同时，也模拟了耳蜗基底膜运动的频率选择特性以及前庭阶内压力。模型计算结果与相关文献的实验数据具有较好的一致性。本研究用真实体积的非螺旋结构简化耳蜗模型，将耳道内的声音激励与基底膜的振动和耳蜗内压力的变化联系起来分析，所建立的模型可为不同条件下的耳声传递模拟和数值仿真研究提供基础的计算分析手段。

本研究建立的耳蜗模型没有考虑耳蜗放大效应等主动机制，也没有考虑耳蜗内流体的黏度影响，这些有待做进一步的工作。

［1］Gan Rongzhu，Feng Bin，Sun QunLi.Three-dimensional finite element modeling of human ear for sound transmission［J］.Ann Biomed Eng，2004，32:847－859.

［2］Gan Rongzhu，Wang Xuelin.Multifield coupled finite element analysis for sound transmission in otitis media with effusion［J］.J Acoust Soc Am，2007，122:3527－3538.

［3］Koike T，Wada H.Modeling of the human middle ear using the finite-element method［J］.J Acoust Soc Am，2002，111:1306－1317.

［4］刘迎曦，李生，孙秀珍.人耳传声数值模型［J］.力学学报，2008，40:107 －113.

［5］姚文娟，李武，付黎杰，等.中耳结构数值模拟及传导振动分析［J］.系统仿真学报，2009，21:651－654.

［6］Bohnke F，Arnold W.3D-finite element model of the human cochlea including fluid－structure couplings［J］.ORL，1999，61:305－310.

［7］Lim KM，Steele CR.A three-dimensional nonlinear active cochlear model analyzed by the WKB-numeric method ［J］.Hear Res，2002，170:190 －205.

［8］Andoh M，Wada H.Prediction of the characteristics of two types of pressure waves in the cochlea:Theoretical considerations［J］.J Acoust Soc Am，2004，116:417－425.

［9］Watts L.The mode-coupling Liouville-Green approximation for a two-dimensional cochlear model［J］.J Acoust Soc Am，2000，108:2266－2271.

［10］Steel CR，Lim K-M.Cochlear model with three dimensional fluid，inner sulcus and food-forward mechanism ［J］.Audiol Neuro-Otol，1999，4:197 －203.

［11］刘后广，塔娜，饶柱石.悬浮振子对中耳声传播特性影响的数值研究 ［J］.力学学报，2010，42:109－114.

［12］Gan Rongzhu，Brian R，Wang Xuelin.Modeling of sound transmission from earcanalto cochlea ［J］. Annals of Biomedical Engineering，2007，35:2180－2195.

［13］WangXuelin， Cheng Tao， Gan Rongzhu. Finite-element analysis of middle-ear pressure effects on static and dynamic behavior of human ear［J］.J Acoust Soc Am，2007，122:906－917.

［14］Igarashi M，Ohashi K，Ishii M.Morphometric comparison of endolymphatic and perilymphatic spaces in human temporal bones［J］.Acta Otolaryngol(stockh)，1986，101:161 －164.

［15］Emadi G，Richter CP，Dallos P.Stiffness of the gerbil basilar membrane:Radial and longitudinal variations ［J］. J Neurophysiol，2003，91:474 －488.

［16］Aibara RJ，Welsh JT，Puria S，et al.Human middle-ear sound transfer function and cochlear impedance［J］.Hear Res，2001，152:100－109.

［17］Puria S，Peake WT，Rosowski JJ.Sound-pressure measurements in the cochlea vestibule of human cadaver ears［J］.J Acoust Soc Am，1997，101:2754－2770.

［18］Hüttenbrink KB，Hudde H.Untersuchungen zur schalleitung durch das rekonstruierte mittelohr mit einem hydrophone［J］.HNO，1994，42:49－57.

［19］NakjimaHH， Dong W， Olson ES， etal. Differential intracochlear sound pressure measurements in normal human temporal bones［J］.JARO，2008，10:23－26.

［20］von Békésy G.Experiments in Hearing ［M］.New York:McGraw-Hill Book Company，Inc，1960.

［21］SteeleCR. Toward three-dimensionalanalysis ofcochlear structure［J］.ORL，1999，61:238 －251.

［22］Greenwood DD. A cochlearfrequency-position function for several species—29 years later［J］.J Acoust Soc Am，1990，87:2592－2605.

［23］Puria S， Allen JB. A parametric study of cochlear input impedance［J］.Acoust Soc Am，1991，89:287－309.

［24］Ren Tianying. Longitudinalpattern ofbasilar membrane vibration in the sensitive cochlea［J］.Proc Natl Acad Sci USA，2002，99:17101－17106.