CN103679701B - 基于支持向量机回归的晶体图像直线轮廓检测方法 - Google Patents

Abstract

本发明公开了一种基于支持向量机回归的晶体图像直线轮廓检测方法，包括采用CCD照相机采集单晶硅生长过程中的边缘轮廓线图像，对该边缘轮廓线图像进行预处理，构建关于直线角度的过完备字典，用LS‑SVR的对偶优化模型求解出与直线角度所对应的稀疏表示系数，求出直线的角度；引入稀疏约束项来优化原始对偶优化模型以提高直线角度的精确度；构建一个直线偏移量的过完备字典，再根据对偶优化模型解出直线的偏移量。本发明能在样本点很少的情况下拟合直线，很好的解决高维数据，在最小二乘支持向量机的原始对偶优化模型中加入了稀疏约束项，进一步提高了算法的鲁棒性，能够准确的估计晶体生长中的中心轴的变化情况，并控制生长中的晶体中心轴位置。

Description

基于支持向量机回归的晶体图像直线轮廓检测方法

技术领域

本发明属于测量技术领域，具体涉及到一种基于支持向量机回归的晶体图像直线轮廓检测方法。

背景技术

在单晶硅生长的直径控制系统中，晶体直径的控制是保证单晶硅能够等径生长的关键步骤，拉单晶的过程始终保持在高温负压的环境中，直径检测必须隔着观察窗在拉晶炉体外部非接触式实现，如图1所示观测模型，拉晶过程中固态晶体与液态溶液的交接处会形成一个明亮的光环，一般通过调整这个光环直径的大小来控制晶体的等径生长。然而，只控制等径生长是不够的，晶柱体可以看作是不同时间段的许多晶体薄片叠加形成的，如图2所示，为了确保晶体呈圆柱型生长，那么除了要求晶体薄片具有相同的直径外，还应该约束所有薄片都具有相同的圆心，如果只控制等径生长就有可能出现等径不等心(每个圆切片不是在一个中心轴上)的情况发生，晶体就有可能弯曲生长，如图3所示模型，严重的情况下，弯曲的单晶会触及炉壁而导致单晶炉坍塌。所以，为了确保安全生产以及生产出高质量的单晶材料，在晶体生长过程中应该同时约束每个薄片等径等轴生长。在拉晶的过程中晶体的两侧直线边缘正好能反映出晶体生长中心轴的控制状态，因此可以通过控制晶体生长边缘的变化趋势来调节拉杆的位置与方向，以便最终生成的晶体接近完美的圆柱形。所以当前晶柱体直线边缘的测量在单晶硅生长控制系统中起着至关重要的作用，检测晶体直线边缘的变化趋势也等价于检测到了晶体中心轴的变化趋势，这样晶体中心轴的估计问题就可以转化为直线检测问题。

对于晶体中心轴控制的先进、可靠的方法是非常重要的。它关系到能否及早发现问题，以降低不必要的损失。直线拟合问题可以描述为：给定一个二值图像，其中有一些像素值为1的点离散的分布在像素值为0背景上，直线检测的目标就是寻找拟合这些共线的像素值为1的点的最佳直线，并且估计出直线的角度以及偏移量。近年来也有人提出几种经典的直线检测算法，《Use of the Hough transform to detect lines and curves inpictures》(Comm.ACM.vol.15,pp.11-15,1972)提到了霍夫曼变换(Hough)算法，Hough变换的基本思想是点线的对偶性，由于图像空间中共线的点对应在参数空间里相交的线，在参数空间中相交于同一个点的所有直线在图像空间里都有共线的点与之对应，Hough变换是对图像中的每个点应用拉东变换(Radon transform)，然后在角度-偏移量二维空间去搜索得到一个最大值点作为最终要估计的直线参数，因为在二维空间里要对所有的参数进行全局搜索，所以这种方法往往是非常耗时的。《Sensor array processing techniques forsuper resolution multi-line-fitting and straight edge detection》(IEEETrans.Image Processing.vol.2,no.4,pp.454-465,1993)和《SLIDE:subspace-basedline detection.IEEE Trans.Pattern Analysis and Machine Intelligence》(PatternAnalysis and Machine Intelligence.vol.16,no.11,pp.1057-1073,1994)涉及基于子空间的直线检测(Subspace-based Line DEtection),它是对图像中的每条直线与一个传感器阵列上的波阵面辐射之间做类比从而估计出直线的参数，这种方法在实际应用中的估计误差比较大，效果也往往不太理想。

《Support vector machines framework for linear signal processing》(Signal Processing.vol.85,no.12,pp.2316-2326,2005)和《Least squares supportvector machine classifiers》(Neural Process.Lett.,vol.9,no.3,pp.293-300,1999)涉及支持向量机(Support Vector Machine)，它是根据统计学理论中的最小化原则提出来的，由有限数据得到判别函数，对独立的测试样本能够得到较小的误差。

发明内容

本发明的目的是提供一种基于支持向量机回归的晶体生长控制方法，解决了现有技术中由于晶体中心轴变化不稳定而导致最终晶体质量不高的问题，对晶体生长的过程进行实时监控，达到了使晶体稳定等径生长的目的。

本发明所采用的技术方案是：一种基于支持向量机回归的晶体图像直线轮廓检测方法，具体按照以下方式实施：

步骤1、首先采用CCD照相机采集单晶硅生长过程中的边缘轮廓线图像，对该边缘轮廓线图像进行预处理，得到用于估计的采样点；

步骤2、基于步骤1中的采样点，构建关于直线角度的过完备字典然后再用最小二乘支持向量机回归LS-SVR模型求解出与直线角度所对应的稀疏表示系数，再由稀疏表示系数求出直线的角度；

步骤3、对步骤2中的LS-SVR的对偶优化模型引入稀疏约束项来优化原始对偶优化模型以提高直线角度的精确度；

步骤4、基于步骤3得到的直线角度构建一个直线偏移量的过完备字典，然后再根据步骤3中的对偶优化模型解出直线的偏移量。

预处理采用基于分块和寻找特征区域的方法，用行扫描线方法对晶体图像的两侧边缘的样本点进行采样，得到用于训练的样本数据,其中，图像预处理得到大小为M×M的二值图像D，灰度为“1”的点表示直线点，而灰度为“0”的点表示背景。

直线的角度的求解过程具体按照以下步骤实施：

CCD照相机采集单晶硅生长过程中的边缘轮廓线图像，边缘轮廓线是一条直线，图像左上角像素视为坐标原点，水平向右的方向为x轴正向，垂直向下方向为y轴正向，则直线上的点的坐标{x,y}，直线的偏移量以倾斜角θ唯一确定，满足下面关系：

$x=ytan\mathrm{\θ}+\tilde{x}---\left(1\right)$

由公式(1)知，估计图像中包含多条直线的偏移量和角度(θ₁,θ₂,…,θ_d)，检测直线实质上是确定直线参数和θ；

定义向量u：

u＝[1 e^jμ e^j2μ…e^j(M-1)μ]^T ，(2)

对二进制图像D进行转换，将二维数据转换为一维数据，即

$\begin{array}{c}z=Du\\ ={\[\underset{i=1}{\overset{d}{\Σ}}{e}^{j\mathrm{\μ}(0\×{\mathrm{tan\θ}}_i+{\tilde{x}}_i)}\underset{i=1}{\overset{d}{\Σ}}{e}^{j\mathrm{\μ}(1\×{\mathrm{tan\θ}}_i+{\tilde{x}}_i)}\mathrm{...}\underset{i=1}{\overset{d}{\Σ}}{e}^{j\mathrm{\μ}\left(\right(M-1)1\×{\mathrm{tan\θ}}_i+{\tilde{x}}_i)}\]}^T\\ =\underset{i=1}{\overset{d}{\Σ}}a\left({\mathrm{\θ}}_i\right)s_i\end{array},---\left(3\right)$

