王振龍 王學林* 胡于進 師 洪 程華茂

1(華中科技大學機械科學與工程學院，武漢 430074)

2(華中科技大學同濟醫(yī)學院附屬協(xié)和醫(yī)院耳鼻咽喉科，武漢 430022)

引言

人耳的結構十分復雜而且精細，從空氣傳導過程來看，外界的聲音信號傳遞到中樞聽覺系統(tǒng)需要經(jīng)過外耳道、中耳、內耳等器官。為了解聲音在人耳內傳遞的機理以及人耳的結構與其功能的關系，已經(jīng)提出了很多中耳和耳蝸傳聲分析模型。有限元數(shù)值方法由于能較準確地模擬復雜的幾何形狀，較容易用統(tǒng)一的算法處理聲音與結構的耦合，因而在耳聲傳遞分析中得到了很大的發(fā)展［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個數(shù)量級。此外，在高頻范圍上，計算的鼓膜和蹬骨底板的位移過小。在建立該模型時還發(fā)現(xiàn)，如采用真實的耳蝸前庭階與鼓階容積，該模型的蹬骨底板位移計算結果在高頻時會變得更小。這表明在建立耳聲傳遞計算模型的時候，中耳和耳蝸模型的匹配還有待深入討論，建立精確的中耳有限元模型就要考慮耳蝸對中耳作用的頻率相關性。同時，合適地選擇耳蝸的流體計算模型，也是改進計算精度的重要研究內容。

本研究利用高解析CT掃描影像數(shù)據(jù)，構建外耳道和中耳的幾何模型和有限元模型，同時建立了真實大小的簡化耳蝸模型;通過外耳道、中耳和被動耳蝸的集成，完成了從外耳道到耳蝸的傳聲模擬，并與相關文獻中中耳轉換功能的實驗數(shù)據(jù)進行了比較。

1 有限元模型

建立的有限元模型包括外耳道、中耳聽骨鏈，中耳韌帶肌腱和簡化耳蝸。

1.1 中耳模型

利用華中科技大學同濟醫(yī)學院附屬協(xié)和醫(yī)院采集的無聽力損傷病史的顳骨螺旋CT掃描影像數(shù)據(jù)集(Siemens Sensation 16)，通過圖像處理軟件和CAD軟件，獲得中耳幾何模型，并通過有限元前處理軟件Hypermesh進行網(wǎng)格劃分，獲得中耳的三維有限元模型。主要建模過程參照文獻［1－2］，只是獲得影像數(shù)據(jù)的途徑不同。

中耳有限元模型的材料參數(shù)，包括鼓膜、砧骨、鐙骨，錘骨、砧骨-錘骨關節(jié)、錘骨柄的密度和楊氏模量，參考文獻［2］和［13］的數(shù)據(jù)。另外，經(jīng)過計算分析過程的反復驗證，最終將砧骨-鐙骨關節(jié)楊氏模量取為5.5 MPa。中耳各個部分的泊松比都取0.3，瑞利(Raleigh)阻尼系數(shù) α 和 β分別假設為 0 s－1和0.75 ×10－4s［12］。中耳結構的韌帶、肌腱以及鼓膜和鐙骨環(huán)形韌帶都假定為彈性材料，材料參數(shù)列于表1。圖1為建立的集成外耳道、中耳和耳蝸的有限元模型。

