利用热成像技术对心率进行无接触检测的研究

2018-02-04 08:26梁智敏肖书明甄庆凯
中国体育科技 2018年1期
关键词:灰度感兴趣均值

梁智敏,陈 骐,肖书明,马 洁,甄庆凯



利用热成像技术对心率进行无接触检测的研究

梁智敏1,陈 骐2,肖书明2,马 洁1,甄庆凯2

1.北京信息科技大学自动化学院,北京100192;2.国家体育总局体育科学研究所,北京100061

目前,在体育领域检测运动员心率的常用方法大多数是接触式检测,而接触式检测会使运动员产生异物感,影响运动员的真实运动状态。利用红外热成像技术,提出了一种非接触式实时测量心率的新方法。对于静态心率,采用图像处理技术对视频帧中感兴趣区域进行分割并追踪,得到视频中感兴趣区域的整体灰度均值变化曲线,对采集到的信号采取小波分析去除信号中的噪声,并进行带通滤波,最后利用小波包重构实现心率的实时检测,取得了较好的效果。对于动态心率变化,采用热红外图像中感兴趣区域的灰度均值的变化来分析动态心率的变化趋势,这两种方法在Bland-Altman分析中取得良好的一致性。

热红外图像;心率测量;滤波;小波包分解与重构

1 引言

运动训练的生理学实质,是通过所施加的训练负荷使运动员身体各器官系统产生适应性变化,从而提高其参与运动竞技的能力。在运动员训练过程中,需要定期或不定期地检测其生理指标,根据运动员的相关指标,调控训练负荷,以改善运动训练的效果。心率作为度量训练负荷的一个重要生理指标,经常用于运动员身体机能的检测和评定。目前,运动员在运动过程中心率的获取一般都是采用心率带、手环等接触式的方式来实现,这种方式可以实现在运动状态中的心率实时采集,但在检测心率的过程中都会与被检测者产生一定的物理接触,且长时间的测量会给被测试者带来不适或异物感,影响运动员的正常运动状态。

基于图像的非接触式检测方法包括可见光成像和热红外成像。与可见光成像相比,热红外成像具有以下几点优势:1)可以发现运动员的潜在身体损伤;2)视频采集过程中受光照变化影响小;3)可以实时、全面的获得人体表面的温度信息。利用热红外图像进行心率检测最早是由美国休斯敦大学的教授Marc Garbey等[12]在2007年提出的,他采用自适应滤波和快速傅里叶变换的方法来实现人体心率的测量,2010年,我国首都医科大学生物医学工程学院的景斌等[2]利用热红外图像并采用功率谱分析的方法也成功的检测出了人体心率,2013年,美国路易斯维尔大学的Travis R.Gault等[11]利用热成像中温度的变化,采用两种不同的方法测试了人在3种不同的状态下的心率,并对两种方法的精确度进行了对比。

由以上文献可见,利用热成像中人体浅表动脉的温度变化来检测人体心率的方法是可行的。本研究利用热红外成像这种非接触式的检测技术并采用新的算法进行心率检测,非接触式的测量方法可以从根本上消除异物感和不适感,新的算法则提高了心率检测的精确度,以及获得更真实有效的心率数据。此外,此类研究获得的结论和方法,对于检测更多其他人体生理参数具有潜在的技术价值。

在基于图像检测心率的数据处理算法中,常用的主要有功率谱分析,傅里叶变换[16],人工智能[1],独立成分分析[6]等,经典的傅里叶变换更适用于平稳随机信号,其特性不随时间变化的信号[8],因此,不适用于动态心率信号的处理。本研究中需要处理的信号是非平稳的随机信号,且信号中含有各种噪声的干扰,所述噪声一般是宽带随机信号,其频带和有用信号的频带存在重叠。目前,在实现非平稳随机信号的消噪中,小波分析被认为是最好的方法[4]。因此,本研究利用小波进行信号的去噪处理,并运用小波包分解与重构对信号进行分析。

2 研究对象与方法

2.1 研究对象

静态心率检测实验中研究对象8人,其中4男、4女,动态心率检测实验中研究对象20人,其中10男、10女,均为在校学生,年龄范围在22~25岁。

2.2 研究方法

