觉得枯燥的话点击网页左下角听几首精选纯音乐吧~

前言

嘛,这是一篇大三上学期《信号处理综合设计》课程的设计论文,设计的内容和论文的结构都是老师指定的。为了水字数写得有些繁琐和生涩,有需要的路人盆友们挑重点看吧~

为什么放在博客里?因为这篇论文写的比较满意丫,居然从很严格的老师那里拿到了还不错的分数(而且想试试博客里LaTeX公式渲染的效果(靠,用KaTeX会蜜汁移位QAQ虽然它加载快(还是用MathJax吧...

强调一下,由于是赶出来的课程论文,不保证所言的严谨性和正确性,仅作学习交流和参考!

(对了,虽然参考了很多CSDN博客的代码,但参考文献不能放网站,只能贴正经论文,现在也不记得当时参考了哪些博客的代码了QAQ)


一、设计要求

1.背景综述

  数字盲水印技术,是一种嵌入式的防伪技术,它是将一些经过加密的防伪信息嵌入到视频、音频和数字图像等数字产品当中,用于数字产品的溯源;也用于区别盗版与正版,保护数字产品版权;也可以作为一种加密技术使用。
  而其中的数字图像不可见水印(盲水印),它不同于普通的数字图像水印,后者使用传统的方式“暴力”地向原始图像中叠加较为明显的水印信息,既影响了原图的信息、影响美观度,又容易二次修改或暴力去除,不能起到很好的防伪作用。而本次设计研究的内容就是使用“盲水印算法”将水印图像“优雅”地嵌入到目标图像中,尽可能地不对原始图像信息造成明显破坏,又能实现隐藏防伪水印。
  盲水印的嵌入基本流程可用下图1表示:
图1-盲水印嵌入流程
  而需要时可以通过复原算法重新将水印内容提取出来,其基本流程如下图2所示:
图2-盲水印提取流程

2.设计内容

  由于盲水印的嵌入和提取算法繁多,因此本次设计主要先是查阅资料,研究不同盲水印嵌入、提取算法的实现方式,对比其特点和优劣,选择其中一种进行详细的算法理论分析和仿真测试。而后通过编写代码,对比计算嵌入前后的峰值信噪比(PSNR)、提取水印前后的误位率(BER)和归一化互相关(Nc),分析其性能,并研究可能影响性能指标的参数。
  最后对嵌入水印的图像实施不同的“攻击”,如图像JPEG压缩、涂鸦或(保持原大小的)剪切、滤波、噪声叠加、锐化等,这些处理可视作对盲水印攻击,降低盲水印的质量。结合前面实现的盲水印嵌入和提取测试、指标,对这些攻击造成的影响进行客观的分析和评价。

二、方案的论述与比较

  总体上来讲,数字图像盲水印的算法可以分为三大类,分别是压缩域、空间域、变换域,总体上他们各自有着不同的特点。

  1. 压缩域:在数字图像进行信息压缩编码的时候嵌入水印信息,要对信息存储和压缩算法进行复杂的变换,总体上设计难度较大,但无论是嵌入前后的不可感知性、提取的前后的鲁棒性,都是比较好的。
  2. 空间域:在图像的像素空间域中嵌入水印信息,其特点是直观、实现较简单,但普遍来讲空域算法不可感知性较差,易察觉与原始图像不同。容易受到空域像素点“攻击”的干扰,因此鲁棒性也较差。
  3. 变换域:在对图像进行变换后,在变换系数上嵌入水印信息,其特点是不可感知性好、鲁棒性强,在合适的嵌入强度下人类肉眼难以察觉。常用的变换域算法有基于离散余弦变换域(DCT)、基于小波变换域(DWT)和基于离散傅里叶变换域(DFT)的算法。
      在空域图像水印算法中,一个典型的、常被人提起的是空域最低有效位 (Last Significant Bit, 简称LSB) 算法。其基本方法是用欲嵌入的秘密信息取代载体图像的最低比特位,原来的图像的高位平面与代表秘密信息的最低平面组成含隐蔽信息的新图像。

  位平面的概念图如下图3,若有一个灰度值为78的像素点,其数据对应二进制为01001110,最低为即最右边的0。
图3-位平面示意图
  可以将一幅图像不同位平面的数据提取出来做对比,会发现位平面越低的数据所贡献的信息量越少,人眼也更不容易察觉到,如下图4所示。
图4-不同位平面信息量
  因此直接对每个像素点的最低位进行信息的嵌入,即可完成水印的嵌入;提取时也是直接将最低位的信息抽取出来,重新组合即可得到水印信息。
  故不难总结出,LSB算法的优点是算法简单、嵌入容量大、不可见性还算好,尤其是提取信息时不需要原始图像,即支持盲检测。当然也可以在最低位替换前将水印信息通过置乱(密钥)等方式加密,提取时再通过密钥还原,增加隐蔽性。
  但是, 该算法实质相当于在图像中添加了高频噪声,因此对图像进行几何变形、信号处理“攻击”如滤波、压缩、加噪声等抵抗能力差, 鲁棒性不好。
  在变换域中,一个常提到的算法是离散余弦变换域(DCT)算法。离散余弦变换是傅立叶变换的一种特殊情况,在傅立叶级数展开式中,如果被展开的函数是实偶函数,那么其傅立叶级数中只包含余弦项,再将其离散化可导出离散余弦变换,因此DCT变换避免了傅里叶变换中的复数运算,它是基于实数的正交变换。
  而DCT变换域数字水印算法的基本原理是将空域图像变换为DCT频域,然后将水印信息嵌入其直流项之中,最后再将频域转换成空域以完成图片的水印的嵌入。
  总的来说,基于DCT变换水印嵌入算法,其特点是算法复杂度中等,同时能有较好的不可见性和鲁棒性。因此后面的详细详细分析与MATLAB仿真验证都会围绕基于DCT变换的水印嵌入算法展开。

三、采用方案的算法分析

1.离散余弦变换(DCT)算法

  鉴于上述提到的,DCT变换域具有优秀特性,因此本组的理论分析研究、MATLAB仿真都将以“基于DCT变换水印嵌入算法”展开。
  对于离散余弦变换(DCT)变换,实质上是在傅立叶变换过程中,如果被展开的函数是实偶函数,那么其傅立叶变换中只包含余弦项,基于这一特点提出的。DCT变换先将函数变成偶函数形式,再对其进行二维离散傅立叶变换,因此DCT变换可以看成是一种简化的傅立叶变换。
  DCT分为一维变换和二维变换,而在图像处理中用到的是二维,其定义式如下:
  若图像为$ f\left(i,j\right) $,则其2DCT变换$ F\left(u,v\right) $为:
$$\mathrm{F}(u, v)=\alpha_{0} \mathrm{c}(u, v) \sum_{x=0}^{N-1} \sum_{y=0}^{N-1} f(x, y) \cos \frac{(2 x+1) u \pi}{2 N} \cos \frac{(2 y+1) v \pi}{2 N}$$
  其中,$u, v=0,1, \ldots, N-1, \alpha_{0}=\frac{2}{N}, c(u, v)=\left\{\begin{array}{cc}1 / 2 & u=v=0 \\ 1 / \sqrt{2} & u v=0 \text { 且 } u \neq v \\ 1 & u v>0\end{array}\right.$

  同样的,2DCT变换有着逆变换2IDCT,其定义式如下:
$$ f\left(x,y\right)=\alpha_1\sum_{u=0}^{N-1}\sum_{v=0}^{N-1}{c\left(u,v\right)F\left(u,v\right)\cos{\frac{\left(2x+1\right)u\pi}{2N}}\cos{\frac{\left(2y+1\right)v\pi}{2N}}} $$
  其中,$u, v=0,1, \ldots, N-1, \alpha_{1}=\frac{2}{N}, c(u, v)=\left\{\begin{array}{cc}1 / 2 & u=v=0 \\ 1 / \sqrt{2} & u v=0 \text { 且 } u \neq v \\ 1 & u v>0\end{array}\right.$

  上述算法可以自己根据算法原理、公式进行编写,也直接可以使用MATLAB工具箱中的函数dct2进行二维DCT变换,调用方法如下:

B = dct2(A)

  上述代码,返回A的二维离散余弦变换给B。其中矩阵B包含离散余弦变换系数B(k1,k2)。
  同样地,也可以使用工具箱中的函数idct2进行二维DCT的反变换,调用的方法如下:

A = idct2(B)

  上述代码,返回B的二维离散余弦变换给A。

2.使用DCT嵌入水印信息

  可以通过上述的代码,很容易地对图像整体进行DCT正变换和逆变换,然而想要嵌入水印,通常是将原始图像划分成无重复的像素块,再对各分块进行DCT变换,对变换后的分块进行水印信息嵌入。当然也有一些学者是计算整个图像的DCT再将水印信息嵌入到中频系数上的,都是利用DCT变换,但实现的途径和得到的效果有所不同。
  我们选择的嵌入水印方案可用下面的流程图5描述:
图5-水印嵌入流程
  假设原始图像是一副512×512的8位灰度图(0~255),水印图像是一个黑白的二进制图像。我们需要先对原始图像进行分块,不妨取分块的大小为8×8,那么最多可以划分64×64个不重复的分块,设块坐标(m,n)。图像分块示意图如下图6所示。
  因此对于需要嵌入的水印图像,其大小不能超过64×64,否则需要裁剪才能嵌入,这会造成信息丢失。当然,可以小于64×64,只不过需要填充像素。
图6-图像分块示意图
  下面将原图像的每个8×8块分别作DCT变换,在此之前我们需要知道每个块起始点点在原始图像中的坐标(x,y)。通过推算不难得出:x=(m-1)×8+1、y=(n-1)×8+1,下图7可以直观地表示。
图7-坐标推算示意图
  由此又可以推出该块中所有像素点的坐标,可以使用数组或范围的形式在MATLAB中使用代码表示:

BLOCK=I(x:x+K-1,y:y+K-1)

  其中,BLOCK为存放分块数据的变量,I为原始图像,K为块大小(本例为8)。
  完成分块的DCT变换后,接着要对各分块嵌入水印信息,由于嵌入的水印是二值黑白图像,故可以通过判断水印像素点的黑或白,选择嵌入不同的强度的信息。可以使用如下的MATLAB代码:

BLOCK(1,1)=BLOCK(1,1)*(1+a*Ins);

  其中Ins是嵌入强度,一般取0.02左右;a是根据对应水印黑白值判断的系数,比如为黑取-1、为白取1,从而实现水印的黑与白嵌入不同强度的信息。
图8-水印信息嵌入示意图
  接着,对嵌入完水印信息的各块再进行DCT反变换,这样就可以还原到图像的空域状态,即得到嵌入完水印的图像。但要注意此时的各像素灰度值会变为浮点型,而图像存储和显示格式不支持浮点型,因此要进行取整再保存。

3.提取还原水印信息

  提取水印可以看作嵌入水印的逆过程,总体上讲它也需要进行分块,而且还需要将原始未嵌入水印图像的分块、嵌入水印后图像的分块作同步比对。将已嵌水印图像分块DCT变换后的直流分量BLOCK2除以对应原始图像分块DCT变换后的直流分量BLOCK1,由于嵌入时的算法会使得这个商在1附近,黑色像素对应商小于1、白色像素对应商大于1。
  为了判断方便,将得到的商-1后直接判断其是否大于0,若小于0则还原为黑像素、否则还原白像素,通过如此往复地对各分块进行运算和比对,就可以得到原始的二值水印图像。提取水印流程图总结如下:
图9-水印提取流程

4.利用Arnold变换置乱水印

  由于上述方法嵌入水印,得到的嵌入水印后图像在一定程度上可以用肉眼、增强对比度的方式查看到水印信息,这不利于水印信息的保密和隐藏。我们可以使用一种可逆的、带有参数的二维图像置乱方法,对水印进行置乱,再进行嵌入处理;验证水印的话,进行水印提取后时用同样参数逆变换去复原水印即可。
  常用的方法是俄国数学家VladimirI Arnold提出的一种变换,被后人称之为Arnold变换,其变换定义式如下:
  设一副M×N的数字图像,则其二维Arnold变换为:
$$\left[\begin{array}{c}x_{n+1} \\ y_{n+1}\end{array}\right]=\left[\begin{array}{cc}1 & b \\ a & a b+1\end{array}\right]\left[\begin{array}{c}x_{n} \\ y_{n}\end{array}\right] \bmod (N)$$
  其中,x、y表示变换前灰度图中像素的位置,$ x_{n+1} $、$ y_{n+1} $表示变换之后的像素位置,a、b为参数,n表示当前变换的次数,N为图像的长或宽(由于该算法只适用于长宽相等的图像,故暂不讨论M不等于N的情况),mod为取余运算。
  同样的,当想要复原置乱的图像时,使用变换矩阵的逆矩阵即可,表达式如下:
$$\left[\begin{array}{c}x_{n+1} \\ y_{n+1}\end{array}\right]=\left[\begin{array}{cc}a b+1 & -b \\ -a & 1\end{array}\right]\left[\begin{array}{c}x_{n} \\ y_{n}\end{array}\right] \bmod (N)$$
  但是要注意使用的参数需要与变换时一致,因此可以理解为使用Arnold变换就是执行了一次简单的水印加密,加密用的参数a、b加上变换(置乱)次数n,也就只有三种可变因素,因此抵抗暴力猜测的能力是有限的。

5.嵌入水印的扩展探讨

  上述分析中只使用了灰度图像作载体、水印图像为黑白二值图像,如果想使用彩色图片作载体,有两种实现方式:①在彩色R、G、B三个通道中任取一个或多个通道进行水印嵌入与提取。②把图像从RGB空间转到YUV空间,提取亮度分量Y,并在分块DCT变换基础上求中频系数绝对值的算术平均值,最终实现水印信息的嵌入与提取。
  对于想要嵌入灰度水印图,一般需要对灰度水印进行预处理,将灰度信息进行映射变换之后再进行嵌入,算法较为复杂,因而不展开讨论。彩色水印的话也可以使用上述的方法进行结合,实现嵌入和提取,但算法难度大、实现后的效果并不如意。

6.嵌入水印后图像评估

  在对嵌入水印引起原图像的失真进行客观评估时,使用峰值信噪比(PSNR)来度量。PSNR的定义如下:
$$ PSNR=10\lg{\frac{\left(2^k-1\right)^2}{MSE}} $$
$$ MSE=\frac{1}{m\times n}\sum_{i=0}^{m-1}\sum_{j=0}^{n-1}\left(X_{ij}-X_{wij}\right)^2 $$
  式中:k为像素灰度的字长,X为原始图像,$ X_{w} $为含水印图像,m、n分别为原始图像X的行数与列数。
  峰值信噪比越大,人们就越难以察觉原始图像嵌入信息的差别,原始图像在嵌入图像后差异就不会那么大,即嵌入水印对原始图像引起的失真就不会很大。在一般情况下,峰值信噪比的值在30dB以上,我们利用人眼就很难将原始图像X与嵌入水印后的图像 $ X_{w} $分辨出来。

7.水印提取质量评估

  提印出来的的水印可能与原始水印不是完全一致,可以通过计算原始水印图像W与提取出来的水印图像W'的误位率BER和归一化互相关Nc来度量。
  误位率BER的定义为:
$$ BER=\frac{Ne}{S\times T} $$
  式中Ne为原始水印图像与提取出来的水印图像的不同像素的个数,S、T分别为水印图像的行数与列数。
  误位率BER的值越接近于0,表明提取出来的水印图像的失真度就非常小,几乎不失真;当BER=0时,W'与W完全一致。
  另一个评估标准是归一化互相关Nc,其定义如下:
$$ Nc=\frac{\sum_{i=1}^{S}\sum_{j=1}^{T}W_{ij}W_{ij}^\prime}{\sum_{i=1}^{S}\sum_{j=1}^{T}\left[W_{ij}\right]^2} $$
  式中$ W_{ij} $为原始水印图像,$ W_{ij}^\prime $为提取出来的水印图像,S、T分别为水印图像的行数与列数。
  或者可以转化为如下的表达式:
$$ NC=\frac{\sum_{m=1}^{M}\sum_{n=1}^{N}{T\left(m,n\right)S^{ij}\left(m,n\right)}}{\sqrt{\sum_{m=1}^{M}\sum_{n=1}^{N}{T^2\left(m,n\right)\sum_{m=1}^{M}\sum_{n=1}^{N}\left(S^{ij}\left(m,n\right)\right)^2}}} $$
  若归一化互相关Nc的值越接近1,说明提取出来的水印图像与原始水印图像的差异就越小;当Nc=1时,W'与W完全一致。

四、算法仿真与结果分析

1.图像整体DCT变换

  我们首先尝试了调用MATLAB工具箱的dct2函数,对目标图像进行整体的DCT变换,代码如下所示:

%--------------图片整体DCT变换----------------
clear;
clc;
img = imread('Peppers.tiff'); % 读取图片
img2gray = rgb2gray(img); % 转为灰度图
img_dct = dct2(img2gray); % 得到2维DCT变换后的图片
% 绘图部分
figure('Name','图片整体DCT变换前后');
subplot(2,2,1);
imshow(img); title('原彩色图');
subplot(2,2,2);
imshow(img2gray); title('原灰度图');
subplot(2,2,3);
imshow(img_dct); title('DCT后图像');
subplot(2,2,4);
imshow(log(abs(img_dct)),[]); title('DCT后图像(取幅度和log)'); % 便于查看高频和低频信息
colormap(gca,jet(64)) %转为色度映射便于直观查看
colorbar

  可以得到如下图所示的结果,不难看出子图4中大部分能量集中在左上角,其对应的是图像的低频部分。
图10-图像整体DCT变换

2.基于DCT的水印无置乱嵌入与提取

  首先尝试了不置乱水印的DCT嵌入,按照之前的原理编写代码如下:

%----------------------------DCT方法水印嵌入与提取--------------------------------
clear;
clc;
%----------------变量定义------------------
M = 512; %原图像长宽
K = 8; %8x8的mask
% 因此水印图片不大于512/8=64(64x64)
N = 64; %水印图像长宽
Ins = 0.02; %嵌入强度
%----------------建立矩阵------------------
I = zeros(M,M); %创建全0的MxM矩阵用于存放图像数据
I_BWM = zeros(M,M);
J = zeros(N,N);
BLOCK = zeros(K,K);
%------------图像读取及预处理---------------
figure('Name','DCT方法水印嵌入与提取');
img = imread('Peppers.tiff'); %读取图片到变量img中
I = rgb2gray(img); %转为灰度图并赋给I
subplot(1,4,1); %显示多幅图,其中的第一个
imshow(I);title('原始灰度图像'); %显示原始图像
 
J = double(rgb2gray(imread('WYU64.bmp'))); %读取水印图片到变量J中
subplot(1,4,2); %显示多幅图,其中的第二个
imshow(J);title('水印图像'); %显示水印图像
%---------------水印嵌入算法----------------
for p=1:N %从1到64循环(mask行扫描)
    for q=1:N %从1到64循环(mask列扫描)
        x=(p-1)*K+1; %mask此次扫描的起始像素点x坐标
        y=(q-1)*K+1; %mask此次扫描的起始像素点y坐标
        BLOCK=I(x:x+K-1,y:y+K-1); %取mask大小(8X8)的像素保存到BLOCK中
        BLOCK=dct2(BLOCK); %对BLOCK(及经过mask取出8x8区域)进行DCT变换
        if J(p,q)==0 %如果水印上该像素值为0(黑)
        a=-1;
        else
        a=1;
        end
        BLOCK(1,1)=BLOCK(1,1)*(1+a*Ins); %目的是取不同嵌入强度嵌入到直流块中
        BLOCK=idct2(BLOCK); %DCT反变换
        I_BWM(x:x+K-1,y:y+K-1)=BLOCK; %反变换后放回图像
    end
end
%----------------嵌入水印后处理----------------
I_BWM_int = uint8(I_BWM); %像素点数据double转int,才能正常显示和保存
subplot(1,4,3); %显示多幅图,其中的第三个
imshow(I_BWM_int);title('嵌入水印后的图像'); %显示嵌入水印后的图像
imwrite(I_BWM_int,'BWMed.bmp','bmp'); %保存嵌入水印后的图像为BWMed.bmp
%-------------提取水印算法(完整版)--------------
for p=1:N
    for q=1:N
        x=(p-1)*K+1;
        y=(q-1)*K+1; %与之前操作相同,移动mask
        BLOCK1=I(x:x+K-1,y:y+K-1); %对原图取此次mask给BLOCK1
        BLOCK2=I_BWM_int(x:x+K-1,y:y+K-1); %对含水印图取此次mask给BLOCK2
        BLOCK1=dct2(BLOCK1); %对BLOCK1进行DCT变换
        BLOCK2=dct2(BLOCK2); %对BLOCK2进行DCT变换
        d=BLOCK2(1,1)/BLOCK1(1,1)-1;
        if d<0
            F(p,q)=0; %小于0则复原该水印像素为黑
        else
            F(p,q)=255; %否则复原该水印像素为白
        end
    end
end
%显示提取的水印
subplot(1,4,4);
imshow(F,[]);
title('提取的水印');

  运行程序可以得到如下图11,这是嵌入强度为0.02时、嵌入块为直流块得到的结果。
图11-强度0.02、直流、无置乱水印嵌入提取结果
  可以由上图11看到,肉眼难以察觉嵌入水印前后图像的区别,得到的提取水印虽然有像素点信息丢失,但结果整体挺好的,这是由于反变换后取整会导致细节信息丢失,从而提取还原水印时,黑白二值的判断会出错。
  除了主观的判断之外,还需要客观的指标定量地去衡量,可以用理论分析中的指标“峰值信噪比(PSNR)”去评估嵌入水印对载体图片的影响(即不可觉察性)、利用指标“误位率(BER)”和“归一化互相关(Nc)”去评估提取水印的质量,还能在水印攻击后用于对比分析该方法的鲁棒性。上述代码如下:

%---------------计算峰值信噪比(PSNR)------------------
%ImageA为原始图像,ImageB为含水印图像
function dPSNR = psnr(ImageA,ImageB)
%保护性判断,图像大小不相等则抛出异常
if (size(ImageA,1) ~= size(ImageB,1)) or (size(ImageA,2) ~= size(ImageB,2))
    warning('图像大小不相等!');
    dPSNR = 0;
    return;
end
ImageA=double(ImageA); %转为double型便于浮点计算
ImageB=double(ImageB);
M = size(ImageA,1); %获取长
N = size(ImageA,2); %获取宽
d = 0 ; %初始值为0
for i = 1:M
    for j = 1:N
        d = d + (ImageA(i,j) - ImageB(i,j)).^2 ; %计算MSE的一部分
    end
end
dPSNR = 10*log10((M*N*max(max(ImageA.^2)))/d) ; %计算PSNR
return
%---------------计算误位率(BER)------------------
%ImageA为原始图像,ImageB为含水印图像
function dBER = ber(ImageA,ImageB)
if (size(ImageA,1) ~= size(ImageB,1)) or (size(ImageA,2) ~= size(ImageB,2))
    warning('图像大小不相等!');
    dBER = 0;
    return;
end
ImageA=double(ImageA); %转为double型便于浮点计算
ImageB=double(ImageB);
M = size(ImageA,1); %获取长
N = size(ImageA,2); %获取宽
error=0;%初始值为0
for i = 1:M
    for j = 1:N
        if(ImageA(i,j) ~= ImageB(i,j)) 
            error = error + 1;
        end
    end
end
dBER=error/(M*N); %根据公式计算误位率BER 
return
%---------------计算归一化相关系数(NC)------------------
%ImageA为原始图像,ImageB为含水印图像
function dNC = nc(ImageA,ImageB)
if (size(ImageA,1) ~= size(ImageB,1)) or (size(ImageA,2) ~= size(ImageB,2))
    warning('图像大小不相等!');
    dNC = 0;
    return;
end
ImageA=double(ImageA); %转为double型便于浮点计算
ImageB=double(ImageB);
M = size(ImageA,1); %获取长
N = size(ImageA,2); %获取宽
d1=0;%初始值为0
d2=0;
d3=0;
for i = 1:M
    for j = 1:N
        d1=d1+ImageA(i,j)*ImageB(i,j) ;
        d2=d2+ImageA(i,j)*ImageA(i,j) ;
        d3=d3+ImageB(i,j)*ImageB(i,j) ;
    end
end
dNC=d1/(sqrt(d2)*sqrt(d3)); %根据公式计算NC
return

  将上述代码添加至单个.m文件中封装成函数,文件名与函数名一致放于工程目录下,在主程序中调用函数,输入需要对比的图片二维数组即可。调用代码如下:

%评价图像失真情况(PSNR)
PSNR = psnr(I,I_BWM_int) %计算PSNR并打印(无分号)
%评价水印提取质量(BER)
BER = ber(J,F) %计算BER并打印(无分号)
%评价水印提取质量(NC)
NC = nc(J,F) %计算NC并打印(无分号)

  可以得到图11结果对应的水印嵌入评估如下所示:
图12-强度0.02、直流、无置乱水印嵌入提取评估
  可以看到PSNR≈38.73dB,在30dB以上,人眼不易察觉嵌入水印前后载体图像的区别。BER≈0.93%、Nc≈0.994,说明提取的水印与原始水印非常接近,出错的像素数较少。

3.水印置乱后的嵌入与提取

  但由于是直接嵌入到直流块中,如果仔细观察图11的嵌入水印后的载体图片,可以隐约观察到水印的像素,再通过增强对比度等方式可以轻易地将水印暴露出来,这不利于加密信息、版权保护。因此利用理论分析中提到的Arnold变换,我们可以进行待嵌水印的置乱、加密。首先根据Arnold变换算法原理将变换过程封装为子函数,代码见下:

%---------------图像Arnold变换置乱算法------------------
%img 灰度图像 a,b为参数 n为变换次数
function arnoldImg = arnold(img,a,b,n)
[h,w] = size(img);
N=h;
arnoldImg = zeros(h,w);
for i=1:n
    for y=1:h
        for x=1:w
            %防止取余过程中出现错误,先把坐标系变换成从0 到 N-1
            xx=mod((x-1)+b*(y-1),N)+1;
            yy=mod(a*(x-1)+(a*b+1)*(y-1),N)+1;  
            arnoldImg(yy,xx)=img(y,x);              
        end
    end
    img=arnoldImg;
end
arnoldImg = uint8(arnoldImg);
end

  同样地,可以由Arnold逆变换编写复原算法如下:

%---------------图像Arnold变换复原算法------------------
%img 灰度图像 a,b为参数 n为变换次数
function img = rearnold(arnoldImg,a,b,n)
[h,w] = size(arnoldImg);
img = zeros(h,w);
N = h;
for i=1:n
    for y=1:h
        for x=1:w
            xx=mod((a*b+1)*(x-1)-b*(y-1),N)+1;
            yy=mod(-a*(x-1)+(y-1),N)+1;
            img(yy,xx)=arnoldImg(y,x);
        end
    end
    arnoldImg=img;
end
img = uint8(img);
end

  将置乱水印的Arnold变换代码插入到主函数中,在水印嵌入前添加如下代码:

J = double(rgb2gray(imread('WYU64.bmp'))); %读取水印图片到变量J中
subplot(2,3,2); %显示多幅图,其中的第二个
imshow(J);title('原始水印图像'); %显示原始水印图像
J_rd = arnold(J,2,3,5); %置乱水印,a=2、b=3、变换次数n=5,赋给J_rd
subplot(2,3,3); %显示多幅图,其中的第三个
imshow(J_rd);title('置乱的水印图像'); %显示置乱的水印图像

  在提取出水印后添加如下代码,恢复置乱:

%显示置乱复原的水印
subplot(2,3,6);
F_rerd = rearnold(F,2,3,5);
imshow(F_rerd,[]);title('复原提取的置乱水印');

  上述代码有些非关键部分省略,运行程序可以得到如下图的结果:
图13-强度0.02、直流、加置乱水印嵌入提取结果
  从主观上讲,由于加入置乱,嵌入水印的载体图片只能在仔细对比的情况下看到一些随机色块,完全无法得知水印的内容,整体观感比之前更好。而且提取到的水印出现错误也更加随机,观感上更容易看出水印信息。客观指标也可以用代码得到如下:
图14-强度0.02、直流、加置乱水印嵌入提取评估
  可以看到,客观指标上相对于图12的无置乱结果,PSNR、BER、Nc指标都要好一点点,但变化不大。
  总的来说,加入水印置乱的效果要更好,也更符合实际应用,我们接下来的参数分析和攻击分析也都会以此为基础进行研究。

五、关键参数的调试与分析

1.嵌入强度对嵌入提取质量的影响

  基于之前带有水印置乱的程序,现只改变嵌入强度,运行多组程序可得到如下表格:
(PS:Markdwon的表格太难用了,就直接放图片了哈qwq)
表格1-不同嵌入强度对指标的影响
  总体上看,随着嵌入强度的增大,对嵌入后载体图片的影响就越大,PSNR也呈现下降趋势,人眼更易察觉图片的异样;但对于水印的提取而言,BER更低、Nc更接近于1,提取的水印更接近于原始图。但是不同的图片,嵌入强度的选取也是不同的,因此选择合适的嵌入强度,对水印嵌入和提取是比较重要的,目前看来还是一开始使用的嵌入强度Ins=0.02比较适合该测试图片。

2.嵌入DCT系数区块不同对效果的影响

  由于上述测试都是对直流块进行嵌入,其水印特征对于空域的反映较为直观、明显,下面尝试嵌入不同的DCT系数区块,DCT系数区块的补充概念如下:
图15-DCT系数的Zig-Zag排列
  对DCT系数区块进行Zig-Zag扫描编号可以得到上图15,按照Zig-Zag扫描编号的顺序,序号从小到大依次是直流分量、交流分量的低频、中频、高频。
  尝试嵌入到中频分量序号25中,即修改代码为嵌入到BLOCK(4,4),由于在该频段下嵌入系数不能太小,否则嵌入不明显故改Ins=5,运行程序查看结果如下:
图16-强度5、交流(4,4)、不置乱水印嵌入提取结果
  仔细观察嵌入水印后图像,没有看到嵌入直流时那种明显色块,但放大看物体边缘可以看到有高频的噪声存在,是这种嵌入方法的特点。因此可以在不加入置乱的情况下,考虑将信息嵌入到交流块中,也能获得比较好的隐蔽性。
图17-嵌入交流分量放大对比图
  主观上看提取的水印,在嵌入系数为5的情况下,能有较为不错的效果,同样地也可以通过客观指标评估,见下图:
图18-强度5、交流(4,4)、不置乱水印嵌入提取评估

3.不同攻击对水印嵌入提取的影响

  下面基于之前嵌入强度0.02、直流、加置乱水印的程序,进行不同“攻击”,综合对比该方法的该参数条件下的鲁棒性。分别对嵌入水印后的图像实施:
①JPEG压缩:将嵌入水印后的图像.bmp使用Windows画图另存为有损压缩格式.jpg,再对压缩后的图像进行水印提取。
②添加白噪声:使用如下MATLAB代码对嵌入水印后的图像叠加白噪声,而后提取水印。

%图像加噪声代码
NoiseP = imread('BWMed.bmp'); %图片载入
Wnoise = 10*randn(size(NoiseP)); %产生白噪声
NoiseP = double(NoiseP) + Wnoise; %白噪声叠加

③涂鸦与剪裁:使用Windows画图工具对嵌入水印后的图像进行涂鸦与剪裁。
④高斯低通滤波:使用以下代码生成高斯低通滤波器模板,并执行滤波。

GSLPF = fspecial('gaussian',[4,4]); %高斯低通滤波模板
LB = imfilter(LB,GSLPF);

⑤均值滤波:使用如下代码生成均值滤波器模板,并执行滤波。

AVER = fspecial('average',[4,4]); %均值滤波模板
LB = imfilter(LB,AVER);

⑥USM锐化滤波:使用如下代码生成USM锐化滤波器模板,并执行滤波。

USM = fspecial('unsharp'); %USM锐化滤波器模板
LB = imfilter(LB,USM);

  根据上述操作或代码,可以得到下表的结果:
表2_1
表2_2
  综上,由于是向直流系数中嵌入水印信息,因此如果攻击的主要对象在低频,或者将高频削到与低频很相近、造成一定混乱,那么影响就会比较大,特别是低通、均值滤波。而对于涂鸦或裁剪而言,因为执行了水印置乱,所以整体上造成的像素破坏是随机的,使得其鲁棒性较好;涂鸦和剪裁所能造成的影响完全取决于其破坏的程度了。
  总体而言,基于DCT变换域的数字图像盲水印算法拥有比较好的不可察觉性和鲁棒性,配合水印置乱算法能实现较好的保密性,十分适合数字图像版权水印、数字图像密文信息传递等使用场景。还可以根据需求,选择不同的DCT系数区块去嵌入,这也可以视为另一层的加密。
  只不过,对于不同图片、不同DCT系数区块,适合的嵌入系数不同,这是要根据需求、实际效果具体实验的。因此这种方法对于通用图片批量水印嵌入可能就不太适合,当然若通过编程修改嵌入系数,对每个图片计算最佳评估指标也是可以的,能够实现辅助选取最佳嵌入系数。

六、参考文献

[1]张晓强,王蒙蒙,朱贵良.图像水印算法研究新进展[J].计算机工程与科学,2012,34(04):17-22.
[2]严正国,雷宇,王浩然,汤英.数字图像水印研究[J].工业控制计算机,2020,33(11):100-102.
[3]胡睿,徐正光.一种基于分块DCT变换和水印置乱的嵌入算法[J].微计算机信息,2005(07):29-31.
[4]张文彬. 基于数字水印的图像版权保护系统的设计与实现[D].中南大学,2010.
[5]王代福,李春杰.一种DCT域的灰度级数字水印算法研究[J].渤海大学学报(自然科学版),2010,31(02):177-182.
[6]李晨,喻枭,田丽华.基于熵编码的JPEG压缩域脆弱图像水印算法[J].计算机应用研究,2019,36(08):2424-2428+2435.
[7]张卫东,寻洋洋,范助文. DCT数字水印算法及Matlab仿真实现[EB/OL]. 北京:中国科技论文在线 [2010-01-15]. http://www.paper.edu.cn/releasepaper/content/201001-651.

尾记

相关的资料,包括这篇论文的pdf文档、MATLAB工程代码、资源文件,我放在链接里:点击跳转
再次强调一下,由于是赶出来的课程论文,不保证上述所言的严谨性和正确性,仅作学习交流和参考!