CN104166994B - 一种基于训练样本优化的骨骼抑制方法 - Google Patents

Abstract

一种基于训练样本优化的骨骼抑制方法，其特征步骤是：①、根据图像灰度与形状信息对肺区进行最优分割；②、使用B样条小波变换的特征和jet进行特征提取，提取有效描述不同尺度骨骼的特征图像是模型建立的基础，特征图像提取包括三阶B样条多尺度小波变换及多尺度N‑jet提取；③、使用Q型因子分析法对样本进行处理；④、使用BP神经网络预测骨骼图像，建立肺部X光图像与DES肋骨图像之间的预测模型；⑤、对正常胸片与预测的骨骼图像做图像减法来预测软组织图像，建立预测模型后，根据肺部x光图像的灰度值分布来预测产生骨结构图像在肺部x光图像中减去预测获得的骨骼图像就可以获得软组织图像。

Description

一种基于训练样本优化的骨骼抑制方法

技术领域

本发明涉及成像技术技术领域，特别是一种采用胸部x光图像经过特征提取到的特征图像与双能减影产生的对应骨骼图像作为训练数据的一种基于训练样本优化的骨骼抑制方法。

背景技术

肺癌是当前对人类健康危害最大的恶性肿瘤之一。因为肺是人体内部器官，多数肺癌在开始的时候只是在身体内悄悄地生长，患者没有任何感觉。当患者因咳嗽、咯血及胸痛等临床症状就诊时，80％的患者已经属于中晚期，其中约75％的患者发现有转移病灶的存在，治疗预后较差。美国国立癌症研究所统计资料显示，如果癌症患者Ⅰ期发现病治疗，其五年生存率可达90％，而诊断治疗发现在Ⅱ期以后，则五年急剧下降至20％以下。因此，肺癌的早期诊断与治疗是提高肺癌患者生存率的关键。在对肺部疾病进行诊断时，医学影像检查应用十分广泛，它是肺癌早期检测的重要手段。随着医学成像技术的飞速发展，形成了x线诊断学、计算机断层扫描成像(CT)、核磁共振成像(MRI)、超声、核医学成像等。x线成像因为其辐射剂量小，价格相对便宜的优势成为胸部疾病的影像诊断的主体，然而，在胸部x光图像中，约有12％-90％的可见肺结节(潜在肺癌病灶)患者会被放射科医生漏诊。其中，82％-95％的漏诊结节位于图像中肋骨或锁骨结构遮挡区域。因此，肺组织和胸部骨骼重叠区域就成为胸部X光图像病变诊断的检查盲区。双能减影(Dual Energy Subtraction,DES)是在数字胸部x线摄影基础上发展的一种较新的成像技术。DES技术根据骨骼与软组织对X线光子能量衰减方式不同，以及不同原子量物质的光电吸收效应差别，利用数字摄影将两种吸收效应的信息进行分离，选择性去除骨骼或软组织的衰减信息，进而获得胸部X光图像、软组织图像及骨骼组织图像。DES技术能够消除骨骼区域对肺部病灶的遮挡，达到提高肺部疾病诊断性能的作用。

虽然DES成像技术能够更清晰的显示肺内病灶，但患者受到的放射线剂量却显著增加。对于成像质量相近的DES与普通胸部X光(Digital radiography,DR)图像，被检者的放射线吸收剂量存在显著差异，且DES的吸收剂量明显高于DR。

国内对胸部X光图像的计算机后处理研究主要集中在图像增强、肺区域分割及肺结节检测等方面，目前研究主要集中在肋骨分割。但是，由于肋骨结构较为复杂，且重叠的部分较多，肋骨的准确分割较为困难，效果并不理想。国外一些学者利用胸部X光图像对骨骼区域的形状及灰度进行估计，以实现骨骼抑制，达到生成无骨骼遮挡的软组织图像的目的。

发明内容

针对现有技术存在的不足，本发明的目的是提供一种采用胸部x光图像经过特征提取得到的特征图像与双能图像产生的对应骨骼图像作为训练数据，通过Q型因子分析对训练样本优化，再利用BP神经网络建立回归模型的一种基于训练样本优化的骨骼抑制方法。

为了实现上述目的，本发明所采用的技术方案为一种基于训练样本优化的骨骼抑制方法，其步骤如下：

步骤一：根据图像灰度与形状信息的对肺区进行分割；

