CN104050671A - 基于灰色关联分析和Zernike矩的油井管套损图像边缘检测方法 - Google Patents

Abstract

基于灰色关联分析和Zernike矩的油井管套损图像边缘检测方法，针对截取到的油井管套损图像，首先采用灰色关联分析算法对图像边缘进行粗定位，然后利用Zernike矩算子实现图像边缘的亚像素级定位，本发明通过将灰色关联分析和Zernike矩相结合的快速边缘检测方法，灰色关联分析能够检测出连续的有效的图像边缘并通过阈值的调节保留大量图像细节，且具有广泛的适用性；Zernike在亚像素边缘检测中，可有效减少亚像素边缘检测所需要的模板数目、降低函数的阶数、增强抗干扰能力，提高边缘定位精度。该方法首先采用灰色关联分析对套损图像边缘检测进行粗定位，再利用Zernike矩算子定位精度高、抗噪性好等特点，将套损图像边缘精确定位到亚像素级。

Description

基于灰色关联分析和Zernike矩的油井管套损图像边缘检测方法

技术领域

本发明涉及图像处理领域，特别涉及基于灰色关联分析和Zernike矩的油井管套损图像边缘检测方法。

背景技术

目前我国各油田套管出现损坏的油气井超过2万余口，每年还在以近千口的趋势增加，严重影响了油田的安全生产和开发效益。油井管套损检测时实施套管修复工作的基础，对准确评价油井管套损程度、合理分析地层应力、套管及时修复至关重要。传统的经典边缘检测算子如Sobel、Priwitt、Canny等均是像素级检测算子，即检测精度至多达到像素级。由于一些外界因素的影响，使得图像在获取和传输过程中容易受到和边缘点频率相近的噪声的干扰，从而导致检测出的图像边缘总存在伪边缘，噪声干扰和边缘精度低等问题，对于实际物体图像的边缘检测未能取得较好的识别效果。

发明内容

为了克服上述现有技术的不足，本发明的目的在于提供了一种基于灰色关联分析和Zernike矩的油井管套损图像边缘检测方法，以套管损坏检测为背景，通过井下电视技术将采集到的油井管内图像传回地面并实时分析套损情况，在对油井管套损图像进行定量分析中，首先需要对套管变形、破裂等位置边缘进行精确定位，具有检测精度高的特点，可以实现套损图像边缘的亚像素级精确定位。

为了达到上述目的，本发明的技术方案是这样实现的：

基于灰色关联分析和Zernike矩的油井管套损图像边缘检测方法，包括如下步骤：

(1)针对经过预处理后的油井管套损图像，利用灰色关联分析算法进行边缘检测，完成对目标边缘的粗定位；

(2)在单位圆内分别计算Zernike矩Z₀₀，Z₁₁，Z₂₀对应的7×7模板：{M₀₀,M₁₁,M₂₀}

(3)利用模板{M₀₀,M₁₁,M₂₀}和第(1)步中用灰色关联分析算法检测出的边缘点的像素进行卷积运算得到图像的n阶m次Zernike矩{Z₀₀,Z₁₁,Z₂₀}；

(4)取一像素点，根据公式(3)计算边缘角度φ，并计算其他3个边缘参数l、k、h；

在确定边缘参数时，需要通过三个Zernike矩，即Z₀₀、Z₁₁和Z₂₀来计算，其各自对应的正交复数多项式V₀₀，V₁₁，V₂₀分别为V₀₀＝1；V₁₁＝x+jy；V₂₀＝2x²+2y²-1，由于旋转前后保持矩不变性，所以有：

Z'₀₀＝Z₀₀；Z′₁₁＝Z₁₁e^jφ；Z'₂₀＝Z₂₀    (1)

由于旋转φ角后的图像关于x轴对称，所以Z′₁₁的虚部为0，则有：

Im[Z′₁₁]＝sinφRe[Z₁₁]-cosφIm[Z₁₁]＝0     (2)

其中，Re[Z₁₁]、Im[Z₁₁]分别为Z₁₁的实部和虚部，其中：

$\mathrm{\φ}={\tan }^{-1}\left[\frac{\mathrm{Im}\left[Z_{11}\right]}{\mathrm{Re}\left[Z_{11}\right]}\right]---\left(3\right)$

因此，可进一步推出图像旋转后的各阶Zernike矩为：

$Z_{00}^{\′}=\mathrm{h\π}+\frac{\mathrm{k\π}}{2}-k{\sin }^{-1}l-\mathrm{kl}\sqrt{1-l^2}---\left(4\right)$

