CN1142511C - 基于自回归滑动平均谱分析的医用红外热图分析方法 - Google Patents

Abstract

本发明属于医用红外热图计算机辅助诊断领域，本发明基于自回归滑动平均谱分析通过对温度变化过程的多幅热图进行现代谱分析，将病变区域对温度变化的响应过程用谱参数表现为谱曲线的形式，将病变区域的温度响应特征以一系列子波的波峰和波谷定量地表现出来。这种特征曲线的表示方法更深刻地揭示了疾病的内在固有特征，并且以更直观的方式显示，有助于对不同疾病特征参数的确定和区分，且减少了诊断的模糊性，提高了诊断的可靠性。

Description

基于自回归滑动平均谱分析的医用红外热图分析方法

技术领域  本发明属于医用红外热图计算机辅助诊断技术领域，特别涉及将自回归滑动平均谱分析方法应用于医用红外热图的分析。

背景技术  自回归滑动平均模型是以传统的傅立叶变换和滤波理论为基础的新型谱分析方法。自现代自回归滑动平均谱产生的二十多年来，在为数众多的领域中都得到了很好的应用，如在工程、医学、气象学、农业学等领域。在医学领域中，应用自回归滑动平均谱对心电图、脑电图、肌电图进行信号分析；在信号传播领域中，用于对语音的识别以及语音的处理与合成；在地震预报方面用作地震信号的建模；在机械***的故障诊断中用来进行机械***的参数识别；在军事信号检测方面，自回归滑动平均谱估计方法也已取得了成熟的应用，如对雷达和声纳信号的检测等等。

自回归滑动平均模型是一种有效的时序分析模型。从谱分析的观点来讲，自回归滑动平均谱分析在记录的时间序列数据较短时，能够产生高分辨率的谱估计，方法本身是数据自适应的，其特例自回归模型和滑动平均模型在进行时序分析的***参数辨识的时候，各有特色。用自回归模型辨识出来的***，是一个全极点的***，在此基础上进行的功率谱估计的波峰较准；而滑动平均模型辨识出来的***是一个全零点的***，在此基础上估计出来的功率谱的波谷较准；自回归滑动平均模型则集中了自回归模型和滑动平均模型的长处，所估计出来的功率谱的各方面的性能都比较好。

时间序列是指按时间顺序排列的一组数据，主要是指一组有序的随机数据。它们是同产生它们的、可称为***的客观事物相关联的。因此，某一时间序列便反映着该***的动态变化过程，就是该***的有关输出或者响应，以及该***同外界相互联系的方式。

医用红外热图是人体体表温度场的反映，是利用红外摄像仪，根据普朗克原理： ${\mathrm{\ω}}_{\mathrm{\λ}}=\frac{C_1}{{\mathrm{\λ}}^5}\frac{1}{\exp (c_2/\mathrm{\λT})-1}$ (ω_λ——对应于波长λ的辐射能量，λ——辐射能波长，T——能量辐射体的绝对温度，C₁C₂——常数)从测量到的人体红外辐射能量计算出来的。根据医学物理学，人体绝大部分代谢都是放热的，并且大部分是通过皮肤辐射出去的，人体出现病变或异常时，首先会在人体体表温度场上有所体现，基于医用红外热图的医学诊断就是分析人体温度场变化的过程。

已有的医用红外热图分析***，主要由热图获取部分(红外摄像仪)和热图分析软件两部分组成。热图获取部分，大多采用的是光机扫描式的红外摄像仪，尽管扫描速度较慢，但其空间分辨率和温度分辨率都较高，已经能够满足医学诊断的要求。热图分析部分，所提供的分析方法，大多较为简单，仅仅局限于分析单幅红外热图，进行一些简单的温度统计(如指定区域的温度直方图)，总体上来说不能够完全分析出医用红外热图中所包含的病理信息，诊断的可靠性不高。之所以出现这样的问题是由于人体体表特征因人而异，差别较大。病灶在热图上的表现也存在有很大的差别，对于一些隐性的、潜在的疾病往往不能用单幅热图简单地分析其热图特征，因而不能对辅助诊断提供有力的分析结果和病变信息。

发明内容  本发明的目的是克服已有技术的不足之处，提供一种现代谱分析法在医用红外热图分析方面的应用，通过多幅红外热图来获取待测人体表面区域对温度变化的响应过程，然后对代表改响应过程的时域温度数据作现代谱参数的估计，从而将体表温度的变化，尤其是将动态过程体表温度的变化转化为谱特征参数，这些参数进而可被提供用于临床诊断的数据。

本发明提出一种基于自回归滑动平均谱分析的医用红外热图分析方法，包括以下步骤：

