CN102903075A - 基于图像特征点全局校正的鲁棒水印方法 - Google Patents

Abstract

本发明公开了一种基于图像特征点全局校正的鲁棒水印方法，主要解决现有水印算法不能有效抵抗常规图像处理和几何攻击的问题。其实现步骤是：(1)利用尺度不变特征变换SIFT法提取原始图像特征点，通过图像分割得到分割区域；(2)结合尺度范围和分割区域构建圆形特征区域；(3)通过抖动量化调制将水印嵌入在特征区域的Zernike矩中；(4)检测时，提取失真图像的SIFT特征点，与原始图像特征点进行匹配，利用RANSAC迭代法对失真图像进行校正；(5)在校正图像的分割区域中构建特征区域，通过抖动量化调制在修改后的Zernike矩中提取水印。本发明具有很好的不可见性，对常规图像处理、几何攻击均具有良好的鲁棒性，可用于互联网上数字作品的版权保护、所有权验证及拷贝控制。

Description

基于图像特征点全局校正的鲁棒水印方法

技术领域

本发明属于信息安全技术领域，具体的说是一种数字图像水印嵌入和盲检测方法，该方法能有效抵抗常规图像处理、几何攻击及组合攻击，可用于互联网上数字作品的版权保护、所有权验证及拷贝控制领域。

背景技术

随着数字技术的不断进步和计算机网络的日益普及，各种形式的多媒体数字作品如图像、视频、音频等纷纷以网络形式发表，多媒体数据正逐渐成为人们获取信息的重要来源。数字化的多媒体数据获取容易、复制简单并且传播迅速，不仅给多媒体信息的存取提供了极大的便利，而且极大提高了信息表达的效率和准确性，但由此引发的盗版问题和版权纷争也成为日益严重的社会问题。任何人都可能在未经信息持有者许可的情况下轻而易举的克隆网络中传播的数字信息或者数字内容并声称自己对原始信息的所有权，甚至伪造他人的数字内容，以期获得非法利益。例如，现代盗版者只需轻点鼠标便可获得原版的复制品，谋取暴利；而一些具有特殊意义的信息，如涉及司法诉讼、政府机要等信息如果遭到恶意攻击和篡改伪造，则会给司法公正和国家安全带来极大的危害。因此如何在利用多媒体信息和计算机网络便利的同时，又可以有效的保护知识产权和保障信息安全已成为了一个亟需解决的现实问题。数字水印技术是实现数字产品版权保护的一种潜在的有效方法，已成为多媒体信息安全领域的一个研究热点，也是信息隐藏研究领域的一个重要分支。数字水印技术一方面弥补了密码技术的缺陷，因为它可以为解密后的数据提供进一步的保护；另一方面，数字水印技术也弥补了数字签名技术的缺陷，因为它可以在原始数据中一次性嵌入大量的秘密信息。

数字水印技术的基本思想是将具有特定意义的标记，即水印，利用数据嵌入的方法隐藏在多媒体数据中，以便保护数字产品的版权、证明产品的真实可靠性、跟踪盗版行为或者提供产品的附加信息。鲁棒的图像水印技术必须具备抵抗多种水印攻击的能力。相对于常规图像处理如噪声、滤波、压缩等，几何攻击如平移、旋转、缩放、剪切，特别是仿射变换、不等比例的缩放、镜面反射等更加难以抵御。几何攻击并不破坏图像水印本身，而是破坏了待检测水印图像与嵌入水印信息之间的同步关系，致使水印检测器无法检测水印信息。

按水印的鲁棒性特点，水印可以分为鲁棒水印和脆弱水印，鲁棒性是大多数数字水印的重要指标，水印鲁棒意味着水印作品能够承受大量的、不同的物理和几何失真。理想情况下，攻击者如果要去掉鲁棒水印必须使水印作品的质量严重下降。而脆弱水印必须对作品的改动很敏感，人们通过脆弱水印的状态判断作品是否被篡改。

基于图像特征的水印方法属于第二代数字水印技术，其基本思想是利用图像中相对稳定的特征点标识出水印嵌入位置，并在与每个特征点对应的局部区域中独立地嵌入水印，检测时仍然利用特征点来定位和检测水印，从而达到抵抗几何攻击的目的。由于图像中提取出的特征点具有一定的稳定性且分布均匀，因此这类方法可以有效的抵抗裁剪攻击。在图像的统计特征中，“矩”具有良好的全局特征表达能力，因此在水印算法中有很好的应用，但因为矩的计算依赖于整幅图像的所有像素，图像如果丢失部分内容必然会导致矩值计算出现很大误差。目前，大部分方法都是结合以上两种方法，即根据图像的特征点选取嵌入区域，在区域中进行矩的计算从而嵌入水印。例如文献Jin-guang Sun,Wei He,“RST Invariant Watermarking Scheme Based on SIFTFeature and Pseudo-Zernike Moment,”IEEE International Symposium on ComputationalIntelligence and Design,vol.2,pp:10-13,2009.先提取图像的特征点，根据特征点构造圆形特征区域，在圆形区域中计算出的矩值上进行水印嵌入和检测。这些方法可以抵抗常规图像处理和一定程度的旋转、缩放等几何攻击，但存在以下问题：几何攻击尤其是复杂的几何攻击会造成图像特征点出现位置的偏移，破坏图像信息与水印信息的同步性，构造的区域内容也会发生改变，圆形特征区域内像素值的改变会导致矩计算出现较大误差，这些问题将大大影响水印检测器的性能，致使水印的检测率较低。

发明内容

本发明的目的在于克服上述已有技术的不足，提供了一种基于特征点匹配的图像矫正方法，将遭受几何攻击后的图像恢复到原始图像状态，以最大程度的恢复图像信息和水印信息的同步性，借助图像分割的信息，提取特征点的位置，在圆形特征区域的Zernike矩值上实现鲁棒水印的嵌入和提取，提高对几何攻击和常规图像处理的鲁棒性。

实现本发明目的的技术方案包括如下步骤：

(1)水印嵌入步骤

(1a)通过密钥Key1生成一个二值的伪随机水印序列b={b₁,b₂,…,b_L}，b_d∈{0,1}，d=1,2，…,L，L是水印序列的位数；

(1b)利用尺度不变特征变换SIFT检测算子，提取原始图像I的SIFT特征点，获得原始图像I的SIFT特征点集F；