表1 有限元模型材料參數(shù)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，其中鼓階進行了半透明處理)，將耳蝸簡化處理為非螺旋狀的、充滿液體的雙腔導管結構，主要包括前庭階、鼓階、基底膜和蝸孔?；啄らL度取32 mm，基底膜在寬度和厚度上的尺寸都是線性變化的，頂端和底端的寬度分別為0.1 mm和0.5 mm，頂端和底端的厚度分別為 7.5和2.5 μm［12］;前庭階和鼓階的截面積自蝸底向蝸頂線性減小，容積分別為48 mm3和47 mm3。根據(jù)文獻［14］的測量數(shù)據(jù)，人耳蝸前庭階與鼓階的總容積約為 75.9 μL［14］。根據(jù)文獻［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，阻尼系數(shù) β均為0.75×10－4。卵圓窗(或者鐙骨環(huán))的楊氏模量為4 MPa，圓窗膜的楊氏模量為0.2 MPa?；啄さ膹椥阅Ａ垦刂啄らL度方向變化［15］，在本研究中，基底膜的彈性模量由底部的25 MPa線性地減小到10 mm處的15 MPa，然后按指數(shù)函數(shù)規(guī)律減小到頂端的1.8 MPa?；啄さ淖枘嵯禂?shù)β是按照雙線性變化的，首先由底端的 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，鼓膜、聽小骨、肌腱和韌帶、關節(jié)軟組織以及耳蝸中的基底膜采用實體單元，流體結構包括外耳道和耳腔內的空氣、耳蝸內的外淋巴液都使用 Fluid30單元模擬。通過定義的流固耦合面，可以將聲壓作為力邊界條件耦合到結構分析中，其中包括:聲音單元與鼓膜、聽小骨以及韌帶表面的耦合，鼓膜兩側分別與中耳腔和外耳道內的流體耦合，卵圓窗和圓窗的兩側分別與耳道和耳蝸內流體的耦合，基底膜與兩側的前庭階和鼓階內的淋巴液的耦合。

在人耳聽力閾0.1～10 kHz范圍內，在外耳道距離鼓膜2 mm處施加90 dB SPL的壓力時，分析中耳在聲音傳播過程中的增益以及耳蝸內基底膜的波動傳播情況。通過分析計算，可以獲得鼓膜、鐙骨足板和基底膜沿長度方向的位移，前庭階以及鼓階內的壓力值，將有限元分析結果與相同條件下的實驗數(shù)據(jù)比較，可以驗證模型的合理性，為模型的進一步利用奠定基礎。

3 計算結果與驗證

3.1 鼓膜振動位移及聲壓增益

圖3(a)所示為模型計算鼓膜臍部(TM)和鐙骨足板(FP)的位移響應，分別對應空心三角形標記的曲線和空心矩形標記的曲線，圖3(b)表示的是對應的相位曲線。

圖(3)中虛線表示的是Gan等對人耳的顳骨進行實驗測量的平均位移曲線［2］;模型的分析結果與Gan的實驗數(shù)據(jù)相比具有較好的一致性，鼓膜和鐙骨足板的有限元分析位移曲線在小于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的聲壓引起的基底膜由蝸底至蝸頂?shù)奈灰魄€如圖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)為用鐙骨足板位移規(guī)范化處理的基底膜沿長度方向振動位移(dB)。頻率在0.5～10 kHz變化時，對應的基底膜位移的峰值位置依次由靠近頂端的位置變化到靠近底端的位置。如圖5(a)所示，7個頻率處的峰值在1.38～129 nm之間變化;相應地，圖5(b)中被動基底膜的最大位移相對于足板位移放大約為25 dB。

為了說明基底膜沿長度方向的頻率選擇特性的模擬結果，將沿基底膜的位移(dBM)相對于鐙骨足板(dFP)的位移進行規(guī)范化，用dBM/dFP表示，并繪于圖6。圖6(a)為不同頻率下基底膜位移相對于鐙骨足板位移的比值，圖6(b)表示與圖6(a)對應的相位角的曲線。從圖6(a)中可明顯觀察到每條曲線都有各自的峰值，而且是有限的并具有單一峰值，基底膜長度方向振動峰值的位置隨著頻率的減小由蝸底部逐漸移動到頂部。隨著頻率的增加，曲線的變化越來越陡，而且峰值區(qū)域越窄。

圖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或者在前庭階靠近頂部的一半時，壓力變化逐漸減小并趨于穩(wěn)定，而且隨著頻率的增加達到穩(wěn)定越快。圖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

