一种基于特征区域几何优化的鲁棒可逆水印嵌入和提取方法

摘要

本发明提供了一种基于特征区域几何优化的鲁棒可逆水印嵌入和提取方法，该水印嵌入方法通过获取原始图像的低频子带、获取初始特征点集合、获取候选特征点集合、获取候选特征区域集合、计算选特征区域关联矩阵、计算候选关联权值向量、获取筛选特征区域集合、获取圆形特征区域集合、得到特征图像、得到密钥图像和获取密钥信息的步骤，克服了现有鲁棒可逆水印方法的缺陷，改善了含水印图像的视觉质量，提高了水印嵌入容量，增强了水印抵抗复杂攻击的鲁棒性，提升了水印的安全性，提高了鲁棒可逆图像水印方法的综合性能。

说明书

技术领域

本发明属于多媒体信息安全领域，涉及一种基于特征区域几何优 化的鲁棒可逆水印嵌入和提取方法，可用于网络环境中数字图像的内 容认证、版权保护。

背景技术

数字图像作为国家基础设施建设的支撑性成果之一，在国民经 济、国防建设中占有十分重要的地位和广泛的应用。随着计算机技术 与数码成像设备的发展与普及，数字图像的安全保护已成为多媒体信 息安全领域亟待解决的重要问题之一。近年来，数字水印技术通过在 宿主图像中嵌入水印的方式为这一问题的解决提供了有效途径，然而 传统方法在水印嵌入过程中常常会对宿主图像造成不可逆失真。即使 这些失真难以被人眼察觉，但却影响了它在医学影像、法庭证据、电 子票据、军事及遥感图像等领域的实际应用。因此，可逆水印技术应 运而生，它利用人类感知与数字图像自身冗余，通过可逆数字嵌入的 方式将水印隐藏到宿主图像中，在水印提取之后能（近）无损地恢复 宿主图像内容，确保其后续应用（如病变检测、分类、目标识别）不 受影响。而且，通过对水印的分析，该技术能够确定版权所有者、认 证图像内容以及跟踪侵权行为，从而为数字图像安全保护提供有力的 技术支撑。同时，由于实际应用中数字图像往往会受到有损压缩、噪 声干扰的影响，嵌入的水印又需要具有抵抗攻击的鲁棒性。这种水印 方案称为鲁棒可逆水印，它以其特有的优势引起了国内外研究者们的 广泛关注。

根据水印嵌入模型的不同，现有的鲁棒可逆图像水印方法可以分 为三类。

第一类是基于直方图旋转的方法，该方法是基于Patchwork理论， 利用数字图像相邻像素间的相关性生成质心向量，并沿不同方向旋转 质心向量来嵌入水印，见文献“De Vleeschouwer C,Delaigle J,and  Macq B.Circular interpretation of bijective transformations in lossless watermarking for media asset management.IEEE Trans. Multimedia,5(1):97-105,2003”。该方法虽然对联合图像专家组 JPEG压缩具有鲁棒性，但由于水印图像存在严重的“椒盐”噪声，

显著降低了含水印图像的视觉质量。

第二类是基于直方图分布约束的方法，见文献“Zou D K,Shi Y  Q,Ni Z C,and Su W.A semi-fragile lossless digital  watermarking scheme based on integer wavelet transform.IEEE  Trans.Circuits and Systems for Video Technology,16(10): 1294-1300,2006”和“Ni Z C,Shi Y Q,Ansari N,Su W,Sun Q B, and Lin X.Robust lossless image data hiding designed for  semi-fragile image authentication.IEEE Trans.Circuits and  Systems for Video Technology,18(4):497-509,2008”。此类方 法根据数字图像的直方图分布，有规则地修改图像的统计特性来实现 水印嵌入与提取。与基于直方图旋转的方法相比，该类方法克服了含 水印图像的“椒盐”噪声，提高了水印图像的视觉质量，但却存在容 量低、可逆性与鲁棒性不稳定的缺陷。

第三类方法是基于广义统计量直方图与聚类的方法，见文献“An  L L,Gao X B,Yuan Y,and Tao D C.Robust lossless data hiding  using clustering and statistical quantity histogram. Neurocomputing,77(1):1-11,2012”和“Li X L,Tao D C,Deng  C and Li J.Robust reversible watermarking via clustering and  enhanced pixel-wise masking.IEEE Transactions on Image  Processing,21(8):3598-3611,2012”。该类方法基于广义统计量 直方图平移与聚类算法实现了水印的无损嵌入与鲁棒提取，增强了水 印抗攻击的鲁棒性。虽然该类方法与前两类方法相比，提升了水印抵 抗JPEG压缩与加性高斯噪声的鲁棒性，但其抵抗复杂攻击的鲁棒性 欠缺，而且在容量、不可感知性及鲁棒性三方面的综合性能有待进一 步提高。

发明内容

本发明目的是增强目前数字图像水印技术对于复杂攻击的鲁棒 性的抵抗性能，提高鲁棒可逆水印方法在容量、不可感知性及鲁棒性 三方面的综合性能，提供一种基于特征区域几何优化的鲁棒可逆水印 嵌入和提取方法。

为达上述目的，本发明技术方案是：一种基于特征区域几何优化 的鲁棒可逆水印嵌入方法，包括如下步骤：

1）、获取原始图像的低频子带：对原始图像I进行三级提升小波 变换，得到其第三级小波分解尺度下的低频子带R_3,a；

2）、获取初始特征点集合：利用多尺度Harris特征检测算子 Harris-Affine在低频子带R_3,a中提取仿射不变特征点，得到初始特 征点集合A：利用多尺度Harris特征检测算子Harris-Affine在低 频子带R_3,a中提取仿射不变特征点，得到初始特征点集合 A＝{A_i,i＝1,2,…m₁}，其中A_i表示第i个特征点，初始特征点集合A的五 个属性包括特征点的横坐标特征点的纵坐标特征点的强度 特征点的特征尺度特征点椭圆参数与m₁表示初 始特征点的个数；