步骤二：采用B样条小波变换和2-jet进行特征提取，提取有效描述不同尺度骨骼的特征图像；

步骤三：使用Q型因子分析法对样本进行处理，其过程如下：

1、对大样本数据进行预处理；

每幅特征图片有p个特征，每个特征得到n个样本，组成数据矩阵X，其表示如下：

用向量表示为

x＝(x₁,x₂,…,x_n) (2)

其中，x_i(i＝1,2,...,n)表示某个特征的第i个样本的值；

2、对数据矩阵X进行按列标准化，采用零均值标准差标准化的方法；

3、求出样本间的相似系数矩阵Q；

记第i个样本和第j个样本的观测p次测试值为n维空间中的两个点，则坐标原点分别到这两个点的n维空间向量间夹角的余弦值为相似系数，其计算公式如下：

所有相似系数构成的矩阵为相似系数矩阵；

4、求相似系数矩阵Q的特征值和相应的特征向量；

假设Q的n个特征值为λ₁,λ₂,.........λ_n，且λ₁≥λ₂≥…≥λ_n，对应的单位正交特征向量为μ₁,μ₂,…,μ_n，在n个样本中选择m个样本(m<n)来代表这组观测样本，n个样本中除m个样本外的样本由选定的m个样本线性组合表示；各个特征值的大小代表各个因子的方差在总方差中所占的比例；η_m为累计贡献率，

累计贡献率越大，说明选出的主因素越拟合于原训练样本的数据，一般累计贡献率为85％-95％；这m个公因子可以看作是m个独立的且不同于原始样本的“典型公共样本”。

5、初始因子载荷矩阵A，其元素a_ij的公式为

μ_ij是λ_j对应的第i个分量，且a_ij表示第i个样本与第j个公因子的相似系数，|a_ij|越大表明第i个样本所具有的载荷越大，第i个样本就越重要，则第i个样本就作为公因子F_j的代表性样本；

6、因子旋转；

采用方差最大正交旋转，所述正交旋转是指坐标轴始终保持90度角旋转，新生成的因子仍然可保持相关性，用公式表示如下

其中：

7、根据因子载荷矩阵G，从原始的样本中按照比例选出具有代表性的样本集。

步骤四：使用BP神经网络预测骨骼图像，建立肺部X光图像的特征图像与DES肋骨图像之间的预测模型；

假设由输入层、隐含层和输出层构成的三层BP神经网络，输入节点x_i，隐含层节点y_j，输出节点z_l；输入节点与隐含层节点间的网络权值为w_ji，隐含层节点与输出节点的网络权值为O_lj，当输出节点的期望输出为t_l时，BP模型的计算公式为：

隐含层检点的输出

其中，θ_j为隐含层的偏置；

输出节点的计算输出

其中，θ_l为输出层的偏置；

输出节点的误差公式

步骤五、对正常胸片与预测的骨骼图像做图像减法来预测软组织图像，建立预测模型后，可以根据肺部x光图像的灰度值分布来预测产生骨结构图像在肺部x光图像中减去预测获得的骨骼图像就可以获得软组织图像。

本发明的优点是通过建立胸部X光图像与对应的双能骨骼图像回归模型这个完全新颖的角度，来获得骨骼图像，然后通过产生的骨骼图像得到软组织。因此，患者在所受X线放射剂量较低的情况下，能够产生高质量无骨骼遮挡的软组织图像。将这些图像应用于肺癌等恶性疾病的早期筛查，能显著提高肺部疾病检测的敏感性与准确率。采用Q型因子分析可以合理的选择样本集，不但减少训练的时间而且提高了神经网络建模的效率。

附图说明

图1是一种基于训练样本优化的骨骼抑制方法的流程图；

具体实施方式

为了进一步了解该基于训练样本优化的骨骼抑制方法，下面结合附图说明如下。

其步骤如下：

步骤一：根据图像灰度与形状信息的对肺区进行最优分割。

先建立肺轮廓的先验模型，再利用灰度与形状相似性信息结合图像特征对肺区域分割。

1：模型初始轮廓位置的确定；

1.1：:标记训练集中每张图像肺边界的边界点；

1.2：训练样本关于一组坐标轴对齐，训练样本包括n张胸部x光图像，对n张胸部x光图像的肺轮廓标定后并将其对齐，其对齐步骤如下：

