一种基于近红外光的膝骨性关节炎病程检测方法与流程

本发明属于医学技术领域，尤其涉及一种基于近红外光的膝骨性关节炎病程检测方法。

背景技术：

膝骨关节炎(osteoarthritis，oa)是一种常见的不可逆慢性关节疾病，发病人群多为60岁以上的老年人。随着社会老龄化，膝骨关节炎的发病率越来越高，据国内统计，65岁以上的老年人发病率为75％，对中老年人的健康和生活质量造成了较大影响。由于膝oa早期的症状不明显，经常被人忽视，而一旦发展到晚期，不仅严重影响患者的身心健康，也带来了巨大的经济负担，因此，寻求侵入小、检测成本低的膝oa早期病程跟踪检测方法，是极有价值的研究课题。膝骨关节炎不仅影响滑膜关节结构的质量，还影响周围组织的功能及感受信号通路，经常表现为关节疼痛，僵硬和运动性能降低。为诊断由膝oa影响的关节中的软骨异常和滑液组成的变化，临床医学采用的成像方法主要包括：放射x光，计算机断层扫描(ct)和核磁共振成像(mri)、骨骼超声波等。其中，x射线由于其分辨率高，易于使用和成本低的特点，成为诊断骨关节疾病被广泛接受的方法。x射线在将骨骼与软组织分化方面具有极高的分辨率，但由于不同软组织之间的成像对比度较低，因此不能很好地区分密度相似的软组织。所以，平片x光片虽能够显示关节间隙狭窄和骨赘形成，却对软骨和液体的变化不敏感，难以捕获早期膝oa的主要特征。ct扫描也常用于诊断oa，由于接受ct检查时会受到一定程度的辐射影响，不宜作为跟踪病程的常规化检测手段，一般只用于膝oa患病晚期检测。mri虽然在使用高对比剂时可以更可靠地检测早期的oa，但该方法不仅昂贵而且耗时，对操作人员的要求也较高。临床医学中，骨骼超声也是评估类风湿关节炎和区域肌肉骨骼病理学的主要方法，但是该方法对边界层的变化敏感，在oa的早期阶段评估中的作用有限。

综上所述，现有技术存在的问题是：(1)会对人体造成辐射影响，无法作为常规体检项目，影响了病情的及早发现与治疗；(2)从图像上无法直接关节软骨和关节腔液体成分的变化，无法发现早期患病信息；(3)进行膝关节成像检查和判定病情时对医疗人员的专业技术要求较高。

技术实现要素：

针对现有技术存在的问题，本发明提供了一种基于近红外光的膝骨性关节炎病程检测方法。

本发明是这样实现的，一种基于近红外光的膝骨性关节炎病程检测方法，所述基于近红外光的膝骨性关节炎病程检测方法，首先采用图像分割处理技术对临床膝关节ct图进行处理，留下肌肉和骨骼组织部分并进行灰度对比加深和边缘提取；随后将处理后的图片按照关节炎患病早期、中期和晚期分成三类；用等量红外光子从同一位置入射上述三类图片，采用蒙特卡洛方法，通过改变关节腔滑液的吸收和散射系数，分别模拟近红外光子在患关节炎早期、中期和晚期的膝关节内部运动轨迹，最后通过高斯函数拟合不同病程下的红外光子出射分布特征；

其中，高斯函数的自变量为出射光子的体素横坐标，因变量为在对应坐标处出射的光子数；最后用统计学方法统计三种病程条件下，运动轨迹通过关节腔并最终出射的有效光子数占所有入射光子数的比率，以及拟合后的高斯函数的对称轴坐标，以此作为判定膝关节炎的患病程度的指标，以范围形式描述。对待检测患者的ct图做相同近红外光子运动轨迹模拟和高斯拟合，计算指标值。若两项指标均落在相同病程范围内，则判定患者的膝关节炎患病程度为对应病程。

进一步，对膝关节ct图片用蒙特卡洛方法模拟红外光子在关节内部的运动轨迹，包括：

临床膝关节原始ct图片处理：

模拟红外光子在膝关节内部的运动轨迹前，采用oa患者的临床检测膝关节图片作为决定光子运动方式的基本信息来源；进行蒙特卡洛模拟时，对骨骼和肌肉组织进行光学特性分析，需剔除毛细血管部分，并且将余下的组织边界清晰化；

近红外光子在膝关节内的运动轨迹模拟：

首先将处理后的临床膝关节ct图片沿x轴和y轴平均分成若干方块，每一个方块称为一个体素，根据每一个体素中的组织类型进行分类，用以计算光子在该体素内的传播方向和步长；