3）、获取候选特征点集合：选取初始特征点集合A中满足特征尺 度向量ξ要求的中间尺度特征点，得到候选特征点集合B：将初始特 征点集合A中每个特征点的特征尺度与特征尺度向量ξ＝[ξ₁,ξ₂]进 行比较，选取满足条件的中间尺度特征点，得到候选特征 点集合B＝{B_i,i＝1,2,…m₂},其中m₂表示候选特征点的个数；

4）、获取候选特征区域集合：利用候选特征点集合B中的每个特 征点，分别构造其仿射协变特征区域，得到候选特征区域集合S：对 候选特征点集合B中的每个特征点B_i，利用公式

$B_i^a{(x-B_i^x)}^2+{2B}_i^b(x-B_i^x)(y-B_i^y)+B_i^c{(y-B_i^y)}^2\le 1$

构造其仿射协变特征区域，得到候选特征区域集合S＝{S_i,i＝1,2,…m₂}， 其中x与y分别表示低频子带R_3,a中任意满足上述不等式约束的元素 横坐标与纵坐标，表示特征点B_i的横坐标，表示特征点B_i的纵坐 标，与表示特征点B_i椭圆参数；

5）、计算候选特征区域关联矩阵：根据候选特征区域集合S，计 算候选特征区域关联矩阵P；

6）、计算候选关联权值向量：利用候选特征区域关联矩阵P，计 算候选关联权值向量L；

7）、获取筛选特征区域集合：从候选关联权值向量L中查找出最 大权值的位置索引k，更新候选特征区域关联矩阵P和候选关联权值 向量L，得到筛选特征区域集合H；

8）、获取圆形特征区域集合：对筛选特征区域集合H进行归一化 处理，得到圆形特征区域集合Q＝{Q_i,i＝1,2,…m₃}，m₃表示筛选特征区 域的个数；

9）、得到特征图像：利用圆形特征区域集合Q对低频子带R_3,a进行 系数置零化处理，得到特征图像C；

10）、得到密钥图像：将灰度水印图像W和特征图像C进行按位 异或运算，得到密钥图像D；

11）、获取密钥信息：利用可逆元胞自动机对密钥图像D与特征 尺度向量ξ进行加密得到密钥信息G。

上述步骤5）具体包括以下子步骤：

5.1)计算候选特征区域集合S中每个候选特征区域S_i对应特征 点的二阶矩阵G_i，i＝1,2,…m₂：

$G_i=\left(\begin{array}{c}{S}_i^a,S_i^b\\ S_i^b,S_i^c\end{array}\right)$

其中，与表示候选特征区域S_i对应特征点的椭圆参数， i＝1,2,…m₂，m₂表示候选特征区域的个数；

5.2)计算候选特征区域集合S中每个候选特征区域S_i的长轴长 度l_i，计算公式为：

E_i＝f(G_i)

$l_i=1/\sqrt{E_i\left(1\right)}$

其中，f(·)计算二阶矩矩阵的特征值矩阵E_i，E_i(1)取E_i中的第一个元 素；

5.3)计算候选特征区域集合S的距离矩阵DIS，其中第i个候选 特征点和第候选j个特征点的距离表示为：

$\mathrm{DIS}(i,j)=\sqrt{{(S_i^x-S_j^x)}^2+{(S_i^y-S_j^y)}^2}$

式中，与分别表示第i个候选特征点和第j个候选特征点的横坐 标，与分别表示第i个候选特征点和第j个候选特征点的纵坐标， i,j＝1,2,…m₂，m₂表示候选特征区域个数；

5.4)计算候选特征区域关联矩阵P，其第(i,j)位置上的元素计算 公式如下：

$P(i,j)=\left(\begin{array}{cc}0& \mathrm{DIS}(i,j)\ge (l_i+l_j)/2\\ 1& \mathrm{DIS}(i,j)<(l_i+l_j)/2\end{array}\right)$

其中，DIS(i,j)表示第i个候选特征点和第j个候选特征点的距离，l_i和 l_j分别表示第i个候选特征点和第j个候选特征点的长轴长度。

上述步骤6）的具体过程是：利用候选特征区域关联矩阵P，计算 候选关联权值向量L，其第i个候选特征点的关联权值计算规则如下：

$L\left(i\right)=\left(\underset{j=1}{\overset{m_2}{\Sigma }}P(i,j)\right)/S_i^u$

其中，表示第i个候选特征点的特征强度，P(i,j)表示候选特征区域 关联矩阵P中第(i,j)位置上的元素，m₂表示候选特征区域的个数。

上述步骤7）具体包括以下步骤：

7.1)从候选关联权值向量L中查找出最大权值的位置索引k，利 用如下公式更新候选特征区域关联矩阵P和候选关联权值向量L；

$L\left(i\right)=L\left(i\right)-P(i,k)\times (1/S_i^u)$

P(i,k)＝0,P(k,i)＝0

其中表示第i个候选特征点的特征强度，L(i)表示第i个候选特征点 的关联权值，P(i,k)与P(k,i)分别表示候选关联矩阵P中第(i,k)与第 (k,i)位置上的元素，i＝1,2,…m₂，m₂表示候选特征区域的个数；

7.2)重复步骤7.1)直到候选关联权值向量L中所有元素都为0， 得到筛选特征区域集合H＝{H_i,i＝1,2,…m₃}，m₃表示筛选特征区域的个 数。

上述步骤11）的具体过程是：

11.1)依次将密钥图像D中每个灰度值转换成8位二进制序列， 其中密钥图像D中第(i,j)位置的灰度值D(i,j)所转换成的二进制序列 表示为$(D_{i,j}^{\left(1\right)},D_{i,j}^{\left(2\right)},...,D_{i,j}^{\left(8\right)}),$这里$D_{i,j}^{\left(\lambda \right)}\in \left\{\mathrm{0,1}\right\},$λ＝1,2,…8，

11.2)按照行扫描顺序将密钥图像D中每个灰度值所转换成的二 进制序列进行连接，得到长度为的密钥图像二进制序 列D_b，表示为：