1.2.1：旋转、缩放和平移每个肺区域形状，使其与训练集中的第一个形状对齐；

1.2.2：根据对齐形状，计算平均形状；

1.2.3：旋转、缩放和平移平均形状使其与第一个形状对齐；

1.2.4：重新将每个形状与当前平均形状对齐；

1.2.5：如果过程收敛或者到指定循环次数，退出；否则转到步:1.2.2；

1.3：建立初始轮廓位置的模型；

训练样本对齐后，利用主成分分析方法找出形状变化的统计信息，据此建立模型；

2、结合灰度信息和形状信息的肺实质分割

在肺分割阶段，同时利用多张特征图像中边界点的灰度与形状信息，使得搜索到的边界灰度、形状信息与训练图像相似。具体步骤如下：

2.1：提取特征图像；

2.1.1：对图像进行高斯平滑处理，抑制图像噪声影响；

2.2：选取边界点的候选点；

对于初始肺边界的每一个点，计算所有特征图像中该点搜索区域内所有像素点的灰度与训练特征图像中相应点灰度的相似程度，选出20个相似程度最大的点，作为该边界点的候选点。相似程度为所有特征图像中该点周围像素点灰度到训练样本特征图像中相应点周围像素点灰度集合的马氏距离h_i。

2.3：使用动态规划进行肺区域分割：

在边界点搜索区域内，像素点的灰度相似性代价为该点周围像素点灰度与训练图像中相应边界点的周围像素点灰度的相似程度h_i；

步骤二、使用B样条小波变换的特征和2-jet进行特征提取；

提取有效描述不同尺度骨骼的特征图像是模型建立的基础，特征图像提取包括三阶B样条多尺度小波变换及多尺度N-jet提取。三阶B样条多尺度小波变换小波变换在空间域和时域中都具有良好的局部变换性质，可以聚焦到对象的任意细节，作为多尺度分析的一种工具。B样条函数随着样条阶数的增加而快速收敛于高斯函数，其一阶导数可以逼近最优边缘检测算子。因此，利用B样条小波进行多尺度边缘增强获得较好的性能。

本专利对胸部x光图像分别做一、二、三尺度小波分解，通过一尺度小波解获得4张图像，1张近似图像，1张水平细节图像，1张垂直细节图像，1张对角线细节图像；通过二尺度小波解得到不同于一尺度的1张近似图像，1张水平细节图像，1张垂直细节图像，1张对角线细节图像；通过三尺度小波解得到不同于二尺度的1张近似图像，1张水平细节图像，1张垂直细节图像，1张对角线细节图像；由于在胸部x光图像中，很少有呈现对角线分布的骨骼结构，一次去掉对角线细节图像，我们得到9张特征图像。

将经过多尺度小波变换后的特征图像再进行高斯滤波局部2-jet的特征提取。根据高斯滤波尺度的不同(sigma＝2,4)，每张小波变换后的图像我们可以得到12张特征图像。

步骤三、使用Q型因子分析法对样本进行处理。

神经网络的泛化能力是指在训练完成以后，输入预测样本时得到正确输出的能力。泛化能力是神经网络的最主要性能之一，如何提高泛化能力是该领域研究者关注的热点。训练样本的选择对网络的泛化能力影响非常大，本专利采用因子分析对大样本数据进行预处理。因子分析可分为R型因子分析和Q型因子分析，R型因子分析的分析对象是变量，Q型因子分析的分析对象是样本。本专利提出的是基于Q型因子分析的训练样本选择，它分析多个样本之间的关系，它可以自动地选择训练样本，而且被选的训练样本具有代表性和典型性，其他的样本都可以用这些潜在的“典型公共样本”的线性组合表示出来，这样就可以充分利用所采集的样本信息。

1、采用三阶B样条多尺度小波变换的特征图像提取对每幅图片p个特征，每个特征就得到n个样本的观测值，组成如下观测数据矩阵

向量表示为x＝(x₁,x₂,…,x_n)

其中，x_i(i＝1,2,...,n)表示某个特征的第i个样本的值。

2、对数据矩阵X进行按列标准化。标准化的目的是为了消除量纲和数量级大小的影响。这里采用的是零均值标准差标准化的方法。

3、求出样本间的相似系数矩阵Q，相似系数是把第i特征和第j特征的两个观测n次测试值，看成是n维空间中的两个点，坐标原点到此两个点的两个n维空间向量间夹角的余弦值，计算公式如下：