1)计算被输入的多幅病变区医用红外热图的温度变化时域数据，生成自回归(n)模型和滑动平均(m)模型的正则方程，其中n和m分别为自回归模型和滑动平均模型的阶次，所说的阶次n和m根据模型适用性检验准则来确定；

2)求解自回归(n)模型的正则方程及代表该病变区所测点的温度变化响应过程的极点参数，将所述参数代入滑动平均(m)模型的正则方程中，求解代表该病变区所测点的温度变化响应过程的零点参数，用该模型估计出来的零、极点参数即表征了所述病变区上的所述点对温度刺激的响应特征；

3)由自回归(n)模型和滑动平均(m)模型生成自回归滑动平均(n，m)模型的正则方程，从中计算代表该病变区所测点的温度变化响应过程的自回归滑动平均(n，m)功率谱曲线；

4)将计算出来的自回归(n)模型的参数同标准病理模板中的自回归(n)模型参数做互相关系数计算，将表征该温度变化过程相关的程度作为对该疾病发生的可能性的辅助诊断指标一；

5)将计算出来的滑动平均(m)模型的参数同标准病理模板中的滑动平均(m)模型参数做互相关系数计算，将表征该温度变化过程相关的程度作为对该疾病发生的可能性的辅助诊断指标二；

6)对指标一和指标二综合考虑从而为临床诊断提供有关断数据。

所说的自回归(n)模型的阶次、滑动平均(m)模型的阶次及自回归滑动平均(n，m)模型的阶次n和m均可为10。

上述的标准病理模板的自回归(n)模型参数和滑动平均(m)模型参数以及自回归滑动平均(n，m)功率谱曲线的计算方法与病变区的所述参数及所述功率谱曲线的计算方法相同。

本发明方法的理论基础在于：

如上所述，现代谱分析法尤其是自回归滑动平均谱分析法已在包括医学信息处理在内的诸多频域得到了应用，这种方法的优点在于能利用有限的时间序列数据来对分析对象建立模型，并且估计出模型的参数，利用参数来揭示被测对象的内在本质。因此，本发明将自回归滑动平均谱分析方法引入到医用红外热图的分析与诊断领域。

任何一部分人体都是一个***，有自己的***参数，发生相同病变的区域的***参数一定有其相似的成分。诊断的过程，就是分析待测区域病理参数，和标准病理模板中的标准病理***参数进行比较的过程。而对于***参数的计算，在工程上已经有了比较成熟的做法，可以对病变区加以外界刺激，分析其响应数据，利用已有的时序分析方法，估算出其***参数(一般为零/极点形的***参数)，并以现代自回归滑动平均谱的形式以曲线的方式显示其变化过程。

本发明所应用的自回归滑动平均谱相关的基础理论说明如下：线性回归模型：

一个随机变量y_t取决于多个随机变量x_1t，x_2t，x_3t...x_nt：

  y_t＝β₁x_1t+β₂x_2t+...β_nx_nt+ε_t   ε_t～NID(0，σ²)其中β_t是模型的参数，ε_t是模型的残差。利用最小二乘法估计β_t： $S=\underset{t=1}{\overset{N}{\mathrm{\Σ}}}{\widehat{\mathrm{\ϵ}}}_1^2$ $=\underset{t=1}{\overset{N}{\mathrm{\Σ}}}\[y_t-({\mathrm{\β}}_1^{^}{x}_{1t}+{\mathrm{\β}}_2^{^}{x}_{2t}+...{\mathrm{\β}}_n^{^}{x}_{\mathrm{nt}}){\]}^2$ 对β_t求导，令其为零，得： $\frac{\∂S}{\∂{\mathrm{\β}}_t}=0---t=\mathrm{0,1,2,3}...N.$ 求解方程组，得到***的线性回归模型： $\left[\begin{array}{c}{y}_1\\ y_2\\ M\\ y_n\end{array}\right]=\left[\begin{array}{c}{\mathrm{\β}}_1\\ {\mathrm{\β}}_2\\ M\\ {\mathrm{\β}}_n\end{array}\right]\left[\begin{array}{cccc}{x}_{11}& x_{21}& \mathrm{\Λ}& x_{n1}\\ x_{21}& x_{22}& \mathrm{\Λ}& x_{n2}\\ \mathrm{\Λ}& \mathrm{\Λ}& \mathrm{\Λ}& \mathrm{\Λ}\\ x_{1n}& x_{2n}& \mathrm{\Λ}& x_{\mathrm{nn}}\end{array}\right]+\left[\begin{array}{c}{\mathrm{\ϵ}}_1\\ {\mathrm{\ϵ}}_2\\ M\\ {\mathrm{\ϵ}}_n\end{array}\right]$

                     y＝βx+ε