$Z_{11}^{\′}=\frac{2k{(1-l^2)}^{\frac{3}{2}}}{3}---\left(5\right)$

$Z_{20}^{\′}=\frac{2k{l(1-l^2)}^{\frac{3}{2}}}{3}---\left(6\right)$

结合式(1)，计算可得到模型中其他参数：

$k=\frac{{3Z}_{11}^{\′}}{2{(1-l^2)}^{\frac{3}{2}}}---\left(7\right)$

$h=\frac{Z_{00}-\frac{\mathrm{k\π}}{2}+k{\sin }^{-1}\left(l\right)+\mathrm{kl}\sqrt{1-l^2}}{\mathrm{\π}}---\left(8\right)$

$l=\frac{Z_{20}}{Z_{11}^{\′}}---\left(9\right)$

(5)设定两阈值l₁和k₁，要求l'＝int(l₁×10+0.5)和k'＝int(k₁×10+0.5),l≤l'∩k≥k'，则该点为边缘点，再根据以下公式计算套损图像边缘点的亚像素坐标；

$\left[\begin{array}{c}{x}_s\\ y_s\end{array}\right]=\left[\begin{array}{c}x\\ y\end{array}\right]+\frac{\mathrm{Nl}}{2}\left[\begin{array}{c}\cos \mathrm{\φ}\\ \sin \mathrm{\φ}\end{array}\right]---\left(10\right)$

其中，(x_s,y_s)是图像边缘的亚像素坐标，(x,y)表示图像的原点坐标。

(6)否则返回第(3)步，取下一像素点继续计算，直到所有边缘像素点都被检测一遍。

本发明通过将灰色关联分析和Zernike矩相结合的快速边缘检测方法，灰色关联分析能够检测出连续的有效的图像边缘并通过阈值的调节保留大量图像细节，且具有广泛的适用性；Zernike在亚像素边缘检测中，可有效减少亚像素边缘检测所需要的模板数目、降低函数的阶数、增强抗干扰能力，提高边缘定位精度。该方法首先采用灰色关联分析对套损图像边缘检测进行粗定位，再利用Zernike矩算子定位精度高、抗噪性好等特点，将套损图像边缘精确定位到亚像素级。

附图说明

图1、图2是本发明实施例2的图像亚像素边缘检测前后示意图。

图3、图4是本发明实施例3的图像亚像素边缘检测前后示意图。

具体实施方式

实施例一

本实例的图像处理可分为两部分，第一部分采用灰色关联分析算法对油井管套损图像进行像素级边缘检测，实现像素级的边缘定位。第二部分采用Zernike矩算子针对上一步中的图像进行二次亚像素级边缘定位。以M×N像素点大小的原始图像为例说明本例的实施步骤：

(1)利用灰色关联分析算法对经过预处理的油井管套损图像进行边缘检测，完成对目标边缘的粗定位；

1、确定参考序列和比较序列

为计算方便，对于M×N大小的图像，取值均为1的3×3模板作为参考序列，即：x₀＝(1,1,1,1,1,1,1,1,1)，比较序列分别由图像中各个像素及其周围的8个相邻位置的像素来组成，即：

x_ij＝(x_i-1,j-1,x_i-1,j,x_i-1,j+1,x_i,j-1,x_i,j,x_i,j+1,x_i+1,j-1,x_i+1,j,x_i+1,j+1)

其中i＝1,2……,M；j＝1,2……,N，当i,j＝1或i＝M，j＝N时，重复其相邻的行或列上相应位置的像素值作为该点的值。为描述方便，使用x_r(k)(r＝1,2,…,M×N；k＝1,2,…,5)表示比较序列x_ij。

2、计算参考序列和各比较序列之间的灰色关联系数ξ_0r

${\mathrm{\ξ}}_{0r}\left(k\right)=\frac{\mathrm{\Δ}\min +\mathrm{\ξ\Δ}\max }{\mathrm{\Δ}0r\left(k\right)+\mathrm{\ξ\Δ}\max }---\left(11\right)$

其中，不妨令Δmin＝min|x₀(k)-x_r(k)|、Δmax＝max|x₀(k)-x_r(k)|、Δ0r(k)＝|x₀(k)-x_r(k)|分别为最小绝对差值、最大绝对差值和绝对差值；ξ为分辨常系数，取ξ＝0.2

3、计算以像素点为中心形成的比较序列和参考序列点之间的灰色关联度：

$R_{\mathrm{or}}=\frac{1}{9}\underset{k=1}{\overset{9}{\mathrm{\Σ}}}{\mathrm{\ξ}}_{0r}\left(k\right)---\left(12\right)$