4、求相似系数矩阵Q的特征值和相应的特征向量，假设Q的n个特征值为λ₁,λ₂,.........λ_n，且λ₁≥λ₂≥…≥λ_n，对应的单位正交特征向量为μ₁,μ₂,…,μ_n。如果要在n个样本中选择m个样本(m<n)来代表这组样本，其他样本由这m个样本线性组合。各个特征值的大小代表各个因子的方差在总方差中所占的比例。η_m为累计贡献率，

累计贡献率越大，说明选出的主因素越拟合于原训练样本的数据，本专利选择的累计贡献率为99％。这m个公因子可以看作是m个独立的且不同于原始观测样本的“典型公共样本”。

5、计算因子载荷矩阵A，其元素a_ij的计算公式为

它表示第i个样本与第j个公因子的相似系数，其绝对值大小决定了样本的相对重要性，进而决定是否要把该样本选择为典型样本。|a_ij|越大就表明第i个样本所具有的载荷越大，第i个样本就越重要，因此第i个样本就作为公因子Fj的代表性样本，这就是训练样本选择的依据。

6、因子旋转

当因子载荷矩阵某一列上的各元素的绝对值差距比较大，并且绝对值大的元素较少时，则该公共因子就容易解释，反之，公共因子就的解释就会比较困难。此时可需要考虑对因子和因子载荷进行旋转，使得旋转后的因子载荷阵的各列元素的绝对值尽可能两极分化即更趋于1或0，这样公共因子就的解释就会比较容易。因子旋转方法有正交旋转和斜交旋转两种，这里我们采用方差最大正交旋转。正交旋转是指坐标轴始终保持90度角旋转，于是新生成的因子仍然可保持相关性。用公式表示如下

方差最大正交旋转的具体步骤如下：

①计算公共因子方差

②将因子载荷矩阵A标准化。

b_ij＝a_ij/h_i (7)

③根据式(6)计算方差V_(k),其中k为迭代的次数。

④任取因子载荷矩阵A中的两列设为a_(ig)和a_(iq)(i＝1,…,n),旋转后的矩阵为

式中，k为迭代的次数，j为样本编号，g(g＝1,…,m-1)和q(q＝g+1,…,m)为因子编号。

⑤将式(8)代入式(5),对求一阶导数,并令其等于零,解出关于角的方程

⑥将上述的旋转过程用下式表示为

B＝AT_gq (10)

式中,T_gq为相应于式(9)的旋转矩阵。

对因子载荷矩阵A中的所有两两成对的因子都进行上述类似的轴旋转计

算。如果有m个主因子,必须对A中所有m列全部配对旋转,共有次,全部旋转完毕算一个循环,此时得到的因子载荷矩阵为

其中,记B₍₁₎为对A进行正交变换C₁而得。经过第一个循环后,可按式(6)计算得V₍₁₎。在第一个循环的基础上从B₍₁₎出发再进行第二个旋转循环,旋转完毕后得B₍₂₎,再按式计算得到V₍₂₎。如此不断重复这个循环,就可得到V值的一个非降序列:

V₍₁₎≤V₍₂₎≤V₍₃₎≤… (12)

由于因子载荷的绝对值不大于1,故这个序列是有上界的,它必然收敛于某一极限V。

⑦确定轴旋转收敛误差：

在旋转过程中,若因子载荷矩阵A的总方差V值的两次计算之差的绝对值小于某个阈值,则停止旋转。这表明找到了旋转轴的最佳位置,否则继续旋转。这个过程可用下式来表示

|V_(k)-V_(k-1)|<ε (13)

式中一般取ε＝10^-7；

⑧计算旋转后的因子载荷矩阵G。

G＝(g_ij),g_ij＝b_ij×h_i (14)

式中，i＝1,…,n,j＝1,…,m。

7、根据上述的因子载荷矩阵G，从原观测样本中选择出具有代表性的样本；把所有选择出的样本作为训练样本进行训练。

步骤四、使用BP神经网络预测骨骼图像，建立肺部X光图像与DES肋骨图像之间的预测模型；

1.基于人工神经网络的骨骼图像预测及软组织图像生成

人工神经网络，也简称为神经网络，是一种应用类似于大脑神经突触连接的结构进行信息处理的数学模型。神经网络所具有的非线性特性、大量的并行分布结构以及学习和归纳能力使其在诸如建模、时间序列分析、模式识别、信号处理以及控制等方面得到了广泛的应用。