11.3)将特征尺度向量ξ＝[ξ₁,ξ₂]转换成8位二进制序列，得到原 始特征尺度向量二进制序列${\xi }_b=({\xi }_1^{\left(1\right)},{\xi }_1^{\left(2\right)},...{\xi }_1^{\left(8\right)},{\xi }_2^{\left(1\right)},{\xi }_2^{\left(2\right)},...{\xi }_2^{\left(8\right)});$

11.4)采用补零方式将原始特征尺度向量二进制序列ξ_b进行零填 充，生成特征尺度向量二进制序列 ${\xi }_b^{\prime }=({\xi }_1^{\left(1\right)},{\xi }_1^{\left(2\right)},...{\xi }_1^{\left(8\right)},{\xi }_2^{\left(1\right)},{\xi }_2^{\left(2\right)},...{\xi }_2^{\left(8\right)},\mathrm{0,0},...0),$它的长度与步骤11.2)得到的 密钥图像二进制序列Db相同，均为

11.5)将步骤11.2)得到的密钥图像二进制序列D_b与步骤11.4) 得到的特征尺度向量二进制序列ξ_b′作为可逆元胞自动机的初始状态， 使用可逆元胞自动机加密方法对ξ_b′进行加密，得到密钥信息G和共享 特征向量Γ，其中可逆元胞自动机规则为41R，迭代次数为10。

上述的鲁棒可逆水印嵌入方法，还包括步骤12）、将密钥信息G 注册到知识产权信息数据库中。

同时，本发明提供了一种基于特征区域几何优化的鲁棒可逆水印 提取方法，包括如下步骤：

a)、利用可逆元胞自动机对密钥信息G进行解密，得到密钥图像 D与特征尺度向量ξ；

b)、将待检测图像I′进行三级提升小波变换，得到其第三级小波 分解尺度下的低频子带R′_3,a：

将大小为M×N的原始图像I′进行三级提升小波变换，得到一组小 波分解子带序列R′＝{R′_t,a,R′_t,h,R′_t,v,R′_t,d}，其中，整数t为分解尺度，1≤t≤3， R′_t,a为第t级小波分解尺度下的低频子带，R′_t,h为第t级小波分解尺度下 的水平子带，R′_t,v为第t级小波分解尺度下的垂直子带，R′_t,d为第t级小 波分解尺度下的对角子带，从子带序列R′中选择大小为的第三级小波分解尺度下的低频子带R′_3,a，其中表示选取比 M/8大的最小整数，表示选取比N/8大的最小整数；

c)、利用多尺度Harris特征检测算子Harris-Affine在低频子 带R′_3,a中提取仿射不变特征点，得到初始特征点集合A′：

利用多尺度Harris特征检测算子Harris-Affine在低频子带R′_3,a中提取仿射不变特征点，得到初始特征点集合A′＝{A_i′,i＝1,2,…m₁}，其 中A_i′表示第i个特征点，它的五个属性包括特征点的横坐标特征 点的纵坐标特征点的强度特征点的尺度特征点椭圆参 数与m1表示初始特征点的个数；

d)、选取初始特征点集合A′中满足特征尺度向量ξ要求的中间尺 度特征点，得到候选特征点集合B′：将初始特征点集合A′中每个特征 点的特征尺度与特征尺度向量ξ＝[ξ₁,ξ₂]进行比较，选取满足 条件的中间尺度特征点，得到候选特征点集合 B′＝{B_i′,i＝1,2,…m₂},其中m₂表示候选特征点的个数。

e)、利用候选特征点集合B′中的每个特征点，分别构造其仿射协 变特征区域，得到候选特征区域集合S′：

对候选特征点集合B′中的每个特征点B_i′，利用公式

$B_i^{\prime a}{(x-B_i^{\prime x})}^2+{2B}_i^{\prime b}(x-B_i^{\prime x})(y-B_i^{\prime y})+B_i^{\prime c}{(y-B_i^{\prime y})}^2\le 1$

构造其仿射协变特征区域，得到候选特征区域集合 S′＝{S_i′,i＝1,2,…m₂}，其中x与y分别表示低频子带R′_3,a中任意满足上述 不等式约束的元素横坐标与纵坐标，表示特征点B_i′的横坐标，表 示特征点B_i′的纵坐标，与表示特征点B_i′椭圆参数；

f)、根据候选特征区域集合S′，计算候选特征区域关联矩阵P′；

g)、利用候选特征区域关联矩阵P′，计算候选关联权值向量L′；

h)、从候选关联权值向量L′中查找出最大权值的位置索引k′，更 新候选特征区域关联矩阵P′和候选关联权值向量L′，得到筛选特征区 域集合H′；

i)、对筛选特征区域集合H′进行归一化处理，得到圆形特征区域 集合Q′＝{Q_i′,i＝1,2,…m₃}，m₃表示筛选特征区域的个数；

j)、利用圆形特征区域集合Q′对低频子带R′_3,a进行系数置零化处 理，得到特征图像C′：

在低频子带R′_3,a中将圆形特征区域集合Q′以外区域的系数值设置 为0，仅保留圆形特征区域中的系数值，进而得到特征图像C′。

k)、将特征图像C′与密钥图像D进行按位异或运算，得到水印图 像W′。

利用如下公式将特征图像C′与密钥图像D进行按位异或运算，得 到水印图像W′：

W′(i,j)＝C′(i,j)∧D(i,j)

其中，D(i,j)、C′(i,j)、W′(i,j)分别为密钥图像D、特征图像C′与水 印图像W′在(i,j)处的元素值，∧表示按位异或运算操作，

上述步骤a）包括以下子步骤：

a1)将密钥信息G和共享特征向量Γ作为可逆元胞自动机的初始 状态，利用可逆元胞自动机解密方法对其进行解密，得到密钥图像二 进制序列D_b和特征尺度向量二进制序列ξ_b′，其中可逆元胞自动机规则 为41R，迭代次数为10。

a2)从特征尺度向量二进制序列ξ_b′中首先提取出第1到第8位二 进制值将其转换成十进制数，得到特征尺度向量ξ中 的第一个元素ξ₁；然后从特征尺度向量二进制序列ξ_b′中提取出第9到 第16位二进制值将其转换成十进制数，得到特征尺 度向量ξ中的第二个元素ξ₂，进而得到特征尺度向量ξ＝[ξ₁,ξ₂]；

