ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

毫米波雷达生命体征检测:从I/Q信号到呼吸心率算法实战

毫米波雷达生命体征检测:从I/Q信号到呼吸心率算法实战 毫米波雷达做心率、呼吸检测这个方向我踩了大半年的坑终于把从射频前端到算法输出的链路完整跑通了。网上能搜到的资料多半是演示视频或者零散的科普真正把为什么能测和怎么实现讲透的文章很少。这篇是系列第一篇先把原理和第一版实例交代清楚后面我会继续写算法优化和实测对比。文章会围绕24GHz多普勒雷达模块展开涉及I/Q信号采集、相位解缠、带通滤波、频谱峰值检测和BNR短时呼吸噪声比这些关键环节适合正在做非接触式生命体征检测、或者打算用毫米波雷达做呼吸心率监测的朋友参考。1. 先说原理雷达测的不是心跳是身体表面的毫米级起伏1.1 多普勒效应与微动信号有一个几乎所有人都会产生的误区觉得毫米波雷达能看见心跳甚至能穿透胸腔看到心脏搏动。实际完全不是这么回事。雷达真正捕捉到的是人体胸腔表面因为呼吸和心跳引起的物理位移。呼吸的时候胸廓起伏明显成年人的呼吸位移通常在1到12毫米之间心跳引起的体表振动就要微弱得多大约只有0.2到0.5毫米。这个数值个体差异很大但量级就是如此。雷达发射电磁波碰到胸腔表面反射回来。胸腔表面朝雷达方向移动时反射波的相位会提前远离雷达时相位会滞后。这个相位变化经过混频解调后就变成了一个低频交流信号。只要我们把这个信号解调出来就能还原出胸腔表面的微动波形进而从波形中分别提取呼吸频率和心跳频率。这里需要理解一个关键点心跳信号不是被雷达直接测到的而是通过测量皮肤表面微小的振动间接得到的。因为心脏泵血时动脉血管会产生搏动这个搏动会传递到身体表面尤其是胸壁、颈部和四肢末端。胸壁上的微振动虽然微弱但24GHz雷达的相位灵敏度足以捕捉到这种亚毫米级位移。1.2 为什么24GHz比更低频段更适合生命体征雷达工作频率的选择直接影响微动检测的灵敏度。相位变化量\Delta \phi和物体位移d之间的关系是\Delta \phi \frac{4\pi \cdot d}{\lambda}其中\lambda是雷达波长。物体位移d相同时波长越短相位变化越大检测灵敏度越高。24GHz雷达的波长约为12.5毫米0.5毫米的心跳位移对应约0.25弧度的相位变化这个量级对于现代ADC和数字信号处理来说已经足够明显。如果换成5.8GHz的微波雷达波长约52毫米同样的位移只产生约0.06弧度的相位变化信噪比就会差很多。那为什么不直接用77GHz毫米波雷达77GHz波长约3.9毫米灵敏度更高但模块价格高一个数量级开发门槛也高。24GHz是成本和灵敏度的均衡点市面上有大量成熟的低功耗模块适合做原型验证。60GHz和77GHz的优势主要体现在FMCW测距和更小的天线尺寸上这在后面的文章里我会详细对比。1.3 一个简单的数学关系相位变化量估算假设目标与雷达的距离为R0胸腔表面在雷达视轴方向上有微动位移x(t)那么接收信号的中频相位可以近似表示为\phi(t) \frac{4\pi (R_0 x(t))}{\lambda}去掉常数项后相位变化量正比于x(t)。所以只要我们准确提取了相位信号就能直接得到胸腔微动的波形。我习惯在代码里用这组换算关系做量级估计24GHz波长12.5mm如果呼吸幅度为5mm相位变化为5.028rad心跳幅度为0.3mm相位变化为0.302rad。这说明呼吸信号的相位幅度可能是心跳的16倍以上。这么悬殊的能量差异意味着如果滤波或者谐波抑制做得不好呼吸信号很容易污染心跳频段。后面我会专门讲这个坑。2. 硬件搭建选模块时最容易看走眼的三个参数2.1 我用的24GHz模块及接口市面上的24GHz微波雷达模块主要分两类多普勒雷达模块和FMCW雷达模块。多普勒模块输出的是I/Q两路基带信号直接反映目标运动的多普勒频率和相位信息FMCW模块内部集成了VCO和混频器能够测量距离但数据接口复杂得多。我选的是24GHz多普勒雷达模块芯片内部集成收发天线和混频器模块外部引出I、Q、GND、VCC5V等引脚。这类模块最大的好处是简单MCU只需要用ADC采集I和Q两路模拟电压即可不需要配置复杂的寄存器。不过简单也意味着代价——它没有距离分辨能力所有距离上的反射信号都会混叠在一起。后面你会看到这个特性会让运动干扰变得非常麻烦。挑选模块时有三个参数容易误导人输出幅度范围有些模块的I/Q输出是1V偏置叠加微弱的交流信号有些是差分输出。在选ADC输入范围时要留足余量避免信号截顶。天线波束角波束角决定了检测区域范围。波束太宽容易把周围环境反射收进来太窄则要求目标必须正对雷达。实测中90度波束角的模块在1米外已经能收到较强的多径反射。接口类型模拟输出和数字接口如SPI/I2C的区别。数字接口模块通常已经内置了信号调理但数据刷新率可能受限需要确认采样率是否满足生命体征信号需求。2.2 电源与天线布局对底噪的影响这是新手最容易忽略的地方。24GHz雷达模块的内部电路对电源纹波非常敏感因为电源纹波会通过混频器的本振通路耦合到基带输出上形成低频杂散。我第一次用开关电源给模块供电时I/Q输出上出现了一个约20Hz的周期性波动实际目标的生命体征信号才0.2到2Hz这个杂散正好落在附近几乎把呼吸信号淹没了。后来改成低压差线性稳压器LDO供电并且用LC滤波把数字地和模拟地分开底噪立刻下降了约12dB。模拟前端的地线布局也很关键雷达模块的接地引脚应该尽量短地连接到主板的模拟地平面避免和数字地环路形成串扰。天线布局方面雷达正前方30度范围内不应该有金属物体和大幅度的平面反射体。我最初把雷达放在桌面金属支架旁边结果相位信号里出现了一个固定频率的振荡分量后来挪到塑料支架上才消失。另外雷达模块和人体之间的距离建议控制在0.3到1.5米。太近时近场效应导致信号饱和太远时回波能量不够微动信号会被热噪声淹没。2.3 数据采集采样率与量化精度怎么定呼吸频率通常在0.1到0.5Hz心率在0.8到2.5Hz但这些信号的谐波成分可能达到5Hz以上。更重要的是时域波形中的谐波对于区分呼吸和心跳有帮助所以采样率不能卡在香农定理的临界点。我最终采用100Hz的采样率这样在2秒的窗口内可以完成一个心跳周期的分析在20秒窗口内做FFT时频率分辨率能达到0.05Hz。ADC位数推荐至少12位。以5V量程为例12位ADC的量化步长约1.22mV如果I/Q输出中0.3mm心跳对应的电压变化只有几十毫伏量化噪声还可接受。实测中心电模拟信号在I/Q输出的幅度约为50-200mV使用12位ADC能得到比较干净的波形。如果模块输出是1V偏置建议外部加一个减法器把直流偏置去掉把交流部分放大后再进ADC这样可以充分利用ADC的量程。3. 从原始I/Q到生命体征波形完整信号链路拆解3.1 去直流与I/Q校正拿到I/Q原始数据后第一步不是急着做FFT而是去直流。由于模块的混频器输出自带直流偏置而且这个偏置会随着温度和环境变化缓慢漂移直接取整个数据段的均值减掉是最简单的做法。但对于实时系统我建议用滑动平均或者一阶高通滤波截止频率设为0.05Hz左右既能去掉直流漂移又不会伤到0.1Hz的呼吸信号。I/Q校正也不能跳过。理论上I和Q应该是等幅、相位差90度的两路信号但实际模块由于射频走线不对称两路增益会有偏差正交相位也不会严格等于90度。如果不做校正相位解调出来的信号会出现镜像频率导致心率频带上出现假峰。最粗暴的校正方法是在雷达前方放置一个振动目标比如扬声器纸盆采集一段已知振动频率的数据然后计算I路和Q路的幅度比和相位差再用一个旋转矩阵把两路信号归一化到标准正交坐标系。3.2 相位提取与相位解缠相位提取很简单就是每个采样点执行atan2(Q, I)。但atan2返回的相位范围是-π到π如果真实相位变化超过这个范围就会发生跳变。呼吸幅度大时相位变化能轻松突破π所以必须做相位解缠。相位解缠的标准做法是比较相邻两个相位采样点的差值如果差值大于π就把后面的相位整体减去2π如果差值小于-π就加上2π。我在代码里这样写for (int i 1; i n; i) { float delta phase[i] - phase[i-1]; while (delta M_PI) delta - 2.0f * M_PI; while (delta -M_PI) delta 2.0f * M_PI; phase[i] phase[i-1] delta; }为什么用while而不是if因为存在极端情况如果某个采样点噪声很大相位可能瞬间跳变超过2π用if处理一次可能不够。while能保证最终相邻差值被约束在-π到π之间。解缠后的相位序列就是胸腔微动的时域波形后续处理都在这个序列上进行。3.3 带通滤波与呼吸/心跳分离在解缠后的相位信号上呼吸和心跳是叠加在一起的。呼吸频率低、幅度大心跳频率高、幅度小。用两个带通滤波器把它们分离呼吸带通0.1Hz - 0.5Hz心跳带通0.8Hz - 2.5Hz我使用四阶Butterworth滤波器离线分析时用零相位滤波filtfilt实时系统则用级联二阶IIR。这里有一个容易被忽略的问题滤波器的群延迟。零相位滤波会翻转两次信号没有相位失真但不能实时输出因果滤波器会有固定延迟对实时心率显示有影响需要根据相位延迟做补偿。如果只是离线分析优先推荐零相位滤波。滤波器的阶数也不是越高越好。太高阶会引入振铃在呼吸波形上表现为脉冲状假象。四到六阶是实际工程中比较合适的范围。3.4 频谱峰值检测与频率估计滤波完成后需要估计呼吸率和心率。最直观的方法是FFT加峰值搜索。以100Hz采样率采集2048点约20.48秒做FFT后频率分辨率为100/20480.0488Hz。然后分别在呼吸带通频段和心跳带通频段内寻找幅度最大的谱线对应的频率乘以60就是每分钟的次数。不过FFT峰值搜索容易受频谱泄漏影响。如果信号频率不是恰好落在FFT的频点上峰值幅度会被稀释到相邻谱线。工程上可以用抛物线插值或Chirp-Z变换提高频率估计精度。更笨但有效的方法是把FFT长度凑到2的幂同时选择帧长使得目标频率尽量落在频点上。对实时系统来说每次处理一帧数据要等20秒才能输出一次结果太慢了。我用的滑动窗口方案是每2秒滑动一次窗口长度20秒保证频率分辨率的同时兼顾实时性。这样刷新率约0.5Hz观察实时波形足够用。3.5 BNR短时呼吸噪声比的计算与作用BNR短时呼吸噪声比是我在调算法时发现很实用的一个质量指标。它衡量的是当前数据段中呼吸信号能量与噪声能量的比值。计算公式我采用BNR 10 * log10( \frac{E_{breath}}{E_{noise}} )其中E_breath是呼吸频带0.1-0.5Hz内的频谱能量E_noise是0.5Hz到1.5Hz频带内的频谱能量。严格来说这个噪声频带里可能包含部分心跳能量但从能量规模上看呼吸信号通常远超心跳所以BNR主要反映的是信号质量而非心跳污染。我把BNR的阈值设为15dB。如果BNR低于15dB说明当前数据段的呼吸信号不可靠可能是人体移动了、目标距离太远或者有突发干扰。在这种情况下我会丢弃这一段的呼吸率和心率估计值而不是强行输出一个可疑的结果。这个机制大大提高了系统的鲁棒性避免了屏幕上心率数值乱蹦的尴尬。4. 第一版实例用C跑通一个最小可用的呼吸心率检测4.1 核心类的骨架设计我第一版算法用C实现代码结构分为采集、预处理、解调、滤波、估计五个模块。下面是核心类的骨架struct RadarConfig { float sampleRate 100.0f; int frameSize 2048; float breathLow 0.1f; float breathHigh 0.5f; float heartLow 0.8f; float heartHigh 2.5f; float bnThreshold 15.0f; }; class VitalSignProcessor { public: void process(const float* iData, const float* qData, int length); float getBreathRate() const { return breathRate_; } float getHeartRate() const { return heartRate_; } float getBnr() const { return bnr_; } private: void removeDcOffset(float* phase); void unwrapPhase(float* phase, int n); void bandpassFilter(float* data, int n, float low, float high); float estimateRate(const float* data, int n, float low, float high); void windowedDft(const float* data, int n, int start, int end); RadarConfig config_; float breathRate_; float heartRate_; float bnr_; };process函数按顺序调用私有方法。每次传入100Hz采样的一批新数据内部维护一个长度为2048的环形缓冲区避免频繁拷贝。在实时运行时这个环形缓冲区的实现非常关键建议用幂等容量的数组加写指针和读指针管理。4.2 关键实现片段解缠、滤波、FFT解缠和滤波的代码片段前面已经给出。FFT部分我用的KissFFT体量小适合嵌入式。初始化一次后续每帧执行一次1024点的复数FFT。注意KissFFT的输入是复数数组需要把滤波后的实数数据放进real部分、imag置零。峰值搜索函数要注意边界处理呼吸频带内0.1Hz可能落在FFT的第2或第3个频点如果直接找全频带最大值可能会把直流泄漏算进去。我在搜索前先排除DC和低于0.05Hz的频点。心率频带搜索同理排除掉低于0.8Hz的所有频点防止呼吸谐波或直流漂移混进来。滤波器的实现可以通过预生成的系数。我在MATLAB里用butter函数设计四阶Butterworth带通滤波器然后导出系数到C数组再用Direct Form II Transposed结构实现。这个结构数值稳定性好适合定点或浮点MCU。void VitalSignProcessor::bandpassFilter(float* data, int n, float low, float high) { // 这里用的是预先生成的二阶节系数 // 以呼吸带通为例一组二阶节的差分方程为 // y[n] b0*x[n] b1*x[n-1] b2*x[n-2] - a1*y[n-1] - a2*y[n-2] // 具体系数在初始化时从配置表中读取 float x10, x20, y10, y20; for (int i0; in; i) { float xn data[i]; float yn b0*xn b1*x1 b2*x2 - a1*y1 - a2*y2; x2 x1; x1 xn; y2 y1; y1 yn; data[i] yn; } }4.3 参数标定与验证方法程序写好以后先别急着拿真人测试。我用一个函数信号发生器叠加两路正弦波一路0.25Hz模拟呼吸一路1.2Hz模拟心跳再加上少量白噪声输入到算法里结果频率估计都非常准确。这个方法可以快速验证FFT和峰值搜索逻辑。随后我用一块小型扬声器纸盆当作模拟胸腔把纸盆贴在雷达前方用信号发生器驱动扬声器发出周期性的振动频率设为0.3Hz和1.5Hz的组合。这时候雷达面对的是真实的多普勒回波可以验证整个射频前端和模拟链路。虽然扬声器纸盆的振动幅度不能完全模拟人体但足以确认相位解调链路是否正常。最后才是人测。人测时要坐在雷达视线前方身体尽量不动。我拿了一台指夹式血氧仪同步记录心率与雷达估算值对比。大约测了10组数据误差在±3次/分以内呼吸率用人工秒表计时对比误差在±1次/分以内。对于第一版原型这个精度可以接受。5. 实测中躲不开的坑与我的处理方式5.1 静止目标为什么还会出现频率漂移我一开始以为只要人坐着不动呼吸和心跳频率就是稳定的但实测时发现心率估计值会时不时上下波动。后来把原始相位波形拉出来看发现信号里存在缓慢的基线漂移。原因是人体不可能完全静止即使呼吸平稳躯干也会有无意识的微小移动频率低于0.1Hz能量却非常大。这个超低频运动在心率带通滤波器里可能表现为零频附近的频谱泄漏影响峰值搜索。处理办法是在带通滤波前先做一阶高通滤波截止频率0.05Hz把低频漂移压掉。如果漂移很剧烈我还会用滑动窗口内的线性拟合并减掉趋势项。实测效果很明显心率估计的方差降了一半以上。5.2 谐波污染呼吸的三次谐波长得就像心跳呼吸波形往往不是纯正弦而是接近带尖峰的形状所以谐波能量很强。呼吸基频0.25Hz三次谐波0.75Hz刚好落在心跳频带下沿如果呼吸基频是0.3Hz三次谐波0.9Hz就完全落在心率频带内了。这时候FFT峰值搜索很容易把呼吸谐波当成心率导致心率翻倍或跳到异常值。我采用的解决办法是在心率频带峰值搜索时除了看绝对峰值还要检查这个峰值频率是否与呼吸基频存在整数倍关系。如果接近呼吸基频的2倍或3倍就把它标记为谐波候选继续搜索第二峰值。这个方法不完美但简单有效。更高级的做法是用自适应噪声对消从原始相位信号中减去重建的呼吸谐波我计划在系列第二篇里详细展开。5.3 运动干扰与多径的鉴别用CW多普勒雷达做生命体征检测最难受的问题就是没有距离分辨力人稍微一动I/Q信号的幅度和相位都会剧烈变化。比如挥手、转身这类大幅度运动相位信号会直接饱和。对付这种干扰我用了一个朴素的策略计算每个滑窗内相位信号的方差和峰峰值如果超过阈值就认为这一帧数据不可靠直接不输出结果等待运动结束。同时结合BNR丢弃低质量数据段。多径反射会造成信号抵消现象在某些距离上人体反射和墙面反射的电磁波相位相反I/Q输出幅度接近零微动信号完全被淹没。这个现象在不同位置反复出现雷达放得越靠近墙角越严重。解决方法是调整雷达位置尽量远离反射面或者用吸波材料/纸壳遮挡非目标方向的反射。如果条件允许换成FMCW雷达可以彻底解决因为FMCW能选择胸腔对应的距离门。5.4 下一步的改进方向第一版跑通以后我觉得最值得投入的方向有三个。第一个是心跳信号的信噪比提升。目前心脏微动信号幅度非常小在部分人群体质上甚至低于ADC噪声。可以考虑改到60GHz频段利用更短的波长获得更大的相位变化或者在接收端加入低噪声放大器和窄带抗混叠滤波器。第二个是呼吸谐波的自适应抵消。前面提到用整数倍关系剔除是权宜之计。更可靠的方案是通过呼吸基频的相位信息重建谐波模型然后从原始信号中减掉。这样对心率频带的保护会更彻底也能在运动伪影中保留更多有效信息。第三个是利用机器学习做运动状态分类。雷达信号中静止、微动、大幅度运动在时频域的特征差异很明显。可以先用传统方法提取BNR、谱熵、峰值形态等特征再用简单的分类器判断当前是否适合估算生命体征。这样系统会更智能不是直接丢弃数据而是选择性地信任不同信号分量。做这个项目最大的体会是雷达生命体征检测本质上是一场信号处理和微弱信号检测的硬仗。硬件选型只是起点后面的每一行算法都是在和噪声、谐波、运动伪影作斗争。我现在回头看第一版代码虽然粗糙但正是通过一遍遍踩坑才真正理解了相位解调为什么必须做、呼吸谐波为什么那么烦人。这个系列后续我会放出更完整的实测数据和优化后的算法框架也欢迎正在做类似项目的朋友一起交流。
返回列表