膝骨关节炎(Osteoarthritis，OA)是一种常见的不可逆慢性关节疾病，发病人群多为60岁以上的老年人。随着社会老龄化，膝骨关节炎的发病率越来越高，据国内统计，65岁以上的老年人发病率为75%，对中老年人的健康和生活质量造成了较大影响^[1]。由于膝OA早期的症状不明显，经常被人忽视，而一旦发展到晚期，不仅严重影响患者的身心健康，也带来了巨大的经济负担，因此，寻求侵入小、检测成本低的膝OA早期病程跟踪检测方法，是极有价值的研究课题。

膝骨关节炎不仅影响滑膜关节结构的质量，还影响周围组织的功能及感受信号通路，经常表现为关节疼痛，僵硬和运动性能降低。为诊断由膝OA影响的关节中的软骨异常和滑液组成的变化，临床医学采用的成像方法主要包括：放射X光，计算机断层扫描(CT)和核磁共振成像(MRI)^[2]、骨骼超声波^[3]等。其中，X射线由于其分辨率高，易于使用和成本低的特点，成为诊断骨关节疾病被广泛接受的方法。X射线在将骨骼与软组织分化方面具有极高的分辨率，但由于不同软组织之间的成像对比度较低，因此不能很好地区分密度相似的软组织。所以，平片X光片虽能够显示关节间隙狭窄和骨赘形成，却对软骨和液体的变化不敏感，难以捕获早期膝OA的主要特征^[4-5]。CT扫描也常用于诊断OA，由于接受CT检查时会受到一定程度的辐射影响，不宜作为跟踪病程的常规化检测手段，一般只用于膝OA患病晚期检测。MRI虽然在使用高对比剂时可以更可靠地检测早期的OA，但该方法不仅昂贵而且耗时，对操作人员的要求也较高。临床医学中，骨骼超声已经成为评估类风湿关节炎和区域肌肉骨骼病理学的主要方法，该方法对边界层的变化敏感，在OA的早期阶段评估中的作用有限^[6]。考虑到膝OA的一旦患病则不可逆转的特性，寻找一种能够检测早期膝OA，并且能作为常规体检项目的无损检测手段成为亟待解决的实际问题。近年来，近红外光谱技术(Near Infrared Reflectance Spectroscopy，NIRS)由于其无创、成本低廉等优点，在人体的研究与应用成为国内外新兴的研究领域，在脑力疲劳度检测、血液含氧量检测等领域纷纷发挥重要作用^[7]。但是，采用NIRS技术进行人体膝OA检测的相关研究依然很少，文献[8]最早提出了通过NIRS进行早期膝OA检测的方法，验证了NIRS在成年家兔膝OA的检测应用上的可行性。本文在该研究的基础上，采用蒙特卡洛(MC, Monte Carlo)方法仿真红外光子在人体膝关节内部的运动轨迹，通过分析比对不同病程下出射红外光子的分布特征，实现对膝OA病程的无损快速判断，探索临床膝OA病程的新型检测手段。

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

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

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

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

步骤1：原始CT图片背景、前景分割。

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

$ {T_0} = n \times \left( {{T_1} + {T_2}} \right) $

(1)

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

$ {T_{k + 1}} = n \times \left( {{T_a} + {T_b}} \right) $

(2)

定义一个极小的正数ε为迭代停止精度标准。若满足T_k+1-T_k < ε，则停止迭代，否则继续迭代。迭代结束后，根据最终阈值将原始CT图片分割为前景和背景部分。背景部分为膝关节外界，显示为纯黑色，前景部分为膝关节内部的骨骼组织，显示为亮白色。图 1展示了经过此步骤处理后的CT图片，其中，(a)图为原始CT图片，(b)图为进行分割处理后的结果。

步骤2：边缘锐化。

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

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

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

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

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

$ s = - \ln \left( {1 - \xi } \right)/{\mu _s} = - \ln \left( \xi \right)/{\mu _s} $

(3)

2) 光子新前进方向d由公式(4)计算。其中，散射角φ由Henyey-Greenstein相位函数获得^{[14, 15]}，方位角θ在(0, 2π)内均匀分布。

$ d = \left[ {\begin{array}{*{20}{c}} {\cos \varphi \sin \theta }&{\sin \varphi \cos \theta } \end{array}} \right] $

(4)

2 膝骨性关节炎病程检测

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

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

文献[16-17]基于图像重建方案，研究并分别给出了膝OA患者关节组织液吸收系数和散射系数随着病情加重的变化规律，进行蒙特卡洛仿真时，根据不同病程分别设置吸收和散射系数，具体数值如表 1所示。