Yt由两部分组成，一部分是确定性部分，完全取决于各自变量的线性组合

是Yt的数学期望；另一部分是随机部分，完全独立于各x_1t，由白噪声ε_t决定，是导致Yt称为随机变量的原因。自回归滑动平均(n，m)模型：

按照多元线性回归的思想，一个自回归滑动平均(n，m)模型可表示为：x_t＝₁x_t-1+₂x_t-2+...+_nx_t-n-θ₁α_t-1-θ₂α_t-2-...-θ_mα_t-m+α_t ${\mathrm{\α}}_t~\mathrm{NID}(0,{\mathrm{\σ}}_{\mathrm{\α}}^2)$ n，m分别表示自回归部分和滑动平均部分的阶次；_i(i＝1，2，...，n)，θ_i(i＝1，2，...，m)为各部分的模型参数；x_i(i＝t-1，t-2，...t-n)为采样的离散时间序列；α_t(i＝t，t-1，...t-m)为具有零均值与方差σ_α ²的白噪声。n，m的阶次是否适用于模型，需进行验证。具体检验过程如下：(1)检测α_t是否是白噪声，检测α_t是否与α_t-1，α_t-2，...无关。计算α_t的自相关系数： ${\mathrm{\ρ}}_{\mathrm{\α},n}=\frac{\underset{t=2}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_t{\mathrm{\α}}_{t-n}}{\underset{t=2}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_t^2}$ 若ρ_α，n-＞0，则认为阶次适用于模型。(2)检验α_t是否与x_t-2，x_t-3...无关。计算α_t与x_t的互相关系数： ${\mathrm{\ρ}}_{\mathrm{\αx},n+m}=\frac{\underset{t=3}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_t{x}_{t-3}}{\sqrt{\left(\underset{t=n+2}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_t^2\right)\left(\underset{t=n+2}{\overset{N}{\mathrm{\Σ}}}{x}_{t-n-1}^2\right)}}$ 若ρ_αx，n+m-＞0，则认为阶次适用。

当模型阶次确定后，自回归滑动平均(n，m)模型亦可表示为： $x\left(t\right)=-\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{a}_{k}x(t-k)+\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{b}_{k}u(t-k)$ 其中a_k＝-_k(k＝1，2，...n)，b_k＝-β_k(k＝1，2，...m)。这个模型由两部分组成，一个为n阶自回归模型 $x\left(t\right)=-\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{a}_{k}x(t-k),$ 另一个为m阶滑动平均模型 $x\left(t\right)=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}b\left(k\right)u(t-k).$ 自回归模型：n阶自回归模型的基础表达式为： $x\left(t\right)=-\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{a}_{k}x(t-k)+u\left(t\right),$ 将方程两边同乘以x(t+m)，求均值，得到： $r_x\left(m\right)=E\left\{x\right(t\left)x\right(t+m\left)\right\}=E\left\{\right[-\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{a}_{k}x(t+m-k)+u(t+m)\left]x\right(t\left)\right\}$ 即 $r_x\left(m\right)=-\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{a}_k{r}_t(t-k)+r_{\mathrm{xt}}\left(m\right)$

u(t)是方差为σ²的白噪声，设此***的单位抽样响应为h(t)，所以： $r_{\mathrm{xu}}\left(m\right)=E\left\{u\right(t+m\left)x\right(t\left)\right\}=E\left\{u\right(t+m\left)\underset{k=0}{\overset{\∞}{\mathrm{\Σ}}}h\right(k\left)u\right(t-k\left)\right\}={\mathrm{\σ}}^2\underset{k=0}{\overset{\∞}{\mathrm{\Σ}}}h\left(k\right)\mathrm{\δ}(m+k)={\mathrm{\σ}}^{2}h(-m)$ 根据Z变换定义， ${\underset{z->\∞}{\mathrm{Lim}}}^{H\left(z\right)=h\left(0\right)},z->\∞$ 的时候h(0)＝1，所以推导其矩阵方程为： $\left[\begin{array}{ccccc}{r}_x^{\·}\left(0\right)& r_x\left(1\right)& r_{xx}\left(2\right)& ...& r_x\left(n\right)\\ r_x\left(1\right)& r_x\left(0\right)& r_x\left(1\right)& ...& r_x(n-1)\\ r_x\left(2\right)& r_x\left(1\right)& r_x\left(0\right)& ...& r_x(n-2)\\ ...& ...& ...& ...& ...\\ r_x\left(n\right)& r_x(n-1)& r_x(n-2)& ...& r_x\left(0\right)\end{array}\right]\left[\begin{array}{c}1\\ a_1\\ a_2\\ ...\\ a_n\end{array}\right]=\begin{array}{c}\left[\begin{array}{c}{\mathrm{\σ}}^2\\ 0\\ 0\\ ...\\ 0\end{array}\right]\end{array}$ 计算此方程，便得到n阶自回归模型的模型参数。再通过下式计算时间序列的功率谱： $p_x\left(e^{\mathrm{jw}}\right)=\frac{{\mathrm{\σ}}^2}{{|1+\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\α}}_k{e}^{-\mathrm{jkw}}|}^2}$ 滑动平均模型：