当关联度R_or大于某一给定的阈值θ时，表明该点与参考数列具有相同的特征，不是边缘点，反之，则是边缘点。

(2)在单位圆内分别计算Zernike矩Z₀₀，Z₁₁，Z₂₀对应的7×7模板：{M₀₀,M₁₁,M₂₀}

M₀₀模板

0	0.0287	0.0686	0.0807	0.0686	0.0287	0
							0.0287	0.0815	0.0816	0.0816	0.0816	0.0815	0.0287
0.0686	0.0816	0.0816	0.0816	0.0816	0.0816	0.0686
							0.0807	0.0816	0.0816	0.0816	0.0816	0.0816	0.0807
0.0686	0.0816	0.0816	0.0816	0.0816	0.0816	0.0686
							0.0287	0.0815	0.0816	0.0816	0.0816	0.0815	0.0287
0	0.0287	0.0686	0.0807	0.0686	0.0287	0

M₁₁实数模板

0	-0.0150	-0.0190	0	0.0190	0.0150	0
							-0.0224	-0.0466	-0.0233	0	0.0233	0.0466	0.0224
-0.0573	-0.0466	-0.0233	0	0.0233	0.0466	0.0573
							-0.0690	-0.0466	-0.0233	0	0.0233	0.0466	0.0690
-0.0573	-0.0466	-0.0233	0	0.0233	0.0466	0.0573
							-0.0224	-0.0466	-0.0233	0	0.0233	0.0466	0.0224

0

-0.0150

-0.0190

0

0.0190

0.0150

0

M₁₁虚数模板

0	-0.0224	-0.0573	-0.0690	-0.0573	-0.0224	0
							-0.0150	-0.0466	-0.0466	-0.0466	-0.0466	-0.0466	-0.0150
-0.0190	-0.0233	-0.0233	-0.0233	-0.0233	-0.0233	-0.0190
							0	0	0	0	0	0	0
0.0190	0.0233	0.0233	0.0233	0.0233	0.0233	0.0190
							0.0150	0.0466	0.0466	0.0466	0.0466	0.0466	0.0150
0	0.0224	0.0573	0.0690	0.0573	0.0224	0

M₂₀模板

0	0.0225	0.0394	0.0396	0.0394	0.0225	0
							0.0225	0.0271	-0.0128	-0.0261	-0.0128	0.0271	0.0225
0.0394	-0.0128	-0.0528	-0.0661	-0.0528	-0.0128	0.0394
							0.0396	-0.0261	-0.0661	-0.0794	-0.0661	-0.0261	0.0396
0.0394	-0.0128	-0.0528	-0.0661	-0.0528	-0.0128	0.0394
							0.0225	0.0271	-0.0128	-0.0261	-0.0128	0.0271	0.0225
0	0.0225	0.0394	0.0396	0.0394	0.0225	0

(3)利用模板{M₀₀,M₁₁,M₂₀}和第(1)步中用灰色关联分析算法检测出的边缘点的像素进行卷积运算得到图像的n阶m次Zernike矩{Z₀₀,Z₁₁,Z₂₀}；

(4)取一像素点，根据以下公式计算边缘角度φ，并计算其他3个边缘参数l、k、h；

在确定边缘参数时，需要通过三个Zernike矩，即Z₀₀、Z₁₁和Z₂₀来计算，其各自对应的正交复数多项式V₀₀，V₁₁，V₂₀分别为V₀₀＝1；V₁₁＝x+jy；V₂₀＝2x²+2y²-1。由于旋转前后保持矩不变性，所以有：

Z'₀₀＝Z₀₀；Z′₁₁＝Z₁₁e^jφ；Z'₂₀＝Z₂₀    (13)

由于旋转φ角后的图像关于x轴对称，所以Z′₁₁的虚部为0，则有：

Im[Z′₁₁]＝sinφRe[Z₁₁]-cosφIm[Z₁₁]＝0      (14)

其中，Re[Z₁₁]、Im[Z₁₁]分别为Z₁₁的实部和虚部，其中：

$\mathrm{\φ}={\tan }^{-1}\left[\frac{\mathrm{Im}\left[Z_{11}\right]}{\mathrm{Re}\left[Z_{11}\right]}\right]---\left(15\right)$

因此，可进一步推出图像旋转后的各阶Zernike矩为：

$Z_{00}^{\′}=\mathrm{h\π}+\frac{\mathrm{k\π}}{2}-k{\sin }^{-1}l-\mathrm{kl}\sqrt{1-l^2}---\left(16\right)$