本研究所用热红外图像的采集设备为上海斌瑞检测技术服务有限公司的IRS S65-A,图像的分辨率为384×288像素,红外波长范围为8~14,测温范围在-20°C~ +350°C,热灵敏度<0.05°C,测温精度在±2°C或±2%,取其大值,由血液流经面部造成的温度变化范围在0.1°C~0.5°C(此范围由视频中的温度变化来获取),在检测精度范围内,满足目标的要求。实验在室内进行,在静止状态下,采集包含受试者面部的视频,每人时长2 min,视频格式为AVI,静态实验时的环境温度为22.3°C。在运动状态下,同样采集包含受试者面部的视频,每人时长16 min,视频格式也为AVI,动态实验时环境温度为23.7°C,运动的过程如下:开始时,保持静止状态1 min,然后在跑台上跑动。为了对不同心率状态进行测试,逐步提升运动负荷,跑台速度按5个台阶逐步递增,分别为:3 km/h、4 km/h、6 km/h、8 km/h、10 km/h,每个速度台阶持续3 min。静态和动态试验均采用德国运动心肺功能测试系统Contex Metalyzer 3B(CPET)同步采集受试者的实时心率数据。

采集实验数据时,如何使热红外图像与心率采集设备实现数据采集的时间同步是一个非常重要的问题。本实验时间同步的方法是为热成像视频设备和心率采集设备设置一个相同的时间起点基准,利用能够产生与环境温度差别较大的物体(例如一个发热设备)作为发令装置,使热成像视频采集设备和心率采集设备在发令装置指挥下,同时启动采集工作。为了进一步提高同步基准点的一致性,再在两种设备采集到的数据启动点之后,延迟2 s,截取数据,再按照3.3.1中所述的方式来同步两种设备采集到的信号。

3 原理及实现

3.1 物理原理

众所周知,自然界中一切温度高于绝对零度(-273 ℃)的物体均不断向空间发射红外辐射,且红外辐射与物体温度存在一定关系。当物体表面具有相同的辐射率时,物体的温度越高,其单位表面积的红外辐射功率越强。红外热像仪的测温原理,即为通过对红外辐射功率的度量来获得物体表面温度的分布情况。

关于人体的皮肤温度和心脏,通常认为存在如下关系:皮肤温度受皮肤血管内的血流量的影响[7],而皮肤血管内的血流量又与心脏活动相关:心脏的各个心房和心室以一定的频率收缩和舒张,把血液压入动脉血管,以及促使血液从静脉血管流回心脏。血液作为一种传热介质,会将人体内核部位的热量传递到人的身体外层(包括皮肤)。理论上来说,血液带来的热量,会导致身体各处皮肤的温度发生变化,1999年,中国中医研究院针灸研究院的张栋等[5]利用红外热像的方法总结了面部皮肤温度与面部血流量的关系,得到结论:在严格控制其他相关条件的情况下,温度与血流量之间有着基本一致的对应关系。因此判断,皮肤温度的变化与心脏的跳动之间存在潜在的相关关系。本研究依据上述的关系,通过对人体皮肤的热红外图像特征进行分析,来研究皮肤的温度数据特征与心率之间的关系。

3.2 静态心率的算法

3.2.1 心率检测算法的整体设计

本实验在采集人体脸部热红外图像的同时,同步采集人体心率数据。首先,在热红外图像中寻找潜在地与人体心率相关性较高的区域作为本次实验的感兴趣区域,然后对采集到的热红外图像中的感兴趣区域进行追踪。由于本实验采集的是人体处于安静状态下的视频,因此,经典的Camshift算法[9]能够满足对感兴趣区域的追踪要求,对追踪到的每个感兴趣区域进行运算,得到感兴趣区域(Region of Interesting,ROI)的灰度均值,然后对该灰度均值进行小波去噪以及带通滤波,最后对得到的滤波结果进行小波包分解和重构,得到对人体心率的估计数据。

3.2.2 视频追踪和数据获取

根据热红外图像中人脸区域亮度很高的特点,对图像进行二值化、腐蚀和图像填充。然后进行人工选择,仅保留人脸和裸露的颈部区域(图1)。

图1 二值化、腐蚀、填充操作后取得的人脸区域示意图

Figure1. Face Area Obtained after Two-valued, Corroded and Filled Operation

