NSCT-GBT-SVD结合特征区域的鲁棒水印算法

张天骐 马焜然 杨宗方 邹 涵

（重庆邮电大学通信与信息工程学院，重庆 400065）

1 引言

在物联网和5G 时代，有关信息安全和版权归属的问题至关重要［1］。数字水印是保证数字产品信息安全和版权的一道重要防线。目前，数字水印的版权保护依然存在许多问题，其中主要的问题之一是数字产品在遭受攻击后如何保证版权的鲁棒性。近年来，鲁棒水印算法层出不穷，许多优秀的算法在遭受常规信号攻击后提取的水印依然具有很高的鲁棒性，但是在面对几何攻击（旋转，缩放，平移和剪切）时提取的水印表现的效果依然不佳。

目前，在解决几何攻击问题方面已有许多方法，其中有基于矩的不变域，Radon 变换校正，图像的归一化等操作。文献［2］通过调制分数阶多通道正交指数矩的幅值将二值水印嵌入，该类算法在面对旋转时表现出较好的抵抗能力，但是却受剪切攻击的限制。文献［3］采用Radon 变换对图像旋转进行校正，在面对旋转攻击时有很好的效果，但是无法很好的应对平移攻击。文献［4］利用非下采样轮廓波变换提取图像低频区域，再对此低频图像做归一化处理。全局的图像归一化仅在面对图像旋转时有很好的效果，在发生其他复杂的几何攻击时无法得到标准的归一化图像，因此抗几何攻击的综合能力较差。

越来越多的文献开始注意图像特征点的重要性，将图像特征点提取与水印算法相结合，可以使水印算法在抵抗几何攻击时表现出很好的效果。文献［5］采用轮廓波变换（Contourlet Transform，CT）分解原始宿主图像，将得到的低频图像采用Blob-Harris 算子提取稳定特征点并构造特征区域，同时结合图像局部归一化嵌入水印。该算法在平移，缩放，剪切体现出很好的性能，但是该算子不够稳定，在发生旋转时会导致特征点检测偏差，导致抗旋转的性能不佳。文献［6］采用Harris-Laplace 算子检测特征点并选择较为稳定的特征点，在特征点之间形成三角形，接着通过仿射变换微调每个三角形将水印信息嵌入。但是该算子提取的特征点稳定性不高，容易受缩放攻击的影响。文献［7］采用尺度不变特征变换（Scale Invariance Feature Transform，SIFT）提取原始图像的局部特征不变区域，然后对该区域应用条带波变换和离散余弦变换算法提取图像的视觉特征向量，最后结合零水印技术完成水印的嵌入和提取，该算法对常见几何攻击具有鲁棒性，但是对剪切表现出弱鲁棒性。文献［8］采用SIFT 算子提取原始图像的特征，选择非重叠并去除边缘附近的特征来确定有效嵌入的合适特征，然后根据特征的坐标和比例系数将水印信息奇偶量化嵌入到圆形区域。该方法在面对几何攻击中的旋转和剪切时表现出很好的效果，但是对于缩放攻击表现不佳。文献［9］采用仿射尺度不变特征变换（Affine-Scale Invariance Feature Transform，ASIFT）方法获得几何不变特征点实现几何变换后的同步，水印信息嵌入在离散小波变换-奇异值分解域。由于采用了特征点匹配和图像校正，在应对大部分几何攻击都有很好的效果，但是在面对剪切时表现不佳。文献［10］利用加速鲁棒局部特征（Speeded-up Robust Features，SURF）算子提取特征点用来几何校正，水印信息被分别嵌入到感兴趣区域和非感兴趣区域。该方法在应对特定角度旋转时有一定的抵抗能力，但是抗几何攻击的综合能力较弱。

为了解决这些问题，本文结合特征点与归一化局部区域嵌入，提出了一种NSCT-GBT［11］-SVD 结合（多尺度多方向结构张量［12］-模糊C 均值聚类［13］）特征区域的鲁棒水印算法。首先对宿主图像进行NSCT 获得低频图像，并采用多尺度多方向结构张量检测方式提取更加稳定的特征点，然后利用模糊C 均值聚类（Fuzzy C-Means algorithm，FCM）将特征点划分为三簇并寻找每个簇中相对稳定、分散且不重叠的特征区域。接着利用图像局部区域归一化操作，寻找水印的嵌入方形区域，将方形区域分块并用GBT 获取分块后的系数矩阵。最后对水印和系数矩阵分别进行奇异值分解，将水印的奇异值嵌入到系数矩阵的奇异值中。

2 NSCT、GBT和特征点检测

2.1 非下采样轮廓波变换