$Z_{11}^{\′}=\frac{2k{(1-l^2)}^{\frac{3}{2}}}{3}---\left(17\right)$

$Z_{20}^{\′}=\frac{2k{l(1-l^2)}^{\frac{3}{2}}}{3}---\left(18\right)$

结合式(13)，计算可得到模型中其他参数：

$k=\frac{{3Z}_{11}^{\′}}{2{(1-l^2)}^{\frac{3}{2}}}---\left(19\right)$

$h=\frac{Z_{00}-\frac{\mathrm{k\π}}{2}+k{\sin }^{-1}\left(l\right)+\mathrm{kl}\sqrt{1-l^2}}{\mathrm{\π}}---\left(20\right)$

$l=\frac{Z_{20}}{Z_{11}^{\′}}---\left(21\right)$

(5)设定两阈值l₁和k₁，要求l'＝int(l₁×10+0.5)和k'＝int(k₁×10+0.5),l≤l'∩k≥k'，则该点为边缘点，再根据公式(22)计算套损图像边缘点的亚像素坐标；

至此，已经计算出图像参数l，k，h，φ，从而推导出图像亚像素边缘的位置为：

$\left[\begin{array}{c}{x}_s\\ y_s\end{array}\right]=\left[\begin{array}{c}x\\ y\end{array}\right]+\frac{\mathrm{Nl}}{2}\left[\begin{array}{c}\cos \mathrm{\φ}\\ \sin \mathrm{\φ}\end{array}\right]---\left(22\right)$

其中，(x_s,y_s)是图像边缘的亚像素坐标，(x,y)表示图像的原点坐标。

(6)若该点不满足阈值条件，则返回第(3)步，取下一像素点继续计算。直到所有边缘像素点检测完毕。

实施例二

以井下电视软件系统在某一井深处的视频截图为例说明本例的实施步骤：

(1)利用灰色关联分析算法对经过预处理的油井管套损图像进行边缘检测，完成对目标边缘的粗定位；

1、确定参考序列和比较序列

为计算方便，取值均为1的3×3模板作为参考序列，即：x₀＝(1,1,1,1,1,1,1,1,1)，比较序列分别由图像中各个像素及其周围的8个相邻位置的像素来组成，即：

x_ij＝(x_i-1,j-1,x_i-1,j,x_i-1,j+1,x_i,j-1,x_i,j,x_i,j+1,x_i+1,j-1,x_i+1,j,x_i+1,j+1)

其中i＝1,2,…,M；j＝1,2,…,N，当i,j＝1或i＝M，j＝N时，重复其相邻的行或列上相应位置的像素值作为该点的值。为描述方便，使用x_r(k)(r＝1,2,…,M×N；k＝1,2,…,5)表示比较序列x_ij。

2、计算参考序列和各比较序列之间的灰色关联系数ξ_0r

${\mathrm{\ξ}}_{0r}\left(k\right)=\frac{\mathrm{\Δ}\min +\mathrm{\ξ\Δ}\max }{\mathrm{\Δ}0r\left(k\right)+\mathrm{\ξ\Δ}\max }---\left(23\right)$

其中，不妨令：Δmin＝min|x₀(k)-x_r(k)|、Δmax＝max|x₀(k)-x_r(k)|、Δ0r(k)＝|x₀(k)-x_r(k)|分别为最小绝对差值、最大绝对差值和绝对差值；ξ为分辨常系数，取ξ＝0.2。

3、计算以各像素点为中心形成的比较序列和参考序列之间的灰色关联度：

$R_{\mathrm{or}}=\frac{1}{9}\underset{k=1}{\overset{9}{\mathrm{\Σ}}}{\mathrm{\ξ}}_{0r}\left(k\right)---\left(24\right)$

当计算出来的关联度R_or大于某一给定的阈值θ时，说明该点与参考数列具有相同的特征，不是边缘点，反之，则是边缘点。本实例中选取阈值θ为0.892。

(2)在单位圆内分别计算Zernike矩Z₀₀，Z₁₁，Z₂₀对应的7×7模板：{M₀₀,M₁₁,M₂₀}

M₀₀模板

0	0.0287	0.0686	0.0807	0.0686	0.0287	0
							0.0287	0.0815	0.0816	0.0816	0.0816	0.0815	0.0287
0.0686	0.0816	0.0816	0.0816	0.0816	0.0816	0.0686
							0.0807	0.0816	0.0816	0.0816	0.0816	0.0816	0.0807
0.0686	0.0816	0.0816	0.0816	0.0816	0.0816	0.0686
							0.0287	0.0815	0.0816	0.0816	0.0816	0.0815	0.0287