使用BP神经网络，它是目前应用最为广泛的神经网络之一。它是一种多层前馈型神经网络，它可以实现从输入到输出的任意非线性映射。其基本思想：学习过程由信号的正向传播与误差反向传播两个过程组成。正向传播时，输入样本从输入层传入，经隐含层逐层处理后传向输出层。若输出层的实际输出与期望输出不符，则转向误差的方向传播阶段。误差的反向传播是将输出误差以某种形式通过隐含层向输入层逐层反传，并将误差分摊给各层的所有单元，从而获得各层单元的误差信号，此误差信号即作为修正各单元权值的依据。这种信号正向传播与误差反向传播的各层权值调整周而复始地进行。权值不断调整的过程也是神经网络的学习训练过程。这个过程一直进行到网络预先设定的学习次数，或满足误差为止。

使用BP神经网络建立肺部X光图像与DES肋骨图像之间的预测模型。BP模型可以实现多层网络的学习的设想，是典型的有教师学习。以简单的3层BP网络来说，它由4个过程组成。对于给定的一组训练模式，不断用一个个训练模式网络重复前向传播和误差反响传播过程，当个训练模式都满足要求时，我们说BP网络已经学习好了。

假设由输入层、隐含层和输出层构成的三层BP神经网络，输入节点x_i，隐含层节点y_j，输出节点z_l。输入节点与隐含层节点间的网络权值为w_ji，隐含层节点与输出节点的网络权值为O_lj，当输出节点的期望输出为t_l时，BP模型的计算公式为：

隐含层检点的输出

其中

输出节点的计算输出

其中

输出节点的误差公式

建立预测模型后，可以根据肺部x光图像的灰度值分布来预测产生骨结构图像。在专利中，在肺部x光图像中减去预测获得的骨骼图像就可以获得软组织图像。

步骤五、对正常胸片与预测的骨骼图像做图像减法来预测软组织图像，建立预测模型后，可以根据肺部x光图像的灰度值分布来预测产生骨结构图像在肺部x光图像中减去预测获得的骨骼图像就可以获得软组织图像。

Claims (5)

1.一种基于训练样本优化的骨骼抑制方法，其步骤如下：

步骤一：根据图像灰度与形状信息对肺区进行分割，即建立肺轮廓的先验模型，利用灰度与形状相似性信息并结合图像特征对肺区域分割；

步骤二：采用B样条小波变换和N-jet进行特征提取，提取有效描述不同尺度骨骼的特征图像；

步骤三：使用Q型因子分析法对样本进行处理，其过程如下：

1、采用Q型因子分析对大样本数据进行预处理，提取特征图像的观测数据矩阵X；

2、对数据矩阵X进行按列标准化，采用零均值标准差标准化的方法；

3、求出样本间的相似系数矩阵Q；

记第i个样本和第j个样本的观测p次测试值为n维空间中的两个点，则坐标原点分别到这两个点的n维空间向量间夹角的余弦值为相似系数，其计算公式如下：

$q_{ij}=\frac{\underset{k=1}{\overset{p}{\&Sigma;}}{x}_{ki}{x}_{kj}}{\sqrt{\underset{k=1}{\overset{p}{\&Sigma;}}{x}_{ki}^2\underset{k=1}{\overset{p}{\&Sigma;}}{x}_{kj}^2}},(i,j=1,2,...,n)$

其中x_ki为第i个样本的第k个特征值，x_kj为第j个样本的第k个特征值；所有相似系数构成的矩阵为相似系数矩阵；

4、求相似系数矩阵Q的特征值和相应的特征向量；

假设Q的n个特征值为λ₁,λ₂,.........λ_n，且λ₁≥λ₂≥…≥λ_n，对应的单位正交特征向量为μ₁,μ₂,…,μ_n，在n个样本中选择m个样本(m<n)来代表这组观测样本，n个样本中除m个样本外的样本由选定的m个样本线性组合表示；各个特征值的大小代表各个因子的方差在总方差中所占的比例；η_m为累计贡献率，

${\mathrm{\&eta;}}_m=\frac{\underset{i=1}{\overset{m}{\&Sigma;}}{\mathrm{\&lambda;}}_i}{\underset{i=1}{\overset{n}{\&Sigma;}}{\mathrm{\&lambda;}}_i}$