首先基于第1节构建的蒙特卡洛程序寻找合适的红外光子入射点。由于膝OA主要造成关节腔内滑液光学特性的改变，为保证出射光子携带尽可能多的关节腔内滑液信息，要求近红外光子运动路径尽量通过关节腔；且需保证出射光子数量足以进行后续的出射分布特征拟合。图 4展示了500个红外光子从图中所示位置入射后，在关节内部的运动轨迹图。图中，黑色部分表示膝关节皮肤的外界，灰色部分表示膝关节内部肌肉组织，白色部分表示膝关节内部骨骼组织，两块白色骨骼组织的缝隙即为关节腔，红色路径表示入射近红外光子的运动轨迹。从图中看出，在膝关节内部运动一定距离后，部分红外光子被吸收，其余主要从两个位置出射。其中，从位置1出射的光子没有通过关节腔，这部分出射光子未携带有效关节腔内滑液信息，属于无效光子；仅有从位置2出射的光子经过关节腔滑液的吸收和散射作用，能够作为后续出射分布特征提取的数据。从图 4的出射结果看，无效光子比重过大，需重新调整入射点和入射角度。

图 5展示了采用同一个病患的早期(a)、中期(b)和晚期(c)膝关节CT图片进行蒙特卡洛模拟的近红外光出射分布情况，各组数据均从同一入射点处打入2.5万个红外光子。经大量实验测试，最终选择正对关节腔的膝关节肌肉与外界分界点作为入射点，并令所有光子从图 5(a)所示的入射点入射。测试结果显示，红外光子从此处入射后，大多数出射(最终到达关节肌肉和外界分界)光子都从关节腔中经过。从图 5的(a), (b), (c)三幅图对比可以看出，随着病情加重，膝关节骨骼组织没有发生明显的变化，但出射的红外光子数却显著减少，这是由于关节滑液浑浊度增强，大量光子在传播过程中被吸收导致。同时，由于散射系数增强，改变了光子运动方向, 光子出射的位置分布情况也略有不同。临床采用近红外光快速检测膝OA病程时，可在入射点处照射红外光，在图 5所示的出射位置附近安装红外光子探测装置。

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

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

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

多项式函数和高斯函数是进行数据拟合最常用的两种函数。记待拟合的自变量矩阵为

$ \mathit{\boldsymbol{X}} = {\left[{\begin{array}{*{20}{c}} {x\left( 1 \right)}&{x\left( 2 \right)}& \cdots &{x\left( N \right)} \end{array}} \right]^{\rm{T}}}, $

因变量矩阵为

$\mathit{\boldsymbol{Y}} = {\left[{\begin{array}{*{20}{c}}{y\left( 1 \right)}&{y\left( 2 \right)}& \cdots &{y\left( N \right)}\end{array}} \right]^{\rm{T}}}, $，N为数据总量。

以三阶多项式函数为例，该函数需寻找参数矩阵$\left[{\begin{array}{*{20}{c}}a&b&c&d\end{array}} \right]$，使其满足公式(5)的条件:

$ \begin{array}{*{20}{c}} {y\left( k \right) = a{x^3}\left( k \right) + b{x^2}\left( k \right) + cx\left( k \right) + d, }\\ {\left( {k = 1, 2, \cdots N} \right)} \end{array} $

(5)

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

$ \begin{array}{*{20}{c}} {y\left( k \right) = {y_{\max }}\exp \left[{{{\left( {\frac{{x\left( k \right)-{x_{\max }}}}{\sigma }} \right)}^2}} \right], }\\ {k = 1, 2, \cdots N} \end{array} $

(6)

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

$ \begin{array}{*{20}{c}} {\ln y\left( k \right) = \ln {y_{\max }}- \frac{{{{\left[{x\left( k \right)-{x_{\max }}} \right]}^2}}}{\sigma } = }\\ {\left( {\ln {y_{\max }} -\frac{{x_{\max }^2}}{\sigma }} \right) + \frac{{2x\left( k \right){x_{\max }}}}{\sigma } -\frac{{{x^2}\left( k \right)}}{\sigma }} \end{array} $

(7)

令：

$ \begin{array}{*{20}{c}} {\ln y\left( k \right) = z\left( k \right), \ln {y_{\max }}-\frac{{x_{\max }^2}}{\sigma } = {b_0}, }\\ {\frac{{2{x_{\max }}}}{\sigma } = {b_1}, -\frac{1}{\sigma } = {b_2}} \end{array} $