选择了人脸和颈部区域之后,根据人脸的宽高比确定人脸的边界(图2)。考虑到心脏跳动与动脉血管中的血流量存在的直接关系,以及动脉血管主要分布于口角两侧至鼻根区的三角区域的事实(图3)。将热红外图像中脸部的下1/3区域设定为所要跟踪的感兴趣区域(图4),由此区域推算得到的心率数据,准确性要高于由脸部其他区域。

图2 人脸区域标记示意图

Figure2. Human Face Area Mark

为了演示本处理方法,以下给出针对一个人脸的热红外视频的处理例子。在该视频中,对感兴趣区域的灰度作均值处理,得到该区域的灰度变化的时间序列。感兴趣区域的灰度均值随时间的变化如图5所示。由彩色图像转化为灰度图像的计算公式如下:

=0.2989×+0.5870×+0.1140×(1)

图3 脸部动脉血管分布图[13]

Figure3. Vascular distribution of facial arteries

图4 选定的脸部感兴趣区域示意图(ROI)

Figure4. Selected Face Area of Interest

图5 时间序列中灰度均值的变化图

Figure5. Change of Gray Mean Value in Time Series

感兴趣区域的灰度均值的计算公式如下:

为了考察感兴趣区域的灰度均值随时间变化的情况,对采集到的原始时间序列进行一阶差分,得到感兴趣区域灰度均值的变化曲线。一阶差分的计算公式如下:

C=(i)-(i-1) (i≥2) (3)

3.2.3 小波去噪

在获取热红外图像的过程中,由于受到检测对象、工作环境、数据采集硬件的影响,图像数据不可避免地存在噪声。小波去噪可以在消除噪声的同时,很好地保留有用信号。小波去噪效果的好坏,与小波基函数和阈值的选取有关。表1所示为几种常用小波基函数的性质对比。

表1 常用的小波基函数及其性质

本研究采集的数据是离散数据。由表1可以排除Morlet、Mexician hat这些无法进行离散小波变换的基函数。由于心率信号受噪声干扰严重,为了能够有效分离信号与噪声,应选取具有一定阶次消失矩的基函数,因此,排除只有一阶消失矩的Haar小波。基于以上分析得出:db小波、sym 小波、coif 小波是适合本研究的小波基函数。从小波函数的分解来看,同一小波族的小波函数其序数越高,分解效果越好,最终选择sym8作为小波基函数。常用的阈值估计方法有4种:固定阈值、自适应阈值、启发式阈值、极大极小阈值。其中,极大极小阈值和自适应阈值规则比较保守,去噪不完全。计算实验的结果说明,固定阈值法的效果较好,因此,本研究选用了固定阈值法。小波分解的层次越高去掉的低频成分越多,去噪效果虽然越显著,但失真度也相应增大。经过多次试验,确定分解层数为2层。为了演示去噪效果,针对一个实际采集的热红外视频的灰度均值利用式(3)所示公式作一阶差分后进行小波去噪(图6)。

图6 小波去噪的计算结果示意图

Figure6. The Computational Results of Wavelet Denoising

根据人的正常心率范围在0.8~3 Hz之间这个事实,认为该频带之外的信号成分为噪声,因此,以该频带为通带,对小波去噪后的信号进行带通滤波,在保留心率信号的同时,消除信号中的噪声成分。针对上述的小波去噪后的信号,进行带通滤波处理(图7)。

图7 带通滤波后的信号示意图

Figure 7. Signal after Bandpass Filtering

3.2.4 小波包分解与重构

小波包分析能够为信号提供一种更精细的分析方法,它将频带进行多层次的划分,对多分辨率分析没有细分的高频部分进行进一步分解,并能根据被分析信号的特征,自适应地选择相应的频带,使之与信号频谱匹配,从而提高了时频分辨率。此处结合图8来说明小波包分解的基本思想。在图8中,S代表被处理的信号,A表示低频,D表示高频,末尾的序号表示分解层数。经小波包分解后每个节点中保存有相应频段的系数,最后利用节点的系数对相应频带的信号进行重构。

图8 小波包分解树示意图

Figure8. Wavelet Packet Decomposition Tree

3.3 动态心率趋势分析

3.3.1 数据获取及处理

要进行动态心率趋势分析,必须采集人体运动状态下的热红外图像,并设定需要追踪的感兴趣区域。热红外图像中脸部的下1/3区域设定为所要跟踪的感兴趣区域(图9)。此处所采用的感兴趣区域跟踪算法为TLD(Tracking-Learning-Detection)跟踪算法[13]。