0	0.0287	0.0686	0.0807	0.0686	0.0287	0

M₁₁实数模板

0	-0.0150	-0.0190	0	0.0190	0.0150	0
							-0.0224	-0.0466	-0.0233	0	0.0233	0.0466	0.0224
-0.0573	-0.0466	-0.0233	0	0.0233	0.0466	0.0573
							-0.0690	-0.0466	-0.0233	0	0.0233	0.0466	0.0690
-0.0573	-0.0466	-0.0233	0	0.0233	0.0466	0.0573
							-0.0224	-0.0466	-0.0233	0	0.0233	0.0466	0.0224
0	-0.0150	-0.0190	0	0.0190	0.0150	0

M₁₁虚数模板

0	-0.0224	-0.0573	-0.0690	-0.0573	-0.0224	0
							-0.0150	-0.0466	-0.0466	-0.0466	-0.0466	-0.0466	-0.0150
-0.0190	-0.0233	-0.0233	-0.0233	-0.0233	-0.0233	-0.0190
							0	0	0	0	0	0	0
0.0190	0.0233	0.0233	0.0233	0.0233	0.0233	0.0190
							0.0150	0.0466	0.0466	0.0466	0.0466	0.0466	0.0150
0	0.0224	0.0573	0.0690	0.0573	0.0224	0

M₂₀模板

0	0.0225	0.0394	0.0396	0.0394	0.0225	0
							0.0225	0.0271	-0.0128	-0.0261	-0.0128	0.0271	0.0225
0.0394	-0.0128	-0.0528	-0.0661	-0.0528	-0.0128	0.0394
							0.0396	-0.0261	-0.0661	-0.0794	-0.0661	-0.0261	0.0396
0.0394	-0.0128	-0.0528	-0.0661	-0.0528	-0.0128	0.0394
							0.0225	0.0271	-0.0128	-0.0261	-0.0128	0.0271	0.0225

0

0.0225

0.0394

0.0396

0.0394

0.0225

0

(3)利用模板{M₀₀,M₁₁,M₂₀}和第(1)步中用灰色关联分析算法检测出的边缘点的像素进行卷积运算得到图像的n阶m次Zernike矩{Z₀₀,Z₁₁,Z₂₀}；

(4)取一像素点，根据以下公式计算边缘角度φ，并计算其他3个边缘参数l、k、h；

在确定边缘参数时，需要通过三个Zernike矩，即Z₀₀、Z₁₁和Z₂₀来计算，其各自对应的正交复数多项式V₀₀，V₁₁，V₂₀分别为V₀₀＝1；V₁₁＝x+jy；V₂₀＝2x²+2y²-1。由于旋转前后保持矩不变性，所以有：

Z'₀₀＝Z₀₀；Z′₁₁＝Z₁₁e^jφ；Z'₂₀＝Z₂₀      (25)

由于旋转φ角后的图像关于x轴对称，所以Z′₁₁的虚部为0，则有：

Im[Z′₁₁]＝sinφRe[Z₁₁]-cosφIm[Z₁₁]＝0    (26)

其中，Re[Z₁₁]、Im[Z₁₁]分别为Z₁₁的实部和虚部，其中：

$\mathrm{\φ}={\tan }^{-1}\left[\frac{\mathrm{Im}\left[Z_{11}\right]}{\mathrm{Re}\left[Z_{11}\right]}\right]---\left(27\right)$

因此，可进一步推出图像旋转后的各阶Zernike矩为：

$Z_{00}^{\′}=\mathrm{h\π}+\frac{\mathrm{k\π}}{2}-k{\sin }^{-1}l-\mathrm{kl}\sqrt{1-l^2}---\left(28\right)$

$Z_{11}^{\′}=\frac{2k{(1-l^2)}^{\frac{3}{2}}}{3}---\left(29\right)$

$Z_{20}^{\′}=\frac{2k{l(1-l^2)}^{\frac{3}{2}}}{3}---\left(30\right)$

结合式(25)，计算可得到模型中其他参数：

$k=\frac{{3Z}_{11}^{\′}}{2{(1-l^2)}^{\frac{3}{2}}}---\left(31\right)$

$h=\frac{Z_{00}-\frac{\mathrm{k\π}}{2}+k{\sin }^{-1}\left(l\right)+\mathrm{kl}\sqrt{1-l^2}}{\mathrm{\π}}---\left(32\right)$

$l=\frac{Z_{20}}{Z_{11}^{\′}}---\left(33\right)$