将所有ct图片分成512*512个体素，根据每一个体素的灰度值，将其标记为外界、肌肉组织以及骨骼组织；红外光子从外界和肌肉组织边界垂直入射后，根据当前所处位置的体素类型，判断光子与膝关节组织的作用方式，按照以下规则进行步长和散射方向抽样的迭代计算，每走一步，光子能量都根据吸收系数被吸收一部分，直至能量降为0,光子消亡；

(1)步长s仅与散射系数μs相关，由公式(1)抽样获得；其中，ξ为(0,1)内均匀分布的随机变量；

s＝-ln(1-ξ)/μs＝-ln(ξ)/μs(1)；

(2)光子新前进方向d由公式(2)计算；其中，散射角由henyey-greenstein相位函数获得，方位角θ在(0，2π)内均匀分布；

进一步，通过高斯函数分析和拟合不同病程下的出射光子分布特征，包括：

基于高斯函数的光子分布特征拟合：

记待拟合的自变量矩阵为x＝[x(1)x(2)…x(n)]^t，因变量矩阵为y＝[y(1)y(2)…y(n)]^t，n为数据总量；

高斯函数由公式(3)描述：

其中，待估参数ymax,xmax,σ分别为高斯曲线的峰值、峰值对应横坐标和半宽。其中，σ与y(k)呈非线性关系，若直接进行求解，采用非线性寻优算法；为简化寻优难度，避免非线性寻优，对公式(3)等号两边同时取以自然对数e(e≈2.718)为底的对数，得：

令：

则公式(4)写为：

经过此处理，将原有的求解ymax,xmax,σ命题转化成了求解b0,b1,b2；由于b0,b1,b2与z(k)呈线性关系，用最小二乘法进行求解；定义待估参数矩阵θ＝[b0b1b2],数据矩阵z和φ分别为：

通过极小化拟合误差求解θ,得其最小二乘解为：

将代入公式(8),得高斯函数的原始待估参数ymax,xmax,σ分别为：

进一步，所述临床膝关节原始ct图片处理，具体包括：

步骤一：原始ct图片背景、前景分割；

首先求出原始ct图中的最大灰度值t1和最小灰度值t2，定义自定义系数n∈[0,1]，设定图像的前景、背景分割阈值t的初值t0：

t0＝n×(t1+t2)(10)；

用k表示迭代次数，根据当前阈值tk将图像分割为前景和背景，分别求出两者的平均灰度值ta和tb，根据公式(11)计算下一次迭代的新阈值tk+1：

tk+1＝n×(ta+tb)(11)；

定义一个极小的正数ε为迭代停止精度标准。若满足tk+1-tk＜ε，则停止迭代，否则继续迭代；迭代结束后，根据最终阈值将原始ct图片分割为前景和背景部分；

步骤二：边缘锐化：

经过步骤一处理后，膝关节骨骼组织提取完毕，将原图中相邻像素点由0变化至非零点的位置标记为边缘，将步骤一获得的分割后图像和锐化后的边缘叠加，并将边缘和骨骼中部用灰色填充。

本发明的另一目的在于提供一种基于近红外光的膝骨性关节炎病程检测系统。

本发明的优点及积极效果为：(1)无损性：相比于x射线和ct检测手段，近红外光子检测是一种无损检测技术，不会给人体带来任何损害，可以作为中老年人常规体检项目，便于及早发现病情；(2)灵敏性：通过捕捉出射红外光子的分布信息，能够发现早期关节炎引起的关节腔滑液的微小病变信息，滑液浓度变化10％以上即可在出射光子的分布特征上体现出区别，相比于x射线，在关节炎的早期判定上更加灵敏有效；(3)病程判定便捷，给出量化指标，通过计算机程序进行数值分析在5分钟内给出病程判断结果，客观性较强，相比于mri和超声波诊断技术，对检测人员的专业水平要求较低；(4)准确率高：本发明的实验结果表明该方法的准确率达到92％以上，在膝骨性关节炎的临床检测应用上具有较大的应用价值。

附图说明

图1是本发明实施例提供的分割前后的膝关节ct图。

图中：a、图为原始ct图；b、图为进行分割处理后的结果图。

图2是本发明实施例提供的处理后的膝关节ct图。

图3是本发明实施例提供的体素化的膝关节ct图。

图4是本发明实施例提供的不同膝oa病程下出射红外光子分布情况高斯拟合图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

