鼻腔气流运动的个体化测量

卢九星 韩德民 张 罗

(首都医科大学附属北京同仁医院耳鼻喉头颈外科北京市耳鼻咽喉科研究所教育部耳鼻咽喉头颈科学重点实验室，北京100730)

鼻腔位于呼吸道的首端，功能有通气、过滤、清洁、加温、加湿、共鸣、反射及嗅觉，无一不依赖于呼吸气流通过鼻腔而加以实现。鼻腔鼻窦结构复杂，在人类进化过程中逐渐形成，具有物种特异性，并非简单的气流通道。以呼吸气流运动为切入点，找出气流运动与鼻腔解剖结构、功能状态之间的规律性关系，进而帮助评估鼻腔生理、病理状态，对上呼吸道的功能判断、疾病预防以及矫正手术的量化分析均有一定价值。

鼻腔空气动力学兴起于20世纪80年代［1］，目前的主要研究方法是数值模拟，利用有限元的方法对上呼吸道进行气流场分析、计算和预测。Tan等［2］根据20例鼻窦计算机断层扫描(computed tomography，CT)图像，三维重建后对气流场进行网格划分，求解纳维叶-斯托克斯(Navier-Stokes)方程，得到中国人正常鼻腔在吸气、呼气2个时相，稳态(气流速率不随时间变化)条件下的气流运动参数。Liu等［3］采用相似的方法，分析鼻中隔偏曲对气流的影响，发现偏曲位置会影响气流，若合并下鼻甲代偿性肥大，偏曲侧之凹面气流会出现明显变化。鼻腔空气动力学在不侵入气道的情况下得到气流场的全部信息，具有任何器械检查不可比拟的优势，但人为控制边界条件，使得模拟结果只具有参考性而不能反映真实情况，例如稳态条件的设定就简化了呼吸过程。且据目前文献［1-3］，往往一个样本都设置相同的边界条件，虽然如此设置可以充分显示不同解剖结构对气流的影响，但却不能个体化反映当时的气流情况。为了逐渐将鼻腔空气动力学应用于临床，将气流运动情况加入评估体系，本研究综合利用CT、鼻阻力仪和鼻声反射仪3种检查手段，进行上呼吸道气流场的数值模拟，力图对受试者进行个体化测量，真实反映气流场情况。

1 资料与方法

1.1 研究对象

受试者10名，汉族，19～21岁，平均年龄(19.5±10.97)岁。其中男性5名，平均(19.6±1.68)岁;女性5名，平均(19.4±2.32)岁。受试者身体质量指数(body mass index，BMI)在18～25之间，平均21.31±0.97，其中男性平均21.67±1.11，女性20.96±0.87。身体健康，五官端正，无先天性颅面部畸形，齿列整齐，无急慢性鼻炎、鼻窦炎、鼻外伤、鼻出血、鼻囊肿、明显鼻中隔偏曲穿孔、鼻塞、鼻异物、鼻畸形等，近3月未患急性鼻科疾病、无鼻腔局部用药史、无鼻堵。试验符合人体试验伦理学标准，受试者严格签署知情同意书。

1.2 医学检查

1)鼻功能检查:鼻阻力仪(NR6，英国GM公司)、鼻声反射仪(A1，英国 GM公司)。测试房间室温20℃ ～25℃，相对湿度50% ～60%，受试者静坐20 min后端坐位进行检查。鼻喷减充血剂后再做1次鼻阻力、鼻声反射检查，每次2项检查在5 min内完成。

2)鼻窦CT扫描:鼻功能检查后30 min内行CT(Brilliance 64，荷兰飞利浦公司)检查。扫描上呼吸道，层厚0.67 mm，软组织窗，向上包括额窦、向下包括甲状软骨，断层图像数据的保存格式为数字影像和通信标准(digital imaging and communications in medicine，Dicom)。

1.3 鼻腔空气动力学分析

将CT断层图像读入Mimics10.0(比利时Materialise集团)软件，建立三维模型，以Stl(Stereo Lithography)文件格式导出(图1)。

图1 Mimics软件界面Fig.1 Mimics software interface