(5)设定两阈值l₁和k₁，要求l'＝int(l₁×10+0.5)和k'＝int(k₁×10+0.5)，l≤l'∩k≥k'，则该点为边缘点，再根据公式(34)计算套损图像边缘点的亚像素坐标。

经过实验确定，本例中可设阈值k₁为0.942，l₁为0.712。

至此，已经计算出图像参数l，k，h，φ，从而推导出图像亚像素边缘的位置为：

$\left[\begin{array}{c}{x}_s\\ y_s\end{array}\right]=\left[\begin{array}{c}x\\ y\end{array}\right]+\frac{\mathrm{Nl}}{2}\left[\begin{array}{c}\cos \mathrm{\φ}\\ \sin \mathrm{\φ}\end{array}\right]---\left(34\right)$

其中，(x_s,y_s)是图像边缘的亚像素坐标，(x,y)表示图1原点坐标。

(6)若该点不满足阈值条件，则返回第(3)步，取下一像素点继续计算。直到所有边缘像素点检测完毕。

最终可定位出套损图像亚像素边缘位置，即得到经本算法检测后的油井管套损图像边缘。如图2所示，在该阈值条件下，从检测出的油井管套损图像边缘可看出，在该位置井深处，套管存在明显变形，并且存在一定的缩径长度。并且在该位置正下方，有一定面积油污存在。为后续的定量解释缩径长度，油污面积大小等做好准备

实施例三

以井下电视软件系统在某一井深处的视频截图为例说明本例的实施步骤：

(1)利用灰色关联分析算法对经过预处理的油井管套损图像进行边缘检测，完成对目标边缘的粗定位；

1、确定参考序列和比较序列

为计算方便，取值均为1的3×3模板作为参考序列，即：x₀＝(1,1,1,1,1,1,1,1,1)，比较序列分别由图像中各个像素及其周围的8个相邻位置的像素来组成，即：

x_ij＝(x_i-1,j-1,x_i-1,j,x_i-1,j+1,x_i,j-1,x_i,j,x_i,j+1,x_i+1,j-1,x_i+1,j,x_i+1,j+1)

其中i＝1,2,…,M；j＝1,2,…,N，当i,j＝1或i＝M，j＝N时，重复其相邻的行或列上相应位置的像素值作为该点的值。为描述方便，使用x_r(k)(r＝1,2,…,M×N；k＝1,2,…,5)表示比较序列x_ij。

2、计算参考序列和各比较序列之间的灰色关联系数ξ_0r

${\mathrm{\ξ}}_{0r}\left(k\right)=\frac{\mathrm{\Δ}\min +\mathrm{\ξ\Δ}\max }{\mathrm{\Δ}0r\left(k\right)+\mathrm{\ξ\Δ}\max }---\left(35\right)$

其中，不妨令：Δmin＝min|x₀(k)-x_r(k)|、Δmax＝max|x₀(k)-x_r(k)|、Δ0r(k)＝|x₀(k)-x_r(k)|分别为最小绝对差值、最大绝对差值和绝对差值；ξ为分辨常系数，ξ为分辨常系数，取ξ＝0.2。

3、计算以各像素点为中心形成的比较序列和参考序列之间的灰色关联度：

$R_{\mathrm{or}}=\frac{1}{9}\underset{k=1}{\overset{9}{\mathrm{\Σ}}}{\mathrm{\ξ}}_{0r}\left(k\right)---\left(36\right)$

当计算出来的关联度R_or大于某一给定的阈值θ时，说明该点与参考数列具有相同的特征，不是边缘点，反之，则是边缘点。本实例中选取阈值θ为0.935。

(2)在单位圆内分别计算Zernike矩Z₀₀，Z₁₁，Z₂₀对应的7×7模板：{M₀₀,M₁₁,M₂₀}

M₀₀模板

0	0.0287	0.0686	0.0807	0.0686	0.0287	0
							0.0287	0.0815	0.0816	0.0816	0.0816	0.0815	0.0287
0.0686	0.0816	0.0816	0.0816	0.0816	0.0816	0.0686
							0.0807	0.0816	0.0816	0.0816	0.0816	0.0816	0.0807
0.0686	0.0816	0.0816	0.0816	0.0816	0.0816	0.0686
							0.0287	0.0815	0.0816	0.0816	0.0816	0.0815	0.0287
0	0.0287	0.0686	0.0807	0.0686	0.0287	0

M₁₁实数模板

0	-0.0150	-0.0190	0	0.0190	0.0150	0
							-0.0224	-0.0466	-0.0233	0	0.0233	0.0466	0.0224
-0.0573	-0.0466	-0.0233	0	0.0233	0.0466	0.0573
							-0.0690	-0.0466	-0.0233	0	0.0233	0.0466	0.0690
