本文介紹了一種界面友好可與醫院聯網的便攜式強迫振蕩呼吸阻抗測試系統。該系統采用STM32產生單頻率或復頻率的正弦振蕩信號，經放大電路放大后驅動揚聲器產生正弦振蕩氣流于受試者口腔中，應用STM32的模數據轉換器（ADC）模塊分別測量由壓力傳感器和流量傳感器得到的模擬電壓信號，再經運算后將得到的參數顯示在觸摸顯示屏上，同時可通過GPRS將數據上傳到上位機進行保存和顯示。最后使用模擬肺和征集志愿者等方式來評價該系統的可靠性。測試結果表明該系統穩定可靠，可以達到監測呼吸阻抗狀況的目的。
引用本文： 周垂柳, 王武, 謝聯昇, 曾碧新. 一種便攜式強迫振蕩呼吸阻抗測試儀的設計與實現. 生物醫學工程學雜志, 2016, 33(5): 951-957. doi: 10.7507/1001-5515.20160153 復制
引言
慢性阻塞性肺疾病(chronic obstructive pulmonary disease，COPD)和哮喘(asthma)等是常見的呼吸系統疾病。COPD是一種慢性肺部疾病，主要表現為呼出氣流受限不完全可逆和呼吸困難[1]。在我國，COPD在15歲以上的人群患病率為3%，全國患者約為2 500萬，在城市的十大疾病死因中居第四位，且患病率與死亡率仍呈上升趨勢[2]。哮喘是一種包含氣道重塑和氣道感染的慢性疾病，并和氣道過于緊縮和過分收縮聯系在一起，我國超過2 500萬人患有此病，每年死于哮喘病的人數高達60萬之多[3-4]。
目前臨床上主要用常規肺功能儀測定上述呼吸系統疾病，但由于這種測定需要被測試者高度配合，所以對嬰幼兒、癱瘓患者、殘疾人和高齡人群實施起來很困難，且檢測時需有專業技師重復細致的指導才能得到準確的結果。而強迫振蕩技術(forced oscillation technique, FOT)在測量中無需被測者高度配合，其基本設想是在呼吸系統中施加一定頻率的振蕩氣流于受試者的口腔并疊加在呼吸氣流之上，隨氣流進入氣道和肺組織，測量經氣道和肺組織吸收并折射后的振蕩壓力和振蕩流量，并由此計算出一系列有用的肺阻抗參數。目前，此測量方法被證明有很好的重復性且整個過程完全無創，適合所有患者，尤其是老人、兒童和重癥患者。使用單頻率模式可追蹤呼吸阻抗隨時間的變換，而復頻率模式可快速檢測各個頻率下的呼吸阻抗值，但國內的研究僅僅只停留在單一模式。鑒于此，本文提出一種集成單頻率模式和復頻率模式的強迫振蕩呼吸阻抗測試儀，在此基礎上增加了無線傳輸功能使測得的數據可以通過通用分組無線服務技術(general packet radio service，GPRS)上傳到上位機，為進一步實現遠程醫療奠定了基礎。
1 系統方案設計
系統由下位機和上位機組成，整體設計框圖如圖 1所示。下位機使用功率為80 W、阻抗為4 Ω的揚聲器、內徑為2.2 cm的通氣管路、內徑為2.1 cm的高慣性管路、流量傳感器和壓力傳感器來構建整個呼吸管路和采集系統。使用ARM STM32F4單片機為核心，利用STM32的模數轉換器(digital to analog converter，DAC)產生特定頻率的振蕩信號，經放大電路放大后激勵揚聲器振動。揚聲器振動箱體內的氣體運動產生振蕩氣流，經出口和連接管進入受試者的口腔并疊加在呼吸氣流之上，隨氣流進入氣道和肺組織，由氣道和肺組織吸收并折射后的壓力和流量信號，經濾波電路后由STM32的模數轉換器(analog to digital converter，ADC)轉換成數字信號。系統在單頻率模式下用移動濾波算法和互相干算法算出每個時刻對應的呼吸阻抗(resistance，Rrs)和呼吸電抗(reactance，Xrs)值并以波形的形式實時顯示在觸摸顯示屏上，檢測結束后顯示整個階段的Rrs和Xrs平均值。系統在復頻率模式下用傅里葉變換將時域信號轉化成頻率信號，利用功率譜方法算出對應頻率點的Rrs和Xrs，并在觸摸顯示屏上顯示出頻譜分析圖、振蕩頻率為5 Hz和20 Hz下的呼吸阻抗R5和R20、振蕩頻率為5 Hz下的呼吸電抗X5以及共振頻率Fres等。利用GPRS模塊與上位機通信，建立連接后將數據上傳到上位機進行進一步的處理。
圖1
系統結構框圖
Figure1.
System scheme
2 電路設計
系統電路主要包括STM32控制模塊、振蕩信號發生模塊、流量和壓力采集模塊、外部靜態隨機存取存儲器(static random access memory, SRAM)模塊、液晶觸摸顯示模塊、GPRS無線傳輸模塊和電源供電模塊。
2.1 STM32控制模塊
本文使用ST公司基于Cortex-M4內核的32位STM32F407ZGT6單片機完成對氣體流動引起的壓力差電信號的采集、模數轉換、數據的分析、數據的無線發送、控制對外部存儲器的讀寫以及液晶觸摸顯示等功能。該單片機最高工作頻率168 MHz，內帶1 MB閃存，具有可變靜態存儲控制器(flexible static memory controller，FSMC)、定時器和3個ADC模塊和2個通道DAC等外設資源。本文直接使用內置的12位DAC產生單頻率和復頻率信號以驅動揚聲器，用內置的2個12位ADC采集氣體壓力和流量信號。
2.2 振蕩信號發生模塊
振蕩信號發生模塊由DAC信號發生接口、信號放大電路、重低音揚聲器、圓錐箱體和氣體管路組成。DAC發生的信號經放大電路后激勵揚聲器產生氣體振蕩流，推動箱體內的氣流進入到管路內，再通過管路進入人體的呼吸道。
2.3 流量和壓力采集模塊
2.3.1 流量傳感器和壓力傳感器的選擇
流量和壓力采集模塊由流量傳感器、壓力傳感器、低通濾波電路和ADC采樣電路組成。流量和壓力傳感器將氣體流量和壓力轉化成模擬電壓信號，經二階低通濾波電路濾除外界干擾后，由ADC采樣電路采集并轉化成可處理的數字信號。
使用的流量傳感器是由矽翔公司生產的FS6022系列一次性呼吸傳感器，可測量滿量程為0~5 L/s，電壓信號輸出范圍為0~5 V，其中0~2.5 V為正向流量，2.5~5.0 V為反向流量，電壓輸出與壓力呈線性關系。
使用的壓力傳感器是由美國SMI公司生產的SM5852壓差式傳感器，該傳感器具有多階壓力修正和溫度補償雙列，確保了采集壓力數據的線性和可靠性。其測量滿量程是0~0.15 PSI，電壓輸出0.50~4.5 V，其中0.5~2.5 V為正壓差，2.5~4.5 V為負壓差。
2.3.2 濾波電路設計
有用的呼吸信號主要集中在1~30 Hz，因此需在流量傳感器和壓力傳感器的輸出端連接低通濾波電路以濾除高頻信號干擾和50 Hz的市電干擾。本文采用的是二階巴特沃思低通濾波電路，使用單電源電平轉換芯片MAX232將+5 V的輸入電壓轉化成10 V電壓供OP07D放大器使用；由于傳感器輸出的電壓大于ADC采樣引腳可接受的電壓范圍(0~3.3 V)，故在其濾波電路后添加分壓電路；最后得到的二階低通濾波電路的截止頻率fc=38 Hz。
2.3.3 ADC采樣電路
STM32自帶8通道、12位精度的ADC轉換電路。ADC轉換的模擬信號包括氣路壓力的模擬電壓信號和氣路流量的模擬電壓信號，經二階低通濾波電路濾波后分別由STM32外部復用端口PC0和PA0采集。采樣頻率由系統滴答定時器(SYSTICK)配置生成，單頻率模式下振蕩周期5 Hz，采樣率100 Hz，復頻率模式下振蕩周期1~20 Hz，采樣率128 Hz。
2.4 外部SRAM
STM32F407ZGT內部具有容量為192 KB的隨機存取存儲器(random access memory，RAM)，但由于強迫振蕩需采集大小16 bit，以10 ms為間隔共32 s的流量和壓力數據，在求解傅里葉變換中還需存儲大量的中間變量，遠超內部自帶的RAM容量，故需外擴SRAM。系統選用ISSI公司生產的16位寬、1 024 K容量的靜態內存芯片IS62WV102416，將該芯片連接到STM32的FSMC接口上，通過對FSMC控制器的設置可以實現在不增加外部器件的情況下對存儲器的直接擴展。
2.5 液晶觸摸顯示模塊
使用瑞佑公司生產的256/64 K色帶有文字和繪圖模式的雙圖層液晶顯示驅動器RA8875，其支持16位的8080微處理機接口的傳輸協議。配置FSMC接口來模擬8080傳輸協議，選擇FSMC_NE4作為片選，再加上RA8875自帶觸摸接口，配置完成即可通過FSMC訪問RA8875相關寄存器來顯示文字、圖像和波形數據并可讀取其觸摸信號。
2.6 GPRS模塊
GPRS是在現有的全球移動通信系統上發展出來的一種新的分組數據承載業務，提供端到端的連接和廣域的無線IP連接。系統采用的GPRS為SIMCom公司的SIM900A模塊，該模塊體積小、功耗低，內嵌TCP/IP協議。STM32處理器與無線模塊的物理接口為RS232，通過發送相應的“AT”指令即可完成對模塊的操作。
3 算法設計
強迫振蕩呼吸阻抗測量原理如圖 2所示，圖中人體呼吸系統總阻抗為Zrs(t)，人體與振蕩源相連管路的呼吸阻抗為Z1，人體與外界相連的呼吸管路的阻抗為Z2。揚聲器產生高頻可控壓力源Pap(tf)由振蕩箱體流經呼吸管路后到達受試者內部，再提取經氣道和肺組織吸收并折射的振蕩壓力和振蕩流量信號分別為p(t, tf)和q(t, tf)。Q(t, tf)為流量傳感器響應的流量，P(t, tf)為壓力傳感器響應的壓力，根據呼吸力學模型可以得到式(1)。
圖2
強迫振蕩呼吸阻抗測量原理圖
Figure2.
Schematic figure of forced oscillation respiratory resistance detection
| $\begin{align} & q\left( t,{{t}_{f}} \right)=\frac{1}{Z\text{rs+}{{Z}_{1}}}{{P}_{\text{ap}}}\left( {{t}_{f}} \right) \\ & p\left( t,{{f}_{f}} \right)=\frac{Z\text{rs}}{Z\text{rs+}{{Z}_{1}}}{{P}_{\text{ap}}}\left( {{t}_{f}} \right) \\ \end{align}$ |
由式(1)可知，雖然q(t, tf)和p(t, tf)都受振蕩源和呼吸系統的影響，但其比值只受呼吸阻抗Zrs的影響，而流量傳感器響應的流量和壓力傳感器響應的壓力為人體的呼吸信號與強迫振蕩信號之和，其表達式如下：
| $\begin{align} & Q\left( t,{{f}_{f}} \right)=Q\left( t \right)+q\left( t,{{t}_{f}} \right) \\ & P\left( t,{{t}_{f}} \right)=P\left( t \right)+p\left( t,{{t}_{f}} \right) \\ \end{align}$ |
為了得到呼吸阻抗Zrs，需要將呼吸信號從中分離出來。
3.1 單頻率算法設計
單頻率檢測是在單頻率振蕩源下分析呼吸阻抗隨時間發生的變化，信號在時域上分析處理，從而可以更直觀地反映呼吸阻力在呼吸周期內發生的變化[5]。計算過程為：首先利用移動平均濾波法消除呼吸成分對呼吸阻抗計算的影響和外界的噪聲干擾，再利用互相干法計算呼吸阻抗、呼吸電抗和呼吸總阻抗參數。
3.1.1 移動平均濾波
移動平均濾波技術是對有限時間內的數據信號進行求和取平均；首先取平均時間窗窗長為一個振蕩周期，其對應的采樣點數為N，其次對振蕩周期內的采樣點進行累加取平均得到ΣP(i)/N，再用該點原始采樣值與該點振蕩周期內的平均值進行差值得到新的數值P(i)-ΣP(i)/N。該數值即為濾除呼吸信號后的強迫振蕩信號，其流程圖如圖 3所示。
圖3
移動平均濾波法原理圖
Figure3.
Schematic figure of moving average filtering method
3.1.2 互相干法
互相干法是以離散序列的離散傅里葉級數為理論基礎，以其基波分量的傅里葉級數近似代替基波與各次諧波的傅里葉級數。利用三角函數形式的傅立葉展開式，因此振蕩壓力信號可表達為：
| $P\left( t \right)={{a}_{0}}+\sum\limits_{k=1}^{\infty }{\left[ {{a}_{k}}\cos \left( \frac{2\pi nt}{T}+{{b}_{k}}\sin \frac{2\pi nt}{T} \right) \right]}$ |
k為k階諧波(k=1, 2, 3, …)，a0為直流分量，ak和bk為實數，與系統對外加信號的響應有關，其值由下式確定：
| $\begin{align} & {{a}_{k}}=\frac{2}{T}\int\limits_{t}^{t+T}{p\left( t \right)}\cos \left( \frac{2\pi kt}{T} \right)\text{d}t \\ & {{b}_{k}}=\frac{2}{T}\int\limits_{t}^{t+T}{p\left( t \right)}\sin \left( \frac{2\pi kt}{T} \right)\text{d}t \\ \end{align}$ |
運用互相干法原理，以每一時刻的基波分量(k=1)的傅里葉級數近似替代基波與各次諧波的傅里葉級數。設一個振蕩周期T內有N個采樣點，故在t時刻有[6]：
| $\begin{align} & {{A}_{t}}=\frac{2}{N}\int\limits_{t-\frac{2}{N}}^{t+\frac{2}{N}}{p\left( t \right)}\sin \left( \frac{2\pi t}{N} \right) \\ & {{B}_{t}}=\frac{2}{N}\int\limits_{t-\frac{2}{N}}^{t+\frac{2}{N}}{p\left( t \right)}\cos \left( \frac{2\pi t}{N} \right) \\ \end{align}$ |
在t時刻的振蕩信號振幅和相位為：
| $\begin{align} & {{\left| p \right|}_{t}}=\sqrt{A_{t}^{2}+B_{t}^{2}} \\ & \theta \left( {{p}_{t}} \right)=\arctan \left( \frac{{{B}_{t}}}{{{A}_{t}}} \right) \\ \end{align}$ |
取移動平均濾波后的采樣點t，并以該點為中心將前后半個振蕩周期內的采樣點數內每點乘以sin(2πt/N)并累加，再除以2/N得到At，按照同樣的思路可以得到該點的Bt，At和Bt的平方并相加取平方根后得到該點的幅值｜p｜t，由arctan(Bt/At)可以計算出該點的相位值，其互相干法流程示意圖如圖 4所示。
圖4
互相干法原理圖
Figure4.
Schematic figure of cross-correlation
按照同樣的方法可以計算出流量振幅｜q｜t和流量相位θ(qt)。再按照公式(7)可計算出相應呼吸總阻抗和相位信息。
| $\begin{align} & {{\left| Z \right|}_{t}}=\frac{{{\left| p \right|}_{t}}}{{{\left| q \right|}_{t}}} \\ & {{\Phi }_{t}}=\theta \left( {{p}_{t}} \right)-\theta \left( {{q}_{t}} \right) \\ \end{align}$ |
由極坐標公式可以得到：
| $\begin{align} & R\text{rs=}{{Z}_{t}}*\cos \left( {{\Phi }_{t}} \right) \\ & X\text{rs=}{{Z}_{t}}*\sin \left( {{\Phi }_{t}} \right) \\ \end{align}$ |
依照上述公式可計算出每個時刻的Rrs和Xrs。
3.2 復頻率算法和設計
頻譜分析法即以橫坐標表示頻率，縱坐標表示對應頻率的函數幅值，可簡便描述各頻率下的頻譜特性。而從時域到頻域，需使用頻譜分析技術即快速傅里葉變換，設一個離散的時域信號x(n)采樣點數為N，離散傅里葉變換X(k)方程為[7]：
| $X\left( k \right)=\frac{1}{N}\sum\limits_{n=0}^{N-1}{x\left( n \right){{e}^{-j2\pi nt/N}}=A\left( k \right)+jB\left( k \right)}$ |
式(9)中，X(k)為頻域值，k為頻率值的索引k=(0, …, N 1)，A(k)為傅里葉變換的實部，B(k)為傅里葉變換的虛部。將采集到的數據進行基4快速傅里葉變換，則得到由各個頻率量的實部和虛部組成的數組，并從數組中按其對應下標序列來提取出各個頻率量的實部和虛部。由于振蕩成分中混入一定的呼吸高頻成分，為避免呼吸高頻成分的影響，本文采用近似替代法，即用該頻率點的往左和往右的第二位信號相加取平均來近似替代在該頻率量的呼吸高頻成分。用頻譜中的該點減去該呼吸高頻成分，便可提取出振蕩波成分，其原理圖如圖 5所示。
圖5
提取振蕩成分原理圖
Figure5.
Schematic figure of extracting oscillatory component
設提取到的壓力和流量的頻譜實部分別為Ap(fk)、Aq(fk)，虛部分別為Bp(fk)、Bq(fk)，則可得到其自功率譜和互功率譜：
| $\begin{align} & {{G}_{\text{pp}}}=\left( {{A}_{\text{p}}}+\text{j}{{B}_{\text{p}}} \right)*\left( {{A}_{\text{p}}}-\text{j}{{B}_{\text{p}}} \right)=A_{\text{p}}^{2}+B_{\text{p}}^{2} \\ & {{G}_{\text{qq}}}=\left( {{A}_{\text{q}}}+\text{j}{{B}_{\text{q}}} \right)*\left( {{A}_{\text{q}}}-\text{j}{{B}_{\text{q}}} \right)=A_{\text{q}}^{2}+B_{\text{q}}^{2} \\ & {{G}_{\text{pq}}}=\left( {{A}_{\text{p}}}{{A}_{\text{q}}}+{{B}_{\text{p}}}{{B}_{\text{q}}} \right)+\text{j}\left( {{A}_{\text{p}}}{{B}_{\text{q}}}-{{B}_{\text{p}}}{{A}_{\text{q}}} \right) \\ \end{align}$ |
上述式中，Gpp為呼吸壓力的自功率譜；Gqq為呼吸流量的自功率譜；Gqp為呼吸壓力與流量的互譜；其相應頻率f下的阻抗的幅值和相位角由下列式子求出：
| $\begin{align} & {{Z}_{f}}=\frac{{{G}_{\text{pp}}}}{\left| {{G}_{\text{qp}}} \right|}=\frac{A_{\text{p}}^{2}+B_{\text{p}}^{2}}{\sqrt{{{\left( {{A}_{\text{p}}}{{A}_{\text{q}}}+{{B}_{\text{p}}}{{B}_{\text{q}}} \right)}^{2}}+{{\left( {{A}_{\text{p}}}{{B}_{\text{q}}}-{{B}_{\text{p}}}{{A}_{\text{q}}} \right)}^{2}}}} \\ & {{\theta }_{f}}=-\arctan \left( \frac{{{A}_{\text{p}}}{{B}_{\text{q}}}-{{B}_{\text{p}}}{{A}_{\text{q}}}}{{{A}_{\text{p}}}{{A}_{\text{q}}}+{{B}_{\text{p}}}{{B}_{\text{q}}}} \right) \\ \end{align}$ |
再由下式可得到相應頻率下的呼吸阻抗和呼吸電抗:
| $\begin{align} & {{R}_{f}}=\left| {{Z}_{f}} \right|\cos {{\theta }_{f}} \\ & {{X}_{f}}=\left| {{Z}_{f}} \right|\sin {{\theta }_{f}} \\ \end{align}$ |
由于測量呼吸系統壓力和流量信號時存在自主呼吸源等無關的噪聲干擾，因此還應當引入頻域相關函數r2(Coherence function)來計算其準確性
| ${{r}^{2}}=\frac{{{\left| {{G}_{\text{qp}}} \right|}^{2}}}{{{G}_{\text{pp}}}{{G}_{\text{qq}}}}$ |
r2是表征線性系統輸入、輸出關系的指標，r2值介于0~1之間，通常r2 > 0.95時所測得參數值才是可靠的。
4 系統性能測試
系統實物如圖 6所示，將管路連接到呼吸管路接口后在顯示屏上觸摸相應的測試選項即可進行測試。為了評價整個系統的可靠性還需進行一些測試驗證，包括機械阻抗模擬測試、模擬肺重復性測試以及初步在體測試。
圖6
系統實物圖
Figure6.
Photo of the system
4.1 機械阻抗測試
為了驗證系統的可靠性，將塑料軟管連接至呼吸管路接口上測定，根據空氣動力學原理，氣流流過圓管時的流體雷諾數(Reynolds number，Re)方程為[8]：
| $Re=\frac{\rho dQ}{\eta S}*{{10}^{3}}$ |
式中d：塑料管直徑(cm)，Q流量(L/s)，S：管截面積(cm2)，ρ：氣體密度(g/cm3)，η：氣體粘滯系數(g/cm)。在測試中我們使用的流量不超過0.2 L/s，選用的管直徑d=1.5 cm，取ρ空氣=1.29×10-3 g/cm3，η=1.81×10-4 g/cm代入上式，得雷諾系數Re < 1997 < Re(K)(對圓管，臨界雷諾數Re(K)一般取2 300)，因此在層流范圍內。由公式L=ρl/r2，R=8ηL/πr4(其中：R為管阻力，L為管慣性，η為粘滯系數，r為管半徑，l為管長度)，可以計算阻抗。但Z與氣體成分、溫度以及管子連接等機械結構有關，絕對值難以與實測值精確比較，而對與同一根塑料管，Z正比于l，故可依此做相對檢驗，測量的不同長度對應的Z值如表 1所示。
表1
不同長度對應的Z值
Table1.
Corresponding Z values of different lengths
用Matlab軟件對呼吸阻抗測試儀測得的Z值對長度l做回歸得：Z=1.093+0.049 4*l, r=0.993。可見阻抗值Z與阻抗管長度成正比，兩者的線性相關性為r=0.993。
4.2 用模擬肺和志愿者測試可靠性
為了驗證系統的重復性，將呼吸機和模擬肺連接至系統，設置呼吸機為CPAP模式，壓力6 cm H2O，呼吸頻率為12次/min，用復頻率方式測試其32 s內的平均阻抗值，重復5次，得到測量值如表 2所示。用統計軟件SPSS 19.0對數據進行相關分析，得到的參數用均值±標準差(x±s)形式，可見基本上滿足一致性要求。
表2
用復頻率方式測試模擬肺呼吸阻抗參數
Table2.
Resistance parameters of simulated lung respiration obtained using complex frequency method
為了驗證系統的有效性，征集2組成年人，第一組為哮喘病患者5例，第二組為正常人5例，分別在本文研制的系統上用單頻率模式下測試3次，數據結果采用均值±標準差形式，如圖 7所示。
圖7
5例哮喘患者和5例正常人的測試結果(x±s)
Figure7.
Test results of 5 patients with asthma and 5 healthy people (x±s)
從圖中可以看出：哮喘患者的Rrs和Xrs明顯高于正常人，對兩組數據進行配對t檢驗，均得到P < 0.01，表明兩組數據差異有統計學意義。
5 總結
本文以STM32單片機為主控芯片，控制呼吸流量信號的采集、處理、顯示和無線傳輸。實驗結果表明儀器不僅可以使用單頻率模式追蹤呼吸阻抗隨時間的變化，也可以使用復頻率模式快速檢測各個頻率上的呼吸阻抗并顯示出頻譜分析圖。儀器提供的UI界面可為用戶提供可視化操作，且測量的肺阻抗重復性高、可靠性好，具有監測肺部狀況的意義。此外該儀器還具有遠程傳輸的功能，可為進一步研究呼吸阻抗測試儀在家用市場上的發展提供參考。
引言
慢性阻塞性肺疾病(chronic obstructive pulmonary disease，COPD)和哮喘(asthma)等是常見的呼吸系統疾病。COPD是一種慢性肺部疾病，主要表現為呼出氣流受限不完全可逆和呼吸困難[1]。在我國，COPD在15歲以上的人群患病率為3%，全國患者約為2 500萬，在城市的十大疾病死因中居第四位，且患病率與死亡率仍呈上升趨勢[2]。哮喘是一種包含氣道重塑和氣道感染的慢性疾病，并和氣道過于緊縮和過分收縮聯系在一起，我國超過2 500萬人患有此病，每年死于哮喘病的人數高達60萬之多[3-4]。
目前臨床上主要用常規肺功能儀測定上述呼吸系統疾病，但由于這種測定需要被測試者高度配合，所以對嬰幼兒、癱瘓患者、殘疾人和高齡人群實施起來很困難，且檢測時需有專業技師重復細致的指導才能得到準確的結果。而強迫振蕩技術(forced oscillation technique, FOT)在測量中無需被測者高度配合，其基本設想是在呼吸系統中施加一定頻率的振蕩氣流于受試者的口腔并疊加在呼吸氣流之上，隨氣流進入氣道和肺組織，測量經氣道和肺組織吸收并折射后的振蕩壓力和振蕩流量，并由此計算出一系列有用的肺阻抗參數。目前，此測量方法被證明有很好的重復性且整個過程完全無創，適合所有患者，尤其是老人、兒童和重癥患者。使用單頻率模式可追蹤呼吸阻抗隨時間的變換，而復頻率模式可快速檢測各個頻率下的呼吸阻抗值，但國內的研究僅僅只停留在單一模式。鑒于此，本文提出一種集成單頻率模式和復頻率模式的強迫振蕩呼吸阻抗測試儀，在此基礎上增加了無線傳輸功能使測得的數據可以通過通用分組無線服務技術(general packet radio service，GPRS)上傳到上位機，為進一步實現遠程醫療奠定了基礎。
1 系統方案設計
系統由下位機和上位機組成，整體設計框圖如圖 1所示。下位機使用功率為80 W、阻抗為4 Ω的揚聲器、內徑為2.2 cm的通氣管路、內徑為2.1 cm的高慣性管路、流量傳感器和壓力傳感器來構建整個呼吸管路和采集系統。使用ARM STM32F4單片機為核心，利用STM32的模數轉換器(digital to analog converter，DAC)產生特定頻率的振蕩信號，經放大電路放大后激勵揚聲器振動。揚聲器振動箱體內的氣體運動產生振蕩氣流，經出口和連接管進入受試者的口腔并疊加在呼吸氣流之上，隨氣流進入氣道和肺組織，由氣道和肺組織吸收并折射后的壓力和流量信號，經濾波電路后由STM32的模數轉換器(analog to digital converter，ADC)轉換成數字信號。系統在單頻率模式下用移動濾波算法和互相干算法算出每個時刻對應的呼吸阻抗(resistance，Rrs)和呼吸電抗(reactance，Xrs)值并以波形的形式實時顯示在觸摸顯示屏上，檢測結束后顯示整個階段的Rrs和Xrs平均值。系統在復頻率模式下用傅里葉變換將時域信號轉化成頻率信號，利用功率譜方法算出對應頻率點的Rrs和Xrs，并在觸摸顯示屏上顯示出頻譜分析圖、振蕩頻率為5 Hz和20 Hz下的呼吸阻抗R5和R20、振蕩頻率為5 Hz下的呼吸電抗X5以及共振頻率Fres等。利用GPRS模塊與上位機通信，建立連接后將數據上傳到上位機進行進一步的處理。
圖1
系統結構框圖
Figure1.
System scheme
2 電路設計
系統電路主要包括STM32控制模塊、振蕩信號發生模塊、流量和壓力采集模塊、外部靜態隨機存取存儲器(static random access memory, SRAM)模塊、液晶觸摸顯示模塊、GPRS無線傳輸模塊和電源供電模塊。
2.1 STM32控制模塊
本文使用ST公司基于Cortex-M4內核的32位STM32F407ZGT6單片機完成對氣體流動引起的壓力差電信號的采集、模數轉換、數據的分析、數據的無線發送、控制對外部存儲器的讀寫以及液晶觸摸顯示等功能。該單片機最高工作頻率168 MHz，內帶1 MB閃存，具有可變靜態存儲控制器(flexible static memory controller，FSMC)、定時器和3個ADC模塊和2個通道DAC等外設資源。本文直接使用內置的12位DAC產生單頻率和復頻率信號以驅動揚聲器，用內置的2個12位ADC采集氣體壓力和流量信號。
2.2 振蕩信號發生模塊
振蕩信號發生模塊由DAC信號發生接口、信號放大電路、重低音揚聲器、圓錐箱體和氣體管路組成。DAC發生的信號經放大電路后激勵揚聲器產生氣體振蕩流，推動箱體內的氣流進入到管路內，再通過管路進入人體的呼吸道。
2.3 流量和壓力采集模塊
2.3.1 流量傳感器和壓力傳感器的選擇
流量和壓力采集模塊由流量傳感器、壓力傳感器、低通濾波電路和ADC采樣電路組成。流量和壓力傳感器將氣體流量和壓力轉化成模擬電壓信號，經二階低通濾波電路濾除外界干擾后，由ADC采樣電路采集并轉化成可處理的數字信號。
使用的流量傳感器是由矽翔公司生產的FS6022系列一次性呼吸傳感器，可測量滿量程為0~5 L/s，電壓信號輸出范圍為0~5 V，其中0~2.5 V為正向流量，2.5~5.0 V為反向流量，電壓輸出與壓力呈線性關系。
使用的壓力傳感器是由美國SMI公司生產的SM5852壓差式傳感器，該傳感器具有多階壓力修正和溫度補償雙列，確保了采集壓力數據的線性和可靠性。其測量滿量程是0~0.15 PSI，電壓輸出0.50~4.5 V，其中0.5~2.5 V為正壓差，2.5~4.5 V為負壓差。
2.3.2 濾波電路設計
有用的呼吸信號主要集中在1~30 Hz，因此需在流量傳感器和壓力傳感器的輸出端連接低通濾波電路以濾除高頻信號干擾和50 Hz的市電干擾。本文采用的是二階巴特沃思低通濾波電路，使用單電源電平轉換芯片MAX232將+5 V的輸入電壓轉化成10 V電壓供OP07D放大器使用；由于傳感器輸出的電壓大于ADC采樣引腳可接受的電壓范圍(0~3.3 V)，故在其濾波電路后添加分壓電路；最后得到的二階低通濾波電路的截止頻率fc=38 Hz。
2.3.3 ADC采樣電路
STM32自帶8通道、12位精度的ADC轉換電路。ADC轉換的模擬信號包括氣路壓力的模擬電壓信號和氣路流量的模擬電壓信號，經二階低通濾波電路濾波后分別由STM32外部復用端口PC0和PA0采集。采樣頻率由系統滴答定時器(SYSTICK)配置生成，單頻率模式下振蕩周期5 Hz，采樣率100 Hz，復頻率模式下振蕩周期1~20 Hz，采樣率128 Hz。
2.4 外部SRAM
STM32F407ZGT內部具有容量為192 KB的隨機存取存儲器(random access memory，RAM)，但由于強迫振蕩需采集大小16 bit，以10 ms為間隔共32 s的流量和壓力數據，在求解傅里葉變換中還需存儲大量的中間變量，遠超內部自帶的RAM容量，故需外擴SRAM。系統選用ISSI公司生產的16位寬、1 024 K容量的靜態內存芯片IS62WV102416，將該芯片連接到STM32的FSMC接口上，通過對FSMC控制器的設置可以實現在不增加外部器件的情況下對存儲器的直接擴展。
2.5 液晶觸摸顯示模塊
使用瑞佑公司生產的256/64 K色帶有文字和繪圖模式的雙圖層液晶顯示驅動器RA8875，其支持16位的8080微處理機接口的傳輸協議。配置FSMC接口來模擬8080傳輸協議，選擇FSMC_NE4作為片選，再加上RA8875自帶觸摸接口，配置完成即可通過FSMC訪問RA8875相關寄存器來顯示文字、圖像和波形數據并可讀取其觸摸信號。
2.6 GPRS模塊
GPRS是在現有的全球移動通信系統上發展出來的一種新的分組數據承載業務，提供端到端的連接和廣域的無線IP連接。系統采用的GPRS為SIMCom公司的SIM900A模塊，該模塊體積小、功耗低，內嵌TCP/IP協議。STM32處理器與無線模塊的物理接口為RS232，通過發送相應的“AT”指令即可完成對模塊的操作。
3 算法設計
強迫振蕩呼吸阻抗測量原理如圖 2所示，圖中人體呼吸系統總阻抗為Zrs(t)，人體與振蕩源相連管路的呼吸阻抗為Z1，人體與外界相連的呼吸管路的阻抗為Z2。揚聲器產生高頻可控壓力源Pap(tf)由振蕩箱體流經呼吸管路后到達受試者內部，再提取經氣道和肺組織吸收并折射的振蕩壓力和振蕩流量信號分別為p(t, tf)和q(t, tf)。Q(t, tf)為流量傳感器響應的流量，P(t, tf)為壓力傳感器響應的壓力，根據呼吸力學模型可以得到式(1)。
圖2
強迫振蕩呼吸阻抗測量原理圖
Figure2.
Schematic figure of forced oscillation respiratory resistance detection
| $\begin{align} & q\left( t,{{t}_{f}} \right)=\frac{1}{Z\text{rs+}{{Z}_{1}}}{{P}_{\text{ap}}}\left( {{t}_{f}} \right) \\ & p\left( t,{{f}_{f}} \right)=\frac{Z\text{rs}}{Z\text{rs+}{{Z}_{1}}}{{P}_{\text{ap}}}\left( {{t}_{f}} \right) \\ \end{align}$ |
由式(1)可知，雖然q(t, tf)和p(t, tf)都受振蕩源和呼吸系統的影響，但其比值只受呼吸阻抗Zrs的影響，而流量傳感器響應的流量和壓力傳感器響應的壓力為人體的呼吸信號與強迫振蕩信號之和，其表達式如下：
| $\begin{align} & Q\left( t,{{f}_{f}} \right)=Q\left( t \right)+q\left( t,{{t}_{f}} \right) \\ & P\left( t,{{t}_{f}} \right)=P\left( t \right)+p\left( t,{{t}_{f}} \right) \\ \end{align}$ |
為了得到呼吸阻抗Zrs，需要將呼吸信號從中分離出來。
3.1 單頻率算法設計
單頻率檢測是在單頻率振蕩源下分析呼吸阻抗隨時間發生的變化，信號在時域上分析處理，從而可以更直觀地反映呼吸阻力在呼吸周期內發生的變化[5]。計算過程為：首先利用移動平均濾波法消除呼吸成分對呼吸阻抗計算的影響和外界的噪聲干擾，再利用互相干法計算呼吸阻抗、呼吸電抗和呼吸總阻抗參數。
3.1.1 移動平均濾波
移動平均濾波技術是對有限時間內的數據信號進行求和取平均；首先取平均時間窗窗長為一個振蕩周期，其對應的采樣點數為N，其次對振蕩周期內的采樣點進行累加取平均得到ΣP(i)/N，再用該點原始采樣值與該點振蕩周期內的平均值進行差值得到新的數值P(i)-ΣP(i)/N。該數值即為濾除呼吸信號后的強迫振蕩信號，其流程圖如圖 3所示。
圖3
移動平均濾波法原理圖
Figure3.
Schematic figure of moving average filtering method
3.1.2 互相干法
互相干法是以離散序列的離散傅里葉級數為理論基礎，以其基波分量的傅里葉級數近似代替基波與各次諧波的傅里葉級數。利用三角函數形式的傅立葉展開式，因此振蕩壓力信號可表達為：
| $P\left( t \right)={{a}_{0}}+\sum\limits_{k=1}^{\infty }{\left[ {{a}_{k}}\cos \left( \frac{2\pi nt}{T}+{{b}_{k}}\sin \frac{2\pi nt}{T} \right) \right]}$ |
k為k階諧波(k=1, 2, 3, …)，a0為直流分量，ak和bk為實數，與系統對外加信號的響應有關，其值由下式確定：
| $\begin{align} & {{a}_{k}}=\frac{2}{T}\int\limits_{t}^{t+T}{p\left( t \right)}\cos \left( \frac{2\pi kt}{T} \right)\text{d}t \\ & {{b}_{k}}=\frac{2}{T}\int\limits_{t}^{t+T}{p\left( t \right)}\sin \left( \frac{2\pi kt}{T} \right)\text{d}t \\ \end{align}$ |
運用互相干法原理，以每一時刻的基波分量(k=1)的傅里葉級數近似替代基波與各次諧波的傅里葉級數。設一個振蕩周期T內有N個采樣點，故在t時刻有[6]：
| $\begin{align} & {{A}_{t}}=\frac{2}{N}\int\limits_{t-\frac{2}{N}}^{t+\frac{2}{N}}{p\left( t \right)}\sin \left( \frac{2\pi t}{N} \right) \\ & {{B}_{t}}=\frac{2}{N}\int\limits_{t-\frac{2}{N}}^{t+\frac{2}{N}}{p\left( t \right)}\cos \left( \frac{2\pi t}{N} \right) \\ \end{align}$ |
在t時刻的振蕩信號振幅和相位為：
| $\begin{align} & {{\left| p \right|}_{t}}=\sqrt{A_{t}^{2}+B_{t}^{2}} \\ & \theta \left( {{p}_{t}} \right)=\arctan \left( \frac{{{B}_{t}}}{{{A}_{t}}} \right) \\ \end{align}$ |
取移動平均濾波后的采樣點t，并以該點為中心將前后半個振蕩周期內的采樣點數內每點乘以sin(2πt/N)并累加，再除以2/N得到At，按照同樣的思路可以得到該點的Bt，At和Bt的平方并相加取平方根后得到該點的幅值｜p｜t，由arctan(Bt/At)可以計算出該點的相位值，其互相干法流程示意圖如圖 4所示。
圖4
互相干法原理圖
Figure4.
Schematic figure of cross-correlation
按照同樣的方法可以計算出流量振幅｜q｜t和流量相位θ(qt)。再按照公式(7)可計算出相應呼吸總阻抗和相位信息。
| $\begin{align} & {{\left| Z \right|}_{t}}=\frac{{{\left| p \right|}_{t}}}{{{\left| q \right|}_{t}}} \\ & {{\Phi }_{t}}=\theta \left( {{p}_{t}} \right)-\theta \left( {{q}_{t}} \right) \\ \end{align}$ |
由極坐標公式可以得到：
| $\begin{align} & R\text{rs=}{{Z}_{t}}*\cos \left( {{\Phi }_{t}} \right) \\ & X\text{rs=}{{Z}_{t}}*\sin \left( {{\Phi }_{t}} \right) \\ \end{align}$ |
依照上述公式可計算出每個時刻的Rrs和Xrs。
3.2 復頻率算法和設計
頻譜分析法即以橫坐標表示頻率，縱坐標表示對應頻率的函數幅值，可簡便描述各頻率下的頻譜特性。而從時域到頻域，需使用頻譜分析技術即快速傅里葉變換，設一個離散的時域信號x(n)采樣點數為N，離散傅里葉變換X(k)方程為[7]：
| $X\left( k \right)=\frac{1}{N}\sum\limits_{n=0}^{N-1}{x\left( n \right){{e}^{-j2\pi nt/N}}=A\left( k \right)+jB\left( k \right)}$ |
式(9)中，X(k)為頻域值，k為頻率值的索引k=(0, …, N 1)，A(k)為傅里葉變換的實部，B(k)為傅里葉變換的虛部。將采集到的數據進行基4快速傅里葉變換，則得到由各個頻率量的實部和虛部組成的數組，并從數組中按其對應下標序列來提取出各個頻率量的實部和虛部。由于振蕩成分中混入一定的呼吸高頻成分，為避免呼吸高頻成分的影響，本文采用近似替代法，即用該頻率點的往左和往右的第二位信號相加取平均來近似替代在該頻率量的呼吸高頻成分。用頻譜中的該點減去該呼吸高頻成分，便可提取出振蕩波成分，其原理圖如圖 5所示。
圖5
提取振蕩成分原理圖
Figure5.
Schematic figure of extracting oscillatory component
設提取到的壓力和流量的頻譜實部分別為Ap(fk)、Aq(fk)，虛部分別為Bp(fk)、Bq(fk)，則可得到其自功率譜和互功率譜：
| $\begin{align} & {{G}_{\text{pp}}}=\left( {{A}_{\text{p}}}+\text{j}{{B}_{\text{p}}} \right)*\left( {{A}_{\text{p}}}-\text{j}{{B}_{\text{p}}} \right)=A_{\text{p}}^{2}+B_{\text{p}}^{2} \\ & {{G}_{\text{qq}}}=\left( {{A}_{\text{q}}}+\text{j}{{B}_{\text{q}}} \right)*\left( {{A}_{\text{q}}}-\text{j}{{B}_{\text{q}}} \right)=A_{\text{q}}^{2}+B_{\text{q}}^{2} \\ & {{G}_{\text{pq}}}=\left( {{A}_{\text{p}}}{{A}_{\text{q}}}+{{B}_{\text{p}}}{{B}_{\text{q}}} \right)+\text{j}\left( {{A}_{\text{p}}}{{B}_{\text{q}}}-{{B}_{\text{p}}}{{A}_{\text{q}}} \right) \\ \end{align}$ |
上述式中，Gpp為呼吸壓力的自功率譜；Gqq為呼吸流量的自功率譜；Gqp為呼吸壓力與流量的互譜；其相應頻率f下的阻抗的幅值和相位角由下列式子求出：
| $\begin{align} & {{Z}_{f}}=\frac{{{G}_{\text{pp}}}}{\left| {{G}_{\text{qp}}} \right|}=\frac{A_{\text{p}}^{2}+B_{\text{p}}^{2}}{\sqrt{{{\left( {{A}_{\text{p}}}{{A}_{\text{q}}}+{{B}_{\text{p}}}{{B}_{\text{q}}} \right)}^{2}}+{{\left( {{A}_{\text{p}}}{{B}_{\text{q}}}-{{B}_{\text{p}}}{{A}_{\text{q}}} \right)}^{2}}}} \\ & {{\theta }_{f}}=-\arctan \left( \frac{{{A}_{\text{p}}}{{B}_{\text{q}}}-{{B}_{\text{p}}}{{A}_{\text{q}}}}{{{A}_{\text{p}}}{{A}_{\text{q}}}+{{B}_{\text{p}}}{{B}_{\text{q}}}} \right) \\ \end{align}$ |
再由下式可得到相應頻率下的呼吸阻抗和呼吸電抗:
| $\begin{align} & {{R}_{f}}=\left| {{Z}_{f}} \right|\cos {{\theta }_{f}} \\ & {{X}_{f}}=\left| {{Z}_{f}} \right|\sin {{\theta }_{f}} \\ \end{align}$ |
由于測量呼吸系統壓力和流量信號時存在自主呼吸源等無關的噪聲干擾，因此還應當引入頻域相關函數r2(Coherence function)來計算其準確性
| ${{r}^{2}}=\frac{{{\left| {{G}_{\text{qp}}} \right|}^{2}}}{{{G}_{\text{pp}}}{{G}_{\text{qq}}}}$ |
r2是表征線性系統輸入、輸出關系的指標，r2值介于0~1之間，通常r2 > 0.95時所測得參數值才是可靠的。
4 系統性能測試
系統實物如圖 6所示，將管路連接到呼吸管路接口后在顯示屏上觸摸相應的測試選項即可進行測試。為了評價整個系統的可靠性還需進行一些測試驗證，包括機械阻抗模擬測試、模擬肺重復性測試以及初步在體測試。
圖6
系統實物圖
Figure6.
Photo of the system
4.1 機械阻抗測試
為了驗證系統的可靠性，將塑料軟管連接至呼吸管路接口上測定，根據空氣動力學原理，氣流流過圓管時的流體雷諾數(Reynolds number，Re)方程為[8]：
| $Re=\frac{\rho dQ}{\eta S}*{{10}^{3}}$ |
式中d：塑料管直徑(cm)，Q流量(L/s)，S：管截面積(cm2)，ρ：氣體密度(g/cm3)，η：氣體粘滯系數(g/cm)。在測試中我們使用的流量不超過0.2 L/s，選用的管直徑d=1.5 cm，取ρ空氣=1.29×10-3 g/cm3，η=1.81×10-4 g/cm代入上式，得雷諾系數Re < 1997 < Re(K)(對圓管，臨界雷諾數Re(K)一般取2 300)，因此在層流范圍內。由公式L=ρl/r2，R=8ηL/πr4(其中：R為管阻力，L為管慣性，η為粘滯系數，r為管半徑，l為管長度)，可以計算阻抗。但Z與氣體成分、溫度以及管子連接等機械結構有關，絕對值難以與實測值精確比較，而對與同一根塑料管，Z正比于l，故可依此做相對檢驗，測量的不同長度對應的Z值如表 1所示。
表1
不同長度對應的Z值
Table1.
Corresponding Z values of different lengths
用Matlab軟件對呼吸阻抗測試儀測得的Z值對長度l做回歸得：Z=1.093+0.049 4*l, r=0.993。可見阻抗值Z與阻抗管長度成正比，兩者的線性相關性為r=0.993。
4.2 用模擬肺和志愿者測試可靠性
為了驗證系統的重復性，將呼吸機和模擬肺連接至系統，設置呼吸機為CPAP模式，壓力6 cm H2O，呼吸頻率為12次/min，用復頻率方式測試其32 s內的平均阻抗值，重復5次，得到測量值如表 2所示。用統計軟件SPSS 19.0對數據進行相關分析，得到的參數用均值±標準差(x±s)形式，可見基本上滿足一致性要求。
表2
用復頻率方式測試模擬肺呼吸阻抗參數
Table2.
Resistance parameters of simulated lung respiration obtained using complex frequency method
為了驗證系統的有效性，征集2組成年人，第一組為哮喘病患者5例，第二組為正常人5例，分別在本文研制的系統上用單頻率模式下測試3次，數據結果采用均值±標準差形式，如圖 7所示。
圖7
5例哮喘患者和5例正常人的測試結果(x±s)
Figure7.
Test results of 5 patients with asthma and 5 healthy people (x±s)
從圖中可以看出：哮喘患者的Rrs和Xrs明顯高于正常人，對兩組數據進行配對t檢驗，均得到P < 0.01，表明兩組數據差異有統計學意義。
5 總結
本文以STM32單片機為主控芯片，控制呼吸流量信號的采集、處理、顯示和無線傳輸。實驗結果表明儀器不僅可以使用單頻率模式追蹤呼吸阻抗隨時間的變化，也可以使用復頻率模式快速檢測各個頻率上的呼吸阻抗并顯示出頻譜分析圖。儀器提供的UI界面可為用戶提供可視化操作，且測量的肺阻抗重復性高、可靠性好，具有監測肺部狀況的意義。此外該儀器還具有遠程傳輸的功能，可為進一步研究呼吸阻抗測試儀在家用市場上的發展提供參考。