NSCT 克服了轮廓波变换的冗余度小、无移位不变性等缺点［14］，因此可以更好地表达图像的二维几何结构，有助于提取更多重要信息。NSCT 既有CT 的优点，又有近似移位不变的特点，能有效解决伪吉布斯效应和轮廓波变换的特征移位问题。由于使用了非下采样，NSCT 减少了由采样操作引起的滤波器失真，同时也确保了分解后各尺度方向的子带和原图大小一致。NSCT 由两部分组成，一部分是非下采样金字塔（Non-Subsampled Pyramid，NSP），该部分保证了分解过程中的多尺度特性。另一部分是非下采样方向滤波器组（Non-Subsampled Directional Filter Bank，NSDFB），该部分保证了分解中的多方向性。NSCT 分解结构与频率划分如图1所示。

图1 NSCT分解结构与频率划分Fig.1 NSCT decomposition structure and frequency division

2.2 基于图的变换

GBT 是一种以信号形式表示图形的新变换。GBT 有助于发现图像中重要的信息，并消除压缩冗余。图定义为G={V，E，s}，其中V和E是图的顶点和边，s表示图像的块。s中的像素与图的顶点相对应，像素之间的关系用邻接矩阵A来表示：

在公式（1）中，aij是G中顶点i和j之间边的权重。图拉普拉斯矩阵L：L=D-A，其中D为以下的对角矩阵：

对于实对称矩阵L，将其进行特征值分解生成一组非负特征值：

在公式（3）中，Λ={λ1，…，λN}是特征值集，V={v1，…，vN}是正交特征向量集。使用正交特征向量去相关s得到与s大小相同的GBT系数矩阵c：

GBT已经应用于音频/视频编码、音频水印等领域。GBT 的能量压缩性能与DCT 相当，但速度更快［11］。

2.3 多尺度多方向结构张量特征点检测

目前基于图像特征点选取的水印算法在数字水印领域应用范围不广，主要是因为提取特征点的算法存在不稳定的因素，特征点不稳定会影响特征区域的选择，对水印算法的抗攻击性能会有直接的影响。大多数特征点检测算法是基于2×2 结构张量的算法，而这些算法是基于Moravec（1979）的理论，即特征点周围的像素强度在各个方向上都变化很大。对于特征点，如果使用两个正交方向导数进行检测，在图像发生旋转时，特征点的水平和垂直方向导数可能会变小，从而导致检测不到旋转之后的特征点［12］。

Zhang 等人（2020）提出了一种多尺度多方向结构张量的特征点检测新方法，该方法是基于各向异性高斯方向导数滤波器来检测特征点。生成多方向各向异性高斯核和多尺度各向异性高斯方向导数，假定多尺度σs，各向异性因子ρ和方向角θk=(k-1)π/d，其中s为尺度等级，ρ取1.5，d为方向数，k=1，2，…，d。多方向各向异性高斯核函数的离散形式gσ，ρ，θ(x，y)和多尺度各向异性高斯方向导数ϕσ，ρ，θ(x，y)分别为：

该方法是在建立在多尺度的情况下进行的，首先建立一个多尺度空间σs(s=0，1，2)，再进行特征点检测：

步骤一：在当前尺度下使用多方向各向异性高斯滤波器平滑输入图像，并推导当前尺度下多方向导数。

步骤二：对于图像的每个像素点，在当前尺度下构造多方向结构张量矩阵M，M的定义如下。

其中，u，v均为6，k为某一方向，d为方向数(d=8)。

步骤三：使用d×d多方向结构张量，在每个尺度上的d个特征值{λ1，λ2，…，λd}形成一个新的特征点量。新的特征点量Es(x，y)定义如下：

其中，λk是多方向结构张量M的特征值，τ为常数(τ=2.22×10-16)，为了避免分母为0。

步骤四：执行非最大值抑制，如果像素对应的特征点量是7×7 窗口内的局部最大值，并且在当前尺度下大于阈值Th(Th=1.0×108)，则将该像素标记为候选特征点。将候选特征点在下一尺度下再次进行特征点检测，即执行步骤一到步骤四。

步骤五：如果候选特征点的对应特征点度量在所有尺度上都大于阈值Th，则将其标记为特征点。本文的σ0，σ1和σ2分别为1.5，3和4.5。

图2（a）是Lena在第一个尺度下的候选特征点，图2（b）是第二个尺度下的候选特征点，图2（c）是最终选择的特征点。

图2 三个尺度情况下提取的特征点Fig.2 Feature points extracted at three scales

3 构造水印嵌入区域

3.1 构造特征区域