(1c)以高斯平滑滤波后的图像I_s的每一个像素为一个点构造网格图，按照距离约束进行图像分割，将平滑图像I_s划分为不同的区域，得到不同分割区域集S={s₁,s₂,…,s_k}，s_t表示一个分割区域，t=1,2,…,k，k是分割区域的个数；

(1d)从分割集S中每一个区域所对应原始图像I的相应区域内，在中频尺度内选择特征强度最大的SIFT特征点，作为圆心构造半径为R的稳定且彼此独立的圆形特征区域，若该区域不存在SIFT特征点，则不选择该分割区域所对应原始图像I的相应区域，从而获得一系列圆形特征区域O={o₁,o₂,…,o_h}，o_l表示一个圆形特征区域，l=1,2，…,h，h是圆形特征区域的个数，h≤k；

(1e)将获得的圆形特征区域四周补0，得到外接方形子图像，计算该外接方形子图像的Zernike矩，并利用抖动量化调制的方法将水印嵌入到筛选出的L个Zernike矩的幅度值中，该L个Zernike矩的位置信息保存为密钥Key2；

(1f)对Zernike矩进行重构，获得含水印的外接方形子图像，并将这些含水印的外接方形子图像去掉周围0值后逐个替换原始圆形特征区域，得到含有水印的图像。

(2)校正受攻击图像步骤：

(2a)利用尺度不变特征变换SIFT检测算子提取受攻击图像I′的SIFT特征点，获得受攻击图像I′的SIFT特征点集F′；

(2b)利用原始图像I的SIFT特征点集F和受攻击图像I′的SIFT特征点集F′，根据距离约束做特征点匹配；

(2c)在匹配好特征点的基础上，利用随机抽样一致RANSAC方法，对匹配好的点进行优化迭代，除去匹配错误的点，计算原始图像I到受攻击图像I′的变换参数T为：

$T=\left[\begin{array}{ccc}{t}_{11}& t_{12}& 0\\ t_{21}& t_{22}& 0\\ t_{31}& t_{32}& 1\end{array}\right],$

式中，t_pq表示待计算的参数，p=1,2,3，q=1,2；

(2d)根据受攻击图像I′和变换参数T，逐像素恢复其到未受攻击时的位置，得到受攻击图像I′校正之后的图像I₁。

(3)水印检测步骤：

(3a)利用尺度不变特征变换SIFT检测算子，提取校正之后图像I₁的尺度不变特征变换SIFT特征点；

(3b)以校正之后图像I₁的每一个像素为一个结点构造网格图，按照距离约束进行图像分割，将校正之后的图像I₁划分为不同区域，得到不同分割区域集

表示一个分割区域，t′=1,2,…,k′，k′是分割区域的个数；

(3c)从分割区域集S′中每一个区域所对应的校正之后图像I₁的相应区域内，在中频尺度内选择特征强度最大的SIFT特征点，作为圆心构造半径为R的稳定且彼此独立的圆形特征区域，若该区域不存在SIFT特征点，则不选择该分割区域所对应的校正之后图像I₁的相应区域，从而获得一系列圆形特征区域

表示一个圆形特征区域，l′=1,2，…,h′，h′是圆形特征区域的个数，h′≤k′；

(3d)对获得的圆形特征区域四周补0，得到外接方形子图像，计算该外接方形子图像的Zernike矩，对Zernike矩进行筛选，得到校正之后图像I₁的Zernike矩集合为：

m≤M_max，n≠4g，g=0,1,2,…，

式中，M_max是最大阶数，表示阶数为m，重复度为n的Zernike矩；

(3e)利用与嵌入水印过程相同的密钥Key2，从中选择L个Zernike矩

用于水印提取，其对应的幅值为

是阶数为m_r，重复度为n_r的Zernike矩，r=1,2,…,L，

是

的幅值；

(3f)通过最小距离解码提取水印

(3g)定义匹配检测阈值X=23，定义x为正确匹配的水印位数，逐位比较原始水印b={b₁,b₂,…,b_L}与提取的水印

得到正确匹配的水印位数x，并将该水印位数x与预先定义的匹配检测阈值X进行比较，判断该圆形特征区域中是否嵌入水印，当x≥X时，则该圆形特征区域嵌入了水印；当x<X时，则该圆形特征区域没有嵌入水印；依次检测校正之后图像I₁的所有圆形特征区域，若检测出大于等于2个的圆形特征区域嵌入了水印，则认为校正之后的图像I₁嵌入了水印，否则认为校正之后的图像I₁没有嵌入水印。

本发明具有以下优点：

(1)本发明由于利用尺度不变特征变换SIFT检测算子和图像分割技术获得一组稳定且彼此独立的圆形特征区域，有效提高了数字水印对几何攻击特别是裁剪攻击的鲁棒性；

(2)本发明由于利用Zernike矩表示圆形特征区域的全局统计特性，克服了因特征点偏移以及插值误差造成的水印检测率低的问题，同时Zernike矩良好的抗噪性、抗压缩能力和旋转不变性，增强了数字水印对常规图像处理和几何攻击的抵抗能力；

(3)本发明由于利用了基于尺度不变特征变换SIFT特征点匹配策略，计算出受攻击图像相对于原始图像的变换矩阵，根据此变换矩阵将受攻击图像矫正至原始图像状态，克服了因图像信息与水印信息的同步性被破坏而造成的水印检测率低的问题，有效提高了数字水印对几何攻击以及组合攻击的鲁棒性。

附图说明

图1是本发明的总体流程图；

图2是本发明中构建的网格示意图；

图3是本发明分割实验的原始图；

图4是本发明分割实验的仿真结果图；

图5是本发明圆形特征区域选择示意图；

图6是本发明中形成的外接方形子图像示意图；

图7是本发明利用抖动量化方式嵌入水印的子流程图；

图8是本发明几何失真Lena图像及校正的仿真结果图；

图9是本发明嵌入水印、全局校正及检测水印的仿真结果图；

图10是本发明中水印对原始图像产生影响的效果示意图。

具体实施方案

参照图1，本发明的实施包括水印嵌入、全局校正和水印检测三个方面。

一.水印嵌入

步骤1，设置密钥Key1，并通过密钥Key1生成一个二值的伪随机水印序列b={b₁,b₂,…,b_L}，b_d∈{0,1}，d=1,2,…,L，L是水印序列的位数。