a3)将步骤a1)得到的密钥图像二进制序列D_b从第一个元素开始 依次按照每8位一组的方式将其转换成十进制数，并按照行扫描的方 式将得到的十进制数序列排列成的矩阵，进而得到密钥 图像D。

上述步骤f）包括以下子步骤：

f1)计算候选特征区域集合S′中每个候选特征区域S_i′对应特征 点的二阶矩矩阵G_i′，i＝1,2,…m₂：

$G_i^{\prime }=\left(\begin{array}{c}{S}_i^{\prime a},S_i^{\prime b}\\ S_i^{\prime b},S_i^{\prime c}\end{array}\right)$

其中，与表示候选特征区域S_i′对应特征点的椭圆参数， i＝1,2,…m₂，m₂表示候选特征区域的个数；

f2)计算候选特征区域集合S′中每个候选特征区域S_i′的长轴长 度l_i′，计算公式为：

E_i′＝f(G_i′)

$l_i^{\prime }=1/\sqrt{E_i^{\prime }\left(1\right)}$

其中，f(·)计算二阶矩矩阵的特征值矩阵E_i′，E_i′(1)取E_i′中的第一个元 素；

f3)计算候选特征区域集合S′的距离矩阵DIS′，其中第i个候选特 征点和第候选j个特征点的距离表示为

${\mathrm{DIS}}^{\prime }(i,j)=\sqrt{{(S_i^{\prime x}-S_j^{\prime x})}^2+{(S_i^{\prime y}-S_j^{\prime y})}^2}$

式中，与分别表示第i个候选特征点和第j个候选特征点的横坐 标，与分别表示第i个候选特征点和第j个候选特征点的纵坐 标，i,j＝1,2,…m₂，m₂表示候选特征区域个数；

f4)计算候选特征区域关联矩阵P′，其第(i,j)位置上的元素计算 公式如下：

$P^{\prime }(i,j)=\left(\begin{array}{cc}0& {\mathrm{DIS}}^{\prime }(i,j)\ge (l_i^{\prime }+l_j^{\prime })/2\\ 1& {\mathrm{DIS}}^{\prime }(i,j)<(l_i^{\prime }+l_j^{\prime })/2\end{array}\right)$

其中，DIS′(i,j)表示第i个候选特征点和第j个候选特征点的距离，l_i′和 l_j′分别表示第i个候选特征点和第j个候选特征点的长轴长度。

上述步骤g）的具体过程是：利用候选特征区域关联矩阵P′，计 算候选关联权值向量L′，其第i个候选特征点的关联权值计算规则如 下：

$L^{\prime }\left(i\right)=\left(\underset{j=1}{\overset{m_2}{\Sigma }}{P}^{\prime }(i,j)\right)/S_i^{\prime u}$

其中，表示第i个候选特征点的特征强度，P′(i,j)表示候选特征区 域关联矩阵P′中第(i,j)位置上的元素，m₂表示候选特征区域的个数。

上述步骤h）包括以下子步骤：

h1)从候选关联权值向量L′中查找出最大权值的位置索引k′，利 用如下公式更新候选特征区域关联矩阵P′和候选关联权值向量L′；

$L^{\prime }\left(i\right)=L^{\prime }\left(i\right)-P^{\prime }(i,k^{\prime })\times (1/S_i^{\prime u})$

P′(i,k′)＝0,P′(k′,i)＝0

其中表示第i个候选特征点的特征强度，L′(i)表示第i个候选特征点 的关联权值，P′(i,k′)与P′(k′,i)分别表示候选关联矩阵P′中第(i,k′)与第 (k′,i)位置上的元素，i＝1,2,…m₂，m₂表示候选特征区域的个数；

h2)重复步骤h1)直到候选关联权值向量L′中所有元素都为0， 得到筛选特征区域集合H′＝{H_i′,i＝1,2,…m₃}，m₃表示筛选特征区域的 个数。

与现有技术相比，本发明具有以下优点：

(1)本发明由于利用零水印机制设计鲁棒可逆水印的嵌入与提 取模型，其“不嵌入性”有效避免了水印嵌入过程中的像素溢出，提 升了含水印图像的视觉质量；

(2)本发明由于将原始图像进行三级提升小波变换，并利用生成 的低频子带进行水印嵌入，大大提升了水印嵌入容量；

(3)本发明由于采用特征区域几何优化的方法进行特征区域的 筛选，降低了现有基于最小支撑树方法的复杂度，提升了水印抵抗复 杂攻击的鲁棒性；

(4)本发明由于采用元胞自动机加密与解密算法进行密钥图像 与特征尺度的加密与解密，有效提升了水印的安全性；

(5)本发明由于改善了含水印图像的视觉质量，提高了水印嵌入 容量，增强了水印抵抗复杂攻击的鲁棒性，提升了水印的安全性，进 而提高了鲁棒可逆水印嵌入与提取方法的综合性能。

具体实施方式

为增强目前数字图像水印技术对于复杂攻击的鲁棒性的抵抗性 能，提高鲁棒可逆水印方法在容量、不可感知性及鲁棒性三方面的综 合性能，本实施例提供了一种基于特征区域几何优化的鲁棒可逆水印 嵌入方法，其实现步骤如下：

步骤1：将原始图像I进行三级提升小波变换，得到其第三级小 波分解尺度下的低频子带R_3,a：

将大小为M×N的原始图像I进行三级提升小波变换，得到一组小 波分解子带序列R＝{R_t,a,R_t,h,R_t,v,R_t,d}，其中，整数t为分解尺度，1≤t≤3， R_t,a为第t级小波分解尺度下的低频子带，R_t,h为第t级小波分解尺度下 的水平子带，R_t,v为第t级小波分解尺度下的垂直子带，R_t,d为第t级小 波分解尺度下的对角子带，从子带序列R中选择大小为的第三级小波分解尺度下的低频子带R_3,a，其中表示选取比 M/8大的最小整数，表示选取比N/8大的最小整数。