稳定特征点选定后，要确保所选特征区域同样具有较好的稳定性。同时也要保证选取的特征区域不能重叠且足够分散，前者是为了避免特征区域重叠造成水印提取效果不佳，后者是为了避免特征区域因为几何攻击而丢失。因此使用模糊C 均值聚类算法来划分特征点进而构造特征区域，从而使构造的特征区域达到稳定、分散且不相互重叠。

特征区域构造具体流程如下：

步骤一：通过多尺度多方向结构张量特征点检测方法提取宿主图像的特征点，具体过程如2.3节。

步骤二：保存最终得到的特征点以及每个点对应的特征点量。

步骤三：用模糊C 均值聚类将提取到的特征点分成三簇。

步骤四：选取每个簇中特征点量最大的特征点，将三个最大特征点量进行从大到小的排序，排序结果决定了特征区域构造的顺序。首先以最大特征点量的特征点为圆心构造一个圆形特征区域。其次选择第二大特征点量的点进行圆形特征区域的构造。在构造的过程若与第一个簇中的圆形区域重叠，则选取该簇中第二大的特征点量，并以该点为圆心构造圆形特征区域。最后在第三个簇中进行圆形特征区域构造，同时也要满足不与第一和第二簇中的特征区域重叠。在整个构造过程中将靠近边缘区域的特征点删除以免构造的圆形区域超出边界。圆形区域的半径要与图像大小有关，以抵抗缩放攻击。假设图像的长和宽分别为w和h，则圆的半径为：

其中，β=22.76。图3（a），（b），（c），（d）和（e）分别是Lena，Baboon，Peppers，Portofino 和House 的特征点分类以及构造的特征区域。

图3 分类+特征区域Fig.3 Classification+Feature region

3.2 构造水印嵌入区域

为了进一步抵挡旋转攻击，本文采用局部特征区域图像归一化技术［4-5］。图像归一化实则矩的归一化，首先寻找图像可以抵抗仿射变换的矩，用不变矩去寻找一组参数，该参数可以消除其他变换对图像的影响。接着利用该参数将原始的图像转换为相应的标准图像，即使在图像发生仿射变换后依然可以变为标准图像。图4（a）为选取Lena 的特征区域再进行四周补零后的结果，该矩形区域的大小为2r×2r。图4（b）是对图4（a）区域进行归一化操作后的结果，归一化后的图像大小为512×512。图4（c）是以圆形特征区域为圆心向四周扩展出一个256×256的方形区域。

图4 水印嵌入区域Fig.4 Watermark embedding area

4 水印嵌入与提取

4.1 水印嵌入

水印嵌入过程如图5所示。

图5 嵌入过程Fig.5 Embedding process

步骤一：对水印预处理。将一幅64×64 大小的灰度图像W作为水印图像并对其进行Arnold 映射变为置乱后图像W′，对其进行奇异值分解：

其中，Sw作为嵌入信息，Uw和Vw作为密钥Key1存储以提取水印信息。

步骤二：提取特征点和构造特征区域。将一幅512×512 大小的灰度图像I作为宿主图像，对宿主图像I进行两层NSCT分解得到低频图像ILL，在低频图像上按照3.1节的步骤构造特征区域，得到三个稳定且不相互重叠的特征区域RLL={ri，i=1，2，3}。

步骤三：构造水印嵌入区域。取出一块特征区域ri，按照3.2 节的方法构造水印嵌入区域，先将特征区域补零，再将补零后的特征区域进行归一化操作，得到归一化后大小为512×512 的特征区域qLL，再选取圆形区域中大小为256×256 的特定方形区域sLL。

步骤四：对sLL进行4×4 分块，得到分块集合OLL={oi，j，1 ≤i≤64，1 ≤j≤64}。对 每个4×4 块oi，j进行GBT 变换，选取变换后的左上角系数，每个分块贡献一个GBT 系数ci，j，最终形成一个64×64的系数集合CGBT={ci，j，1 ≤i≤64，1 ≤j≤64}。对系数集合CGBT进行奇异值分解：

步骤五：将水印奇异值Sw按加性嵌入的方式嵌入到系数矩阵的奇异值Sc中，嵌入后的奇异值为S′，嵌入方式为：

其中，k为嵌入强度，本文取0.25。

步骤七：对其他两个水印嵌入区域执行步骤三到步骤六操作，同时将每个区域形成系数矩阵的奇异值Sc保存为S={Sc，c=1，2，3}，用于提取水印，这里设为密钥Key2。将重构后的低频图像与原始图像高频子带重构为带水印图像I′。

4.2 水印提取

水印提取过程如图6所示。

图6 提取过程Fig.6 Extracting process