m阶滑动平均模型的基本表达式： $x\left(t\right)=u\left(t\right)+\underset{k=1}{\overset{m}{\mathrm{\Σ}}}b\left(k\right)u(t-k)$ $H\left(z\right)=1+\underset{k=1}{\overset{m}{\mathrm{\Σ}}}b\left(k\right)z^{-k}$ $r_x\left(t\right)=E\left\{\right[u(t+n)+\underset{k=t}{\overset{m}{\mathrm{\Σ}}}b\left(k\right)u(t+n)\left]x\right(n\left)\right\}=\underset{k=0}{\overset{m}{\mathrm{\Σ}}}b\left(k\right)r_{\mathrm{xt}}(t-k)$

b(0)＝1。因为：

r_x(t-k)＝E{x(n)u(n+t-k)}＝σ²h(k-t)

h(i)＝b(i)，i＝0，1，2，...，m便可得到滑动平均模型的正则方程：

计算此方程，可得到m阶滑动平均模型的模型参数。并通过下式计算其功率谱： $p_x\left(e^{\mathrm{jw}}\right)={\mathrm{\σ}}^2\underset{k=-m}{\overset{m}{\mathrm{\Σ}}}{r}_x\left(k\right)e^{-\mathrm{jwk}}$

综合运用自回归模型和滑动平均模型的算法，得到自回归滑动平均(n，m)模型的正则方程为：通过求解上式方程，可得到模型的n阶自回归模型参数和m阶滑动平均模型参数，进而可计算出此序列的功率谱为： $P_x\left(e^{\mathrm{jw}}\right)=\frac{{\mathrm{\σ}}^2{|1+\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{b}_k{e}^{-\mathrm{jwk}}|}^2}{{|1+\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{a}_k{e}^{-\mathrm{jwk}}|}^2}$

本发明的突出贡献是：

首先将自回归滑动平均谱分析方法引入到医用红外热图的分析与诊断领域，提供一种现代谱分析法在医用红外热图分析方面的应用，通过多幅红外热图来获取待测人体表面区域对温度变化的响应过程，然后对代表改响应过程的时域温度数据作现代谱参数的估计，从而将体表温度的变化，尤其是将动态过程体表温度的变化转化为谱特征参数，这些参数进而可被提供用于临床诊断的数据。

已有的热图分析方法都是基于对单幅热图的分析，利用病变区域和正常区域的温度分布的不对称性及温度差来进行疾病诊断。这种单一的作为某种病症诊断的温度差数据由于病人生理条件的差异及病变发生部位的不同往往难于确定，尚无统一的标准，因此常借助于医生的临床经验，从而更多地引入了诊断的主观性和不精确性。

与已有技术相反，本发明的方法通过对温度变化过程的多幅热图进行现代谱分析，将被测体表区域对温度变化的响应过程用谱参数表现为谱曲线的形式，将病变区域的温度响应特征以一系列子波的波峰和波谷定量地表现出来。这种特征曲线的表示方法更全面、更深刻地揭示了疾病的内在固有特征，并且以更直观的方式显示，有助于对不同疾病特征参数的确定和区分，也有助于减少采用单一温度差标准进行诊断的模糊性，从而提高了诊断的可靠性。

本发明的特点是：

(1)采用自回归滑动平均模型进行***特征参数估计来获取人体体表的温度变化特征参数。

(2)利用有限的时间序列数据来对分析对象建立模型，并且估计出模型的参数，利用参数来揭示被测对象的内在本质。

(3)采用现代自回归滑动平均谱参数估计的方法来分析红外热图。可将体表温度参数的变化转化为功率谱特征曲线，从而可为诊断提供更为准确的病理数据。

(4)采用现代自回归滑动平均谱参数估计的方法可更充分地揭示病变区域对温度变化的响应过程及其内在的病理信息，从而有可能对某些隐性疾病的诊断提供有用的数据。