其中，μ是一个常量或者变量，标量向量

$a\left({\mathrm{\θ}}_i\right)={\left[\begin{array}{ccccc}1& e^{{\mathrm{j\μ}\tan \mathrm{\θ}}_i}& e^{j2\mathrm{\μ}\tan {\mathrm{\θ}}_i}& ...& e^{j(M-1)\mathrm{\μ}{\tan \mathrm{\θ}}_i}\end{array}\right]}^T,---\left(4\right)$

经过上述转换，二进制图像D转换为由M个传感器组成的均匀线性阵列的虚拟快拍信号z,其中第i个信号幅度为s_i，第i个信号即第i条直线为平面波，平滑移动到图像第1列像素的左侧由实心原点表示的虚拟均匀线性阵列上，其中，向量u中的μ满足μ|tanθ_i|≤π条件；

为了从单快拍z中估计θ_i，建立起来一个关于角度估计的稀疏表示模型，由于这d个信号的入射角，其中d个信号的入射角即为d个直线的倾斜角，取值范围为：[-90°,90°]，为此，将[-90°,90°]细分为N_θ个均匀空间，基于上述均匀间隔化得到的角度集基于公式(4)描述的信号导向矢量结构以及这N_θ个可能入射方向，构造如下的虚拟导向矢量过完备字典：

$A=\left[\begin{array}{cccc}a\left({\widetilde{\mathrm{\θ}}}_1\right)& a\left({\widetilde{\mathrm{\θ}}}_2\right)& \mathrm{...}& a\left({\widetilde{\mathrm{\θ}}}_{N_{\mathrm{\θ}}}\right)\end{array}\right]\∈C^{M\×N_{\mathrm{\θ}}},---\left(5\right)$

该过完备字典由N_θ列组成，每列具有和公式(4)相同的结构，这样虚拟均匀线性阵列接受信号z可借助这N_θ列重新表示为：

z＝As+n (6)

其中与N_θ个潜在的入射角对应的信号幅度为n为对应的噪声，当s的第i个元素对应于第d个入射角之一时，则取非零值，否则，其取值为0；

公式(6)可以写成另外一种实数形式：

$\tilde{z}=\tilde{A}\tilde{s}+\tilde{n}---\left(7\right)$

其中

其中real(·)和imag(·)分别表示·的实部与虚部，事实上，公式(7)可以看作是一个虚拟线性系统或是一个线性回归问题，其中A(i,:)表示虚拟输入，z(i)表示相关输出i＝1,2,…,2M，SVR来解下面的线性回归问题得到稀疏向量

$\begin{array}{cc}\underset{\tilde{s},b}{\min }& \frac{1}{2}{\tilde{s}}^{T}s+\frac{C}{2}\underset{i=1}{\overset{2M}{\Σ}}{e}^2\left(i\right)\end{array}$

$\begin{array}{cc}s.t.& \tilde{z}\left(i\right)=\tilde{A}(i,:)\tilde{s}+b+e\left(i\right),i=1,2,\mathrm{...},2M,\end{array}---\left(8\right)$

其中b偏移量，e(i)为误差项，引入拉格朗日乘子α_i,i＝1,2,…,2M,公式(8)的对偶形式为

$\begin{array}{c}L(\tilde{s},b,\tilde{e},\mathrm{\α})=\frac{1}{2}{\tilde{s}}^T\tilde{s}+\frac{C}{2}\underset{i=1}{\overset{2M}{\Σ}}{\tilde{e}}^2\left(i\right)\\ +\underset{i=1}{\overset{2M}{\Σ}}{\mathrm{\α}}_i\{\tilde{z}\left(i\right)-\tilde{A}(i,:)\tilde{s}-b-\tilde{e}\left(i\right)\}\end{array}---\left(9\right)$

优化条件如下：