步骤四：由密钥Key2提供的原始系数矩阵奇异值Sc进行水印系数提取，提取规则如下：

步骤五：对其他两个特征区域做同样步骤二到四操作，最终会得到三个对角矩阵。分别对三个对角矩阵上执行逆SVD：

其中，Uw和Vw由密钥Key1 提供。最终会得到三个置乱水印，进行Arnold 逆映射，选择效果最好的水印为提取出的水印。

综上所述，水印系统的整个流程分为两个环节：水印嵌入和水印提取环节。在水印嵌入环节中，需要进行两级NSCT 的分解与合并，以及多尺度多方向结构张量特征点检测，因此水印嵌入过程用时26.2 s，将未受攻击的图像提取出来需要进行两级的NSCT 分解和多尺度多方向结构张量特征点检测，因此水印提取需要17.4 s，整个算法的流程需要43.6 s，该执行时间可能随着计算机配置的变化而变化。

5 仿真结果与分析

本文实验在Windows10 操作系统上进行，采用的仿真工具为Matlab2018a。实验采用的宿主图像分别为512×512 标准灰度图像Lena、Baboon、Peppers、Portofino 和House，水印图像采用64×64 标准灰度图像Peugeot。水印算法通常通过两个重要参数进行评估：不可见性和鲁棒性。不可见性是指与原始图像相比，水印图像既不会扭曲也不会降低水印图像的感知质量。鲁棒性意味着仍然可以从被攻击后的图像中提取出水印。

5.1 不可见性分析

对不可见性的分析需要满足人类视觉主观上的判断和数据客观的评价。本文采用峰值信噪比（Peak Signal To Noise Ratio，PSNR）来表示嵌入水印后图像不可见性的客观评价。对于灰度图像，PSNR的定义式如下：

其中，m和n分别是宿主图像的长度和宽度。

图7 中（a），（b），（c），（d）和（e）分别为原始宿主图像Lena，Baboon，Peppers，Portofino 和House。（f），（g），（h），（i）和（j）分别为嵌入水印后的图像。在主观视觉层面上，嵌入水印后的图像与原始宿主图像几乎没有任何差别。在客观的数据层面，PSNR 大于40 dB 便有很好的不可见性，而本文五幅图像嵌入后的PSNR 分别为45.903 dB、46.171 dB、44.757 dB、46.187 dB 和46.261 dB。图8 中（a），（b），（c），（d），（e）和（f）分别是原始水印，Lena，Baboon，Peppers，Portofino 和House 未受攻击所提取出的水印，提取出的水印归一化相关系数（Normalized Correlation Coefficient，NC）值分别为0.999，0.998，0.999，0.998和0.997。由于本文在构造特征区域的过程中，采用了图像的局部归一化的操作，而该操作会将原始局部区域放大后嵌入水印再进行缩小，因此会丢失一部分的水印信息，导致提取水印的NC 值无法达到1，也是造成图8 中视觉效果有一定差异的原因。但五幅嵌入水印后的宿主图像PSNR 均大于44 dB，且提取的水印NC 均大于0.997，说明了本文算法满足不可见性的要求。

图8 原始水印与未受攻击提取的水印Fig.8 The original watermark and the watermark extracted without attack

5.2 鲁棒性分析

鲁棒性表示了算法的抗攻击性能，本文采用归一化相关系数作为评价算法的鲁棒性的标准。NC值越接近1，提取出来的水印效果越好，算法的鲁棒性越好。NC的定义式为：

其中，w和h分别是水印图像的长度和宽度，W*和W分别是提取的水印和原始水印。为了验证本文算法的鲁棒性，分别从常规攻击和几何攻击进行实验。

本文用于评估水印鲁棒性的常规攻击包括JPEG压缩因子为70、50和30；均值为0，方差分别为0.01 和0.02 的椒盐噪声、斑纹噪声和高斯噪声；窗口分别为2×2 和3×3 的中值滤波、均值滤波和高斯滤波以及直方图均衡。水印的抗常规信号攻击性能的实验结果见表1。

表1 水印对常规信号攻击的鲁棒性Tab.1 Robustness of watermark against conventional signal attacks

由表1 可知，本文算法在面对常规信号攻击（JPEG压缩，噪声攻击，滤波攻击）都有很好的效果，最大NC 值均在0.96 以上。特别是在面对JPEG 压缩时，提出的水印NC均在0.99以上，因为水印最终是嵌入在GBT 变换后的抗压缩位置上，拥有很好的抗压缩能力。但是在面对直方图均衡时表现出来的效果较差，因为该攻击大幅度的改变宿主图像的像素值，导致提出的水印效果较差。