附图简要说明：

图1为本发明的一种实施例的热图自回归滑动平均谱估计流程图。

具体实施方式

以下结合实施例及附图来详细说明本发明的方法及特点。

在本发明所述的方法的一个实施例中，采用时间序列分析手段，引入外界刺激--温度刺激，并对温度刺激求取其响应，并对响应数据作现代功率谱参数的估计，从而将体表温度的变化，尤其是动态过程体表温度的变化转化为谱特征参数。这些特征参数进而可被用来与标准(健康人)的病理参数模板进行比较，从而可为医务工作者的临床诊断提供必要的病理分析数据。

为获取本发明所需的温度分析数据，采取以下的步骤来摄取人体红外热图：

(1)患者在室温为25摄氏度的房间内静坐5分钟。

(2)选择可疑的病变区域，用45摄氏度的热水袋热敷两分钟，使该区域局部升温。

(3)撤去热水袋，用医用红外热像仪(型号：DH98)连续拍摄64张热图，记录此可疑病变区域的温度衰减过程。

(4)自动采集64张图中可疑病变区的同一点，形成分析需要的温度衰减曲线的时域数据。

图1为本发明的一实施例的热图自回归滑动平均谱估计流程图，所采用的自回归滑动平均模型为一10阶自回归10阶滑动平均模型。首先计算被输入的所述时域数据的互相关函数，生成自回归(10)和滑动平均(10)的正则方程。求解自回归(10)的正则方程及病变区的极点参数。将所述参数代入滑动平均(10)的正则方程中，求解病变区的零点参数。根据求出的病变区的零极点参数，计算病变区的自回归滑动平均(10，10)功率谱曲线。将计算出来的自回归(10)模型的参数同标准病理模板中的自回归(10)参数做互相关系数计算，得到辅助诊断指标一。将计算出的滑动平均(10)模型的参数同标准病理模板中的滑动平均(10)参数做互相关系数的计算，得到辅助诊断指标二。对指标一和指标二综合考虑从而为临床诊断提供有关数据。其中标准病理模板的自回归(10)参数和滑动平均(10)参数以及自回归滑动平均(10，10)功率谱曲线的计算方法与病变区的所述参数及所述功率谱曲线的计算方法相同。

Claims (2)

1、一种基于自回归滑动平均谱分析的医用红外热图分析方法，包括：对待测人体表面区域进行温度刺激；获取该体表区域温变过程的多幅红外热图；对该多幅红外热图的温度数据进行谱分析以获取表征该体表区域温度特征的参数，用作为临床诊断数据；其特征在于，

从多幅红外热图中获取待测人体表面区域所对应的动态过程体表温度的时域温度数据；

利用该有限的时间序列温度数据来对分析对象建立模型，对该时域温度数据作现代谱参数估计；

将该动态过程体表温度的变化转化为谱特征参数，得到病变区域对温度变化的频率响应过程；

将温度的变化的谱特征参数同标准病理模板的谱特征参数做相关计算，其结果提供用于临床诊断的数据；

所述的现代谱参数估计为自回归滑动平均谱参数估计，包括以下步骤：

1)计算被输入的多幅病变区医用红外热图的温度变化时域数据，生成自回归模型和滑动平均模型的正则方程，根据模型适用性检验准则来确定自回归模型和滑动平均模型的阶次；

2)求解自回归模型的正则方程及代表该病变区所测点的温度变化响应过程的极点参数，将所述参数代入滑动平均模型的正则方程中，求解代表该病变区所测点的温度变化响应过程的零点参数，用该模型估计出来的零、极点参数即表征了所述病变区上的所述点对温度刺激的响应特征；

3)由自回归模型和滑动平均模型生成自回归滑动平均模型的正则方程，从中计算代表该病变区所测点的温度变化响应过程的自回归滑动平均功率谱曲线；

4)将计算出来的自回归模型的参数同标准病理模板中的自回归模型参数做互相关系数计算，将该相关的程度作为对该疾病发生的可能性的辅助诊断指标一；

5)将计算出来的滑动平均模型的参数同标准病理模板中的滑动平均模型参数做互相关系数计算，将该相关的程度作为对该疾病发生的可能性的辅助诊断指标二；

6)综合作为辅助诊断指标一的相关系数和作为辅助诊断指标二的相关系数，得到总的加权相关系数，该系数的大小即代表了该疾病产生的可能性的大小。

2、如权利要求1所述的基于自回归滑动平均谱分析的医用红外热图分析方法，其特征在于，所说的自回归模型的阶次、滑动平均模型的阶次及自回归滑动平均模型的阶次均为10。