图9 动态数据采集过程中感兴趣区域示意图

Figure9. The Area of Interest during Dynamic Data Collection

在实验过程中,受试者在跑台上跑动,红外热像仪集中拍摄人脸部视频,同时采用运动心肺功能测试仪同步采集运动心率。运动心肺功能测试仪每间隔5 s采集1次心率数据,而热红外图像的帧率为每秒7帧,为了使二者同步,对采集的热红外视频进行重采样(每35帧采样一次),然后计算重采样的热红外图像中的感兴趣区域的灰度均值。

4 结果及分析

4.1 静态心率检测结果

4.1.1 小波包重构结果

针对上述带通滤波后的信号(图7)进行小波包的分解与重构,分解层数为2层,小波基函数选择db8,熵的类型选择shannon,重构后第2层每个节点的信号分量如图10所示。分析计算结果,认为只有节点0的信号分量与心跳有较明显的相关,而其他的3个节点的信号分量均为噪声。并且,认为节点0的信号分量中的每一个局部峰对应于一次心跳。

图10 小波重构后第二层各节点的信号示意图

Figure10. The Signal of the Second Layer after the Wavelet Packet Reconstructed

4.1.2 心率值的计算方法及结果

上述心率信号的提取过程是对人脸图像中的感兴趣区域进行处理。为了得到实时性的心率信号估计,需要对每一次得到的感兴趣区域的相应数值进行处理,利用公式(7)估计受试者每分钟的心率:

=**60/(7)

其中表示当前的帧数,表示经小波包重构后信号的波峰数,表示视频的采样频率。

图11 选择不同感兴趣区域时的结果对比图

Figure11. Comparison on the Results of Different Areas of Interest

根据以上算法,得到了上述一个热红外视频中,选定人脸的不同区域作为感兴趣区域时的心率信号估计(图11)。在图11中,蓝色曲线为以人脸的下部1/3区域为感兴趣区域得到的心率估计信号,绿色曲线为以额头为感兴趣区域得到的心率估计信号。在本例中,实验者1的平均心率为58,由图可以看出,蓝色曲线的误差比绿色曲线的误差小。因此,图11的结果证明,以人脸的下部1/3区域(面动脉即分布在此区域)作为感兴趣区域,能够得到更准确的心率估计,也验证了3.2.2节的有关假定。

对于8位静态受试者,心率实时结果如图12所示,蓝色曲线为本算法得出的实时心率值,红色曲线表示运动心肺功能测试仪实时采集到的心率均值:

Figure12. Real-time Calculation Results of Heart Rate

采集8名静止仰卧状态下的受试者的面部热红外视频,设定感兴趣区域并追踪,并根据前述算法计算出安静心率;以本研究方法和以运动心肺功能测试仪得到的心率数据,分别计算它们在测试时间内的均值,并将得到的均值进行比较(表2)。

表2 本研究算法的心率计算结果和运动心肺功能测试仪实测的心率数据

4.2 动态心率趋势分析结果

首先分析其中10位男性受试者的动态心率变化趋势,由于实验者5的心率数据丢失,因此,在图13中列出了9位受试者的心率变化趋势。从图13中可以看出,男性在运动过程中面部皮肤的温度变化较大,9名实验者在运动负荷不断增大的情况下,灰度均值也随之不断增大,在心率小于130次/min时,根据热红外图像计算得到的灰度均值曲线中的台阶变化,与运动心肺功能测试仪实测心率的台阶变化基本趋势一致。但当心率超过130次/min时,计算心率曲线不再随着实测心率的上升而上升,反而出现下降趋势。这可能是由于人体在运动过程中温度不断升高,人体表面的发射率(又称为黑度,反映物体向外发射辐射的能力)会随温度的变化而变化,此时热成像设备会根据现有的温度输出信号标定曲线作为计算温度的原始数据,从而对输出视频的灰度进行调节修正,因此导致灰度均值在心率大于130次/min以后逐渐下降。

Figure13. Comparison of the Gray Mean of Thermal Imaging with the Real Heart Rate in 10 Male Subjects

Figure14. Contrast between the Gray Mean of Thermal Imaging and the Change of Real Heart Rate in Female Subjects