常规的计算机断层扫描检测手段由于会带来辐射影响，无法作为常规体检项目，影响了病情的及早发现与治疗。

下面结合附图对本发明的应用原理作详细描述。

本发明实施例提供的基于近红外光的膝骨性关节炎病程检测方法包括：结合临床膝关节ct图片用蒙特卡洛方法模拟红外光子在关节内部的运动轨迹；通过高斯函数分析和拟合不同病程下的出射光子分布特征；以有效光子出射率和拟合函数对称轴位置作为判定的指标。

下面结合具体实施例对本发明的应用原理作进一步描述。

本发明实施例提供的基于近红外光的膝骨性关节炎病程检测方法中，

1、基于蒙特卡洛的近红外光子运动轨迹模拟：

近红外光谱技术因其分析速度快、对检测样本具有无损性等优点，被广泛地用于生物医学领域。近红外光波段的光学窗口是600-1300nm,在该波段内，人体内含有的水对近红外光的吸收很低，若用近红外光波照射膝关节，光子能够在膝关节组织内传输较远的距离，并经过一系列碰撞、散射和反射作用后，携带深层关节内部信息溢出组织边界。根据近红外光这一特性，本发明首先采用蒙特卡洛方法，模拟近红外光子从膝关节外部入射后在关节内部的运动轨迹。蒙特卡洛方法是一种通过采集大量随机样本模拟某一物理过程并得到其统计规律的方法，不需要对组织的几何参数作限制。采用该方法，通过设置吸收系数，散射系数等参数模拟膝关节组织模型，方便灵活，无需从数学方程或表达式出发，当实验次数足够多时，能够有效反映光子在膝关节内部的运动规律。将一束射入膝关节的近红外光看成大量单个光子依次以相同入射点、相同角度进入关节内部，结合临床膝关节ct图片中的骨头和肌肉组织分布情况，根据光子所遇组织的光学特性，迭代产生膝关节组织中单个光子的行走步长和方向，追踪每一个光子的行进轨迹直到其能量在传播过程中衰减至零。

1.1临床膝关节原始ct图片处理

模拟红外光子在膝关节内部的运动轨迹前，为准确描述膝关节内部组织结构，本发明采用oa患者的临床检测膝关节图片作为决定光子运动方式的基本信息来源。在膝关节原始ct图中，包含了骨骼、肌肉组织以及细小的毛细血管。进行蒙特卡洛模拟时，仅针对骨骼和肌肉组织进行光学特性分析，故需剔除毛细血管部分，并且将余下的组织边界清晰化。处理步骤如下：

步骤1：原始ct图片背景、前景分割；

首先求出原始ct图中的最大灰度值t1和最小灰度值t2，定义自定义系数n∈[0,1]，设定图像的前景、背景分割阈值t的初值t0：

t0＝n×(t1+t2)(1)；

用k表示迭代次数，根据当前阈值tk将图像分割为前景和背景，分别求出两者的平均灰度值ta和tb，根据公式(2)计算下一次迭代的新阈值tk+1：

tk+1＝n×(ta+tb)(2)；

定义一个极小的正数ε为迭代停止精度标准。若满足tk+1-tk＜ε，则停止迭代，否则继续迭代。迭代结束后，根据最终阈值将原始ct图片分割为前景和背景部分。背景部分为膝关节外界，显示为纯黑色，前景部分为膝关节内部的骨骼组织，显示为亮白色。

图1展示了经过此步骤处理后的ct图片，其中，a图为原始ct图片，b图为进行分割处理后的结果。

步骤2：边缘锐化。

经过步骤1处理后，膝关节骨骼组织提取完毕，但膝关节边缘和内部肌肉组织随之消失。由于原始图中，膝关节外侧为纯黑色，内侧为灰色，对比明显，将原图中相邻像素点由0变化至非零点的位置标记为边缘，将步骤1获得的分割后图像和锐化后的边缘叠加，并将边缘和骨骼中部用灰色填充，表示膝关节肌肉组织。图2为最终处理结果。

1.2近红外光子在膝关节内的运动轨迹模拟：

膝关节组织内部结构复杂，组织中含有骨头、肌肉等成分。各成分之间的光学性质具有很大的差异。红外光子进入膝关节内部后，其传播方向、步长将随着所遇组织的不同而发生改变。入射瞬间，光子遇到空气与肌肉组织边界，折射率发生改变，在分界面同时发生透射和全反射；进入膝关节内部后，若在肌肉组织中传播，则按照随机概率被散射粒子散射或被吸收粒子吸收；若遇到骨组织，则发生全反射；同时，随着运动距离的增大，光子能量逐步被吸收，步长逐渐减小，直至消亡。研究显示，当膝关节发生病变时，内部组织的光学性质也随之发生改变。随着关节病情的加重，骨组织的形态不会发生明显改变，但关节腔内积水增多，伴随着粘多糖的消失，关节腔处的滑液浑浊度增加，导致组织的吸收系数，散射系数变大。因此，可以通过研究膝关节的光学特性，作为组织病变的早期诊断依据。