步骤2：利用多尺度Harris特征检测算子Harris-Affine在低频 子带R_3,a中提取仿射不变特征点，得到初始特征点集合A；

利用多尺度Harris特征检测算子Harris-Affine在低频子带R_3,a中提取仿射不变特征点，得到初始特征点集合A＝{A_i,i＝1,2,…m₁}，其 中A_i表示第i个特征点，它的五个属性包括特征点的横坐标特征 点的纵坐标特征点的强度特征点的尺度特征点椭圆参 数与m₁表示初始特征点的个数。

步骤3：选取初始特征点集合A中满足特征尺度向量ξ要求的中间 尺度特征点，得到候选特征点集合B；

将初始特征点集合A中每个特征点的特征尺度与特征尺度向量 ξ＝[ξ₁,ξ₂]进行比较，选取满足条件的中间尺度特征点，得到 候选特征点集合B＝{B_i,i＝1,2,…m₂},其中m₂表示候选特征点的个数。

步骤4：利用候选特征点集合B中的每个特征点，分别构造其仿 射协变特征区域，得到候选特征区域集合S；

对候选特征点集合B中的每个特征点B_i，利用公式

$B_i^a{(x-B_i^x)}^2+{2B}_i^b(x-B_i^x)(y-B_i^y)+B_i^c{(y-B_i^y)}^2\le 1$

构造其仿射协变特征区域，得到候选特征区域集合S＝{S_i,i＝1,2,…m₂}， 其中x与y分别表示低频子带R_3,a中任意满足上述不等式约束的元素 横坐标与纵坐标，表示特征点B_i的横坐标，表示特征点B_i的纵坐 标，与表示特征点B_i椭圆参数。

步骤5：根据候选特征区域集合S，计算候选特征区域关联矩阵 P；

5.1)计算候选特征区域集合S中每个候选特征区域S_i对应特征 点的二阶矩矩阵G_i，i＝1,2,…m₂：

$G_i=\left(\begin{array}{c}{S}_i^a,S_i^b\\ S_i^b,S_i^c\end{array}\right)$

其中，与表示候选特征区域S_i对应特征点的椭圆参数， i＝1,2,…m₂，m₂表示候选特征区域的个数；

5.2)计算候选特征区域集合S中每个候选特征区域S_i的长轴长 度l_i，计算公式为：

E_i＝f(G_i)

$l_i=1/\sqrt{E_i\left(1\right)}$

其中，f(·)计算二阶矩矩阵的特征值矩阵E_i，E_i(1)取E_i中的第一个元 素；

5.3)计算候选特征区域集合S的距离矩阵DIS，其中第i个候选 特征点和第候选j个特征点的距离表示为

$\mathrm{DIS}(i,j)=\sqrt{{(S_i^x-S_j^x)}^2+{(S_i^y-S_j^y)}^2}$

式中，与分别表示第i个候选特征点和第j个候选特征点的横坐 标，与分别表示第i个候选特征点和第j个候选特征点的纵坐标， i,j＝1,2,…m₂，m₂表示候选特征区域个数；

5.4)计算候选特征区域关联矩阵P，其第(i,j)位置上的元素计算 公式如下：

$P(i,j)=\left(\begin{array}{cc}0& \mathrm{DIS}(i,j)\ge (l_i+l_j)/2\\ 1& \mathrm{DIS}(i,j)<(l_i+l_j)/2\end{array}\right)$

其中，DIS(i,j)表示第i个候选特征点和第j个候选特征点的距离，l_i和 l_j分别表示第i个候选特征点和第j个候选特征点的长轴长度。

步骤6：利用候选特征区域关联矩阵P，计算候选关联权值向量 L；

利用候选特征区域关联矩阵P，计算候选关联权值向量L，其第i 个候选特征点的关联权值计算规则如下：

$L\left(i\right)=\left(\underset{j=1}{\overset{m_2}{\Sigma }}P(i,j)\right)/S_i^u$

其中，表示第i个候选特征点的特征强度，P(i,j)表示候选特征区域 关联矩阵P中第(i,j)位置上的元素，m₂表示候选特征区域的个数；

步骤7：从候选关联权值向量L中查找出最大权值的位置索引k， 更新候选特征区域关联矩阵P和候选关联权值向量L，得到筛选特征 区域集合H；

7.1)从候选关联权值向量L中查找出最大权值的位置索引k，利 用如下公式更新候选特征区域关联矩阵P和候选关联权值向量L；

$L\left(i\right)=L\left(i\right)-P(i,k)\times (1/S_i^u)$

P(i,k)＝0,P(k,i)＝0

其中表示第i个候选特征点的特征强度，L(i)表示第i个候选特征点 的关联权值，P(i,k)与P(k,i)分别表示候选关联矩阵P中第(i,k)与第 (k,i)位置上的元素，i＝1,2,…m₂，m₂表示候选特征区域的个数；

7.2)重复步骤7.1)直到候选关联权值向量L中所有元素都为0， 得到筛选特征区域集合H＝{H_i,i＝1,2,…m₃}，m₃表示筛选特征区域的个 数。

步骤8：对筛选特征区域集合H进行归一化处理，得到圆形特征 区域集合Q；

利用文献“BaumbergA.Reliablefeaturematchingacross  widelyseparatedviews.In:ProceedingsofInternational  ConferenceonComputerVisionandPatternRecognition.Hilton  Head,USA:IEEE,2000.774-781”中的方法对筛选特征区域集合H 中的每个筛选特征区域利用公式进行归一化处理，得 到圆形特征区域集合Q＝{Q_i,i＝1,2,…m₃}，其中，G_i表示筛选特征区域 集合H中第i个筛选特征区域H_i对应特征点的二阶矩矩阵，x和x^*分 别表示第i个筛选特征区域H_i中任一元素映射前后的横坐标，y和y^*分别表示第i个筛选特征区域H_i中任一元素映射前后的纵坐标，上标 T表示矩阵转置运算符，m₃表示筛选特征区域的个数。

步骤9：利用圆形特征区域集合Q对低频子带R_3,a进行系数置零 化处理，得到特征图像C；

在低频子带R_3,a中将圆形特征区域集合Q以外区域的系数值设置 为0，仅保留圆形特征区域中的系数值，进而得到特征图像C。