对10名女性受试者进行分析,由于女性实验者10采集的是侧面的热红外视频,无法提取本研究中感兴趣区域的灰度均值,因此不做分析。对于剩余的9名女性受试者动态心率变化趋势做出分析(图14)。可以看出,女性在运动过程当中脸部的温度变化不大,实验者1、3、6、9同样出现了运动负荷增大,灰度均值反而下降的趋势,并且男性和女性随运动负荷增加,灰度曲线变化的趋势也明显不同:女性受试者在增加运动负荷时,灰度均值曲线都没有太大的变化,也没有明显的随心率不断变化的趋势,这一点与男性灰度均值曲线变化有明显的差别。

图13中男性实验者10的两组数据中心率变化的趋势与灰度均值曲线的变化趋势基本一致,只在心率大于130次/min时有差别,因此,进一步对这两组数据进行统计学分析,其相关性达到了0.857 62,根据图13中男性实验者10的数据,对两种方法的一致性进行Bland-Altman分析,相关结果如图15所示。Bland-Altman分析的基本思想是计算出两种测量结果的一致性界限,并用图形的方法直观地反映这个一致性界限。在图15中,上下两条水平虚线代表95%一致性界限的上下限,中间的虚线代表差值的平均数,两种测量方法的一致性程度越高,代表差值均数的虚线越接近均数0。图中所示两种测量方法的结果只有3个位于95%的界限以外,其余结果均位于界限内。

图15 Bland-Altman分析示意图

Figure15. Bland-Altman Analysis

5 总结

由静态实验的结果可以看出,本研究采用的小波包分解与重构的算法具有较高的时-频分辨率,可以计算出人体静态的实时心率,且具有较好的效果。在可见光领域,普通的成像受光照变化影响很大,而红外热成像受光照变化影响小,因而有更广泛的适用性[14]。本研究采用的相关方法,实现了对心率的非接触实时检测,并且相关系数达到了0.873。在信号采集的过程中,由于人的面部血管的分布具有唯一性[15],而在本实验中对所有对象采用的是同一个感兴趣区域,这可能是造成测量误差的一个原因。也有其他相关研究利用红外序列对心率进行了检测,虽然相关性达到了0.957[2],但其要求受试者保持静止状态,对实验对象的约束较多,且没有实现心率检测的实时性。另外,受目前热成像硬件设备限制,仅能够使用已经拍摄好的视频进行追踪和计算,尚无法实现利用佩戴式心率表、可穿戴设备,以及基于普通CCD的实时心率数据反馈[3]。

由动态的实验结果可以看出,利用动态热红外图像的灰度均值与真实心率的对比,热红外图像的灰度均值能很好地反映出运动过程心率的变化趋势,另外,对于男性实验者10运用这两种方法得到的数据进行Bland-Altman分析的实验结果均落在了95%的置信区间内,因此可以认为,在评价心率变化趋势时,利用热红外成像灰度均值和运用运动心肺功能测试仪这两种方法具有较好的一致性。

另外,本研究提出的小波包分解与重构算法同样适用于动态心率的测量,原因是动态心率信号为非平稳随机信号,并且频率是随时间变化而逐渐增大的。小波包分析有很好的时频特性,在分析动态心率时可以分析随时间变化的心跳频率变化,但计算心率值的算法与静态心率的算法却有所不同,这是因为在运动过程中,心率的频率是不断变化的。因此,不可以运用本文4.1.2中计算波峰数的方法来计算人体运动状态的实时心率,具体的动态心率算法还有待进一步研究。

目前,利用红外热成像视频来检测动态心率为以非接触式的方法检测心率提供了新的思路和数据支持,也是以非接触的方式检测心率的最终目标,并且在众多的实验中,有一位实验者的温度变化和心率变化具有很好的一致性,后续研究应考虑各相关因素对皮肤温度变化的影响,建立皮肤温度变化模型,利用抗运动干扰的算法对心率信号进行提取,这也是本研究团队下一步的工作方向。

[1] 蔡成贤,王伟.基于人工智能的心率检测算法[J].中国医疗器械杂志,2010,34(1):13.

[2] 景斌,李海云.基于红外序列图像的心率无损检测方法研究[J].中国生物医学工程学报,2010,(6):943-946.

[3] 刘蕾.基于视频的实时心率测量算法研究[D].无锡:江南大学, 2014.