常见的几何攻击有旋转，剪切，缩放和平移，这些攻击会使宿主图像像素位置和内容发生改变，导致水印提取位置发生偏差而无法正确的提取出水印信息，是许多水印算法所忽略的。为了进一步验证本文算法的鲁棒性，本文对已选取的五幅图像进行旋转，剪切，缩放和平移，以及与常规信号组合攻击来测试其鲁棒性。水印算法的抗几何攻击和抗组合攻击的性能结果如表2所示。图9中（a）～（l）和图10（a）～（l）是选取部分攻击和提取出的水印效果。

表2 水印对几何攻击和组合攻击的鲁棒性Tab.2 The robustness of watermarking against geometric attacks and combinatorial attacks

图9 遭受几何攻击的图像Fig.9 Image under geometric attack

图10 不同几何攻击下提出的水印Fig.10 Extracted watermarks under different geometric attacks

从表2 和图10 可知，本文算法在面对旋转时，即使发生大幅度旋转，提出的水印NC 也在0.98 以上。在面对剪切时，只有Lena这幅宿主图像在左边剪切1/2 时，提取出来的效果略差于其他情况，因为Lena 这幅宿主图像的特征区域位置较左边，在发生左边大面积剪切时，对Lena的三个特征区域有一定的影响。在面对平移时，水印提取的效果最稳定，提取的水印NC值均在0.98以上。对于三者的两两组合攻击，提取的水印NC也在0.95以上。

在面对缩放时，提取的水印无法达到0.99的效果。原因在于特征半径的选取与图像的边界有关，而利用本文设置的半径因子和图像大小去求取半径时，所求出的半径会有一定的误差，这个误差将会影响特征区域的大小，从而导致提取出的水印NC 值略小于其他强度的攻击所提取的效果，但在面对单个缩放攻击时提取的水印NC 也有0.96 以上。对于旋转、平移、剪切和缩放任意两者的组合攻击，提取的水印NC均在0.94以上。

在面临常规信号和几何信号的组合攻击时，提取的水印效果略差于单个类型的攻击。从图9（c）可以看出，该组合攻击对图片改变很大，提取出的水印也较差于其他情况。这是因为缩放会影响特征区域的半径，方差为0.01的高斯噪声会对宿主图像有一定影响，导致提取出的水印NC 值在0.96 左右。但对于常规信号攻击与几何攻击的组合，提取的水印NC值均在0.95以上。以上分析体现了本文算法对于几何攻击以及组合攻击的鲁棒性。

5.3 对比实验

为了进一步验证本文算法的鲁棒性，将文献［5］，文献［9］以及文献［15］与本文算法进行对比实验，宿主图像选用标准灰度图像Lena，大小为512×512，水印选择64×64标准灰度图像。对嵌入后的图像进行几何攻击和非几何攻击，所提取出的水印效果如表3。

表3 鲁棒性实验对比Tab.3 Comparison of robustness experiments

文献［5］与本文是特征区域嵌入型，文献［9］与文献［15］是特征点匹配校正型。根据表3的结果，本文算法在面对常规信号攻击时，抗噪声性能明显优于文献［5］和文献［15］，与文献［9］略微持平，但本文算法在抵抗滤波攻击时略差于文献［9］。在面对几何攻击时，仅在缩放的情况下略低于其他算法的性能，但在面对旋转，剪切和平移时均强于文献［5］、文献［9］和文献［15］。本文算法在抵抗几何攻击综合能力方面上优于其他三种算法。综上分析可知，本文在抵抗常规信号攻击和几何攻击时，更具有鲁棒性。

6 结论

本文采用具有更加精准的特征点检测方式进行多尺度多方向特征点检测，保证了图像在发生几何形变后能够再次精准的定位到选取的特征点。为了能够应对剪切和平移攻击，本文采用FCM 的方式划分特征点，在每个特征点簇中选取分散、不重叠的稳定特征区域。为了更好的应对旋转和缩放攻击，本文对局部区域进行归一化，并将嵌入区域的半径与图像大小联系。水印信息的嵌入是在NSCT 低频域上进行的，保证了水印的稳定性，同时加入了GBT 变换，保证了水印的抗压缩能力。实验表明，本文算法在水印嵌入后PSNR 可以达到44 dB以上。在面对噪声和JPEG 常规信号攻击时，提取的水印NC值均在0.98以上。在面对旋转，平移，剪切，缩放和组合的几何攻击时，提取的水印NC 均在0.94以上，特别是旋转，平移这两种攻击时，提出的水印NC可达0.98以上。体现了本文算法的抗几何攻击的鲁棒性。