步骤10：将灰度水印图像W和特征图像C进行按位异或运算， 得到密钥图像D；

选取大小为的灰度水印图像W，并利用如下公式将 其与特征图像C进行按位异或运算，得到密钥图像D：

D(i,j)＝C(i,j)∧W(i,j)

其中，D(i,j)、C(i,j)、W(i,j)分别为密钥图像D、特征图像C与水印 图像W在(i,j)处的元素值，∧表示按位异或运算操作，

步骤11：利用可逆元胞自动机对密钥图像D与特征尺度向量ξ进 行加密得到密钥信息G：

11.1)依次将密钥图像D中每个灰度值转换成8位二进制序列， 其中密钥图像D中第(i,j)位置的灰度值D(i,j)所转换成的二进制序列 表示为$(D_{i,j}^{\left(1\right)},D_{i,j}^{\left(2\right)},...,D_{i,j}^{\left(8\right)}),$这里$D_{i,j}^{\left(\lambda \right)}\in \left\{\mathrm{0,1}\right\},$λ＝1,2,…8，

11.2)按照行扫描顺序将密钥图像D中每个灰度值所转换成的二 进制序列进行连接，得到长度为的密钥图像二进制序 列D_b，表示为：

11.3)将特征尺度向量ξ＝[ξ₁,ξ₂]转换成8位二进制序列，得到原 始特征尺度向量二进制序列${\xi }_b=({\xi }_1^{\left(1\right)},{\xi }_1^{\left(2\right)},...{\xi }_1^{\left(8\right)},{\xi }_2^{\left(1\right)},{\xi }_2^{\left(2\right)},...{\xi }_2^{\left(8\right)});$

11.4)采用补零方式将原始特征尺度向量二进制序列ξ_b进行零填 充，生成特征尺度向量二进制序列 ${\xi }_b^{\prime }=({\xi }_1^{\left(1\right)},{\xi }_1^{\left(2\right)},...{\xi }_1^{\left(8\right)},{\xi }_2^{\left(1\right)},{\xi }_2^{\left(2\right)},...{\xi }_2^{\left(8\right)},\mathrm{0,0},...0),$它的长度与步骤11.2)得到的 密钥图像二进制序列D_b相同，均为

11.5)将步骤11.2)得到的密钥图像二进制序列D_b与步骤11.4) 得到的特征尺度向量二进制序列ξ_b′作为可逆元胞自动机的初始状态， 利用文献“冀峰,安玲玲,邓成,高新波.基于多重元胞自动机的图 像水印加密算法.自动化学报,38(11):1824-1830,2012”中所使 用的可逆元胞自动机加密方法对其进行加密，得到密钥信息G和共享 特征向量Γ，其中可逆元胞自动机规则为41R，迭代次数为10。

通过上述步骤1～步骤11可以实现水印的嵌入，得到密钥信息 G；然后将其注册到知识产权(Intellectual Property Right,简称 IPR)信息数据库中来保护版权。

同时，本实施例提供了一种基于特征区域几何优化的鲁棒可逆水 印提取方法，其实现步骤如下：

步骤A：利用可逆元胞自动机对密钥信息G进行解密，得到密钥 图像D与特征尺度向量ξ；

A1)将密钥信息G和共享特征向量Γ作为可逆元胞自动机的初始 状态，利用文献“冀峰,安玲玲,邓成,高新波.基于多重元胞自动 机的图像水印加密算法.自动化学报,38(11):1824-1830,2012” 中所使用的可逆元胞自动机解密方法对其进行解密，得到密钥图像二 进制序列D_b和特征尺度向量二进制序列ξ_b′，其中可逆元胞自动机规则 为41R，迭代次数为10。

A2)从特征尺度向量二进制序列ξ_b′中首先提取出第1到第8位二 进制值将其转换成十进制数，得到特征尺度向量ξ中 的第一个元素ξ₁；然后从特征尺度向量二进制序列ξ_b′中提取出第9到 第16位二进制值将其转换成十进制数，得到特征尺 度向量ξ中的第二个元素ξ₂，进而得到特征尺度向量ξ＝[ξ₁,ξ₂]；

A3)将步骤A1)得到的密钥图像二进制序列D_b从第一个元素开始 依次按照每8位一组的方式将其转换成十进制数，并按照行扫描的方 式将得到的十进制数序列排列成的矩阵，进而得到密钥 图像D；

步骤B：将待检测图像I′进行三级提升小波变换，得到其第三级 小波分解尺度下的低频子带R′_3,a；

将大小为M×N的原始图像I′进行三级提升小波变换，得到一组小 波分解子带序列R′＝{R′_t,a,R′_t,h,R′_t,v,R′_t,d}，其中，整数t为分解尺度，1≤t≤3， R′_t,a为第t级小波分解尺度下的低频子带，R′_t,h为第t级小波分解尺度下 的水平子带，R′_t,v为第t级小波分解尺度下的垂直子带，R′_t,d为第t级小 波分解尺度下的对角子带，从子带序列R′中选择大小为 的第三级小波分解尺度下的低频子带R′_3,a，其中表示选取比 M/8大的最小整数，表示选取比N/8大的最小整数。

步骤C：利用多尺度Harris特征检测算子Harris-Affine在低 频子带R′_3,a中提取仿射不变特征点，得到初始特征点集合A′；

利用多尺度Harris特征检测算子Harris-Affine在低频子带R′_3,a中提取仿射不变特征点，得到初始特征点集合A′＝{A_i′,i＝1,2,…m₁}，其 中A_i′表示第i个特征点，它的五个属性包括特征点的横坐标特征 点的纵坐标特征点的强度特征点的尺度特征点椭圆参 数与m₁表示初始特征点的个数。

步骤D：选取初始特征点集合A′中满足特征尺度向量ξ要求的中 间尺度特征点，得到候选特征点集合B′；

将初始特征点集合A′中每个特征点的特征尺度与步骤A2)得到 的特征尺度向量ξ＝[ξ₁,ξ₂]进行比较，选取满足条件的中间尺 度特征点，得到候选特征点集合B′＝{B_i′,i＝1,2,…m₂},其中m₂表示候选 特征点的个数。

