（电路与系统专业论文）骨传导超声助听的特性分析与系统设计.pdf

摘要 发现骨传导人可以感知超声以来，人们通过人体实验的方法，研究和分析了 骨传导超声的特性。目前，骨传导超声已应用于医学领域。骨传导超声助听器就 是其中一个非常重要的应用，它可以帮助一些利用传统助听器无法恢复听力的感 音神经性耳聋患者和某些严重耳聋患者在一定程度上恢复听力。但是，骨传导超 声助听器的性能并不理想。同时，由于超声感知机理尚不明确，要提出直接有效 地提高助听器性能的方法还比较困难。 本文主要的研究对象就是骨传导超声助听器，通过研究骨传导超声听力的特 性，进一步提高骨传导超声助听器的性能。本文主要从计算机仿真、硬件系统设 计和心理声学实验三个方面展开工作。在计算机仿真中，模拟超声在人脑中传播 的情形，利用时域有限差分方法计算声场分布。通过分析声场，与传统的骨传导 声音形成的声场作比较，得到了骨传导超声听力的新特性。然后根据骨传导超声 助听器的原理，设计了基于数字信号处理器的骨传导超声助听器系统。利用该助 听器设备，进行了几组心理声学实验。通过实验，更直观地了解并研究了骨传导 听力的特性，总结出一些有利于提高骨传导超声助听器性能的方法。 仿真实验结果表明：双耳骨传导超声听觉没有传统双耳听觉的定位能力，而 单耳骨传导超声可以同时激励双耳耳蜗，传输效率比骨传导声音高。因此，单耳 骨传导超声助听的效果可能比双耳骨传导超声助听的效果更为理想。同时，利用 骨传导超声助听器进行的心理声学实验的结果表明：大多数实验参与者均能感知 到超声产生的高频音调，部分实验者通过骨传导超声助听器可以进行音频识别。 通过实验，本文还预测了骨传导超声的解调方式，进一步估计了骨传导超声助听 中载波频率的最优选择。这些结果将有助于骨传导超声助听性能的进一步提高。 关键词：骨传导：超声：助听器；时域有限差分方法；数字信号处理器；心理声 学实验 a b s t r a c t p e o p l eh a sb e g a nt os t u d yc h a r a c t e r i s t i c so ft h eb o n e c o n d u c t e du l t r a s o u n d ( b c u ) w i t l lh u m a nh e a r i n gt e s t s s i n c ei tw a sf o u n dt h a tp e o p l ec a l lp e r c e i v et h e u l t r a s o u n dt h r o u g ht h eb o n ec o n d u c t i o n t h eb c uh a ss e v e r a lm e d i c a la p p l i c a t i o n s ， a m o n gw h i c ht h eb o n ec o n d u c t e du l t r a s o n i ch e a r i n ga i d ( b c u h a ) i st h em o s t i m p o r t a n t o n e 1 1 1 eb c u h ac a n h e l p t h o s e h e a r i n gi m p a i r e dp e r s o n s w i t h s e n s o r i n e u r a lh e a r i n gl o s sa n ds o m ep r o f o u n dh e a t i n gl o s sw h oc a nn o tr e s t o r et h e i r h e a r i n g 、 ，i t t lt h et r a d i t i o n a lh e a r i n ga i d h o w e v e r , t h ep e r f o r m a n c eo ft h i sn e w h e a t i n g a i di sn o tq u i t ef a v o r a b l e m e a n w h i l e ，a st h e u n d e r l y i n gp e r c e p t i o n m e c h a n i s mo ft h eb c ui si n k n o w r ly e t ，i ti sd i f f i c u l tt oi m p r o v et h eh e a r i n ga i d d i r e c t l ya n de f f e c t i v e l yf r o mi t sp r i n c i p l e t h i sd i s s e r t a t i o ni sa i m e da ti m p r o v i n gt h ep e r f o r m a n c eo ft h eb c u h ab y s t u d y i n gc h a r a c t e r i s t i c so ft h eb c u t h es t u d yc a nb ed i v i d e di n t ot h r e ea s p e c t s ：t h e c o m p u t e rm o d e l i n g ，t h eh a r d w a r ed e s i g na n dt h eh u m a np s y c h o a c o u s t i c a lt e s t i nt h e c o m p u t e rm o d e l i n g ，t h ef i n i t ed i f f e r e n c et i m ed o m a i n ( f d t d ) m e t h o di su s e dt o c a l c u l a t et h ea c o u s t i cf i e l di nt h eh u m a nh e a dm o d e lt os i m u l a t et h et r a n s m i s s i o no f t h e u l t r a s o u n d b yc o m p a r i n g w i t ht h eb o n ec o n d u c t e da u d i b l e s o u n d ，n e w c h a r a c t e r i s t i c so ft h eb c ua r eo b t a i n e d t h e n ，g u i d e db yt h ep r i n c i p l eo ft h eb c u h a ， ab c u h ad e v i c eb a s e do nt h ed i g i t a ls i g n a lp r o c e s s o r ( d s p ) s y s t e mi sd e v e l o p e d w i t l lt h i sd e v i c e ，af e ws e t so fp s y c h o a c o u s t i c a lt e s t sa r ec o n d u c t e dt os t u d ym o r e a b o u tc h a r a c t e r i s t i c so ft h eb c u s e v e r a lm e t h o d sh a v e b e e np r o p o s e dt oi m p r o v et h e p e r f o r m a n c eo ft h eb c u h a s i m u l a t i o nr e s u l t ss h o wt h a tt h eb i l a t e r a lb c uh e a r i n gd o e s n th a v et h es p a c e p e r c e p t i o na b i l i t ya st h ec o m m o nb i n a u r a lh e a t i n g ，w h i l et h eu n i l a t e r a lb c us t i m u l u s c a l ls t i m u l a t eb o t hc o c h l e a e t h es t u d yi n d i c a t e st h a tt h eu n i l a t e r a lb c u h am a y p e r f o r mb e t t e rt h a nt h eb i l a t e r a lo n e m e a n w h i l e ，r e s u l t si np s y c h o a c o u s t i c a lt e s t s s h o wt h a tm o s ts u b j e c t sc a np e r c e i v et h eh i g h - f r e q u e n c yp i t c ho ft h eu l t r a s o u n da n d p a r t o fs u b j e c t sc a nd i s c r i m i n a t ed i f f e r e n ta u d i b l e s i g n a l sm o d u l a t e do n t ot h e u l t r a s o n i cc a r r i e r i na d d i t i o n ，t h i sd i s s e r t a t i o na l s oe s t i m a t e st h ed e m o d u l a t i o nt y p eo f t h eb c ua n dt h eb e s tc h o i c ef o rt h eu l t r a s o n i cc a r r i e rf r e q u e n c y i ti sh o p e dt h a tt h e s e r e s u l t sc a l lh e l pt of u l t h e ri m p r o v et h ep e r f o r m a n c eo ft h eb c u h a k e y w o r d s ：b o n ec o n d u c t i o n ；u l t r a s o u n d ；h e a r i n ga i d ；f d t d ；d s p ； p s y c h o a c o u s t i c a lt e s t 1 1 研究背景 第一章绪论 1 1 1 骨传导超声听觉简介 上世纪初，人们一直认为人耳可听到的声音最高频率为2 4k h z 。然而，1 9 5 0 年r p u m p h r e y 首次报导人通过骨传导可以感知到频率高达1 2 0k h z 的超声【i 】。 从此，关于骨传导超声的研究拉开了序幕。 骨传导超声在医学中的应用主要在耳疾的诊断和治疗方面，如耳鸣的治疗【2 j 、 突发性耳聋的诊断、美尼尔症的诊断【3 】和噪声性听障的诊断【4 1 。其中，一个非常 重要的应用是骨传导超声应用于助听技术。骨传导超声助听技术起源于1 9 9 1 年 m l l e n h a r d t 等的研究垆j 。他们首次报导，人不但可以在一定频率范围内区分超 声的频率差异，而且还可以识别调制在超声上的语音信号。同时，他们对老年听 力受损和某些严重耳聋患者做同样的实验，结果表明：很多听力受损的患者，特 别是通过传统的助听器仍无法听到声音的感音神经性耳聋患者，通过骨传导超声 技术可以感知到声音，并具有一定的语音识别能力。在此基础上，2 0 0 4 年，日 本产业综合研究所( a d v a n c e di n d u s t r i a ls c i e n c ea n dt e c h n o l o g y , a i s t ) 根据骨传导 超声助听的原理及特性研发出骨传导超声助听器( b o n e c o n d u c t e du l t r a s o u n d h e a r i n g a i d ，b c u h a ) t 。实验结果表明：利用该助听器，4 0 的严重耳聋患者可 以感知到声音，1 7 的严重耳聋患者可以识别语音。骨传导超声助听器的使用比 人工耳蜗移植方便很多，可以使患者免受精神上和肉体上的痛苦：而且由于超声 感知机理与传统助听器助听原理不同，对于一些利用传统助听器无法改善听力的 患者，骨传导超声助听器可达到较好的效果。因此，骨传导超声助听器具有较大 的应用价值。 目前，尽管骨传导超声助听技术的研究已取得了一定的成果，人们对骨传导 超声的特性有了一定的了解，但对于其超声感知机理仍不明确，有待进一步研究。 1 1 2 国内外研究现状 自从上世纪中叶发现人可以感知超声这一现象后，人们就一直在研究超声感 知的机理。人们认为，由于声特性阻抗不匹配，超声肯定不是通过中耳传导，而 只能通过骨传导【l 】。但单单这一理论还不足以解释其感知机理。对超声感知特性 的研究有利于人们进一步揭示超声感知的机理。因此，目前有很多人从不同的角 度研究骨传导超声的特性。 超声感测区域 1 9 5 0 年，r p u m p h r e y 就提出如果只是通过气传导，人只能听到频率在1 6 5 k h z 以下的声音，但是如果将传感器固定在耳乳突骨或是颞叶区域，人可以感知 到频率高达1 0 0k h z 的超声。实验表明：耳乳突骨和颞骨是超声感测的较理想区 域【l 】o 超声语音感知 1 9 9 1 年，m l l e n h a r d t 等通过在听力正常人体上的实验发现：人不但可以在 一定频率范围内识别超声的频率差异，而且还可以识别调制在超声上的语音信号 【5 】。他们进一步对老年听力受损和严重耳聋患者做了同样的实验，结果表明：这 些人对超声频率的识别也有较高的准确率。于是他们估计骨传导超声可能对耳聋 患者复聪有所帮助。这一研究成果开辟了将骨传导超声技术运用于人体助听的未 来。 1 9 9 8 年，h h o s o i 等用频率为1k h z 和2k h z 的气传导纯音以及将1k h z 和 2k h z 的纯音通过调幅方式调制到超声载波上作为激励，将它们听觉诱发的n l m 的等效电流偶极子( e c d ) 叠加到由磁共振成像( m 酣) 所获得的解剖结构图像上， 以准确地标记出e c d 的位置。结果表明：气传导纯音听觉诱发的n l m 的e c d 和骨传导纯音调制超声听觉诱发的n l m 的e c d 都位于听觉皮层，而且两者位置 十分靠近。结果表明：骨传导超声激励刺激的是人脑的听觉皮层1 7 1 。 感知频率范围和听阈值 1 9 6 3 年，j f c o r s o 通过人体实验测定了频率在5 1 0 0l d - i z 范围的骨传导声波 的听阈值。结果如图1 1 所示。 暑 e 鲁 蓄宗 o 卜 幂 0i o i l 2 0 3 0 4 8 5 56 7 7 7 ，5 频s f k i - i = ) 图i 1 高音频及超声频率范围内的骨传导声波的听力阈值曲线 2 从图1 1 可以看出：1 0 2 0k h z 的频率范围内声压有一个明显的增长，斜率近 似为5 0d b o c t a v e ( 分贝倍频程) ，而2 0 1 0 0k h z 斜率较为平缓，约1 5d b o c t a v e ( 分 贝倍频程) 。j f c o r s o 分析：频率低于2 0k h z ，声波存在两种传导途径：气传导 和骨传导，并以气传导为主。频率超过2 0k h z ，声波只能由骨传导。以2 0k h z 为分界的曲线斜率的显著差异表明：气传导和骨传导可能涉及的是两种不同的听 觉机制，若是同一种听觉机制则可能运作方式不一样( 8 】。 超声听力动态范围 2 0 0 3 年，t n i s h i m u r a l 等测量了超声听力的动态范刮9 1 。心理声学上动态范围 指的是引起人不舒适响度值( u n c o m f o r t a b l el e v e l ，u c l ) 与可承受的响度阈值之 间的差异。正常的听力动态范围比较宽( 约大于1 0 0d b ) ，而他们所测得的骨传导 超声的动态范围为1 3 2 3d b ，比正常气传导高频声音的动态范围窄。因此，如果 要研制骨传导超声助听器，需要对助听器的输出进行压缩。 超声的非线性解调 2 0 0 5 年k f u j i m o t o 通过测量纯音双边带调制超声的频率阈值( d i f f e r e n c e l i m e n sf o rf r e q u e n c y d l f ) ，发现调制纯音频率在l - 4k h z 的范围内d l f 和气传导 纯音一样较小( 低于1 ) ，但调制纯音频率低于0 2 5k h z 或高于6l d - i z 时，d l f 就 会增大。同时他们还分别测量了以单频超声作载波和叠加一定频带噪声的超声作 载波调制纯音后的d l f 。实验结果是：单频超声载波调制信号下的d l f 与气传导 纯音的d l f 具有相似的特性，这表明解调过程中包含有传统气传导的频率处理过 程。而对于叠加一定频带噪声的超声作载波调制下的信号，如果单靠线性解调是 很难将窄带的音频信息解调出来的，而一阶非线性解调就可以实现。此外，他们 还对比了双边带载波传输调幅( a md s b t c ) 信号和双边带抑制载波调幅信号 ( a md s b s c ) 信号的d l f ，结论为：超声a md s b t c 信号产生的音调对应的是 调制的气传导纯音的音调。而超声a md s b s c 信号产生的音调对应的是两倍频 率的气传导纯音音调。这些实验结果可以说明：骨传导超声的解调过程是非线性 过程【1 们。 骨传导超声的感知机理 关于骨传导超声的感知机理，目前存在很多猜测。1 9 8 5 年，k o h y a m a 等测 量了由骨传导超声引发豚鼠的皮层脑电图，可以看出它与直接对耳蜗神经加电激 励时的皮层脑电图相似。因此他们认为对骨传导超声作出响应的区域应该在耳蜗 里。他们同时还研究了注射了卡那徽素引起耳蜗外毛细胞受损的豚鼠的皮层脑电 图。将注射了卡那徽素的豚鼠和正常的豚鼠对比，其骨传导超声诱发的动作电位 并无明显差异。于是他们估计超声感知是通过内毛血细胞起作用的，而与外毛细 胞无关【】。以前有人认为：耳蜗受损患者的听力动态范围比听力正常的人要窄， 其中一个原因是由于外毛细胞的缺失【1 2 1 。使用人工耳蜗的患者的听力动态范围 也比较窄，这是由于内毛细胞或耳蜗神经是直接受到激励的【1 3 】。2 0 0 3 年t n i s h i m u r a l 澳1 量的动态范围显示骨传导超声的动态范围也较窄，且超声掩蔽音引 起的频率掩蔽范围较大。由此推断超声感知可能与外毛细胞无关。人耳耳蜗的特 征频率是不可能包含到高达1 2 0k h z 的高频声波的。超声可以被感知是因为它通 过引起基底膜的振动来刺激耳蜗底转的内毛细胞。如果超声的强度低，由于耳蜗 上的外毛细胞对超声的响应并不理想，骨传导超声无法激励内毛细胞，所以无法 感知。随着超声激励的强度增加，超声激励导致基底膜在耳蜗底转内强烈振动， 从而足够刺激内毛细胞产生感应。这些结果表明：骨传导超声刺激的是耳蜗底转 的内毛细胞【9 】。这一点与k o h y a m a 在豚鼠身上做实验得出的结论相符。 1 9 9 1 年，m l l e n h a r d t 等提出人对语音调制超声具有感知能力的同时，也提 出一种假设，他们认为球囊可能是超声感知器。球囊是一个对重力和直线加速度 做出响应的耳石器官，在耳蜗受损时，球囊负责声音传导。球囊和耳蜗之间存在 相互作用，而且球囊中的毛细胞长有一些纤毛，这些纤毛会对2 0 1 0 0k h z 的频率 振动做出响应【5 】。但是2 0 0 5 年k f u j i m o t o 提出：如果是球囊起作用的话，超声的 解调过程应该是线性过程，而这与后来提出的超声非线性解调是相矛盾【1 0 1 。此 外，在2 0 0 5 年k f u j i m o t o 的实验中，实验者可以同时感知到解调信号的音调和由 超声载波引起的高频音调( 1 0 1 4l ( i - t z ) 。他们证明通过一阶非线性解调可将超声调 制的语音信号解调成可听信号。但是他们仍没办法确认非线性解调的具体类型 t o l o 目前，尽管人们对于骨传导超声特性已有了一定的认识，但对于骨传导超声 的感知机理至今仍无定论，还有待更深入的研究。 1 2 本论文主要工作 本论文研究的主要目的是分析骨传导超声的特性，研究超声感知的机理，进 4 一步提高骨传导超声助听器的性能，使其在助听技术方面得到更广泛的应用。 本论文从以下三个方面展开研究： 1 ) 骨传导超声的计算机仿真。 利用数值计算方法，计算超声激励下人脑内的声场分布，研究骨传导超声在 人脑中的传播特性。 2 ) 骨传导超声助听器的系统设计。 根据骨传导超声助听的原理，设计基于d s p 平台的骨传导超声助听器设备， 为进行心理声学实验研究做准备。 3 ) 骨传导超声的实验研究 利用骨传导超声助听器进行心理声学实验，一方面验证计算机仿真的结果， 另一方面进一步发掘骨传导超声的特性，对骨传导超声助听器的设计提出改进方 案。 本论文的主要创新点在于通过计算机仿真实验和心理声学实验发现了骨传 导超声的一些新特性，将有助于对超声感知机理的进一步探究。具体表现为： 1 ) 利用时域有限差分方法分别计算骨传导超声激励和骨传导声音激励在人 脑内形成的声场。通过分析人耳耳蜗处的声压响应，发现骨传导超声可同时刺激 双耳耳蜗，且传输能力高于骨传导声音。同时，双耳激励的仿真实验结果表明： 双耳骨传导超声没有普通双耳听觉的定位能力。因此，单耳骨传导超声助听的效 果可能更为理想。实验进一步模拟了骨传导超声助听器的工作过程，分析得到了 骨传导超声助听的最佳载波频率。 2 ) 利用基于d s p 的骨传导超声助听器进行心理声学实验，实验结果表明： 通过该设备可进行一定的音频识别，这在国内的相关报导中尚属首次。并且，通 过实验预测了骨传导超声的解调方式，提出了更有利于提高助听效果的设计方 案。 第二章基于时域有限差分方法的超声骨传 导助听特性研究 第一章中介绍了骨传导超声的相关背景和研究现状。目前的大部分研究都是 基于人体实验结果的一些分析。本章将介绍一种利用计算机仿真技术研究骨传导 超声特性的方法。计算机仿真具有高效、安全、受环境条件约束较少及参数改变 方便等优点。仿真模型一般是对实际系统某些属性的逼近，利用这样的模型进行 实验，有助于我们对实际系统的认识，帮助我们解决实际问题。 本章利用时域有限差分方法，建立声波在人脑中传播的模型。通过分析声场， 得到骨传导超声的一些特性，从而指导骨传导超声助听系统的设计和改进。 2 1 时域有限差分方法运用于声场计算 2 1 1 时域有限差分方法介绍 时域有限差分( f i n i t ed i f f e r e n c et i m ed o m a i n ，f d t d ) 方法是19 6 6 年k s y e e 提出的电磁场数值计算方法【1 4 1 。f d t d 方法是求解麦克斯韦微分方程的直接时域 方法。对电磁场的e 、h 分量在空间和时间上采取交替抽样的离散方式，每一个 e ( 或h ) 场分量周围有四个h ( 或e ) 场分量环绕，应用这种离散方式将麦克斯韦旋 度方程转化为一组差分方程，并在时间轴上逐步推进地求解空间电磁场。这一方 法可以用于处理复杂形状目标和非均匀介质物体的电磁散射和辐射等问题。 近年来，越来越多的研究表明：f d t d 方法同样适用于声场的计算，有助于 人们了解复杂声场中的反射、散射和干涉等特性【1 5 椰1 。下面介绍f d t d 方法用 于声场数值计算的过程。 小振幅声波传播的波动方程为： 劾 1 一口 = v ，一一p( 2 1 ) - - j 戗kk a 哥 1 一 - = = 一一vp ( 2 2 ) 、， o t p 其中p 是声压，哥是质点振动速度，a 是衰减系数，r 是介质的压缩系数 ( k = l 2 ) ，p 是介质密度，c 是介质中声速。 在直角坐标系中( 2 1 ) 和( 2 2 ) 可写为： 6 望= 一上( 等+ 誓+ 誓) 一竺p ( 2 3 ) o tr 、a xa 1 ， a z 7 盯 、。7 ( 2 4 ) ( 2 5 ) ( 2 6 ) 声压分量和质点振速分量均采用间隔采样的方式，在时间和空间域中的离散 点用以下符号表示： p 一( f ，k ) = p i x ，i a y ，k a z ；n a t 】 以+ ( f 一五1 ，舭) = v 朋一三) 缸，脚，脸；( 玎+ 三) & 】 ( f ，_ 一三1 ，栌v y t i a x , ( j 一三1 ) 每，；( 肼三) 刎 ( 2 7 ) ( 2 8 ) ( 2 9 ) ( f ，舭一主) = 1 ，： i a x , j a y , ( 七一) a z ；( 斛2 ) a t ( 2 1 0 ) 其中i 、j f 、k 表示三维空间坐标，a x 、a y 、a z 分别是工、y 和z 轴方向上的空间 离散间隔：刀是时间序号，a t 为时间离散间隔。 在f d t d 中声压和质点振速各分量的空间分布如图2 1 所示。 歹一 图2 1f d t d 的离散网格 f d t d 方法中称这种离散网格单元为y e e 元胞。其特点剐1 8 】： 声压和振速分量在空间上交替抽样。每个坐标平面上的声压分量的四周由振 速分量环绕，振速分量四周由声压分量环绕。 每个场分量，自身相距一个空间步长，声压与振速相距半个空间步长。声压 7 p x p一吵p一叨印一趴望砂印一跳 1一户lp_一p 一 一 一 = = = 饥百盟甜毗百 取刀时刻、磁振速取刀+ 去时刻。 声压的以+ l 时刻的值由，l 时刻值得到，振速刀+ i 1 时刻的值由疗一i 1 时刻的值 得到。 y e e 元胞内媒质只取一种。 这种声场各分量的空间相对位置适合于声波波动方程的差分计算，能够恰当 地描述声场的传播特性。此外，声压和振速在时间上交替抽样，使声波波动方程 离散以后构成显式差分方程，从而可以在时间上迭代求解，而不需要进行矩阵求 逆运算。因而，给定相应声场问题的初始值，利用f d t d 方法就可以逐步推进 地求得以后各个时刻空间声场的分布。 利用f d t d 方法将声波波动方程在时间和空间上离散化，从而得到差分迭代 方程： p 打+ 。( f ，_ ，七) = ( i - a - - a t ) p 月( f ，七) 一 笪【：；兰! ! 二圭：兰：! ! 二兰兰! ! 二主：兰：竺+ fr y ，2 ( + 主，七) 一屹2 ( 7 二主，j ，七三。 + 匕n + - - ：( f 至i - - 4 - ，工七) = 以 ( j 一三，七) 一 一一 l 匕2 ( f一，j ，七) = 屹2 ( j 一去，七) 一 去_ 七) _ p n ( f - u 】 + 吾( 瓦一吾，七) = 蟛一圭( 一吾，七) 一 面a t 叭f ，尼) 一朋，一1 ，七) 】 y 二+ 三( f ，k 一三1 ) = y 二一圭( _ ，七一丢) 一 y 二+ i ( f ，一i ) = y 二一i ( f ，_ ，七一去) 一 祟 p 一( f ，七) 一p 一( f ，k 一1 ) 】 五【p ”【2 ，- ，纠一川，_ ，一1 川 ( 2 1 1 ) ( 2 1 2 ) ( 2 1 3 ) ( 2 1 4 ) 2 1 2 数值稳定条件 f d t d 方法是以一组有限差分方程来代替连续的声波波动方程，即以差分方 程组的解来代替原偏微分方程组的解。只有离散后差分方程组的解是收敛和稳定 的，这种代替才有意义。收敛性是指当离散间隔趋于零时，差分方程的解在空间 任意一点和任意时刻都一致趋于原方程的解。稳定性是指寻求一种离散间隔所满 足的条件，在此条件下差分方程的数值解与原方程的严格解之间的差为有界【1 9 】。 要使f d t d 计算结果是收敛且稳定的，时间步长，与空间步长( 缸、) ，、止) 间应当满足c o u r a n t 稳定性条件【1 9 1 ，即： f s ，= = = = = = = ；：一 c 撇( ) 2 + ( ) 2 + ( 勉y 气、 其中c 懈为计算区域内包含的各种介质的声速最大值。 另外，由于对波动方程作差分近似，这种离散处理会导致波的数值色散，即 差分近似后波数七和角频率国之间不再是简单线性关系。这种色散与离散间隔 ( 缸、知、蚴有关。但由于介质本身是无色散的，差分近似的结果会引入误差。 为了减小这种数值色散带来的误差，f d t d 空间离散间隔也要满足一定的要求。 假设空间域上x 、y 和z 三个方向采用等间隔离散a x = a y = a z = - c ；，所关心频段的频 率上限为矗“，对应波长为a m i n ，则应满足【1 9 】： 6s 九咖| n q 1 6 ) 通常脸l o 。 2 1 3 吸收边界条件 通常f d t d 方法需要模拟的声场都是无限开域的。但由于f d t d 计算时，任 一时刻空间内每个网格的声场分量都要存贮下来供下一时间步计算之需。模拟的 声场越大，所需的存储空间越大。由于计算机容量的限制，f d t d 计算只能在有 限区域进行。为了能模拟无限开域的声波散射过程，在计算区域的截断边界处作 特殊处理，使得向边界面行进的波在边界处保持外向行进的特征，无明显的反射 现象，并且不会使内部空间的场产生畸变。这种边界条件称为吸收边界条件。 2 1 3 1 几类吸收边界条件 吸收边界条件从开始简单的插值边界，到后来的m u r 吸收边界【2 0 1 、廖氏吸收 边界条件2 1 1 ，再到近来广泛采用的完全匹配层( p e r f e c tm a t c hl a y e r , p m l ) 吸收边 界 2 2 1 。 下面简单介绍几种应用较为广泛的吸收边界条件。 9 1 1m u r 吸收边界条件2 0 】 设一维波动方程为： 磐o x 专誓o t = o ( 2 1 7 ) z 1 ，z z 、 将它分解为2 个单向波方程： r 警七鲁= 0 一 亿埘8r 二+ 一二= i z 1 ) i 缸 1 ，研 、7 1 塑一三鲤：o 、o xvo t 其解为： 蝴m ( 姓州，剖圳赶了1 掣l ( 2 1 9 ) 上式表达的是沿x 方向传播的波，该波无反射波。因此，只要在截断边界x = o 处满足( 2 1 8 ) 就不会发生射。将此条件用于y e e 元胞中，可得： q ) - ( 1 ) 一伊”( o ) = ! 兰【缈一+ 1 ( o ) 一9 ”( o ) 】 ( 2 2 0 ) 即： 伊州( o ) ：缈一( o ) ( 1 一祟) + 尝伊一( 1 ) ( 2 2 1 ) x x 当缸= v a t 时，有。 伊肿1 ( o ) = 9 一( 1 ) ( 2 2 2 ) 式( 2 2 2 ) 表明：时l 时刻的- - - - 0 处的场值等于刀时刻= 1 处的场值，此条件满足一 个时间步向边界x = o 处移动了一个网格，好像边界不存在，这就是一阶m u r 条 件。为了取得更高的精度，m u r 条件也有类似的二阶近似。 m u r 吸收边界条件的优势在于它实施简单，吸收效果较好。但是其问题在于 直角坐标系下采用y e e 网格划分，在角区域存在较大误差，而二阶近似尽管精度 较高，但编程复杂，且对三维的情况还可能出现结果发散的现象。 2 ) 廖氏吸收边界条件2 1 】 廖氏吸收边界条件是利用牛顿后向差分多项式在时间域和空间域对波函数 进行外插的结果。它比同阶的m u r 吸收边界条件反射小约一个数量级( 2 0d b ) ， 并且矩形计算区域的角点处也易于实现。其各阶吸收边界条件可用统一的公式表 示。由于其计算和推导过程较为繁琐，这里只给出廖氏吸收边界条件最终的公式 表达： 缈( x ，f + ，) = ( 一1 ) 川c 加( x - j v a t ，t - ( j - - 1 ) a t ) ( 2 2 3 ) 其中c ；5 为组合数，表示廖氏吸收边界条件的阶数。 3 ) p m l 吸收边界条件【2 2 l 1 0 前面介绍的几种吸收边界条件在f d t d 计算区域外边界存在o 5 0 o - - 5 的数 值反射。为了更好地模拟实际情况，更加有效的吸收边界条件仍是研究热点。1 9 9 4 年j p b e r e n g e r 提出了p m l 吸收边界条件来吸收外向电磁波。它将电磁场分量 在吸收边界区分裂，并能分别对各个分裂的场分量赋以不同的损耗。它相当于在 f d t d 计算区域的边界设置一种非物理的吸收媒质，以使任意入射角和任意频率 入射的平面波，投射到真空介质表面的反射系数的理论值都为0 。据报道，p m l 的反射系数是前述方法二阶或三阶吸收边界条件的1 3 0 0 0 ，总的网格噪声能量 是以前方法的1 1 07 。本论文采用了p m l 吸收边界条件，将其稍做变换运用于声 场的计算。 2 1 3 2p m l 的提出 卜回冗歹卜缁j p b e r e n g e r 针对电慨饭堤出的p m l 的于匪导过程。以二维t m 坡 为例。在笛卡尔坐标系中t m 波有三个场分量：凤、毋、匠，假设媒质的电导率 和磁导率分别为盯和矿，媒质的介电常数为国，磁导系数为肋。麦克斯韦方程组 可以写成： 岛鲁+ 葩：= 等一等 亿2 4 ， 。等矗h 工= 等 亿2 5 ) o 瓦王+ 仃工= 言 ( 2 鳓警“嘭= 一鲁 亿2 6 ， 如果满足条件： 旦：! ( 2 2 7 ) 一= 一 ， 6 0 , u o 即媒质的阻抗和真空的阻抗相等，当平面波垂直入射到真空介质表面时就不会 发生反射现象。 在p m l 媒质中，在x 和y 方向分别设置吸收媒质以便对风、马进行吸收， 所以相应地把磁场分量匠分解成艮和两个分量。于是，p m l 媒质中就有四 个电磁场分量，即凤、凰、坛、岛，相应的麦克斯韦方程组就变成： 鲁帆耻警 ( 2 2 8 ) 占。鲁一驴一警 ( 2 2 9 ) 。警“即一掣 ( 2 3 0 ) o 言+ 仃y 月l 一 ( 2 , a 0 譬+ o - , 一ny = 掣 ( 2 3 1 ) +y = j 上( 2 3 1 ) 式中的参数( q ，仃：，1 7 ，仃：) 分别是各向同性的电导率和磁导率。 p m l 在y e e 单元网格空间的构造如图2 2 所示。中间的真空区域为计算区域， 计算区域外包围着p m l 媒质，p m l 外包围的是理想导体。计算区域的左、右边 界，吸收材料是匹配的p m l ( o r , ，盯：，0 ，0 ) 媒质，它能让外行波无反射地通过自由 空间p m l 的分界面a b 和c d 。同样，计算区域上、下边界，吸收材料是匹配 的p m l ( 0 ，0 ，o r ，仃：) 媒质。四个角点，采用p m l ( o r , ，盯：，盯，盯：) 媒质。 图2 2p m l 媒质的设置 如p m l 厚度为a 距离交界面，处外行波的幅度可表示为： 少( ，) = ( o ) g - ( a 湖训知c ) 7( 2 3 2 ) 其中，e y g j , b 行波相对于交界面的入射角，提a ；或西。穿过p m l 后，波将被外 围的理想导体反射回来再次穿过p m l ，回到中间的真空区域。反射系数为： r ( o ) = e - 2 ( 口5 p 7 舻) 艿 ( 2 3 3 ) 当日= 叹时，反射系数为l ，此时电磁波会平行此面传播，最终被垂直此 面的p m l 媒质吸收。 由( 2 3 3 ) 可见，反射系数是甜的函数。因此，理论上若给定p m l 层衰减值 大小，可以尽可能地减小p m l 层的厚度以减小计算区域的大小。但是，实际上， 电导率跳变太大，在计算时会导致数值反射。因此，实际计算中p m l 的厚度须 1 2 取几个网格的厚度，且导电率从自由空间p m l 分界面的0 逐渐增加到最外层 边界的默。设距离p m l 分界面，处的导电率为a ( 一，则分界面反射系数为： r ( 乡) = 2 - 2 ( a c o s o 舻) 如办 嘶) = 仃一( 詈) 疗 将( 2 3 5 ) 代x ( 2 3 4 ) 得： ( 2 3 4 ) ( 2 3 5 ) r ( o ) = e - 2 3 口( 肘1 ) e o c ( 2 3 6 ) 由( 2 3 5 ) 式可看出，r = 0 时，萨0 ，对应于真空媒质：当严耐，o = a m 缸，对 应于p m l 的最外层。所以，先确定戤，以便确定电导率取多大时在最外层上 的反射量是可以容忍的。 在p m l 媒质层中，j p b e r e n g e r 建议采用指数差分格式。利用上述的嘲设 置，得到p m l 中的差分迭代方程： 珧胪p 吒“归黜+ 器 乒“2 ( f + l 2 歹) 一厶陀( f 一1 2 , j ) l( 2 3 7 ) 黜力_ e - 口, o , j ) t u s o 黜一岩 乒2 ( i , j + l 2 ) 一厶2 ( f ，一1 2 ) ( 2 3 8 ) 研川2 ( f ，j + 1 2 ) = e e - a y ( i , j + 1 1 2 ) a t 鳓r r 善n 州2 ( f ，j + 1 2 ) 1。一( f ，j + 1 1 2 ) a t l l 一号长而 磁( j + 1 ) + 鬈( ，一，+ 1 ) 一磁( f ，) 一鬈( 7 ，一，) 】( 2 3 9 ) 日：+ 1 7 2 ( f + 1 2 ，) = e - o ：件2 & 7 鳓h ；一1 7 2 ( f + 1 2 ，) + 葫丽郴三( 1 ，卅e 黔1 ，沪e 2 ( “) 一日( 例( 2 4 0 ) 类比于二维情况，三维p m l 中将6 个场分量分解成1 2 个分量。处理同二维， 此处不再详细介绍。 p m l 虽然是针对电磁波提出的，但是类推，它同样也适用于声场的计算。将 电磁场中的参变量做如下替换【2 3 1 ： e z ，一日j ，片j ，占，a ，o r ，矿 _ p ，v 工，r ，p ，口，口 替换后，j p b e r e n g e r 的p m l 吸收边界条件就可直接运用于声场。 1 3 棚_ e - a , ( i , j ) a , r p 抑一器 【+ 1 趁( i + l 2 , j ) 一佗( f 一1 2 , j ) l ( f - - e - a y ( t , j ) a t l a p y ( “) 一等 x r u n + l 坨( i , j + l 2 ) 一哆2 ( f ，j - l 2 ) v ，1 2 ( f 十1 2 ，n = e ( “u z , j ) ，p o + 1 1 2 , ) ) v 7 1 2 ( i + 1 2 ，) ( 2 4 1 ) ( 2 4 2 ) 一三三j 石j 了万 p ：( + 1 ，) + p ；( + 1 ，j ) 一p ：( ，) 一p ；( ，) 】( 2 4 3 ) 1 ，2 ( f ，j + 1 2 ) = p 一们p m 址7 p 1 ，2 ( f ，+ 1 2 ) 1 一d - a y ( 1 ，j + 1 1 2 ) a t l p ( 1 ，j + 1 1 2 ) 一三i r _ j 虿j i _ 矿【p ：( ，_ ，+ 1 ) + p ；( ，_ ，+ 1 ) 一p ：( ，) 一p ；( ，_ ，) 】( 2 4 4 ) 根据上节的推导，在声场的p m l 计算中，满足条件竺：生时，可使声波在 r p 边界处无反射。在电磁场的p m l 推导中，导电率瘌磁导率，从自由空间p m l 分界面的0 逐渐增加到最外层边界的戤和，m 默。在声场中，口和口存在着同 样的规律。设置m 层p m l 。 图2 2 中的l 和r 区域，口，和口：为0 ，吼和从内到外逐渐增加，按下式 设置： 口j ( f ，力= 口撇，( 旦尘! ) 2 , i = 1 ，肌一1 ，= m ，一m 一1 ( 2 4 5 ) 口：( f ，j f ) = 口二( 竺些) 2 r ，f ：0 ，朋一1 ，_ ，：朋，n ，一朋一1 ( 2 4 6 ) b 和t 区域，口，和口：为0 ，口，和口：从内到外逐渐增加 口，( f ，) = 口懈( ! l 二2 ) 2 ，= l ，朋一1 ，i ：肌，n x 一所一1 ( 2 4 7 ) 口：( f ，) = 口0 ( 旦业) 2 t = o ，肌一l ，汪棚，工一肌一1 ( 2 4 8 ) 四个角点处： 口，( ，) = 口一( 掣) 2 ，f = l ，朋，_ ，= l ，朋 m + 1 啪胪口一( 等) 2 ，= l ，棚= l ，册 口“_ ) = 口二( 型导) 2 ，污1 棚，_ ：l ，棚 历+ l 1 4 ( 2 4 9 ) ( 2 5 0 ) ( 2 5 1 ) a j ( ? ，) = 口：。( 竺己寺 ) 2 ，f = l ，竹，= 1 ，聊