[4] 张吉先,钟秋海,戴亚平.小波门限消噪法应用中分解层数及阈值的确定[J].中国电机工程学报,2004,24(2):118-122.

[5] 张栋, 薛立功, 魏正岫, 等.面部皮肤温度与面部血流量关系的对照观察[J].生物医学工程学杂志,1999,(1):81-85.

[6] 周秦武,隋芳芳,白平,等.嵌入式无接触视频心率检测方法[J].西安交通大学学报,2013,12:55-60.

[7] 朱大年 ,王庭槐. 生理学[M].北京:人民卫生出版社, 2013:85-95,255-260.

[8] CHEKMENEV S Y, FARAG AA, MILLER W M,. Multiresolution approach for noncontact measurements of arterial pulse using thermal imaging[M]//augmented vision perception in Infared,2009:87-112.

[9] DONGARRA J. THe top 10 algorithms[J]. IEEE Comput Sci Eng,2000,2(1):22-23.

[10] DORLING KINDERSLEY. The Concise Human Body Book[M]. Dorling Kindersley, 2009:146-147.http://joyousoasis.com/wp-content/uploads/2013/02/CardiovascularSystem.gif.

[11] GAULT T R, FARAG A A.A fully automatic method to extract the heart-rate from thermal video[J]. IEEE Conf Comput Vis Pattern Recognit Workshops,2013,13(4):336-341.

[12] GARBEY M, SUN N,MERLA A,. Contact-free measurement of cardia-c pulse based on the analysis of thermal imagery[J]. IEEE Trans Biomed Engi,2007,54(8):1418-26.

[13] KALAL Z, MIKOLAJCZYK K, MATAS J. Tracking-learning-detection[J]. IEEE Trans Pattern Anal Mach Intell, 2012,34(7):1409-1422.

[14] SOCOLINSKY D A, WOLFF L B, NEUHEISEL J D. Illumination invariant face recognition using thermal infrared imagery [C]// In: Proceedings of the 2001 IEEE Computer Society Conference on Computer Vision and Pattern Recognition. Kauai:IEEE,2001:527-534.

[15] YANGJ, LIUL, JIANGT. A modified gabor filter design method for fingerprint image enhancement[J]. Pattern Recogn Lett2003, 24(12):1805-1817.

[16] ZENG W, ZHANG Q. Infrared video based non-invasive heart rate measurement[C]//IEEE Conference on Robotics and Biomimetics, Zhuhai, China, 2015:6-9.

Research on Non-contact Heart Rate Detection Using Thermal Imaging

LIANG Zhi-min1,CHEN Qi2,XIAO Shu-ming2,MA Jie1,Zhen Qing-kai2

1.Beijing Information Science & Technology University, Beijing 100192, China; 2. China Institute of Sport Science, Beijing 100061, China.

At present, the commonly used method to detect heart rate in the field of sports is contact detection, contact detection will affect the status of athlete. Therefore, a new method for non-contact measuring changes in heart rate in real time with infrared thermal imaging technology is proposed. For static heart rate, we using image processing technology to segmentation and tracking for region of interest in video frames, and get overall gray-scale average curve in video region of interest, and using wavelet to eliminate noise in the collected signals and band-pass filters, then wavelet packets reconstruction to achieve real-time heart rate detection, and achieved very good results. For dynamic changes in heart rate, we using gray-scale of ROI in thermal video to analyze dynamic change trend of heart rate, the two methods get good consistecy in the analysis of Bland-Altman.

G804.49

1002-9826(2018)01-0136-10

10.16470/j.csst.201801019

2016-09-22;

2017-10-12

国家体育总局2014年度重点研究领域课题(2014B028);国家体育总局体育科学研究所基本科研业务费项目(基本1434)。

梁智敏,女,在读硕士研究生,主要研究方向为控制工程-图像处理,Email:809916630@qq.com。

陈骐,男,副研究员,主要研究方向为体育工程学,E-mail: chenqi@ciss.cn。

猜你喜欢
灰度感兴趣均值
航空滤光片阵列多光谱图像条带灰度调整算法
采用改进导重法的拓扑结构灰度单元过滤技术
对自己感兴趣
天津港智慧工作平台灰度发布系统和流程设计
Arduino小车巡线程序的灰度阈值优化方案
均值—方差分析及CAPM模型的运用
均值—方差分析及CAPM模型的运用
均值不等式的小应用
编读往来
现在是几点