采用Icem-cfd(Ansys12.1，美国Ansys公司)软件划分流场网格。读入Stl模型文件，在面部前方设置外流场，与呼吸面罩近似。整个上呼吸道流场划分非结构性网格，八叉树(Octree/Robust)方法，呼吸道网格最大边长1 mm，平滑网格至98%网格质量大于0.6 (图2)。

图2 Icem-cfd软件界面Fig.2 Icem-cfd software interface

网格文件读入Fluent(Ansys12.1，美国Ansys公司)软件。视气流为连续介质。外流场设为入口，取值标准大气压;呼吸道末端为出口，瞬态条件，代入气流速率曲线方程(图3);鼻腔壁无滑移边界条件。采用剪切压力传输(shear stress transport，SST)K-ω 模型［4］计算呼吸气流的三维流动，收敛残差0.001。

连续方程

动量方程

湍流能量方程

能量方程

标准K-ε方程转换为

SST 由此可得，φ1 为 K-ω，φ2 为 K-ε

图3 呼吸模式图Fig.3 Patterning breathing，hypotheses for asymmetric sine curve

运算结果进入Cfd-post(Ansys12.1，美国Ansys公司)软件行三维可视化处理，得到1个完整呼吸过程中，气流在整个上呼吸道中的流动物理参数(气压、速率、壁面剪切力)，及在整体和局部的详细分布图像。1个呼吸周期中，近似平均取9个时间点，顺序说明呼吸过程中气流的变化(图3);每例呼吸道，从鼻孔前至鼻咽，依据解剖特征取9个截面，说明流经时的气流状态(图4)。

1.4 统计学方法

应用医学统计分析软件SPSS 17.0(美国SPSS公司)对气流运动参数进行统计。计量数据用均数±标准差(±s)表示，配对数据比较用配对t检验，数值模拟参数和仪器实际测量值比较采用t检验。以P＜0.05为差异具有统计学意义。

图4 上呼吸道截面示意图Fig.4 The schematic section of upper respiratory tract

2 结果

2.1 验证有限元模型

利用鼻声反射仪测量受试者的鼻腔最小截面积和其距前鼻孔的距离，0～5 cm、2～5 cm、0～7 cm的鼻腔容积［6］，同一受试者两侧鼻腔形态差异无统计学意义(P＞0.05)。同时测量有限元模型的鼻阈截面积(图4，C截面)、C截面距前鼻孔的距离，距前鼻孔5 cm、2～5 cm之间、7 cm的模型体积。经比较两者之间的差异无统计学意义(均P＞0.05)(表1)。

表1 受试者的鼻腔形态数据Tab.1 The data of volunteers’nasal cavity

2.2 验证气流计算结果

利用鼻阻力仪，测量男性、女性受试者在鼻压差150 Pa、75 Pa时的鼻腔单侧阻力值。同一受试者两侧鼻腔阻力之间的差异无统计学意义(P＞0.05)。鼻阻力取决于鼻腔的气流量和气流产生的压力，鼻阻力=压力差/流量。

在鼻前庭处(图4，B截面)、后鼻孔处(图4，H截面)取1个呼吸周期的压力值，制成曲线，求出经鼻压差，再与流量相比得到鼻阻力值。计算得到的鼻阻力值和实际测量得到的鼻阻力值相比较，2者间的差异无统计学意义(P＞0.05)(表2)。

2.3 气流的运动情况

呼吸过程中，气流主要流经总鼻道和中、下鼻道，以层流为主，平静呼吸时未形成明显涡流，在呼吸道边缘，如鼻前庭顶、底、侧部，嗅裂、鼻底、鼻腔后上部、鼻甲表面，咽鼓管圆枕前、后等多处散在有小的反向流动或涡流形成，部位、大小不定。

表2 受试者鼻腔阻力Tab.2 Volunteers’nasal resistance (Pa·cm-3·s-1)