(8)

则公式(7)可写为：

$ \begin{array}{*{20}{c}} {z\left( k \right) = {b_0} + {b_1}x\left( k \right) + {b_2}{x^2}\left( k \right) = }\\ {\left[{\begin{array}{*{20}{c}} 1&{x\left( k \right)}&{{x^2}\left( k \right)} \end{array}} \right]\left[{\begin{array}{*{20}{c}} {{b_0}}\\ {{b_1}}\\ {{b_2}} \end{array}} \right]} \end{array} $

(9)

经过此处理，将原有的求解y_max, x_max, σ命题转化成了求解b₀, b₁, b₂。由于b₀, b₁, b₂与z(k)呈线性关系，可以用最小二乘法进行求解。定义待估参数矩阵$\mathit{\boldsymbol{\theta }} = \left[{\begin{array}{*{20}{c}}{{b_0}}&{{b_1}}&{{b_2}}\end{array}} \right]$，数据

矩阵Z和Φ分别为：

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{Z}} = {{\left[{\begin{array}{*{20}{c}} {z\left( 1 \right)}&{z\left( 2 \right)}& \cdots &{z\left( N \right)} \end{array}} \right]}^{\rm{T}}}, }\\ {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} = \left[{\begin{array}{*{20}{c}} 1&{x\left( 1 \right)}&{{x^2}\left( 1 \right)}\\ 1&{x\left( 2 \right)}&{{x^2}\left( 2 \right)}\\ \vdots&\vdots&\vdots \\ 1&{x\left( N \right)}&{{x^2}\left( N \right)} \end{array}} \right]} \end{array} $

(10)

通过极小化拟合误差${\varepsilon _{OE}} = \sum\limits_{k = 1}^N {{{\left[{z\left( k \right)-\hat z\left( k \right)} \right]}^2}} $求解θ，得其最小二乘解${\mathit{\boldsymbol{\hat \theta }}}$为：

$ \mathit{\boldsymbol{\hat \theta }} = {\left[{{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}} \right]^{ -1}}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\rm{T}}}\mathit{\boldsymbol{Z}} $

(11)

将${\mathit{\boldsymbol{\hat \theta }}}$代入公式(8)，得高斯函数的原始待估参数y_max, x_max, σ分别为：

$ \sigma =-\frac{1}{{{b_2}}};\;\;\;{x_{\max }} =-\frac{{{b_1}}}{{2{b_2}}};\;\;\;{y_{\max }} = \exp \left( {{b_0}-\frac{{b_1^2}}{{4{b_2}}}} \right) $

(12)

分别对2.1节中提取的膝OA患病早期、中期和晚期情况下的出射红外光子坐标进行高斯函数数据拟合。经过对大量测试数据对比分析得到，早期、中期和晚期病程下，出射红外光子在体素化后的CT图片x轴上的分布基本一致，但在CT图片y轴上的分布却有明显差异，同时，随着病程的发展，出射光子的数量也急剧减少。故进行红外光子分布特征的高斯拟合时，令自变量(横轴)为出射光子纵坐标，因变量(纵轴)为对应自变量坐标下的出射光子数。图 6展示了采用从图 5提取的光子分布信息进行高斯拟合的结果。从图中可以看出，随着膝OA病情的加重，出射光子数显著减少，与关节腔滑液吸收系数增大，在传输过程中被吸收的光子数增多的理论吻合；同时，出射光子数最多处所对应的光子纵坐标(即高斯函数的对称轴)也发生了改变，随着病情加重，该值逐渐向右偏移。

2.3 膝OA病程判定与检验

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

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

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

3 结论

1) 模拟近红外光子从膝关节CT图片上正对关节腔的位置入射，能够使大多数光子的运行路径经过关节腔，从而携带尽可能多的关节腔滑液病变信息出射，提高病程检测的准确度。

2) 随着膝OA患病程度的加重，骨组织的形态不会发生明显改变，但关节腔内滑液浑浊度增加，组织的吸收系数，散射系数随之变大，以等量近红外光子从同一位置入射膝关节，出射的红外光子数也随之减少，分布规律也发生改变。

3) 采用高斯函数对不同膝OA病程下的出射红外光子分布特征进行拟合，可知病情越重，高斯函数的对称轴坐标值越大。仿真实验结果显示，同时以有效近红外光子出射率和高斯函数对称轴坐标作为判定膝OA病程的两项指标，判定准确率达92%以上。