步骤E：利用候选特征点集合B′中的每个特征点，分别构造其仿 射协变特征区域，得到候选特征区域集合S′；

对候选特征点集合B′中的每个特征点B_i′，利用公式

$B_i^{\prime a}{(x-B_i^{\prime x})}^2+{2B}_i^{\prime b}(x-B_i^{\prime x})(y-B_i^{\prime y})+B_i^{\prime c}{(y-B_i^{\prime y})}^2\le 1$

构造其仿射协变特征区域，得到候选特征区域集合 S′＝{S_i′,i＝1,2,…m₂}，其中x与y分别表示低频子带R′_3,a中任意满足上述 不等式约束的元素横坐标与纵坐标，表示特征点B_i′的横坐标，表 示特征点B_i′的纵坐标，与表示特征点B_i′椭圆参数。

步骤F：根据候选特征区域集合S′，计算候选特征区域关联矩阵 P′；

F1)计算候选特征区域集合S′中每个候选特征区域S_i′对应特征 点的二阶矩矩阵G_i′，i＝1,2,…m₂：

$G_i^{\prime }=\left(\begin{array}{c}{S}_i^{\prime a},S_i^{\prime b}\\ S_i^{\prime b},S_i^{\prime c}\end{array}\right)$

其中，与表示候选特征区域S_i′对应特征点的椭圆参数， i＝1,2,…m₂，m₂表示候选特征区域的个数；

F2)计算候选特征区域集合S′中每个候选特征区域S_i′的长轴长 度l_i′，计算公式为：

E_i′＝f(G_i′)

$l_i^{\prime }=1/\sqrt{E_i^{\prime }\left(1\right)}$

其中，f(·)计算二阶矩矩阵的特征值矩阵E_i′，E_i′(1)取E_i′中的第一个元 素；

F3)计算候选特征区域集合S′的距离矩阵DIS′，其中第i个候选特 征点和第候选j个特征点的距离表示为

${\mathrm{DIS}}^{\prime }(i,j)=\sqrt{{(S_i^{\prime x}-S_j^{\prime x})}^2+{(S_i^{\prime y}-S_j^{\prime y})}^2}$

式中，与分别表示第i个候选特征点和第j个候选特征点的横坐 标，与分别表示第i个候选特征点和第j个候选特征点的纵坐 标，i,j＝1,2,…m₂，m₂表示候选特征区域个数；

F4)计算候选特征区域关联矩阵P′，其第(i,j)位置上的元素计算 公式如下：

$P^{\prime }(i,j)=\left(\begin{array}{cc}0& {\mathrm{DIS}}^{\prime }(i,j)\ge (l_i^{\prime }+l_j^{\prime })/2\\ 1& {\mathrm{DIS}}^{\prime }(i,j)<(l_i^{\prime }+l_j^{\prime })/2\end{array}\right)$

其中，DIS′(i,j)表示第i个候选特征点和第j个候选特征点的距离，l_i′和 l_j′分别表示第i个候选特征点和第j个候选特征点的长轴长度。

步骤G：利用候选特征区域关联矩阵P′，计算候选关联权值向量 L′；

利用候选特征区域关联矩阵P′，计算候选关联权值向量L′，其第 i个候选特征点的关联权值计算规则如下：

$L^{\prime }\left(i\right)=\left(\underset{j=1}{\overset{m_2}{\Sigma }}{P}^{\prime }(i,j)\right)/S_i^{\prime u}$

其中，表示第i个候选特征点的特征强度，P′(i,j)表示候选特征区 域关联矩阵P′中第(i,j)位置上的元素，m₂表示候选特征区域的个数；

步骤H：从候选关联权值向量L′中查找出最大权值的位置索引 k′，更新候选特征区域关联矩阵P′和候选关联权值向量L′，得到筛选 特征区域集合H′；

H1)从候选关联权值向量L′中查找出最大权值的位置索引k′，利 用如下公式更新候选特征区域关联矩阵P′和候选关联权值向量L′；

$L^{\prime }\left(i\right)=L^{\prime }\left(i\right)-P^{\prime }(i,k^{\prime })\times (1/S_i^{\prime u})$

P′(i,k′)＝0,P′(k′,i)＝0

其中表示第i个候选特征点的特征强度，L′(i)表示第i个候选特征点 的关联权值，P′(i,k′)与P′(k′,i)分别表示候选关联矩阵P′中第(i,k′)与第 (k′,i)位置上的元素，i＝1,2,…m₂，m₂表示候选特征区域的个数；

H2)重复步骤H1)直到候选关联权值向量L′中所有元素都为0， 得到筛选特征区域集合H′＝{H_i′,i＝1,2,…m₃}，m₃表示筛选特征区域的 个数。

步骤I：对筛选特征区域集合H′进行归一化处理，得到圆形特 征区域集合Q′；

利用文献“Baumberg A.Reliable feature matching across  widely separated views.In:Proceedings of International  Conference on Computer Vision and Pattern Recognition.Hilton  Head,USA:IEEE,2000.774-781”中的方法对筛选特征区域集合H′ 中的每个筛选特征区域利用公式进行归一化处理， 得到圆形特征区域集合Q′＝{Q_i′,i＝1,2,…m₃}，其中，G_i′表示筛选特征区 域集合H′中第i个筛选特征区域H_i′对应特征点的二阶矩矩阵，x和x^*分别表示第i个筛选特征区域H_i′中任一元素映射前后的横坐标，y和 y^*分别表示第i个筛选特征区域H_i′中任一元素映射前后的纵坐标，上 标T表示矩阵转置运算符，m₃表示筛选特征区域的个数。

步骤J：利用圆形特征区域集合Q′对低频子带R′_3,a进行系数置零 化处理，得到特征图像C′；

在低频子带R′_3，a中将圆形特征区域集合Q′以外区域的系数值设置 为0，仅保留圆形特征区域中的系数值，进而得到特征图像C′。

步骤K：将特征图像C′与密钥图像D进行按位异或运算，得到水 印图像W′。