为方便进行红外光子的运动轨迹模拟，如图3所示，首先将处理后的临床膝关节ct图片沿x轴和y轴平均分成若干方块，每一个方块称为一个体素，根据每一个体素中的组织类型进行分类，用以计算光子在该体素内的传播方向和步长。体素越小,越能够保证每一个体素内的组织类型唯一，从而提高模拟精度。将所有ct图片分成512*512个体素，根据每一个体素的灰度值，将其标记为外界(黑色)、肌肉组织(灰色)以及骨骼组织(白色)。用体素的编号作为ct图片上的横纵坐标，即横纵坐标范围均为[0,512]。红外光子从外界和肌肉组织边界垂直入射后，根据当前所处位置的体素类型，判断光子与膝关节组织的作用方式，按照如下规则进行步长和散射方向抽样的迭代计算，每走一步，光子能量都根据吸收系数被吸收一部分，直至能量降为0,光子消亡。

(1)步长s仅与散射系数μs相关，由公式(3)抽样获得。其中，ξ为(0,1)内均匀分布的随机变量。

s＝-ln(1-ξ)/μs＝-ln(ξ)/μs(3)；

(2)光子新前进方向d由公式(4)计算。其中，散射角由henyey-greenstein相位函数获得，方位角θ在(0，2π)内均匀分布。

2膝骨性关节炎病程检测：

随着膝oa患病程度的增加，膝关节组织的光学特性随之发生改变，这一特性使通过分析近红外光射入膝关节后的出射分布特征来判断膝oa病程成为可能。为使验证效果更具真实性，以下所有实验均采用项目合作方--厦门大学附属中山医院提供的临床医学膝关节ct图片进行测试。实验样本共包含300位病患的膝关节ct图片，其中，100位患者的膝oa病程属于早期，100位处于中期，100位处于晚期。在三个病程中各选择50份ct图作为实验用例，用于执行蒙特卡洛仿真程序并提取出射红外光子的分布特征，余下作为验证数据，检验实验结果的准确性。

2.1出射红外光子分布特征提取

现有技术中，基于图像重建方法，分别给出了膝oa患者关节组织液吸收系数和散射系数随着病情加重的变化规律，进行蒙特卡洛仿真时，根据不同病程分别设置吸收和散射系数。

首先基于构建的蒙特卡洛程序寻找合适的红外光子入射点。由于膝oa主要造成关节腔内滑液光学特性的改变，为保证出射光子携带尽可能多的关节腔内滑液信息，要求近红外光子运动路径尽量通过关节腔；

为方便后续进行光子分布特征拟合，首先需处理和提取不同膝oa病程下的数字化出射红外光子信息。构建的蒙特卡洛程序中已保存所有光子从入射至消亡每一次运动的步长、方向以及横、纵坐标。提取符合如下2个条件的出射红外光子坐标：(1)光子行进路径通过关节腔；(2)光子最终位置的横纵坐标处于膝关节肌肉组织和外界分界处。

对所有测试数据进行相同实验，采用相同入射光子数量和相同的入射点、入射角度，记录每一组实验的出射红外光子坐标。

2.2基于高斯函数的光子分布特征拟合

多项式函数和高斯函数是进行数据拟合最常用的两种函数。记待拟合的自变量矩阵为x＝[x(1)x(2)…x(n)]^t，因变量矩阵为y＝[y(1)y(2)…y(n)]^t，n为数据总量。

以三阶多项式函数为例，该函数需寻找参数矩阵[abcd],使其满足公式(5)的条件:

y(k)＝ax³(k)+bx²(k)+cx(k)+d,(k＝1,2,...n)(5)；

三阶多项式函数有4个待定参数，各参数均与y(k)呈线性关系，可以用最小二乘法极小化拟合误差直接求解，无需采用非线性寻优算法。虽然参数求解方法简便，但本文旨在通过红外光子的出射分布特征快速判定膝oa病程，若分布特征由4个甚至更多待定参数共同描述，将增大病程判断的界定难度。前期实验结果显示，不同病程下的出射光子分布均呈现正态分布特征，仅在对称轴坐标上表现出差异，高斯函数图像也符合此特点，且该函数仅由1个待估参数来描述对称轴位置，能够方便地用于病程界定，在此采用高斯函数对光子分布特征进行拟合。高斯函数可由公式(6)描述：