步骤2：尺度不变特征变换SIFT检测算子利用局部图像特征提取原始图像I的特征点，并描述每个特征点的属性，即位置、尺度和方向，得到SIFT的特征点集F。

2.1)检测尺度空间极值

由不同尺度的高斯核与原始图像I的卷积，得到不同尺度的图像，表示为：

L(x,y,σ)=G(x,y,σ)*I(x,y)

L(x,y,kσ)=G(x,y,kσ)*I(x,y)

式中，G(x,y,σ)表示高斯核函数，σ及kσ表示高斯核函数不同的尺度信息，I(x,y)是原始图像I在第y行、第x列的像素点的像素值，L(x,y,σ)表示尺度信息为σ的高斯核函数与原始图像卷积得到的图像，L(x,y,kσ)表示尺度信息为kσ的高斯核函数与原始图像卷积得到的图像；

原始图像I的每一个像素与同一尺度的临近8个像素以及相邻的上下两个不同尺度对应位置的9×2个像素进行像素值大小的比较，检测该像素的局部极值，确定极值点的位置和所在的尺度，即特征点的位置和所在的尺度；

2.2)消除非稳定的特征点：

用高斯差DoG算子定义为两个不同尺度高斯核的差分与原始图像I的卷积，表示为：

D(x,y,σ)=(G(x,y,kσ)-G(x,y,σ))*I(x,y)=L(x,y,kσ)-L(x,y,σ)，

式中，D(x,y,σ)表示高斯差DoG算子；

在特征点的位置和所在的尺度确定的基础上，用2×2Hessian矩阵H计算每个特征点的稳定性，根据稳定性量度标准来剔除非稳定的特征点，2×2Hessian矩阵H表示为：

$H=\left[\begin{array}{cc}{D}_{\mathrm{xx}}& D_{\mathrm{xy}}\\ D_{\mathrm{xy}}& D_{\mathrm{yy}}\end{array}\right],$

式中，D_xx表示D(x,y,σ)在x处二阶导数，D_yy表示D(x,y,σ)在y处二阶导数，D_xy表示D(x,y,σ)先在x处求一阶偏导数，再在y处求二阶偏导数；

特征点的稳定性量度标准为：

其中r是2×2Hessian矩阵H的最大特征值和最小特征值的比率，用来控制特征点的稳定性。

2.3)指定特征点的主方向：

为了达到图像旋转不变性，计算每个特征点的梯度方向θ，对每个特征点指定一个主方向，表示为：

$\mathrm{\&theta;}={\tan }^{-1}\left(\frac{L_{x,y+1}-L_{x,y-1}}{L_{x+1,y}-L_{x-1,y}}\right),$

式中，tan^-1是反正切操作，L_x+1，y表示尺度图像L(x,y,kσ)在第y行、第x+1列的像素值，L_x-1，y表示尺度图像L(x,y,kσ)在第y行、第x-1列的像素值，L_x，y+1表示尺度图像L(x,y,kσ)在第y+1行、第x列的的像素值，L_x，y-1表示尺度图像L(x,y,kσ)在第y-1行、第x列的像素值；

在特征点为中心的邻域窗口内采样，并用直方图统计特征点所在邻域像素的梯度方向，直方图的峰值处代表了该特征点所在邻域梯度的主方向，作为该特征点的主方向，主方向的峰值代表了该特征点的幅值。

1.4)生成SIFT描述子：

将坐标轴旋转为特征点的主方向，取特征点为中心的4×4的子块，每个像素点定义8个方向的向量信息，因此一个特征点就可以得到128个方向描述符，即为特征描述符，然后将这个1×128的特征描述符归一化，以使其对亮度变化不敏感，归一化后的特征描述符成为特征描述子。

在确定了特征点的位置、尺度、方向及特征描述子后，就得到了原始图像I的特征点集F。

步骤3：对原始图像I进行高斯平滑预处理，去除图像的噪声，得到高斯平滑滤波后的图像I_s：

$I_s=G(x,y)\&CircleTimes;I(x,y),$

式中，

表示线性卷积操作，I(x,y)是原始图像I在第y行、第x列的像素点的像素值，G(x,y)是高斯平滑函数，表示为：

$G(x,y)=\frac{1}{2{\mathrm{\&pi;\&sigma;}}^2}{e}^{-\frac{x^2+y^2}{2{\mathrm{\&sigma;}}^2}},$

式中，σ表示方差。

步骤4，以高斯平滑滤波后的图像I_s的每一个像素为一个点构造网格图，即将平滑滤波后的图像I_s的一个像素点p_i对应一个结点v_i∈V，V是结点的集合，每一个结点v_i连接其邻域的一个结点v_j，构成边(v_i,v_j)∈E，E是边的集合，从而将平滑滤波后的图像I_s构造成由结点集合V及边集合E组成的网格图G=(V,E)，如图2所示。

步骤5，按照距离约束进行图像分割，将平滑滤波后图像I_s划分为不同的区域：

5.1)将网格图G中边的权值定义为该边所连接的两个结点，两个像素点的像素值之差的绝对值表示为：

w((v_i,v_j))=|I_s(p_i)-I_s(p_j)|，

式中，v_i是像素点p_i对应的结点，v_j是像素点p_j对应的结点，(v_i,v_j)表示连接v_i和v_j的边，w((v_i,v_j))表示边(v_i,v_j)的权值，I_s(p_i)表示平滑滤波后的图像I_s中像素点p_i的像素值，I_s(p_j)表示平滑后的图像I_s中像素点p_j的像素值；

5.2)在定义了网格图G中边的权值的基础上，将集合E中的所有边，按照权值的升序进行排列，得到集合π=(e₁,e₂,…,e_M)，e_f表示集合E中的一条边，f=1,2,…,M，M是边的数量；

5.3)初始化分割状态，将每一个点v_i单独作为一个分割区域，获得初始的分割集S⁰={v₁,v₂,…,v_N}，v_i表示集合V中的一个结点，i＝1,2,…,N，N是平滑滤波后图像I_s中像素点的个数；

5.4)对于分割集S^y-1，y＝1,2,…,N，按照以下方法构造新的分割集S^y：

假设

和

分别是分割集S^y-1中的两个不同区域，

包含点v_i，

包含点v_j，e_y=(v_i,v_j)是连接v_i和v_j的边，w(e_y)是边e_y的权值，如果并且

则合并