将呼吸气流在各截面处的压强、速率依时间顺序绘成曲线，详见图5。A～Ⅰ分别对应鼻孔前、鼻前庭、鼻阈至鼻咽中部9个截面(图5)。A为气流进入、呼出鼻孔的压强速率变化。B至G为左右鼻腔内气流的压强速率变化。整个呼吸周期内，鼻前庭(B)的气流压力都为负压。固有鼻腔(C～G)内的气流压强在吸入时相为负;呼出时相为正，且形成明显的平台期。呼吸道内的气流随着深入体内而逐渐平缓，速率变小。

呼吸周期内Ⅰ～Ⅸ时刻，气流在呼吸道各截面处的压强速率变化，详见图6。Ⅴ是吸呼转换时相，之前是吸气时相，之后是呼气时相。吸气时相(Ⅰ～Ⅳ)，气流速率进入鼻腔后逐渐减小;压力为负，接近体表而不断变小。呼气时相，固有鼻腔内的气流压力为正。

图5 一个呼吸周期内，呼吸道A～I截面处的压强、速率曲线Fig.5 One respiratory cycle，the pressure，velocity curve in the respiratory tract section A～Ⅰ

3 讨论

鼻腔空气动力学可以排除仪器检查中的各种不确定干扰因素，捕捉现有检查手段难以观测到的现象;其价值突出表现在能描述呼吸道内任意部位的气流运动。因此，验证数值模拟的真实性、结果的可信性极其必要，否则鼻腔空气动力学无异于美轮美奂的空中楼阁。

图6 Ⅰ～Ⅸ时刻，气流流经呼吸道A～I截面时的压强、速率曲线Fig.6 Time pointⅠ ～Ⅸ，the pressure，velocity curve in the respiratory tract section A ～ Ⅰ

本研究以Dicom连续数据集为图像基础，采用人-机交互的方式将上呼吸道和面部从背景图像中分割出来。图像分割技术是三维重建的基石［7］，目前的技术只能逼近而无法完全反映真实的解剖结构，对特定的解剖结构，计算机自动分割也很难达到视觉思维的准确性。因此，全自动图像分割很难实现，人-机合作十分必要，主要体现在自动分割算法的选择和结果的修正上。图像分割的准确性与操作者密切相关，结果具有不可重复性。本研究为保持一致性，尽量采用计算机自动分割。Mimics软件生成的三维模型是用三角面片(Stl格式)来描述的，一例呼吸道三维模型的三角面片数量平均2×105左右，局部细节的修整可以通过删改三角面片来完成，改动误差在10-6m2。

在Icem-cfd软件中对三维模型进行网格划分，四面体网格最大边长设为1 mm，能够捕捉呼吸道大于1 mm的解剖特征。通过鼻声反射仪的检测数据，验证模型的形态准确性。本研究得到的有限元模型在形态上能较真实地反映上呼吸道的解剖结构特征，与实际的结构轮廓有很好的相关性，应用此模型对呼吸气流进行数值模拟，其结果能够较好的符合实际。

鼻孔入口是一个曲面，很难进行准确的图像分割，外鼻及周围面部形态不仅影响呼吸气流进入鼻孔的方式，还会影响气流在呼吸道内的运动情况［8］。因此，本研究在鼻孔前加设外流场，形似呼吸面罩，直径9 cm左右，边界条件取一个大气压。如此设置，呼吸气流呈自然状态进入鼻孔;因受试者鼻周形态不同，而对气流的不同影响也都会被保留。

本研究比较由鼻阻力仪实际测得的鼻阻力值和由数据模拟结果计算得来的鼻阻力值，来验证数值模拟结果是否合理。男女受试者测定鼻阻力值时，佩戴面罩，呼吸力度较弱，潮气量较小，得到的气流平均流速也较小，由此计算出的气流速率曲线方程偏于正常值范围下游。据文献［4］，实际测量的呼吸速率方程曲线近似矩形，呼气时相长于吸气时相。类似研究中，呼吸模式多设定为两时相相等或不相等的三角形［9］、正弦曲线［10］。本研究将呼吸曲线设置为呼吸时相/吸气时相=1.2/1，曲线下面积相等的半正弦曲线(图3)，较为符合真实情况。