其中，待估参数ymax,xmax,σ分别为高斯曲线的峰值、峰值对应横坐标和半宽。其中，σ与y(k)呈非线性关系，若直接进行求解，需采用非线性寻优算法。为简化寻优难度，避免非线性寻优，对公式(6)等号两边同时取以自然对数e(e≈2.718)为底的对数，得：

令：

则公式(7)可写为：

经过此处理，将原有的求解ymax,xmax,σ命题转化成了求解b0,b1,b2。由于b0,b1,b2与z(k)呈线性关系，可以用最小二乘法进行求解。定义待估参数矩阵θ＝[b0b1b2],数据矩阵z和φ分别为：

通过极小化拟合误差求解θ,得其最小二乘解为：

将代入公式(8),得高斯函数的原始待估参数ymax,xmax,σ分别为：

分别对提取的膝oa患病早期、中期和晚期情况下的出射红外光子坐标进行高斯函数数据拟合。经过对大量测试数据对比分析得到，早期、中期和晚期病程下，出射红外光子在体素化后的ct图片x轴上的分布基本一致，但在ct图片y轴上的分布却有明显差异，同时，随着病程的发展，出射光子的数量也急剧减少。故进行红外光子分布特征的高斯拟合时，令自变量(横轴)为出射光子纵坐标，因变量(纵轴)为对应自变量坐标下的出射光子数。

图4展示了采用从提取的光子分布信息进行高斯拟合的结果。从图中可以看出，随着膝oa病情的加重，出射光子数显著减少，与关节腔滑液吸收系数增大，在传输过程中被吸收的光子数增多的理论吻合；同时，出射光子数最多处所对应的光子纵坐标(即高斯函数的对称轴)也发生了改变，随着病情加重，该值逐渐向右偏移。

2.3膝oa病程判定与检验

针对前文提到的用于实验的处于膝oa早期、中期和晚期各50例，共150例原始ct图片进行出射近红外光分布特征的高斯拟合。各病程下出射光子分布特征如表1所示。表1中，不同病程下的有效光子出射率和高斯函数参数xmax均以数值范围表示，表示所有测试用例中，此两项指标的最小值和最大值。

表1膝oa不同病程下的出射红外光子分布特征

以表1列出的膝oa各病程出射光子分布特征作为依据，对剩余150例ct图片同样进行出射近红外光子分布的高斯拟合并记录其有效光子出射率和高斯函数参数xmax。若这两个指标均处于表1中同一个病程范围内，则判定该样本病患处于该病程。若两个指标分属不同病程，考虑到不同病程间的有效光子出射率差异明显，则以其作为判断依据。表2基于近红外光检测的膝oa病程判定准确率。

表2基于近红外光检测的膝oa病程判定准确率

从表2中可以看出，根据本文提出的方法进行的膝oa病程检测，准确率达到了92％上，能够有效地用于临床膝oa的无损快速检测中，该方法由于采用对人体无损的近红外光进行检测，可以作为中老年人日常体检的检查项目，对膝oa的早期发现、及时治疗干预具有较大的临床应用前景。

3、结论

针对传统膝骨性关节炎ct检测手段容易给患者带来辐射影响，无法作为日常体检项目持续跟踪病情，导致未进行及时治疗而加重病情的问题，本发明提出了基于近红外光的膝oa病程快速检测方法，采用近红外光从特定位置照射膝关节，通过检测和分析携带关节腔滑液信息的出射光子分布特征直接判定患者所处膝oa病程。本发明基于临床患者膝关节的大量ct检测图片，以蒙特卡洛法分别模拟近红外光子患病早期、中期和晚期的膝关节内部的运动路径，并用高斯函数拟合出射光子分布特征，通过有效光子出射率以及拟合函数对称轴坐标两项指标判定患者的患病程度。通过对大量实验数据分析，验证该方法的判定准确率抵92％以上。若将该方法用于临床诊断，由于是非侵入式检测手段，不会对患者造成明显影响，可以作为中老年人常规体检项目，以便及时发现病情进行干预治疗。本发明采用的是蒙特卡洛方法进行红外光子运动路径仿真，后续分析将着力于通过在患者膝关节上放置红外光子探测装置，直接获取光子出射信息，进一步检验本文方法的准确性和有效性。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。