和

获得新的分割集S^y，否则S^y=S^y-1，其中

定义为最小类间差，表示为：

$\mathrm{MInt}(C_i^{y-1},C_j^{y-1})=\min (\mathrm{Int}\left(C_i^{y-1}\right)+\mathrm{\&tau;}\left(C_i^{y-1}\right),\mathrm{Int}\left(C_j^{y-1}\right)+\mathrm{\&tau;}\left(C_j^{y-1}\right))$

式中，min(·,·)是取小操作，即在

与

中取较小值，定义为区域

的最大类内差，其中

是由

中所有的点构成的最小支撑树，e^y-1是由

中所有的点构成的最小支撑树中的边，

定义为区域

的阈值函数，

表示

中点的个数，

定义为的最大类内差，其中

是由

中所有的点构成的最小支撑树，e^y-1是由

中所有的点构成的最小支撑树中的边，

定义为区域

的阈值函数，表示

中点的个数，K表示分割集S^y-1的每个区域至少包含像素点的个数；

5.5)对于分割集S^y-1，y=1,2,…,N，循环计算5.4)，返回S^y，直到分割状态不再改变，得到平滑滤波后的图像I_s的最后的不同分割区域集S={s₁,s₂,…,s_k}，s_t表示一个分割区域，t=1,2，…,k，k是分割区域的个数。

对于如图3所示的原始图像：其中图3(a)是原始Lene图像，图3(b)是原始Peppers图像，经过分割后得到如图4所示的仿真结果图，其中，图4(a)是对图3(a)原始Lene图像分割后的图像，图4(b)是对图3(b)原始Peppers图像分割后的图像。

步骤6：构造局部圆形特征区域。

从分割集S中每一个区域所对应原始图像I的相应区域内，在中频尺度内选择特征强度最大的SIFT特征点，即主方向幅值最大的SIFT特征点，作为圆心构造半径为R的稳定且彼此独立的圆形特征区域，若该区域不存在SIFT特征点，则不选择该分割区域所对应原始图像I的相应区域，从而获得一系列圆形特征区域O={o₁,o₂,…,o_h}，o_l表示一个圆形特征区域，l=1,2,…,h，h是圆形特征区域的个数，h≤k。

上述步骤2至步骤6的过程如图5所示，先对原始Lene图像进行特征提取和图像分割，再经过构造局部圆形特征区域，得到原始Lene图像中稳定且彼此独立的圆形特征区域。

步骤7：构造圆形特征区域的外接方形子图像，并计算该外接方形子图像的Zernike矩。

如图6所示，首先对图6(a)所示的圆形特征区域建立外接正方形，如图6(b)所示，再在圆形特征区域与外接正方形的间隙补0，得到如图6(c)所示外接方形子图像，然后计算该外接方形子图像的Zernike矩。

步骤8：筛选出符合条件的Zernike矩，并利用抖动量化调制的方法将水印嵌入到筛选出的L个Zernike矩的幅度值中。

8.1)计算外接方形子图像的Zernike矩：

由Zernike矩的相关理论得知，部分Zernike矩存在微小的计算误差，必须合理选择Zernike矩，Zernike矩的选择应该考虑以下两个方面：首先，定义最大阶数M_max=30，选择阶数较低的Zernike矩，因为当阶数高于数值M_max时，Zernike矩的计算将不再准确，再次，重复度为n=4g，g=0,1,2,…的Zernike矩存在微小计算误差，所以不选择这些Zernike矩，则建立合理的Zernike矩集合S_Zernike为：

S_Zernike={Z_mn}，m≤M_max，n≠4g，g=0,1，2,…，

式中，Z_mn表示阶数为m，重复度为n的Zernike矩；

从Zernike矩集合S_Zernike中随机选择L个Zernike矩

用于水印嵌入，设其对应的Zernike矩幅值为该L个Zernike矩的位置信息保存为密钥Key2，其中，

是阶数为m_r，重复度为n_r的Zernike矩，r=1,2,…,L，是

对应的幅值，L≤W，W是集合S_Zernike中的元素个数；

8.2)利用步骤1产生的水印序列b={b₁,b₂,…,b_L}中的每位水印b_r，r=1,2,…,L，来量化相应的Zernike矩幅值

实现水印的嵌入，量化公式为：

$A_{m_r{n}_r}^{\&prime;}=\left[\frac{A_{m_r{n}_r}-d_r\left(b_r\right)}{\mathrm{\&Delta;}}\right]\mathrm{\&Delta;}+d_r\left(b_r\right),$

式中，是对应量化后的Zernike矩幅值，[·]是四舍五入操作，Δ是量化步长，d_r(·)是第r个抖动函数，且满足d_r(1)=Δ/2+d_r(0)；矢量(d₁(0)，…,d_L(0))通过密钥Key3产生，且在区间[0,1]上服从均匀分布；

在量化Zernike矩幅值

时，如果n_r≠0，则要同时应用水印b_r来量化它的共轭矩的幅值，以保证它们有相同的幅值；

8.3)通过L个量化后的Zernike矩的幅度值获得L个修改后的Zernike矩，其表示为：

$Z_{m_r{n}_r}^{\&prime;}=\frac{A_{m_r{n}_r}^{\&prime;}}{A_{m_r{n}_r}}{Z}_{m_r{n}_r},$ r=1,…,L，

式中，

为Z中的第r个Zernike矩的幅度值，

为量化后的第r个Zernike矩的幅度值，

为Z中的第r个Zernike矩，为修改后的第r个Zernike矩。

步骤9：对Zernike矩进行重构，获得含水印的局部方形子图像，将其四周去“0”，得到含水印的局部圆形特征区域，将含水印的圆形特征区域逐一替换原始的圆形特征区域，获得含水印的图像。

9.1)含水印的局部方形子图像由两部分合并而成：第一部分是对未被选择的Zernike矩重构得到的局部方形子图像f_rest(x,y)：

f_rest(x,y)=f_o(x,y)-f_Z(x,y)，

式中，f_o(x,y)是原始局部方形子图像，f_Z(x,y)是Z中待修改的L个Zernike矩通过重构得到的部分子图像，且有：

$f_Z(x,y)=\underset{r=1}{\overset{L}{\mathrm{\&Sigma;}}}{Z}_{m_r{n}_r}{V}_{m_r{n}_r}(x,y)+Z_{m_r,-n_r}{V}_{m_r,-n_r}(x,y),$