耳蝸的輸入阻抗是中耳計算中不易準確決定的參數(shù)，而前庭階的截面積是耳蝸輸入阻抗最重要的影響因素［23］。本研究采用人耳蝸尺寸的模型，得到的中耳壓力增益與相關文獻的實驗結果基本一致。在圖4中，計算結果與Aibrara和Nakajima的實驗數(shù)據(jù)曲線相比較，在0.4～10 kHz的頻率范圍內很好符合，因此可以較準確地體現(xiàn)中耳的聲壓增益。我們注意到，Nakajima等的實驗采用了最新的高精度光纖壓力傳感器［19］，測量中對標本的破壞較少，因而有較高的可靠性。在圖4中，也畫出了采用文獻［12］的模型，取正常液體體模量的計算結果，其低頻范圍的計算結果偏低。本模型在低頻計算上有了較大的改進。對中耳壓力增益模擬，只有在中耳和耳蝸集成模型下才能計算得到。目前普遍采用的中耳聲音傳遞計算有限元模型將耳蝸簡化為集中質量和彈簧阻尼單元，不可能分析以中耳壓力增益為指標的中耳轉換功能，因此本研究建立的模型使對中耳轉換功能的模擬更加深入。同時，部分頻率的中耳壓力增益和相位的計算結果與實驗結果有一定的差別，在其他類型的模型上也有相似的問題，按照文獻［23］對耳蝸輸入阻抗的模型分析，相位計算結果的誤差可能是由于采用單調減少的前庭階截面模型所造成的。

從圖6(a)中，可明顯觀察到每頻率對應的基底膜位移曲線都有各自的峰值，基底膜長度方向振動峰值的位置隨著頻率的減小由蝸底部逐漸移動到頂部。隨著頻率的增加，曲線的變化越來越陡，而且峰值區(qū)域越窄。這種變化與von Békésy的實驗測量［20］是一致的，與其他模型分析推斷［21－22］以及對動物的耳蝸實驗觀察［24］相類似。

耳蝸基底膜位置－頻率關系圖定義了沿耳蝸長度方向基底膜速度峰值的位置，它是輸入頻率的函數(shù)，通過對比分析，可以看出模型結果與實驗結果及其他經(jīng)驗模型結果一致?；啄さ膭偠葘啄の恢?頻率關系和基底膜位移曲線的形態(tài)有重要的影響。在耳蝸振動解析計算模型中，有很多基底膜剛度模型，但這些剛度模型在有限元數(shù)值計算中無法直接應用。在有限元數(shù)值計算中，必須確定基底膜彈性模量沿耳蝸長度方向的變化，在缺乏實驗數(shù)據(jù)的情況下，基底膜彈性模量的確定非常困難，只有通過基底膜位置-頻率的關系和基底膜位移曲線的形態(tài)近似估計，因此這是一個非常耗時的計算。由于采用了不同類型的流體單元和不同密度的網(wǎng)格，文獻［12］采用的基底膜彈性模量用于本研究的模型，得不到合理的基底膜位移曲線的形態(tài)，也得不出正確的基底膜位置-頻率關系圖。通過模型的計算分析，基底膜剛度、阻尼等的變化會影響響應峰值的位置以及響應幅值的大小，但剛度對基底膜位置-頻率關系的影響較為顯著。

5 結論

本研究建立了包括外耳道、中耳、中耳腔和耳蝸結構的有限元集成模型，運用聲音-結構耦合動力學分析，模擬中耳聲音傳遞功能，計算鼓膜和鐙骨足板對耳道內壓力的位移響應、中耳的聲壓增益。同時，也模擬了耳蝸基底膜運動的頻率選擇特性以及前庭階內壓力。模型計算結果與相關文獻的實驗數(shù)據(jù)具有較好的一致性。本研究用真實體積的非螺旋結構簡化耳蝸模型，將耳道內的聲音激勵與基底膜的振動和耳蝸內壓力的變化聯(lián)系起來分析，所建立的模型可為不同條件下的耳聲傳遞模擬和數(shù)值仿真研究提供基礎的計算分析手段。

本研究建立的耳蝸模型沒有考慮耳蝸放大效應等主動機制，也沒有考慮耳蝸內流體的黏度影響，這些有待做進一步的工作。

［1］Gan Rongzhu，F(xiàn)eng 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］劉迎曦，李生，孫秀珍.人耳傳聲數(shù)值模型［J］.力學學報，2008，40:107 －113.

［5］姚文娟，李武，付黎杰，等.中耳結構數(shù)值模擬及傳導振動分析［J］.系統(tǒng)仿真學報，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］劉后廣，塔娜，饒柱石.懸浮振子對中耳聲傳播特性影響的數(shù)值研究 ［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.