-0.0573	-0.0466	-0.0233	0	0.0233	0.0466	0.0573
							-0.0224	-0.0466	-0.0233	0	0.0233	0.0466	0.0224
0	-0.0150	-0.0190	0	0.0190	0.0150	0

M₁₁虚数模板

0	-0.0224	-0.0573	-0.0690	-0.0573	-0.0224	0
							-0.0150	-0.0466	-0.0466	-0.0466	-0.0466	-0.0466	-0.0150
-0.0190	-0.0233	-0.0233	-0.0233	-0.0233	-0.0233	-0.0190
							0	0	0	0	0	0	0
0.0190	0.0233	0.0233	0.0233	0.0233	0.0233	0.0190
							0.0150	0.0466	0.0466	0.0466	0.0466	0.0466	0.0150
0	0.0224	0.0573	0.0690	0.0573	0.0224	0

M₂₀模板

0	0.0225	0.0394	0.0396	0.0394	0.0225	0
							0.0225	0.0271	-0.0128	-0.0261	-0.0128	0.0271	0.0225
0.0394	-0.0128	-0.0528	-0.0661	-0.0528	-0.0128	0.0394
							0.0396	-0.0261	-0.0661	-0.0794	-0.0661	-0.0261	0.0396
0.0394	-0.0128	-0.0528	-0.0661	-0.0528	-0.0128	0.0394
							0.0225	0.0271	-0.0128	-0.0261	-0.0128	0.0271	0.0225
0	0.0225	0.0394	0.0396	0.0394	0.0225	0

(3)利用模板{M₀₀,M₁₁,M₂₀}和第(1)步中用灰色关联分析算法检测出的边缘点的像素进行卷积运算得到图像的n阶m次Zernike矩{Z₀₀,Z₁₁,Z₂₀}；

(4)取一像素点，根据以下公式计算边缘角度φ，并计算其他3个边缘参数l、k、h；

在确定边缘参数时，需要通过三个Zernike矩，即Z₀₀、Z₁₁和Z₂₀来计算，其各自对应的正交复数多项式V₀₀，V₁₁，V₂₀分别为V₀₀＝1；V₁₁＝x+jy；V₂₀＝2x²+2y²-1。由于旋转前后保持矩不变性，所以有：

Z'₀₀＝Z₀₀；Z′₁₁＝Z₁₁e^jf；Z'₂₀＝Z₂₀     (37)

由于旋转φ角后的图像关于x轴对称，所以Z′₁₁的虚部为0，则有：

Im[Z′₁₁]＝sinφRe[Z₁₁]-cosφIm[Z₁₁]＝0     (38)

其中，Re[Z₁₁]、Im[Z₁₁]分别为Z₁₁的实部和虚部，其中：

$\mathrm{\φ}={\tan }^{-1}\left[\frac{\mathrm{Im}\left[Z_{11}\right]}{\mathrm{Re}\left[Z_{11}\right]}\right]---\left(39\right)$

因此，可进一步推出图像旋转后的各阶Zernike矩为：

$Z_{00}^{\′}=\mathrm{h\π}+\frac{\mathrm{k\π}}{2}-k{\sin }^{-1}l-\mathrm{kl}\sqrt{1-l^2}---\left(40\right)$

$Z_{11}^{\′}=\frac{2k{(1-l^2)}^{\frac{3}{2}}}{3}---\left(41\right)$

$Z_{20}^{\′}=\frac{2k{l(1-l^2)}^{\frac{3}{2}}}{3}---\left(42\right)$

结合式(37)，计算可得到模型中其他参数：

$k=\frac{{3Z}_{11}^{\′}}{2{(1-l^2)}^{\frac{3}{2}}}---\left(43\right)$

$h=\frac{Z_{00}-\frac{\mathrm{k\π}}{2}-k{\sin }^{-1}\left(l\right)+\mathrm{kl}\sqrt{1-l^2}}{\mathrm{\π}}---\left(44\right)$

$l=\frac{Z_{20}}{Z_{11}^{\′}}---\left(45\right)$

(5)设定两阈值l₁和k₁，要求l'＝int(l₁×10+0.5)和k'＝int(k₁×10+0.5),l≤l'∩k≥k'，则该点为边缘点，再根据公式(46)计算套损图像边缘点的亚像素坐标。

经过实验确定，本例中可设阈值k₁为1.581，l₁为0.942。

至此，已经计算出图像参数l，k，h，φ，从而推导出图像亚像素边缘的位置为：