5、初始因子载荷矩阵A；

其中μ_ij是λ_j对应的第i个分量，且

a_ij表示第i个样本与第j个公因子的相似系数；

6、因子旋转；采用方差最大正交旋转，所述正交旋转是指坐标轴始终保持90度角旋转，得到因子载荷矩阵G，其表示为：

G＝(g_ij),g_ij＝b_ij×h_i

式中，i＝1,…,n,j＝1,…,m；b_ij＝a_ij/h_i；

7、根据因子载荷矩阵G，从原始的样本中按照比例选出具有代表性的样本集；

步骤四：使用BP神经网络预测骨骼图像，建立肺部X光图像的特征图像与DES肋骨图像之间的预测模型；

假设由输入层、隐含层和输出层构成的三层BP神经网络，输入节点x_i，隐含层节点y_j，输出节点z_l；输入节点与隐含层节点间的网络权值为w_ji，隐含层节点与输出节点的网络权值为O_lj，当输出节点的期望输出为t_l时，BP模型的计算公式为：

隐含层检点的输出

$y_j=f(\underset{i}{\&Sigma;}{w}_{ji}{x}_i-{\mathrm{\&theta;}}_j)=f\left({\mathrm{net}}_j\right)---\left(7\right)$

其中，θ_j为隐含层的偏置；

输出节点的计算输出

$z_l=f(\underset{j}{\&Sigma;}{o}_{lj}{y}_j-{\mathrm{\&theta;}}_l)=f\left({\mathrm{net}}_l\right)$

其中，θ_l为输出层的偏置；

输出节点的误差公式

$E=\frac{1}{2}\underset{l}{\&Sigma;}{(t_l-z_l)}^2=\frac{1}{2}\underset{l}{\&Sigma;}{(t_l-f(\underset{j}{\&Sigma;}{o}_{lj}{y}_j-{\mathrm{\&theta;}}_l\left)\right)}^2=\frac{1}{2}\underset{l}{\&Sigma;}{(t_l-f(\underset{j}{\&Sigma;}{o}_{lj}f\left(\underset{i}{\&Sigma;}{w}_{ji}{x}_i-{\mathrm{\&theta;}}_j\right)-{\mathrm{\&theta;}}_l\left)\right)}^2$

步骤五、对正常胸片与预测的骨骼图像做图像减法来预测软组织图像，建立预测模型后，可以根据肺部x光图像的灰度值分布来预测产生骨结构图像在肺部x光图像中减去预测获得的骨骼图像就可以获得软组织图像。

2.根据权利要求1所述的基于训练样本优化的骨骼抑制方法，其特征在于，训练样本包括n张胸部x光图像。

3.根据权利要求1所述的基于训练样本优化的骨骼抑制方法，其特征在于，所述B样条小波变换提取为三阶B样条多尺度小波变换提取，N-jet提取为多尺度2-jet提取。

4.根据权利要求1所述的基于训练样本优化的骨骼抑制方法，其特征在于，所述特征图像变换提取过程为对胸部x光图像分别做一、二、三尺度小波分解，通过一尺度小波解获得4张图像，1张近似图像，1张水平细节图像，1张垂直细节图像，1张对角线细节图像；通过二尺度小波解得到不同于一尺度的1张近似图像，1张水平细节图像，1张垂直细节图像，1张对角线细节图像；通过三尺度小波解得到不同于二尺度的1张近似图像，1张水平细节图像，1张垂直细节图像，1张对角线细节图像；去掉对角线细节图像，得到9张特征图像；将得到的9张特征图像再进行高斯滤波局部N-jet的特征提取；根据高斯滤波尺度的不同，每张小波变换后的图像再得到12张特征图像。

5.根据权利要求1所述的基于训练样本优化的骨骼抑制方法，其特征在于，所述方差最大正交旋转的公式如下：

$V=\underset{j=1}{\overset{j=m}{\&Sigma;}}\left\{\&lsqb;n\underset{i=1}{\overset{n}{\&Sigma;}}{\left(\frac{a_{ij}^2}{h_i^2}\right)}^2-{\left(\underset{i=1}{\overset{n}{\&Sigma;}}\frac{a_{ij}^2}{h_i^2}\right)}^2\&rsqb;/n^2\right\}=\max$

其中：a_ij表示第i个样本与第j个公因子的相似系数。