式中，

是求和操作，

是阶数为m_r，重复度为n_r的Zernike矩，

是阶数为m_r，重复度为n_r的正交函数，

是的共轭矩，

是阶数为m_r，重复度为-n_r的正交函数，x,y表示(x,y)处的像素点位于原始图像的第y行、第x列；

第二部分是已经修改Zernike矩的重构部分子图像f_Z′(x,y)：

$f_{Z^{\&prime;}}(x,y)=\underset{r=1}{\overset{L}{\mathrm{\&Sigma;}}}{Z}_{m_r{n}_r}^{\&prime;}{V}_{m_r{n}_r}(x,y)+Z_{m_r,-n_r}^{\&prime;}{V}_{m_r,-n_r}(x,y),$

式中，

为对应

修改后的Zernike矩，

是

的共轭矩；

9.2)将第一部分正方形子图像块f_rest(x,y)和第二部分正方形子图像块f_Z′(x,y)进行合并，即得到含水印的局部方形子图像f′(x,y)：

f′(x,y)=f_rest(x,y)+f_Z′(x,y)；

9.3)将所有含水印的外接方形子图像块替换原始外接方形子图像块，并将含有水印的外接方形子图像块去“0”得到含水印的圆形特征区域，替换所有的原始圆形特征区域后即可得到含水印图像。

上述步骤8至步骤9的过程如图7所示，先计算该外接方形子图像的Zernike矩，筛选出符合条件的Zernike矩，并利用抖动量化调制的方法将水印嵌入到筛选出的L个Zernike矩的幅度值中，再合并对未被选择的Zernike矩重构得到的局部方形子图像及对修改Zernike矩重构得到部分子图像，得到获得含水印的局部方形子图像。

二.全局校正

步骤10：利用尺度不变特征变换SIFT检测算子提取受攻击图像I′的SIFT特征点，获得受攻击图像I′的SIFT特征点集F′。

步骤11：利用原始图像I的SIFT特征点集F和受攻击图像I′的SIFT特征点集F′，根据距离约束做特征点匹配。

11.1)对于原始图像I特征点集F中的任意一个SIFT特征点，计算其与受攻击图像I′特征点集F′中所有SIFT特征点的欧氏距离；

11.2)设置比例阈值μ=0.6，在欧氏距离中找出距该点的最近距离和次近距离，当最近距离除以次近距离小于比例阈值μ时，则距该点最近距离对应的点就是与该点相匹配的点。

步骤12：在匹配好特征点的基础上，利用随机抽样一致RANSAC方法，对匹配好的点进行优化迭代，除去匹配错误的点，计算原始图像I到受攻击图像I′的变换参数T为：

$T=\left[\begin{array}{ccc}{t}_{11}& t_{12}& 0\\ t_{21}& t_{22}& 0\\ t_{31}& t_{32}& 1\end{array}\right],$

式中，t_pq表示待计算的参数，p=1,2,3，q=1,2。

所属的随机抽样一致RANSAC方法，是文献M.A.Fischler and R.C.Bolles,“Random sample consensus:a paradigm for model fitting with applications to imageanalysis and automated car-tography”，Comm.ACM,vol.24,no.6,pp.381-395,Jun.1981.中所述的方法。

步骤13：根据受攻击图像I′和变换参数T，对受攻击图像I′逐像素恢复到未受攻击时的位置，得到受攻击图像I′校正之后的图像I₁，即：初始化一个零矩阵I_v，大小与原始图像I相同，对于I_v中的一个点(x,y)，做如下变换：

$\left[\begin{array}{c}{x}^{\&prime;}\\ y^{\&prime;}\\ 1\end{array}\right]=T\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]=\left[\begin{array}{ccc}{t}_{11}& t_{12}& 1\\ t_{21}& t_{22}& 1\\ t_{31}& t_{32}& 0\end{array}\right]\left[\begin{array}{c}x\\ y\\ 1\end{array}\right],$

计算出x′和y′的值，对受攻击图像I′采用双线性差值的方法逐个进行像素值的恢复，具体方法如下：

式中，

代表向上取整操作，

代表向下取整操作，

表示受攻击图像I′中位于第

行，第

列的像素值，

表示受攻击图像I′中位于第

行，第

列的像素值，

表示受攻击图像I′中位于第

行，第

列的像素值，表示受攻击图像I′中位于第

行，第

列的像素值，I_v1(x,y)和I_v2(x,y)是中间操作符，将受攻击图像I′中所有的点逐一进行如上操作，所有的点遍历结束后，即可得到矫正之后的图像I₁=I_v。

图8所示为几何失真Lena图像及全局校正后的仿真结果图，其中图8(a)为全局仿射变换后的Lena图像，图8(b)为不等比例缩放后的Lena图像，图8(c)为镜面反射后的Lena图像，图8(d)为旋转15度后的Lena图像，经过全局校正之后得到如图8(e)、图8(f)、图8(g)及图8(h)所示的仿真结果图，即图8(e)为图8(a)全局校正之后的图像，图8(f)为图8(b)全局校正之后的图像，图8(g)为图8(c)全局校正之后的图像，图8(h)为图8(d)全局校正之后的图像。

三.水印检测

步骤14：利用尺度不变特征变换SIFT检测算子，提取校正之后图像I₁的尺度不变特征变换SIFT特征点。

由于利用SIFT检测算子直接从受攻击图像I′中提取特征点，可能会和嵌入水印时从原始图像I中提取出的特征点并非完全一致，甚至会有很大的差异，故本发明利用尺度不变特征变换SIFT检测算子，提取校正之后图像I1的尺度不变特征变换SIFT特征点

步骤15：以校正之后图像I₁的每一个像素为一个结点构造网格图，按照距离约束进行图像分割，将校正之后的图像I₁划分为不同的区域，得到不同分割区域集

表示一个分割区域，t′=1,2，…,k′，k′是分割区域的个数。

步骤16：从分割区域集S′中每一个区域所对应的校正之后图像I₁的相应区域内，在中频尺度内选择特征强度最大的SIFT特征点，作为圆心构造半径为R的稳定且彼此独立的圆形特征区域，若该区域不存在SIFT特征点，则不选择该分割区域所对应的校正之后图像I₁的相应区域，从而获得一系列圆形特征区域

表示一个圆形特征区域，l′=1,2,…,h′，h′是圆形特征区域的个数，h′≤k′。