利用如下公式将特征图像C′与步骤A3)得到的密钥图像D进行按 位异或运算，得到水印图像W′：

W′(i,j)＝C′(i,j)∧D(i,j)

其中，D(i,j)、C′(i,j)、W′(i,j)分别为密钥图像D、特征图像C′与水 印图像W′在(i,j)处的元素值，∧表示按位异或运算操作，

通过上述步骤A～步骤K可以实现水印提取，从待检测图像I′中提 取出水印图像W′。

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

1.实验条件与实验说明

实现本发明的软件环境为美国Mathworks公司开发的MATLAB  R2013a，实验中的灰度图像选自CVG-UGR图像数据库的10张 512×512×8图像，它们的名称分别为13(house).png、33b.png、 46b.png、57b.png、66b.png、97b.png、cmpndd.png、colomtn.png、 malight.png和porthead.png。将本发明方法分别与Gao等人提出的 基于最小支撑树的特征区域选择方法（见文献“Gao X B,DengC,Li  X L,Tao D C.Geometric distortion insensitive image  watermarking in affine covariant regions.IEEE Transactions on  Systems,Man,and Cybernetics,Part C40(3):278-286,2010”） 和基于广义统计量直方图与聚类的鲁棒可逆水印方法进行实验对比， 实验中采用“5/3filter coefficients”小波变换将输入图像分成 一系列低频和高频子带。本发明的一些记号为：将Gao等人提出的基 于最小支撑树的特征区域选择方法记为MST，将基于广义统计量直方 图与聚类的鲁棒可逆水印方法记为WSQH-SC，将本发明方法记为 GPFR。

2.实验内容

实验1：鲁棒性实验

本发明进行鲁棒性实验的具体过程是：在水印嵌入过程中，首先 将测试图像进行三级提升小波变换，得到其第三级小波分解尺度下的 低频子带R_3,a，然后分别用MST和本发明两种方法构造圆形特征区域， 再利用本发明方法中的水印嵌入方法生成密钥信息；在水印提取过程 中，将待检测图像进行三级提升小波变换，得到其第三级小波分解尺 度下的低频子带R′_3,a，然后分别用MST和本发明两种方法构造圆形特 征区域，再利用本发明方法中的水印提取方法得到水印图像W′。实 验中特征尺度向量ξ＝[1,2]，MST方法的阈值取3。

本发明在JPEG压缩、JPEG2000压缩、加性高斯噪声、剪切、平 移五种攻击下测试两种方法的鲁棒性，实验中JPEG压缩的质量因子 为20，JPEG2000压缩的压缩率为0.2，加性高斯噪声的均值为0.02、 方差为0.05，剪切是去除图像中间400*400的区域，平移的大小是 向右向下平移20个像素。本发明采用相似度ρ作为评判依据，测试 两种方法的鲁棒性，ρ的计算公式为：

其中W(i,j)与W′(i,j)分别表示原始的水印图像W与提取的水印W′在 (i,j)处的元素值，ρ的取值越大，鲁棒性越强，反之亦然。

表1给出了五种攻击下MST与本发明方法的鲁棒性对比，其中相 似度的结果为所有测试图像的平均结果。由表1结果可见，本发明方 法的鲁棒性要高于MST方法。

表1.不同方法的鲁棒性

方法/攻击 JPEG JPEG2000 高斯噪声 剪切 平移 MST 0.9594 0.9576 0.9584 0.9540 0.9579 GPFR 0.9617 0.9614 0.9613 0.9587 0.9626

实验2：时间实验

本发明采用MATLAB中的tic与toc命令计算本发明方法与MST 方法在进行特征区域优化过程中的时间，表2给出了所有测试图像的 平均时间对比结果值。由表2结果可见，本发明方法的时间复杂度要 低于MST方法。

表2.不同方法的时间对比

方法 MST GPFR 时间(秒) 0.0644 0.0273

实验3：综合性能实验

将本发明方法与WSQH-SC方法进行实验对比，实验中分别用这两 种方法进行水印嵌入得到含水印图像，然后对生成的含水印图像分别 进行JPEG压缩、JPEG2000压缩、加性高斯噪声攻击得到降质的含水 印图像，最后再利用这些方法从降质的含水印图像中提取水印进行容 量、视觉质量和鲁棒性的对比。实验中，WSQH-SC方法的块大小为8×8， 水印嵌入强度为16。

本发明以原始图像中最多能嵌入的水印位数作为评判依据，测试 两种方法的容量。嵌入的水印位数越多，容量越大，反之亦然。同时， 以客观指标峰值信噪比PSNR作为评判依据，测试两种方法在原始图 像以最大容量嵌入情况下含水印图像的视觉质量，其中PSNR表示为

式中，M×N是原始图像大小，I(i,j)是原始图像在第i行第j列的像素 值，I′(i,j)是含水印图像在第i行第j列的像素值。PSNR越大，视觉质 量越好，反之亦然；当PSNR等于INF时，表示含水印图像与原始图 像完全一样，视觉质量最好。

表3给出了本发明方法与WSQH-SC方法在容量与视觉质量的综合 性能对比。在鲁棒性方面，WSQH-SC方法只能分别抵抗质量因子为70 的JPEG压缩、压缩率为1.2的JPEG2000压缩与均值为0、方差为0.01 加性高斯噪声，但本发明方法还能抵抗几何攻击，而且抗攻击的鲁棒 性更强。同时，本发明方法利用元胞自动机加密和解密算法对密钥图 像与特征尺度的加密与解密，提升了水印的安全性。

表3.不同方法的综合性能

方法/测试依据 容量 视觉质量 WSQH-SC 1006 36.6 GPFR 4096 INF

综上，本发明克服了现有鲁棒可逆水印方法的缺陷，改善了含水 印图像的视觉质量，提高了水印嵌入容量，增强了水印抵抗复杂攻击 的鲁棒性，提升了水印的安全性，提高了鲁棒可逆图像水印方法的综 合性能。

以上例举仅仅是对本发明的举例说明，并不构成对本发明的保护 范围的限制，凡是与本发明相同或相似的设计均属于本发明的保护范 围之内。