$\left[\begin{array}{c}{x}_s\\ y_s\end{array}\right]=\left[\begin{array}{c}x\\ y\end{array}\right]+\frac{\mathrm{Nl}}{2}\left[\begin{array}{c}\cos \mathrm{\φ}\\ \sin \mathrm{\φ}\end{array}\right]---\left(46\right)$

其中，(x_s,y_s)是图像边缘的亚像素坐标，(x,y)表示图3原点坐标。

(6)若该点不满足阈值条件，则返回第(3)步，取下一像素点继续计算，直到所有边缘像素点检测完毕。

最终可定位出套损图像亚像素边缘位置，即得到经本算法检测后的油井管套损图像边缘。如图4所示，在该阈值条件下，检测出的边缘更加详细，连续性也更好一些，相比较实施例二中的检测结果，可以更精确的解释井下套损情况，为下一步的定量解释做好准备工作。

Claims (1)

1.基于灰色关联分析和Zernike矩的油井管套损图像边缘检测方法，包括如下步骤：

(1)针对经过预处理后的油井管套损图像，利用灰色关联分析算法进行边缘检测，完成对目标边缘的粗定位；

(2)在单位圆内分别计算Zernike矩Z₀₀，Z₁₁，Z₂₀对应的7×7模板：{M₀₀,M₁₁,M₂₀}；

(3)利用模板{M₀₀,M₁₁,M₂₀}和第(1)步中用灰色关联分析算法检测出的边缘点的像素进行卷积运算得到图像的n阶m次Zernike矩{Z₀₀,Z₁₁,Z₂₀}；

(4)取一像素点，根据公式(3)计算边缘角度φ，并计算其他3个边缘参数l、k、h；

在确定边缘参数时，需要通过三个Zernike矩，即Z₀₀、Z₁₁和Z₂₀来计算，其各自对应的正交复数多项式V₀₀，V₁₁，V₂₀分别为V₀₀＝1；V₁₁＝x+jy；V₂₀＝2x²+2y²-1，由于旋转前后保持矩不变性，所以有：

Z'₀₀＝Z₀₀；Z′₁₁＝Z₁₁e^jφ；Z'₂₀＝Z₂₀       (1)

由于旋转φ角后的图像关于x轴对称，所以Z′₁₁的虚部为0，则有：

Im[Z′₁₁]＝sinφRe[Z₁₁]-cosφIm[Z₁₁]＝0     (2)

其中，Re[Z₁₁]、Im[Z₁₁]分别为Z₁₁的实部和虚部，其中：

$\mathrm{\φ}={\tan }^{-1}\left[\frac{\mathrm{Im}\left[Z_{11}\right]}{\mathrm{Re}\left[Z_{11}\right]}\right]---\left(3\right)$

因此，可进一步推出图像旋转后的各阶Zernike矩为：

$Z_{00}^{\′}=\mathrm{h\π}+\frac{\mathrm{k\π}}{2}-k{\sin }^{-1}l-\mathrm{kl}\sqrt{1-l^2}---\left(4\right)$

$Z_{11}^{\′}=\frac{2k{(1-l^2)}^{\frac{3}{2}}}{3}---\left(5\right)$

$Z_{20}^{\′}=\frac{2k{l(1-l^2)}^{\frac{3}{2}}}{3}---\left(6\right)$

结合式(1)，计算可得到模型中其他参数：

$k=\frac{{3Z}_{11}^{\′}}{2{(1-l^2)}^{\frac{3}{2}}}---\left(7\right)$

$h=\frac{Z_{00}-\frac{\mathrm{k\π}}{2}+k{\sin }^{-1}\left(l\right)+\mathrm{kl}\sqrt{1-l^2}}{\mathrm{\π}}---\left(8\right)$

$l=\frac{Z_{20}}{Z_{11}^{\′}}---\left(9\right)$

(5)设定两阈值l₁和k₁，要求l'＝int(l₁×10+0.5)和k'＝int(k₁×10+0.5),l≤l'∩k≥k'，则该点为边缘点，再根据以下公式计算套损图像边缘点的亚像素坐标：

$\left[\begin{array}{c}{x}_s\\ y_s\end{array}\right]=\left[\begin{array}{c}x\\ y\end{array}\right]+\frac{\mathrm{Nl}}{2}\left[\begin{array}{c}\cos \mathrm{\φ}\\ \sin \mathrm{\φ}\end{array}\right]---\left(10\right)$

其中，(x_s,y_s)是图像边缘的亚像素坐标，(x,y)表示图像的原点坐标；

(6)否则返回第(3)步，取下一像素点继续计算，直到所有边缘像素点都被检测一遍。