模型和数值模拟得到验证后，就可系统分析呼吸气流的运动情况。本研究选择的是正常鼻腔，对称的结构产生相似的气流，两侧的压力、速率曲线基本伴行。气流进入鼻腔，至鼻阈处压力陡然变大，占压力变化值的55% ～63%;速率在鼻阈处也达到峰值;呼气时，气流压力下降至鼻阈后略有抬升，在鼻孔前接近大气压;气流速率平稳增高，在鼻阈处有一坡值，随后下降出鼻腔。呼吸周期内，压力、速率曲线明显成双峰状，吸气时相压力、速率曲线的峰值和坡度大于呼气时相，说明呼气气流变化小于吸气时相，吸气时气流变化更剧烈。

本研究得到的呼吸道模型在形态上能真实的反映上呼吸道的解剖特征，应用此模型对上呼吸道流场进行数值模拟，结果较好的符合实际。

鼻腔空气动力学是生物力学的一个分支，主要研究呼吸过程中，鼻腔在同气体作相对运动情况下的受力特性，气体流动规律和伴随发生的物理化学变化;是在医学和流体力学的基础上，随着计算流体动力学的发展而成长起来的一个学科。其指导思想贯穿着医学唯物观:“没有功能的结构是没有意义的，没有结构的功能是不存在的”。鼻腔结构是鼻腔功能的物质基础，鼻腔功能是鼻腔结构的外在表现。结构变异首先导致鼻腔气流变化，功能随之改变;反之，气流变化会影响鼻腔功能的实现，同时也作用于鼻腔结构，并在时间维度上逐渐积累，最终导致结构的改变［11］。结构和功能两者矛盾运动，由量变达到质变，最终产生病理症状。“结构-功能-症状”三者的辩证关系体现了鼻腔疾病从发生到发展、量变到质变、微观到宏观的因果关系。呼吸气流是切入点，充分研究气流的运动规律，就可掌握鼻腔“结构-功能-症状”的辩证关系。

［1］Leong S C，Chen X B，Lee H P，et al.A review of the implications of computational fluid dynamic studies on nasal airflow and physiology［J］.Rhinology，2010，48(2):139-145.

［2］Tan J，Han D，Wang J，et al.Numerical simulation of normal nasal cavity airflow in Chinese adult:a computational flow dynamics model［J］.Eur Arch Otorhinolaryngol，2012，269(3):881-889.

［3］Liu T，Han D，Wang J，et al.Effects of septal deviation on the airflow characteristics:using computational fluid dynamics models［J］.Acta Otolaryngol，2012，132(3):290-298.

［4］Liu Y，Edgar A M，Matida，et al.Numerical simulation of aerosol deposition in a 3-D human nasal cavity using RANS，RANS/EIM，and LES［J］.J Aerosol Sci，2007，38:683-700.

［5］严碧涯，张有为，郁秀芳，等.我国健康成人呼吸气流速率曲线图的研究［J］.北京医学，1980，2(5):257-260.

［6］Hilberg O，Pedersen O F.Acoustic rhinometry:recommendations for technical specifications and standard operating procedures［J］.Rhinol Suppl，2000，16:3-17.

［7］苏秀云，刘蜀彬.Mimics 软件临床应用 计算机辅助外科入门技术［M］.北京:人民军医出版社，2011:73.

［8］Taylor D J，Doorly D J，Schroter R C.Inflow boundary profile prescription for numerical simulation of nasal airflow［J］.J R Soc Interface，2010，7(44):515-527.

［9］Zhang J，Liu Y X，Sun X Z，et al.Computational fluid dynamics simulations of respiratory airflow in human nasal cavity and its characteristic dimension study［J］.Acta Mech Sin，2008，24(2):223-228.

［10］Ishikawa S，Nakayama T，Watanabe M，et al.Visualization of flow resistance in physiological nasal respiration:analysis of velocity and vorticities using numerical simulation［J］.Arch Otolaryngol Head Neck Surg，2006，132(11):1203-1209.

［11］里晓红，杨克林.慢性鼻窦炎60例变应性因素分析［J］.临床误诊误治，2010，23(11):1028-1029.