步骤17：对获得的圆形特征区域四周补0，得到外接方形子图像，计算该外接方形子图像的Zernike矩，对Zernike矩进行筛选，得到校正之后图像I₁的Zernike矩集合

为：

m≤M_max，n≠4g，g=0,1，2,…，

式中，M_max是最大阶数，

表示阶数为m，重复度为n的Zernike矩。

步骤18：利用与嵌入水印过程相同的密钥Key2，从

中选择L个Zernike矩

用于水印提取，该L个Zernike矩Z′对应的幅值为

是阶数为m_r，重复度为n_r的Zernike矩，r=1,2,…,L，

是

的幅值。

步骤19：通过最小距离解码提取水印

19.1)采用抖动量化调制方法，对L个Zernike矩幅值

进行量化，量化公式为：

${\left(A_{m_r{n}_r}^{\&prime;}\right)}_Q=\left[\frac{A_{m_r{n}_r}^{\&prime;}-d_r\left(Q\right)}{\mathrm{\&Delta;}}\right]\mathrm{\&Delta;}+d_r\left(Q\right),$ Q=0,1，r=1,…,L，

式中，[·]是四舍五入操作，Δ是量化步长，d_r(Q)是第r个抖动函数，满足d_r(1)=Δ/2+d_r(0)，抖动矢量(d₁(0),…,d_L(0))是利用与嵌入水印过程相同的密钥Key3产生的，

是

采用抖动函数d_r(Q)量化后的Zernike矩幅值；

通过对每一个Zernike矩的幅值进行量化，得到两组量化式

和

r=1,…,L；

19.2)将上述

与

之间的距离定义为：

将上述与

之间的距离定义为：

r=1,…,L，通过比较dis0与dis1的大小，提取L位水印信息提取公式为：

$b_r^{\&prime;}=\underset{Q\&Element;\left\{\mathrm{0,1}\right\}}{\arg \min }{({\left(A_{m_r{n}_r}^{\&prime;}\right)}_Q-A_{m_r{n}_r}^{\&prime;})}^2,$ r=1，…,L，

式中，

是第r个Zernike矩的幅值，argmin操作表示为：如果dis0<dis1，则

否则，r=1,…,L，L是Zernike矩的个数。

步骤20：定义每个圆形特征区域的虚警概率P_fp≤10^-4，确定匹配检测阈值X=23，逐位比较原始水印b={b₁,b₂,…,b_L}与提取的水印

得到正确匹配的水印位数x，并将该水印位数x与预先确定的匹配检测阈值X进行比较，判断该圆形特征区域中是否嵌入水印，当x≥X时，则该圆形特征区域嵌入了水印；当x<X时，则该圆形特征区域没有嵌入水印。

步骤21：依次检测校正之后图像I₁的所有圆形特征区域，定义m是检测出嵌入了水印的圆形区域的个数，定义校正之后图像I₁的虚警概率P_FP≤10^-5，确定m取2，故检测出大于等于2个的圆形特征区域嵌入了水印，则认为校正之后的图像I₁嵌入了水印，否则认为校正之后的图像I₁没有嵌入水印。

上述步骤1至步骤21的仿真结果如图9所示，其中图9(a)是原始Lena图像嵌入水印的圆形区域示意图，图9(b)是对含水印图像进行镜面反射后的示意图，图9(c)是图9(b)经过全局校正之后的图像，图9(d)是从校正之后的图像中检测出嵌入了水印的圆形区域示意图。

本发明的优点可通过以下仿真实验进一步说明：

本发明在大量标准灰度图像上进行了测试实验，其中包括基准测试图像Lena、Peppers，以不可见性和鲁棒性作为本发明性能优劣的评测依据。

(1)不可见性

本发明以客观指标峰值信噪比PSNR作为评价不可见性的依据。本发明中的PSNR值主要受三个因素的影响：当水印长度和圆形区域半径固定时，抖动调制中的量化步长Δ影响PSNR值，Δ越大，PSNR值越小；当Δ和圆形区域半径固定时，水印位数越长，PSNR值越小；当水印长度和Δ固定时，圆形区域半径越大，PSNR值越小。

本发明的实验中，量化步长Δ=5，水印长度L=30，圆形区域的半径作为一个关键因子来提高水印系统的安全性，使PSNR值大于50dB。如图10所示的水印对原始图像产生影响的效果示意图，其中图10(a)为原始Lena图像，图10(b)为含水印的Lena图像，图10(c)为原始Peppers图像，图10(d)为含水印的Peppers图像，图10说明了本发明具有很好的不可见性。

(2)鲁棒性

利用评测工具Stirmark 4.0对含水印图像进行一系列攻击实验，以测试本发明的鲁棒性。表1和表2分别给出了本发明对含水印图像进行常规图像处理和几何攻击后的水印检测结果。以检测率作为评价本发明鲁棒性的指标，表格中的分数即为检测率，检测率DR为：