$\left\{\begin{array}{c}\frac{\∂L}{\∂\tilde{s}}=0\→\tilde{s}=\underset{i=1}{\overset{2M}{\Σ}}{\mathrm{\α}}_i{\tilde{A}}^T(i,:)\\ \frac{\∂L}{\∂b}=0\→\underset{i=1}{\overset{2M}{\Σ}}{\mathrm{\α}}_i=0,\\ \frac{\∂L}{\∂\tilde{e}\left(i\right)}=0\→{\mathrm{\α}}_i=C\tilde{e}\left(i\right),i=1,\mathrm{...2}M,\\ \frac{\∂L}{\∂{\mathrm{\α}}_i}=0\→\tilde{z}\left(i\right)-\tilde{A}(i,:)\tilde{s}-b-\tilde{e}\left(i\right)=0,i=1,\mathrm{...2}M\end{array}\right.---\left(10\right),$

根据KKT(Karush-Kuhn-Tucker)条件，公式(9)中描述的问题简化成如下线性系统

$\left[\begin{array}{cc}Q+\frac{1}{C}{I}_{2M}& 1_{2M\×1}\\ 1_{1\×2M}& 0\end{array}\right]\left[\begin{array}{c}\mathrm{\α}\\ b\end{array}\right]=\left[\begin{array}{c}\tilde{z}\\ 0\end{array}\right]---\left(11\right)$

其中α＝[α₁ α₂…α_2M]^T，1_2M×1表示一个2M×1维向量，每个元素都是1，只要α和b由式(11)解出来，向量由如下的公式计算出来：

$\tilde{s}=\underset{i=1}{\overset{2M}{\Σ}}{\mathrm{\α}}_i{\tilde{A}}^T(i,:)={\tilde{A}}^T\mathrm{\α}---\left(12\right)$

基于计算出阵列信号的虚拟谱其中第i个元素可以表示为

$p_{spec}\left(i\right)=\sqrt{{\tilde{s}}^2\left(i\right)+{\tilde{s}}^2(i+N_{\mathrm{\θ}})}---\left(13\right)$

然后对虚拟谱中的元素从大到小进行降序排序，定义sort为降序排序操作，Angle＝sort{Pspec}，谱线中d个最高峰的位置对应估计的角度值

$({\widehat{\mathrm{\θ}}}_1,{\widehat{\mathrm{\θ}}}_2,\mathrm{...},{\widehat{\mathrm{\θ}}}_d)=Angle(1,2,\mathrm{...}d)---\left(14\right).$

LS-SVR的对偶优化模型引入稀疏约束项来优化原始对偶优化模型以提高直线角度的精确度具体按照如下步骤实施：

公式(11)等价于如下所示的最小化问题：

$\begin{array}{cc}\underset{\mathrm{\α},b}{\min }& \left|\right|\tilde{z}-\left[\begin{array}{cc}Q+\frac{1}{C}{I}_{2M}& 1_{2M\×1}\end{array}\right]\left[\begin{array}{c}\mathrm{\α}\\ b\end{array}\right]|{|}_2^2\\ s.t.& 1_{1\×2M}\mathrm{\α}=0\end{array}---\left(15\right)$

在公式(15)引入下面关于的稀疏约束项：

$\begin{array}{cc}\underset{\mathrm{\α},b}{\min }& \left|\right|\tilde{z}-\left[\begin{array}{cc}Q+\frac{1}{C}{I}_{2M}& 1_{2M\×1}\end{array}\right]\left[\begin{array}{c}\mathrm{\α}\\ b\end{array}\right]|{|}_2^2\\ & +\mathrm{\λ}\underset{i=1}{\overset{N_{\mathrm{\θ}}}{\Σ}}\sqrt{{\left({\tilde{A}}^T\right(:,i\left)\mathrm{\α}\right)}^2+{\left({\tilde{A}}^T\right(:,i+N_{\mathrm{\θ}}\left)\mathrm{\α}\right)}^2}\\ s.t.& 1_{1\×2M}\mathrm{\α}=0\end{array}---\left(16\right)$

其中参数λ是稀疏向量和误差范数之间的折中；

引入N_θ+1个新变量把公式(16)改写成另外一种形式，

$\begin{array}{c}\underset{\mathrm{\η},b,\mathrm{\α},{\mathrm{\γ}}_1,\mathrm{...},{\mathrm{\γ}}_{N_{\mathrm{\θ}}}}{\min }\mathrm{\η}+\mathrm{\λ}\underset{i=1}{\overset{N_{\mathrm{\θ}}}{\Σ}}{\mathrm{\γ}}_i\\ \begin{array}{cc}s.t.& \left|\right|\tilde{z}-\left[\begin{array}{cc}Q+\frac{1}{C}{I}_{2M}& 1_{2M\×1}\end{array}\right]\left[\begin{array}{c}\mathrm{\α}\\ b\end{array}\right]|{|}_2^2\≤\mathrm{\η},\end{array}\\ \left|\right|\left[\begin{array}{c}{\tilde{A}}^T(:,i)\\ {\tilde{A}}^T(:,i+N_{\mathrm{\θ}})\end{array}\right]\mathrm{\α}\left|\right|\≤{\mathrm{\γ}}_i,i=1,2,\mathrm{...}{N}_{\mathrm{\θ}},\\ 1_{1\×2M}\mathrm{\α}=0\end{array}---\left(17\right)$

将公式(17)推导成一个标准的二阶锥规划问题，由Strum开发的用于处理对称锥优化问题的MATLAB工具来解决。

偏移量的计算具体按照如下步骤实施：

假设与角度所对应的潜在的的偏移量是基于和生成一个大小为M×M的二值图像B_i,j，其中只包含一条偏移量和角度分别为和的直线；

基于B_i,j，(i,j)∈[(1,1),…,(1,N₁),(2,1),…,(2,N₂),…,(d,N_d)]，构建另一个有个原子的偏移量过完备字典B，即

$B=\left[\begin{array}{cccc}vec\left(B_{1,1}\right)& vec\left(B_{1,2}\right)& \mathrm{...}& vec\left(B_{d,N_d}\right)\end{array}\right]\∈R^{M^2\×\underset{i=1}{\overset{d}{\Σ}}{N}_i}---\left(18\right)$

其中vec(B_i,j)代表B_i,j的向量形式；

显然，vec(D)表示成vec(B_i,j)的线性组合，

vec(D)＝Bv+n (19)

其中是v是稀疏的线性组合系数，n是噪声，由于即B的行数非常大，定义一个随机映射矩阵去乘以公式(19)的两边来得到一个等价的方程：

$\tilde{D}=\tilde{B}v+\tilde{n},---\left(20\right)$

其中，q＜＜M²,

应用(8)-(12)及(15)-(17)将公式(20)化成如公式(17)所示的二阶锥规划问题，再由Strum开发的用于处理对称锥优化问题的MATLAB工具来解出最优解v，然后对v的绝对值进行降序排序，即Offset＝sort{abs(v)}，v中绝对值最大的d个峰值就代表检测出来的偏移量即：

$({\tilde{x}}_1,{\tilde{x}}_2,\mathrm{...},{\tilde{x}}_d)=Offset(1,2,\mathrm{...},d)---\left(21\right);$

将公式(17)推导成一个标准的二阶锥规划问题，由Strum开发的用于处理对称锥优化问题的MATLAB工具来解决，具体按照以下步骤实施：

在SeDuMi中，标准的优化问题形式定义为：

$\begin{array}{ccc}\underset{z}{\max }& p^{T}h& \\ Subjectto& r_j-q_j^{T}h\∈{\mathrm{SOC}}^{g_j\×l}& j=1,2,\mathrm{...}J\end{array}---\left(22\right)$

其中的p和r_j是任意的向量，q_j是任意的矩阵，h中包含有期望优化的变量，J是二阶锥约束的个数，g_j维的约束定义为：

${\mathrm{SOC}}^{g_j\×l}=\left\{\right|\left|\mathrm{\ϵ}\right||\≤\widetilde{\mathrm{\ϵ}}\}---\left(23\right)$

在这里是g_j维向量中的第一个元素，ε是g_j-1维的向量，包含了中的其他元素，根据式(17)定义各个变量如下所示：

$\begin{array}{c}\underset{\mathrm{\η},b,\mathrm{\α},{\mathrm{\γ}}_1,\mathrm{...},{\mathrm{\γ}}_{N_{\mathrm{\θ}}}}{\min }\mathrm{\η}+\mathrm{\λ}\underset{i=1}{\overset{N_{\mathrm{\θ}}}{\Σ}}{\mathrm{\γ}}_i\\ \begin{array}{cc}s.t.& \left|\right|\tilde{z}-\left[\begin{array}{cc}Q+\frac{1}{C}{I}_{2M}& 1_{2M\×1}\end{array}\right]\left[\begin{array}{c}\mathrm{\α}\\ b\end{array}\right]|{|}_2^2\≤\mathrm{\η},\end{array}\\ \left|\right|\left[\begin{array}{c}{\tilde{A}}^T(:,i)\\ {\tilde{A}}^T(:,i+N_{\mathrm{\θ}})\end{array}\right]\mathrm{\α}\left|\right|\≤{\mathrm{\γ}}_i,i=1,2,\mathrm{...}{N}_{\mathrm{\θ}},\\ 1_{1\×2M}\mathrm{\α}=0\end{array}---\left(17\right)$

其中的定义：

$p=\left[\begin{array}{cccc}-1& 0_{1\×2M}& 0& -{\mathrm{\λI}}_{1\×N_{\mathrm{\θ}}}\end{array}\right],---\left(24.2\right)$

$r_1={\left[\begin{array}{cccc}\frac{1}{2}& \frac{1}{2}& 0& \tilde{z}\end{array}\right]}^T,---\left(24.3\right)$

r₂＝[0 0]^T， (24.5)

r_2+i＝[0 0]^T，i＝1,2,…,N_θ， (24.7)

将式(17)转换成如式(22)的标准的SOCP形式，利用SOCP优化工具包得到很好的解决，

如果解出来公式(17)、(24.1)至(24.8)是角度的计算方法，计算出的拉格朗日乘子为α_角度，与角度所对应的稀疏频谱是直线对应的角度公式(14)计算出来。

本发明的有益效果是，能在样本点很少的情况下进行很好的直线拟合，在解决高维数据时发挥了很大的优势，所以能够准确的估计晶体生长中的中心轴的变化情况，使生长中的晶体中心轴位置得到很好的控制，这是其他方法所不能达到的。

附图说明

图1是本发明方法用于孔径图像提取的CCD照相机系统安装图；

图2是本发明方法中晶体等径生长模型；

图3是本发明方法中晶体弯曲生长模型；

图4是本发明方法中用于描述算法的直线模型；

图5是本发明方法利用CCD照相机获得的单晶图像；

图6是本发明方法采用图像预处理技术得到的二值图像；

图7是本发明方法估计的晶体图像直线边缘所对应的角度频谱图；

图8是本发明方法估计的晶体图像直线边缘所对应的偏移量频谱图；

图9是晶体图像直线边缘的实际拟合效果图；

图10是本发明方法小角度仿真实验图像；

图11是本发明方法小角度实验所对应的角度频谱图；

图12是本发明方法小角度实验所对应的偏移量频谱图。

具体实施方式

下面结合附图和理论推导对本发明方法进行详细说明。

本发明提出了一种基于支持向量机回归的晶体图像直线轮廓检测方法，该方法将晶体中心轴控制问题转化成直线检测问题，在采取数字图像处理技术提取晶体生长边缘特征点的基础上，首先将直线检测问题划分为两个步骤，即角度估计问题与偏移量估计问题。在角度估计步骤中，首先针对这些特征采样点形成均匀线性阵列的虚拟快拍数据，接着，基于可能的虚拟入射角(本发明中采取[-90°,90°])构造的相应的远场虚拟阵列流行，然后将直线的倾斜角估计问题转换为均匀线性阵列单快拍信号入射角的估计问题。然后分别针对直线的角度估计问题和偏移量估计问题建立支持向量机回归模型，并且引入了稀疏表示约束项确保估计参数的准确性，具体按照以下方式实施：

步骤1、首先采用CCD照相机采集单晶硅生长过程中的边缘轮廓线图像，对该边缘轮廓线图像进行预处理，得到用于估计的采样点；

具体按照以下方式实施：

采用CCD照相机采集单晶硅生长过程中的晶体图像，然后对图像进行预处理得到用于估计的采样点，预处理采用基于分块和寻找特征区域的方法，用行扫描线方法对晶体图像的两侧边缘的样本点进行采样，得到用于训练的样本数据，其中，为了便于描述，假设图像预处理得到大小为M×M的二值图像D，灰度为“1”的点表示直线点，而灰度为“0”的点表示背景。

步骤2、基于步骤1中的采样点，针对直线的标准模型构建关于直线角度的过完备字典然后再用最小二乘支持向量机回归(Least squares support vector regression，LS-SVR)的对偶优化模型求解出与直线角度所对应的稀疏表示系数，再由稀疏表示系数求出直线的角度，具体按照以下步骤实施：

CCD照相机采集单晶硅生长过程中的边缘轮廓线图像，边缘轮廓线是一条直线，以图4所示的单条直线为例，直线的偏移量以倾斜角θ唯一确定，图像左上角像素视为坐标原点，水平向右的方向为x轴正向，垂直向下方向为y轴正向，则直线上的点的坐标{x,y}，满足下面关系：

$x=ytan\mathrm{\θ}+\tilde{x}---\left(1\right)$

由公式(1)知，估计图像中包含多条直线的偏移量和角度(θ₁,θ₂,…,θ_d)，检测直线实质上是确定直线参数和θ；

基于上述二值图像中的直线形成均匀线性阵列的虚拟快拍数据，接着，基于可能的虚拟入射角构造相应的远场虚拟阵列模型，将直线的倾斜角估计问题转换为均匀线性阵列单快拍信号入射角的估计问题，基于这样的虚拟阵列流行以及潜在的入射角范围建立一个过完备字典形成一个基于稀疏表示的回归问题,本发明考虑的回归模型是最小二乘支持向量机模型，最终求出稀疏表示系数，由稀疏表示系数求出信号的入射角度。

定义向量u：

u＝[1 e^jμ e^j2μ…e^j(M-1)μ]^T ，(2)

对二进制图像D进行转换，将二维数据转换为一维数据，即

$\begin{array}{c}z=Du\\ ={\[\underset{i=1}{\overset{d}{\Σ}}{e}^{j\mathrm{\μ}(0\×{\mathrm{tan\θ}}_i+{\tilde{x}}_i)}\underset{i=1}{\overset{d}{\Σ}}{e}^{j\mathrm{\μ}(1\×{\mathrm{tan\θ}}_i+{\tilde{x}}_i)}\mathrm{...}\underset{i=1}{\overset{d}{\Σ}}{e}^{j\mathrm{\μ}\left(\right(M-1)1\×{\mathrm{tan\θ}}_i+{\tilde{x}}_i)}\]}^T\\ =\underset{i=1}{\overset{d}{\Σ}}a\left({\mathrm{\θ}}_i\right)s_i\end{array},---\left(3\right)$

其中，μ是一个常量或者变量，标量向量

$a\left({\mathrm{\θ}}_i\right)={\left[\begin{array}{ccccc}1& e^{{\mathrm{j\μ}\tan \mathrm{\θ}}_i}& e^{j2\mathrm{\μ}\tan {\mathrm{\θ}}_i}& ...& e^{j(M-1)\mathrm{\μ}{\tan \mathrm{\θ}}_i}\end{array}\right]}^T,---\left(4\right)$

经过上述转换，二进制图像D转换为由M个传感器组成的均匀线性阵列的虚拟快拍信号z,其中第i个信号幅度为s_i，第i个信号(第i条直线)类似于平面波，辐射到(平滑移动到)图像第1列像素的左侧由实心原点表示的虚拟均匀线性阵列上。其中，向量u中的μ满足μ|tanθ_i|≤π条件；

为了从单快拍z中估计θ_i，我们建立起来一个关于角度估计的稀疏表示模型，由于这d个信号的入射角(即d个直线的倾斜角)，取值范围为：[-90°,90°]，为此，我们考虑将[-90°,90°]细分为N_θ个均匀空间，基于上述均匀间隔化得到的角度集基于公式(4)描述的信号导向矢量结构以及这N_θ个可能入射方向，构造如下的虚拟导向矢量过完备字典：

$A=\left[\begin{array}{cccc}a\left({\widetilde{\mathrm{\θ}}}_1\right)& a\left({\widetilde{\mathrm{\θ}}}_2\right)& \mathrm{...}& a\left({\widetilde{\mathrm{\θ}}}_{N_{\mathrm{\θ}}}\right)\end{array}\right]\∈C^{M\×N_{\mathrm{\θ}}},---\left(5\right)$

该过完备字典由N_θ列组成，每列具有和公式(4)相同的结构，这样虚拟均匀线性阵列接受信号z可借助这N_θ列重新表示为：

z＝As+n (6)

其中与N_θ个潜在的入射角对应的信号幅度为n为对应的噪声，当s的第i个元素对应于第d个入射角之一时，则取非零值，否则，其取值为0；

公式(6)可以写成另外一种实数形式：

$\tilde{z}=\tilde{A}\tilde{s}+\tilde{n}---\left(7\right)$

其中

其中real(·)和imag(·)分别表示·的实部与虚部。事实上，公式(7)可以看作是一个虚拟线性系统或是一个线性回归问题，其中A(i,:)表示虚拟输入，z(i)表示相关输出i＝1,2,…,2M，由于支持SVR可以找到具有最优的最大化间隔超平面的参数，在本发明中，我们应用SVR来解下面的线性回归问题得到稀疏向量

$\begin{array}{cc}\underset{\tilde{s},b}{\min }& \frac{1}{2}{\tilde{s}}^{T}s+\frac{C}{2}\underset{i=1}{\overset{2M}{\Σ}}{e}^2\left(i\right)\end{array}$

$\begin{array}{cc}s.t.& \tilde{z}\left(i\right)=\tilde{A}(i,:)\tilde{s}+b+e\left(i\right),i=1,2,\mathrm{...},2M,\end{array}---\left(8\right)$

其中b偏移量(不同于直线的偏移量)，e(i)为误差项。在这里我们引入拉格朗日乘子α_i,i＝1,2,…,2M,公式(8)的对偶形式为

$\begin{array}{c}L(\tilde{s},b,\tilde{e},\mathrm{\α})=\frac{1}{2}{\tilde{s}}^T\tilde{s}+\frac{C}{2}\underset{i=1}{\overset{2M}{\Σ}}{\tilde{e}}^2\left(i\right)\\ +\underset{i=1}{\overset{2M}{\Σ}}{\mathrm{\α}}_i\{\tilde{z}\left(i\right)-\tilde{A}(i,:)\tilde{s}-b-\tilde{e}\left(i\right)\}\end{array}---\left(9\right)$

优化条件如下：

$\left\{\begin{array}{c}\frac{\∂L}{\∂\tilde{s}}=0\→\tilde{s}=\underset{i=1}{\overset{2M}{\Σ}}{\mathrm{\α}}_i{\tilde{A}}^T(i,:)\\ \frac{\∂L}{\∂b}=0\→\underset{i=1}{\overset{2M}{\Σ}}{\mathrm{\α}}_i=0,\\ \frac{\∂L}{\∂\tilde{e}\left(i\right)}=0\→{\mathrm{\α}}_i=C\tilde{e}\left(i\right),i=1,\mathrm{...2}M,\\ \frac{\∂L}{\∂{\mathrm{\α}}_i}=0\→\tilde{z}\left(i\right)-\tilde{A}(i,:)\tilde{s}-b-\tilde{e}\left(i\right)=0,i=1,\mathrm{...2}M\end{array}\right.---\left(10\right),$

根据KKT(Karush-Kuhn-Tucker)条件，公式(9)中描述的问题简化成如下线性系统

$\left[\begin{array}{cc}Q+\frac{1}{C}{I}_{2M}& 1_{2M\×1}\\ 1_{1\×2M}& 0\end{array}\right]\left[\begin{array}{c}\mathrm{\α}\\ b\end{array}\right]=\left[\begin{array}{c}\tilde{z}\\ 0\end{array}\right]---\left(11\right)$

其中α＝[α₁ α₂…α_2M]^T，1_2M×1表示一个2M×1维向量，每个元素都是1，只要α和b由式(11)解出来，向量由如下的公式计算出来：

$\tilde{s}=\underset{i=1}{\overset{2M}{\Σ}}{\mathrm{\α}}_i{\tilde{A}}^T(i,:)={\tilde{A}}^T\mathrm{\α}---\left(12\right)$

基于我们计算出阵列信号的虚拟谱其中第i个元素表示为

$p_{spec}\left(i\right)=\sqrt{{\tilde{s}}^2\left(i\right)+{\tilde{s}}^2(i+N_{\mathrm{\θ}})}---\left(13\right)$

然后对虚拟谱中的元素从大到小进行降序排序。定义sort为降序排序操作，Angle＝sort{Pspec}，谱线中d个最高峰的位置对应估计的角度值

$({\widehat{\mathrm{\θ}}}_1,{\widehat{\mathrm{\θ}}}_2,\mathrm{...},{\widehat{\mathrm{\θ}}}_d)=Angle(1,2,\mathrm{...}d)---\left(14\right)$

从上式可以看出，如果权重系数是稀疏的，我们很容易可以从(14)式中得到角度的估计值。然而，按照上述方法步骤并不能得到满意的直线角度估计值，这是因为最小二乘支持向量机强调得到的是稀疏的拉格朗日乘子，并没有得到我们期望的权重向量的稀疏。而且，实际图像中存在较大的噪声，这些都有可能进一步对估计结果造成影响，因此我们对目标函数做进一步优化。

步骤3、对步骤2中的LS-SVR的对偶优化模型引入稀疏约束项来优化原始对偶优化模型以提高直线角度的精确度；

为了得到稀疏向量我们给传统的SVR模型引入关于的稀疏项来构建一个与原始的关于α稀疏的模型截然不同的新模型，

由于是满秩的，公式(11)中的α和b都与唯一的解，因此在一个线性约束下，公式(11)等价于如下所示的最小化问题：

$\begin{array}{cc}\underset{\mathrm{\α},b}{\min }& \left|\right|\tilde{z}-\left[\begin{array}{cc}Q+\frac{1}{C}{I}_{2M}& 1_{2M\×1}\end{array}\right]\left[\begin{array}{c}\mathrm{\α}\\ b\end{array}\right]|{|}_2^2\\ s.t.& 1_{1\×2M}\mathrm{\α}=0\end{array}---\left(15\right)$

为了确保是稀疏的，给公式(15)引入下面关于的稀疏约束项：

$\begin{array}{cc}\underset{\mathrm{\α},b}{\min }& \left|\right|\tilde{z}-\left[\begin{array}{cc}Q+\frac{1}{C}{I}_{2M}& 1_{2M\×1}\end{array}\right]\left[\begin{array}{c}\mathrm{\α}\\ b\end{array}\right]|{|}_2^2\\ & +\mathrm{\λ}\underset{i=1}{\overset{N_{\mathrm{\θ}}}{\Σ}}\sqrt{{\left({\tilde{A}}^T\right(:,i\left)\mathrm{\α}\right)}^2+{\left({\tilde{A}}^T\right(:,i+N_{\mathrm{\θ}}\left)\mathrm{\α}\right)}^2}\\ s.t.& 1_{1\×2M}\mathrm{\α}=0\end{array}---\left(16\right)$

其中参数λ是稀疏向量和误差范数之间的折中。

引入N_θ+1个新变量把公式(16)改写成另外一种形式，

$\begin{array}{c}\underset{\mathrm{\η},b,\mathrm{\α},{\mathrm{\γ}}_1,\mathrm{...},{\mathrm{\γ}}_{N_{\mathrm{\θ}}}}{\min }\mathrm{\η}+\mathrm{\λ}\underset{i=1}{\overset{N_{\mathrm{\θ}}}{\Σ}}{\mathrm{\γ}}_i\\ \begin{array}{cc}s.t.& \left|\right|\tilde{z}-\left[\begin{array}{cc}Q+\frac{1}{C}{I}_{2M}& 1_{2M\×1}\end{array}\right]\left[\begin{array}{c}\mathrm{\α}\\ b\end{array}\right]|{|}_2^2\≤\mathrm{\η},\end{array}\\ \left|\right|\left[\begin{array}{c}{\tilde{A}}^T(:,i)\\ {\tilde{A}}^T(:,i+N_{\mathrm{\θ}})\end{array}\right]\mathrm{\α}\left|\right|\≤{\mathrm{\γ}}_i,i=1,2,\mathrm{...}{N}_{\mathrm{\θ}},\\ 1_{1\×2M}\mathrm{\α}=0\end{array}---\left(17\right)$

公式(17)是一个二阶锥规划问题，可以由Strum开发的用于处理对称锥优化问题的MATLAB工具来解决，求出稀疏表示系数后再由公式(13)和(14)求出直线的角度，此时的角度是比较精确的。

这样得到的权重向量是稀疏的，很好的提取了角度特征，用式(13)-(14)就可以得到准确的角度估计值。

步骤4、基于步骤3得到的直线角度构建一个直线偏移量的过完备字典，然后再根据步骤3中的对偶优化模型解出直线的偏移量

假设，与角度所对应的潜在的的偏移量可以是 基于和我们可以生成一个大小为M×M的二值图像B_i,j，其中只包含一条偏移量和角度分别为和的直线。

基于B_i,j，(i,j)∈[(1,1),…,(1,N₁),(2,1),…,(2,N₂),…,(d,N_d)]我们构建另一个有个原子的过完备字典B，即

$B=\left[\begin{array}{cccc}vec\left(B_{1,1}\right)& vec\left(B_{1,2}\right)& \mathrm{...}& vec\left(B_{d,N_d}\right)\end{array}\right]\∈R^{M^2\×\underset{i=1}{\overset{d}{\Σ}}{N}_i}---\left(18\right)$

其中vec(B_i,j)代表B_i,j的向量形式。

显然，vec(D)表示成vec(B_i,j)的线性组合，

vec(D)＝Bv+n (19)

其中是v是稀疏的线性组合系数，n是噪声，由于即B的行数非常大，为了避免“维数灾难”，应用压缩感知理论来对其进行降维，即定义一个随机映射矩阵去乘以公式(19)的两边来得到一个等价的方程：

$\tilde{D}=\tilde{B}v+\tilde{n},---\left(20\right)$

其中，q＜＜M²,

公式(20)和公式(7)的求解方法是类似的，应用(8)-(12)及(15)-(17)将公式(20)化成如公式(17)所示的二阶锥规划问题，再由Strum开发的用于处理对称锥优化问题的MATLAB工具来解出最优解v，然后对v的绝对值进行降序排序，即Offset＝sort{abs(v)}，v中绝对值最大的d个峰值就代表检测出来的偏移量即：

$({\tilde{x}}_1,{\tilde{x}}_2,\mathrm{...},{\tilde{x}}_d)=Offset(1,2,\mathrm{...},d)---\left(21\right).$

根据步骤3和步骤4计算出来的角度和偏移量拟合出晶体的直线轮廓边缘，在晶体生长的过程中，根据所得到的直线的角度及偏移量来调节拉杆的位置与方向以控制晶体的生长趋势。

实施例

为了证明本发明方法的有效性，我们用霍夫变换(HT)方法和基于子空间的(SLIDE)方法来和本发明的方法做比较，下面结合具体实例以及图表进行详细说明。

第一步：单晶图像采集及图像二值化

首先采用CCD照相机采集单晶硅生长过程中的晶柱体图像，如图5所示，采用图像增强技术增加图像的有效信息，然后采取阈值分割来确定单晶图像直线边缘的所在的区域，明亮光圈两侧最高点处就是孔径和晶体直线边缘的分界处。所以控制晶体生长的中心轴也就等于估计分界点上侧直线的参数。对这部分图像用Prewitt边缘检测算子得到二值图像，然后应用行扫描线技术得到两侧边缘的采样点，得到用于训练的样本数据为图6所示。

第二步：直线角度估计。

本发明将直线参数估计问题分为两个阶段进行，即第一阶段是角度估计，然后基于估计的角度来估计直线的偏移量。在角度估计阶段，首先基于上述二值图像中的直线形成均匀线性阵列的虚拟快拍数据(如公式(3)所示)，将直线的倾斜角估计问题转换为均匀线性阵列单快拍信号入射角的估计问题，基于这样的虚拟阵列流行以及潜在的入射角范围建立一个过完备字典形成一个基于稀疏表示的回归问题(如公式(7)所示)。本发明采取最小二乘支持向量机(LS-SVR)的对偶优化模型(式(11))。

第三步：估计直线的角度以及对应的偏移量

针对LS-SVR原始模型的缺点以及实际图像中存在的噪声会对估计结果产生较大干扰的问题，我们给传统的中引入权重向量的稀疏约束项(不同于传统的拉格朗日乘子稀疏)来提高估计的准确性。并用二阶锥规划(Second Order Cone Programming，SOCP)优化工具对公式(17)和(20)推导出来的二阶锥规划问题求解最优解。SeDuMi是Strum开发的用于处理对称锥优化问题的MATLAB工具箱，能够用来求解二阶锥和线性约束下的凸优化问题，在SeDuMi中，标准的优化问题形式定义为：

$\begin{array}{ccc}\underset{z}{\max }& p^{T}h& \\ s.t.& r_j-q_j^{T}h\∈{\mathrm{SOC}}^{g_j\×l}& j=1,2,\mathrm{...}J\end{array}---\left(22\right)$

其中的p和r_j是任意的向量，q_j是任意的矩阵，h中包含有期望优化的变量，J是二阶锥约束的个数，g_j维的约束定义为：

${\mathrm{SOC}}^{g_j\×l}=\left\{\right|\left|\mathrm{\ϵ}\right||\≤\widetilde{\mathrm{\ϵ}}\}---\left(23\right)$

在这里是g_j维向量中的第一个元素，ε是g_j-1维的向量，包含了中的其他元素，根据式(17)并引入变量η：

$\begin{array}{c}\underset{\mathrm{\η},b,\mathrm{\α},{\mathrm{\γ}}_1,\mathrm{...},{\mathrm{\γ}}_{N_{\mathrm{\θ}}}}{\min }\mathrm{\η}+\mathrm{\λ}\underset{i=1}{\overset{N_{\mathrm{\θ}}}{\Σ}}{\mathrm{\γ}}_i\\ \begin{array}{cc}s.t.& \left|\right|\tilde{z}-\left[\begin{array}{cc}Q+\frac{1}{C}{I}_{2M}& 1_{2M\×1}\end{array}\right]\left[\begin{array}{c}\mathrm{\α}\\ b\end{array}\right]|{|}_2^2\≤\mathrm{\η},\end{array}\\ \left|\right|\left[\begin{array}{c}{\tilde{A}}^T(:,i)\\ {\tilde{A}}^T(:,i+N_{\mathrm{\θ}})\end{array}\right]\mathrm{\α}\left|\right|\≤{\mathrm{\γ}}_i,i=1,2,\mathrm{...}{N}_{\mathrm{\θ}},\\ 1_{1\×2M}\mathrm{\α}=0\end{array}---\left(17\right)$

其中的定义：

$p=\left[\begin{array}{cccc}-1& 0_{1\×2M}& 0& -{\mathrm{\λI}}_{1\×N_{\mathrm{\θ}}}\end{array}\right],---\left(24.2\right)$

$r_1={\left[\begin{array}{cccc}\frac{1}{2}& \frac{1}{2}& 0& \tilde{z}\end{array}\right]}^T,---\left(24.3\right)$

r₂＝[0 0]^T， (24.5)

r_2+i＝[0 0]^T，i＝1,2,…,N_θ， (24.7)

将式(17)转换成如式(22)的标准的SOCP形式，利用SOCP优化工具包得到很好的解决。

如果解出来公式(17)和(24.1)至(24.8)是角度的计算方法，计算出的拉格朗日乘子为α_角度，与角度所对应的稀疏频谱是对公式(20)应用式(8)-(12)及(15)-(17)然后化成SOCP标准形式解出偏移量的拉格朗日乘子α_偏移量，所以与偏移量所对应的稀疏系数解是直线对应的角度和偏移量估计值可以分别由公式(14)和公式(21)计算出来。不同方法对于单晶图像直线边缘的估计值如表1所示，本发明方法对应的角度和偏移量的谱图如图7和图8所示,晶体图像直线边缘的实际拟合效果如图9所示。

表1不同方法对于单晶图像直线边缘的估计值

小角度直线估计：

本发明还增加了关于小角度测量的仿真实验，我们假设我们的二值图像是如图10所示，两条直线的参数分别是(θ₁＝-43°,)和(θ₂＝-45°,)用本发明方法估计得到的结果如表2所示，对应的频谱估计结果如图11和图12所示。可以看出，本发明方法能够很好的估计出直线的参数，这进一步说明本发明方法在直线角度和偏移量差别很小的情况下也能够很好区分，也证明了本发明方法的准确性。

表2：小角度情况下不同方法的估计结果

Claims (6)

1.一种基于支持向量机回归的晶体图像直线轮廓检测方法，其特征在于，具体按照以下方式实施：

步骤1、首先采用CCD照相机采集单晶硅生长过程中的边缘轮廓线图像，对该边缘轮廓线图像进行预处理，得到用于估计的采样点；

步骤2、基于步骤1中的采样点，构建关于直线角度的过完备字典然后再用最小二乘支持向量机回归LS-SVR模型求解出与直线角度所对应的稀疏表示系数，再由稀疏表示系数求出直线的角度；

步骤3、对步骤2中的LS-SVR的对偶优化模型引入稀疏约束项来优化原始对偶优化模型以提高直线角度的精确度；

步骤4、基于步骤3得到的直线角度构建一个直线偏移量的过完备字典，然后再根据步骤3中的对偶优化模型解出直线的偏移量

2.根据权利要求1所述的基于支持向量机回归的晶体图像直线轮廓检测方法，其特征在于：所述的步骤1预处理采用基于分块和寻找特征区域的方法，用行扫描线方法对晶体图像的两侧边缘的样本点进行采样，得到用于训练的样本数据,其中，图像预处理得到大小为M×M的二值图像D，灰度为“1”的点表示直线点，而灰度为“0”的点表示背景。

3.根据权利要求1或2所述的基于支持向量机回归的晶体图像直线轮廓检测方法，其特征在于：所述的步骤2直线的角度的求解过程具体按照以下步骤实施：

CCD照相机采集单晶硅生长过程中的边缘轮廓线图像，边缘轮廓线是一条直线，图像左上角像素视为坐标原点，水平向右的方向为x轴正向，垂直向下方向为y轴正向，则直线上的点的坐标{x,y}，直线的偏移量以倾斜角θ唯一确定，满足下面关系：

$x=ytan\mathrm{\θ}+\tilde{x}---\left(1\right)$

由公式(1)知，估计图像中包含多条直线的偏移量和角度(θ₁,θ₂,…,θ_d)，检测直线实质上是确定直线参数和θ；

定义向量u：

u＝[1 e^jμ e^j2μ…e^j(M-1)μ]^T ，(2)

对二进制图像D进行转换，将二维数据转换为一维数据，即

$\begin{array}{c}z=Du\\ ={\[\underset{i=1}{\overset{d}{\Σ}}{e}^{j\mathrm{\μ}(0\×{\mathrm{tan\θ}}_i+{\tilde{x}}_i)}\underset{i=1}{\overset{d}{\Σ}}{e}^{j\mathrm{\μ}(1\×{\mathrm{tan\θ}}_i+{\tilde{x}}_i)}\mathrm{...}\underset{i=1}{\overset{d}{\Σ}}{e}^{j\mathrm{\μ}\left(\right(M-1)\×{\mathrm{tan\θ}}_i+{\tilde{x}}_i)}\]}^T\\ =\underset{i=1}{\overset{d}{\Σ}}a\left({\mathrm{\θ}}_i\right)s_i\end{array},---\left(3\right)$

其中，μ是一个常量或者变量，标量向量

经过上述转换，二进制图像D转换为由M个传感器组成的均匀线性阵列的虚拟快拍信号z，其中第i个信号幅度为s_i，第i个信号即第i条直线为平面波，平滑移动到图像第1列像素的左侧由实心原点表示的虚拟均匀线性阵列上，其中，向量u中的μ满足μ|tanθ_i|≤π条件；

为了从单快拍z中估计θ_i，我们建立起来一个关于角度估计的稀疏表示模型，由于这d个信号的入射角，取值范围为：[-90°,90°]，为此，我们考虑将[-90°,90°]细分为N_θ个均匀空间，基于上述均匀间隔化得到的角度集基于公式(4)描述的信号导向矢量结构以及这N_θ个可能入射方向，构造如下的虚拟导向矢量过完备字典：

$A=\[a\left({\widetilde{\mathrm{\θ}}}_1\right)a\left({\widetilde{\mathrm{\θ}}}_2\right)...a\left({\widetilde{\mathrm{\θ}}}_{\mathrm{\θ}}\right)\]\∈C^{M\×N_{\mathrm{\θ}}},---\left(5\right)$

该过完备字典由N_θ列组成，每列具有和公式(4)相同的结构，这样虚拟均匀线性阵列接受信号z可借助这N_θ列重新表示为：

z＝As+n (6)

其中与N_θ个潜在的入射角对应的信号幅度为n为对应的噪声，当s的第i个元素对应于第d个入射角之一时，则取非零值，否则，其取值为0；

公式(6)可以写成另外一种实数形式：

$\tilde{z}=\tilde{A}\tilde{s}+\tilde{n}---\left(7\right)$

其中

其中real(·)和imag(·)分别表示·的实部与虚部，事实上，公式(7)可以看作是一个虚拟线性系统或是一个线性回归问题，其中A(i,:)表示虚拟输入，z(i)表示相关输出i＝1,2,…,2M，SVR来解下面的线性回归问题得到稀疏向量

$\begin{array}{cc}\underset{\tilde{s},b}{min}& \frac{1}{2}{\tilde{s}}^{T}s+\frac{C}{2}\underset{i=1}{\overset{2M}{\Σ}}{e}^2\left(i\right)\end{array}$

$\begin{array}{cc}s.t.& \tilde{z}\left(i\right)=\tilde{A}(i,:)\tilde{s}+b+e\left(i\right),i=1,2,...,2M,\end{array}---\left(8\right)$

其中b偏移量，e(i)为误差项，引入拉格朗日乘子α_i,i＝1,2,…,2M,公式(8)的对偶形式为

$\begin{array}{c}L(\tilde{s},b,\tilde{e},\mathrm{\α})=\frac{1}{2}{\tilde{s}}^T\tilde{s}+\frac{C}{2}\underset{i=1}{\overset{2M}{\Σ}}{\tilde{e}}^2\left(i\right)\\ +\underset{i=1}{\overset{2M}{\Σ}}{\mathrm{\α}}_i\{\tilde{z}\left(i\right)-\tilde{A}(i,:)\tilde{s}-b-\tilde{e}\left(i\right)\}\end{array}---\left(9\right)$

优化条件如下：

$\left\{\begin{array}{c}\frac{\∂L}{\∂\tilde{s}}=0\→\tilde{s}=\underset{i=1}{\overset{2M}{\Σ}}{\mathrm{\α}}_i{\tilde{A}}^T(i,:)\\ \frac{\∂L}{\∂b}=0\→\underset{i=1}{\overset{2M}{\Σ}}{\mathrm{\α}}_i=0,\\ \frac{\∂L}{\∂\tilde{e}\left(i\right)}=0\→{\mathrm{\α}}_i=C\tilde{e}\left(i\right),i=1,...2M,\\ \frac{\∂L}{\∂{\mathrm{\α}}_i}=0\→\tilde{z}\left(i\right)-\tilde{A}(i,:)\tilde{s}-b-\tilde{e}\left(i\right)=0,i=1,...2M\end{array}\right.---\left(10\right),$

根据KKT(Karush-Kuhn-Tucker)条件，公式(9)中描述的问题简化成如下线性系统

$\left[\begin{array}{cc}Q+\frac{1}{C}{I}_{2M}& 1_{2M\×1}\\ 1_{1\×2M}& 0\end{array}\right]\left[\begin{array}{c}\mathrm{\α}\\ b\end{array}\right]=\left[\begin{array}{c}\tilde{z}\\ 0\end{array}\right]---\left(11\right)$

其中α＝[α₁ α₂ …α_2M]^T，1_2M×1表示一个2M×1维向量，每个元素都是1，只要α和b由式(11)解出来，向量由如下的公式计算出来：

$\tilde{s}=\underset{i=1}{\overset{2M}{\Σ}}{\mathrm{\α}}_i{\tilde{A}}^T(i,:)={\tilde{A}}^T\mathrm{\α}---\left(12\right)$

基于计算出阵列信号的虚拟谱其中第i个元素可以表示为

$p_{spec}\left(i\right)=\sqrt{{\tilde{s}}^2\left(i\right)+{\tilde{s}}^2(i+N_{\mathrm{\θ}})}---\left(13\right)$

然后对虚拟谱中的元素从大到小进行降序排序，定义sort为降序排序操作，Angle＝sort{Pspe，c}谱线中d个最高峰的位置对应估计的角度值

$({\widehat{\mathrm{\θ}}}_1,{\widehat{\mathrm{\θ}}}_2,...,{\widehat{\mathrm{\θ}}}_d)=Angle(1,2,...d)---\left(14\right).$

4.根据权利要求3所述的基于支持向量机回归的晶体图像直线轮廓检测方法，其特征在于：所述的步骤3LS-SVR的对偶优化模型引入稀疏约束项来优化原始对偶优化模型以提高直线角度的精确度具体按照如下步骤实施：

公式(11)等价于如下所示的最小化问题：

$\begin{array}{cc}\underset{\mathrm{\α},b}{min}& \left|\right|\tilde{z}-\[\begin{array}{cc}Q+\frac{1}{C}{I}_{2M}& 1_{2M\×1}\end{array}\]\left[\begin{array}{c}\mathrm{\α}\\ b\end{array}\right]|{|}_2^2\end{array}---\left(15\right)$

s.t.1_1×2Mα＝0

在公式(15)引入下面关于的稀疏约束项：

$\begin{array}{cc}\underset{\mathrm{\α},b}{\min }& \left|\right|\tilde{z}-\[\begin{array}{cc}Q+\frac{1}{C}{I}_{2M}& 1_{2M\×1}\end{array}\]\left[\begin{array}{c}\mathrm{\α}\\ b\end{array}\right]|{|}_2^2\\ & +\mathrm{\λ}\underset{i=1}{\overset{N_{\mathrm{\θ}}}{\Σ}}\sqrt{{\left({\tilde{A}}^T\right(:,i\left)\mathrm{\α}\right)}^2+{\left({\tilde{A}}^T\right(:,i+N_{\mathrm{\θ}}\left)\mathrm{\α}\right)}^2}\\ s.t.& 1_{1\×2M}\mathrm{\α}=0\end{array}---\left(16\right)$

其中参数λ是稀疏向量和误差范数之间的折中；

引入N_θ+1个新变量把公式(16)改写成另外一种形式，

$\underset{\mathrm{\η},b,\mathrm{\α},{\mathrm{\γ}}_1,...,{\mathrm{\γ}}_{N_{\mathrm{\θ}}}}{min}\mathrm{\η}+\mathrm{\λ}\underset{i=1}{\overset{N_{\mathrm{\θ}}}{\Σ}}{\mathrm{\γ}}_i$

$\begin{array}{cc}s.t.& \left|\right|\tilde{z}-\[\begin{array}{cc}Q+\frac{1}{C}{I}_{2M}& 1_{2M\×1}\end{array}\]\left[\begin{array}{c}\mathrm{\α}\\ b\end{array}\right]|{|}_2^2\≤\mathrm{\η},\end{array}$

$\left|\right|\left[\begin{array}{c}{\tilde{A}}^T(:,i)\\ {\tilde{A}}^T(:,i+N_{\mathrm{\θ}})\end{array}\right]\mathrm{\α}\left|\right|\≤{\mathrm{\γ}}_i,i=1,2,...N_{\mathrm{\θ}},---\left(17\right)$

1_1×2Mα＝0

将公式(17)推导成一个标准的二阶锥规划问题，由Strum开发的用于处理对称锥优化问题的MATLAB工具来解决。

5.根据权利要求4所述的基于支持向量机回归的晶体图像直线轮廓检测方法，其特征在于：所述的步骤4偏移量的计算具体按照如下步骤实施：

假设与角度所对应的潜在的的偏移量是i＝1,2,…d，基于和生成一个大小为M×M的二值图像B_i,j，其中只包含一条偏移量和角度分别为和的直线；

基于B_i,j，(i,j)∈[(1,1),…,(1,N₁),(2,1),…,(2,N₂),…,(d,N_d)]，构建另一个有个原子的偏移量过完备字典B，即

$B=\[vec\left(B_{1,1}\right)vec\left(B_{1,2}\right)...vec\left(B_{d,N_d}\right)\]\∈R^{M^2\×\underset{i=1}{\overset{d}{\Σ}}{N}_i}---\left(18\right)$

其中vec(B_i,j)代表B_i,j的向量形式；

显然，vec(D)表示成vec(B_i,j)的线性组合，

vec(D)＝Bv+n (19)

其中是v是稀疏的线性组合系数，n是噪声，由于即B的行数非常大，定义一个随机映射矩阵去乘以公式(19)的两边来得到一个等价的方程：

$\tilde{D}=\tilde{B}v+\tilde{n},---\left(20\right)$

其中，

应用(8)-(12)及(15)-(17)将公式(20)化成如公式(17)所示的二阶锥规划问题，再由Strum开发的用于处理对称锥优化问题的MATLAB工具来解出最优解v，然后对v的绝对值进行降序排序，即Offset＝sort{abs(v)}，v中绝对值最大的d个峰值就代表检测出来的偏移量即：

$({\tilde{x}}_1,{\tilde{x}}_2,...,{\tilde{x}}_d)=Offset(1,2,...,d)---\left(21\right).$

6.根据权利要求4或5所述的基于支持向量机回归的晶体图像直线轮廓检测方法，其特征在于：所述的将公式(17)推导成一个标准的二阶锥规划问题，由Strum开发的用于处理对称锥优化问题的MATLAB工具来解决，具体按照以下步骤实施：

在SeDuMi中，标准的优化问题形式定义为：

$\begin{array}{cccc}\underset{z}{\max }& p^{T}h& & \\ Subject& to& r_j-q_j^{T}h\∈{\mathrm{SOC}}^{g_j\×l}& j=1,2,\mathrm{...}J\end{array}---\left(22\right)$

其中的p和r_j是任意的向量，q_j是任意的矩阵，h中包含有期望优化的变量，J是二阶锥约束的个数，g_j维的约束定义为：

${\mathrm{SOC}}^{g_j\×l}=\left\{\right|\left|\mathrm{\ϵ}\right||\≤\widetilde{\mathrm{\ϵ}}\}---\left(23\right)$

在这里是g_j维向量中的第一个元素，ε是g_j-1维的向量，包含了中的其他元素，根据式(17)定义各个变量如下所示：

$\underset{\mathrm{\η},b,\mathrm{\α},{\mathrm{\γ}}_1,...,{\mathrm{\γ}}_{N_{\mathrm{\θ}}}}{min}\mathrm{\η}+\mathrm{\λ}\underset{i=1}{\overset{N_{\mathrm{\θ}}}{\Σ}}{\mathrm{\γ}}_i$

$\begin{array}{cc}s.t.& \left|\right|\tilde{z}-\[\begin{array}{cc}Q+\frac{1}{C}{I}_{2M}& 1_{2M\×1}\end{array}\]\left[\begin{array}{c}\mathrm{\α}\\ b\end{array}\right]|{|}_2^2\≤\mathrm{\η},\end{array}$

$\left|\right|\left[\begin{array}{c}{\tilde{A}}^T(:,i)\\ {\tilde{A}}^T(:,i+N_{\mathrm{\θ}})\end{array}\right]\mathrm{\α}\left|\right|\≤{\mathrm{\γ}}_i,i=1,2,...N_{\mathrm{\θ}},---\left(17\right)$

1_1×2Mα＝0

其中的定义：

$p=\[\begin{array}{cccc}-1& 0_{1\×2M}& 0& -{\mathrm{\λI}}_{1\×N_{\mathrm{\θ}}}\end{array}\],---\left(24.2\right)$

$r_1={\left[\begin{array}{cccc}\frac{1}{2}& \frac{1}{2}& 0& \tilde{z}\end{array}\right]}^T,---\left(24.3\right)$

r₂＝[0 0]^T， (24.5)

r_2+i＝[0 0]^T，i＝1,2,…,N_θ， (24.7)

将式(17)转换成如式(22)的标准的SOCP形式，利用SOCP优化工具包得到很好的解决，

如果解出来公式(17)、(24.1)至(24.8)是角度的计算方法，计算出的拉格朗日乘子为α_角度，与角度所对应的稀疏频谱是直线对应的角度公式(14)计算出来。