$\mathrm{DR}=\frac{\#\mathrm{hitregions}}{\#\mathrm{hostregions}}$

式中，分数的分子#hit regions表示从校正之后的图像中正确检测到嵌入了水印的圆形区域个数，分母#host regions表示原始图像中嵌入水印的区域个数，即9个。

表1水印抵抗常规图像处理操作的检测率

从表1可以看出，本发明能有效抵抗常规图像处理。这是因为：(1)提取的特征点具有很好的稳定性，在经过常规图像处理后，水印的嵌入位置几乎没有发生改变，保证了水印信息的正确提取；(2)局部Zernike矩本身对加噪和JPEG压缩具有很好的抵抗力。

表2水印抵抗几何攻击的检测率

从表2可以看出，本方法对几何攻击，包括仿射变换、不等比例缩放、镜面发射等复杂的几何攻击都具有较好的鲁棒性，这是因为本发明采用了基于特征点匹配的全局校正方法，克服了因图像信息与水印信息的同步性被破坏而造成的水印检测率低的问题。

综上，本发明提高了对几何失真图像校正恢复的能力和效果，提高了数字水印针对常规图像处理和复杂几何攻击的鲁棒性。

Claims (7)

1.一种基于图像特征点全局校正的鲁棒水印方法，包括：

(1)水印嵌入步骤：

(1a)通过密钥Key1生成一个二值的伪随机水印序列b={b₁,b₂,…,b_L}，b_d∈{0,1}，d=1,2,…,L，L是水印序列的位数；

(1b)利用尺度不变特征变换SIFT检测算子，提取原始图像I的SIFT特征点，获得原始图像I的SIFT特征点集F；

(1c)以高斯平滑滤波后的图像I_s的每一个像素为一个点构造网格图，按照距离约束进行图像分割，将平滑图像I_s划分为不同的区域，得到不同分割区域集S={s₁,s₂,…,s_k}，s_t表示一个分割区域，t=1,2,…,k，k是分割区域的个数；

(1d)从分割集S中每一个区域所对应原始图像I的相应区域内，在中频尺度范围内选择特征强度最大的SIFT特征点，作为圆心构造半径为R的稳定且彼此独立的圆形特征区域，若该区域不存在SIFT特征点，则不选择该分割区域所对应原始图像I的相应区域，从而获得一系列圆形特征区域O={o₁,o₂,…,o_h}，o_l表示一个圆形特征区域，l=1,2，…,h，h是圆形特征区域的个数，h≤k；

(1e)将获得的圆形特征区域四周补0，得到外接方形子图像，计算该外接方形子图像的Zernike矩，并利用抖动量化调制的方法将水印嵌入到筛选出的L个Zernike矩的幅度值中，该L个Zernike矩的位置信息保存为密钥Key2；

(1f)对Zernike矩进行重构，获得含水印的外接方形子图像，并将这些含水印的外接方形子图像去掉周围0值后逐个替换原始圆形特征区域，得到含有水印的图像；

(2)校正受攻击图像步骤：

(2a)利用尺度不变特征变换SIFT检测算子提取受攻击图像I′的SIFT特征点，获得受攻击图像I′的SIFT特征点集F′；

(2b)利用原始图像I的SIFT特征点集F和受攻击图像I′的SIFT特征点集F′，根据距离约束做特征点匹配；

(2c)在匹配好特征点的基础上，利用随机抽样一致RANSAC方法，对匹配好的点进行优化迭代，除去匹配错误的点，计算原始图像I到受攻击图像I′的变换参数T为：

$T=\left[\begin{array}{ccc}{t}_{11}& t_{12}& 0\\ t_{21}& t_{22}& 0\\ t_{31}& t_{32}& 1\end{array}\right],$

式中，t_pq表示待计算的参数，p=1,2,3，q=1,2；

(2d)根据受攻击图像I′和变换参数T，逐像素恢复其到未受攻击时的位置，得到受攻击图像I′校正之后的图像I₁；

(3)水印检测步骤：

(3a)利用尺度不变特征变换SIFT检测算子，提取校正之后图像I₁的尺度不变特征变换SIFT特征点；

(3b)以校正之后图像I₁的每一个像素为一个结点构造网格图，按照距离约束进行图像分割，将校正之后的图像I₁划分为不同区域，得到不同分割区域集

表示一个分割区域，t′=1,2，…,k′，k′是分割区域的个数；

(3c)从分割区域集S′中每一个区域所对应的校正之后图像I₁的相应区域内，在中频尺度内选择特征强度最大的SIFT特征点，作为圆心构造半径为R的稳定且彼此独立的圆形特征区域，若该区域不存在SIFT特征点，则不选择该分割区域所对应的校正之后图像I₁的相应区域，从而获得一系列圆形特征区域

表示一个圆形特征区域，l′=1,2，…,h′，h′是圆形特征区域的个数，h′≤k′；

(3d)对获得的圆形特征区域四周补0，得到外接方形子图像，计算该外接方形子图像的Zernike矩，对Zernike矩进行筛选，得到校正之后图像I₁的Zernike矩集合

为：

m≤M_max，n≠4g，g=0,1,2,…，

式中，M_max是最大阶数，

表示阶数为m，重复度为n的Zernike矩；

(3e)利用与嵌入水印过程相同的密钥Key2，从

中选择L个Zernike矩

用于水印提取，其对应的幅值为

是阶数为m_r，重复度为n_r的Zernike矩，r=1,2,…,L，

是

的幅值；

(3f)通过最小距离解码提取水印

(3g)定义匹配检测阈值X=23，定义x为正确匹配的水印位数，逐位比较原始水印b={b₁,b₂,…,b_L}与提取的水印

得到正确匹配的水印位数x，并将该水印位数x与预先定义的匹配检测阈值X进行比较，判断该圆形特征区域中是否嵌入水印，当x≥X时，则该圆形特征区域嵌入了水印；当x<X时，则该圆形特征区域没有嵌入水印；依次检测校正之后图像I₁的所有圆形特征区域，若检测出大于等于2个的圆形特征区域嵌入了水印，则认为校正之后的图像I₁嵌入了水印，否则认为校正之后的图像I₁没有嵌入水印。

2.根据权利要求1所述的基于图像特征点全局校正的鲁棒水印方法，其中步骤(1c)所述的以高斯平滑滤波后的图像I_s的每一个像素为一个点构造网格图，按如下步骤进行：

(1c1)对原始图像I进行高斯平滑预处理，去除图像的噪声，得到高斯平滑滤波后的图像I_s：

$I_s=G(x,y)\&CircleTimes;I(x,y)$

式中，

表示线性卷积操作，I(x,y)是原始图像I在第y行、第x列的像素点的像素值，G(x，y)是高斯平滑函数，表示为：

$G(x,y)=\frac{1}{2{\mathrm{\&pi;\&sigma;}}^2}{e}^{-\frac{x^2+y^2}{2{\mathrm{\&sigma;}}^2}},$

式中，σ表示方差；

(1c2)将平滑后的图像I_s的一个像素点p_i对应一个结点v_i∈V，V是结点的集合，每一个结点v_i连接其邻域的一个结点v_j，构成边(v_i,v_j)∈E，E是边的集合，从而将平滑后的图像I_s构造成由结点集合V及边集合E组成的网格图G=(V,E)。

3.根据权利要求1所述的基于图像特征点全局校正的鲁棒水印方法，其中步骤(1c)所述的按照距离约束进行图像分割，按如下步骤进行：

(1c3)将网格图G中边的权值定义为该边所连接的两个结点，两个像素点的像素值之差的绝对值，表示为：

w((v_i,v_j))=|I_s(p_i)-I_s(p_j)|，

式中，v_i是像素点p_i对应的结点，v_j是像素点p_j对应的结点，(v_i,v_j)表示连接v_i和v_j的边，w((v_i,v_j))表示边(v_i,v_j)的权值，I_s(p_i)表示平滑后的图像I_s中像素点p_i的像素值，I_s(p_j)表示平滑后的图像I_s中像素点p_j的像素值；

(1c4)在定义了网格图G中边的权值的基础上，对平滑后的图像I_s进行分割，当由边(v_i,v_j)连接的两个结点v_i和v_j位于两个不同区域内，并且w((v_i,v_j))的值小于定义的阈值时，将这两个不同区域合并，否则保持原分割状态，依次进行，直到分割状态不再改变，从而得到平滑后的图像I_s的分割区域集S={s₁,s₂,…,s_k}。

4.根据权利要求1所述的基于图像特征点全局校正的鲁棒水印方法，其中步骤(1e)所述的利用抖动量化调制的方法将水印嵌入到筛选出的L个Zernike矩的幅度值中，按如下步骤进行：

(1e1)定义最大阶数M_max，选择阶数小于等于M_max，且重复度为n的Zernike矩，建立Zernike矩集合S_Zernike为：

S_Zernike={Z_mn}，m≤M_max，n≠4g，g=0,1,2,…，

式中，Z_mn表示阶数为m，重复度为n的Zernike矩；

(1e2)从集合S_Zernike中随机选择L个Zernike矩

用于水印嵌入，其对应的Zernike矩幅值为

是阶数为m_r，重复度为n_r的Zernike矩，r=1,2,…,L，

是

对应的幅值，L≤W，W是集合S_Zernike中的元素个数；

(1e3)利用水印序列b={b₁,b₂,…,b_L}中的每位水印b_r，r=1,2,…,L，来量化相应的Zernike矩幅值

实现水印的嵌入，量化公式为：

$A_{m_r{n}_r}^{\&prime;}=\left[\frac{A_{m_r{n}_r}-d_r\left(b_r\right)}{\mathrm{\&Delta;}}\right]\mathrm{\&Delta;}+d_r\left(b_r\right),$

式中，

是对应量化后的Zernike矩幅值，[·]是四舍五入操作，Δ是量化步长，d_r(·)是第r个抖动函数，且满足d_r(1)=Δ2+d_r(0)；矢量(d₁(0)，…,d_L(0))通过密钥Key3产生，且在区间[0,1]上服从均匀分布；

在量化Zernike矩幅值

时，如果n_r≠0，则要同时应用水印b_r来量化它的共轭矩的幅值。

5.根据权利要求1所述的基于图像特征点全局校正的鲁棒水印方法，其中步骤(1f)所述的对Zernike矩进行重构，获得含水印的外接方形子图像，按如下步骤进行：

(1f1)通过L个量化后的Zernike矩的幅度值获得L个修改后的Zernike矩，其表示为：

$Z_{m_r{n}_r}^{\&prime;}=\frac{A_{m_r{n}_r}^{\&prime;}}{A_{m_r{n}_r}}{Z}_{m_r{n}_r},$ r=1,…，L，

式中，

为Z中的第r个Zernike矩的幅度值，

为量化后的第r个Zernike矩的幅度值，

为Z中的第r个Zernike矩，

为修改后的第r个Zernike矩；

(1f2)利用未被选择的Zernike矩，重构得到第一组正方形子图像块f_rest(x,y)：

f_rest(x,y)=f_o(x,y)-f_Z(x,y)，

式中，f_o(x,y)是原始局部方形子图像，f_Z(x,y)是Z中待修改的L个Zernike矩通过重构得到的部分子图像，x,y表示(x,y)处的像素点位于原始图像的第y行、第x列；

(1f3)利用L个修改后的Zernike矩通过重构得到第二组方形子图像块f_Z′(x,y)；

(1f4)将第一组正方形子图像块f_rest(x,y)和第二组正方形子图像块f_Z′(x,y)进行合并，得到含水印的正方形图像块f′(x,y)：

f′(x,y)=f_rest(x,y)+f_Z′(x,y)。

6.根据权利要求1所述的基于图像特征点全局校正的鲁棒水印方法，其中步骤(2a)所述的利用原始图像I的SIFT特征点集F和受攻击图像I′的SIFT特征点集F′，根据距离约束做特征点匹配，按如下步骤进行：

(2a1)对于原始图像I特征点集F中的任意一个SIFT特征点，计算其与受攻击图像I′特征点集F′中所有SIFT特征点的欧式距离；

(2a2)设置比例阈值μ=0.6，在欧氏距离中找出距该点的最近距离和次近距离，当最近距离除以次近距离小于μ时，则距该点最近距离对应的点就是与该点相匹配的点。

7.根据权利要求1所述的基于图像特征点全局校正的鲁棒水印方法，其中步骤(3f)所述的通过最小距离解码提取水印

按如下步骤进行：

(3f1)采用抖动量化调制方法，对选择出的L个Zernike矩的幅值

进行量化，量化公式为：

${\left(A_{m_r{n}_r}^{\&prime;}\right)}_Q=\left[\frac{A_{m_r{n}_r}^{\&prime;}-d_r\left(Q\right)}{\mathrm{\&Delta;}}\right]\mathrm{\&Delta;}+d_r\left(Q\right),$ Q=0,1，r=1，…,L，

式中，[·]是四舍五入操作，Δ是量化步长，d_r(Q)是第r个抖动函数，且满足d_r(1)=Δ/2+d_r(0)，抖动矢量(d₁(0),…,d_L(0))是利用与嵌入水印过程相同的密钥Key3产生的，

是

采用抖动函数d_r(Q)量化后的Zernike矩幅值；

通过对每一个Zernike矩的幅值进行量化，得到两组量化式

和

r=1，…,L；

(3f2)将上述与

之间的距离定义为：

将上述

与

之间的距离定义为：

r=1,…,L，通过比较dis0与dis1的大小，提取L位水印信息

提取公式为：

$b_r^{\&prime;}=\underset{Q\&Element;\left\{\mathrm{0,1}\right\}}{\arg \min }{({\left(A_{m_r{n}_r}^{\&prime;}\right)}_Q-A_{m_r{n}_r}^{\&prime;})}^2,$ r=1,…,L，

式中，

是第r个Zernike矩的幅值，argmin操作表示为：如果dis0<dis1，则

否则，r=1,…,L，L是Zernike矩的个数。

Images