135.基于矩阵奇异值分解的水印算法

基于矩阵奇异值分解的水印算法

一.实验目的

了解基于矩阵奇异值分解的图像数字水印技术，掌握基于矩阵奇异值分解的图像水印算法原理，设计并实现一种基于矩阵奇异值分解的数字水印算法。

二.实验条件

(1) Windows 10或7操作系统；

(2) MATLAB 2014b；

(3)图像文件

三.实验原理

1.矩阵的奇异值分解（SVD）与图像矩阵的能量

矩阵的奇异值分解变换是一种正交变换，它可以将矩阵对角化。我们知道任何一个矩阵都有它的奇异值分解，对于奇异值分解可用下面的定理来描述。

定理

设A是一个秩为 r的 矩阵，则存在正交矩阵 和 ，使得 ，其中 ，这里 ，是矩阵对应的正特征值。称 为A的奇异值分解， （ ）称为A 的奇异值。

矩阵的 F（Frobenius）范数定义为，  由于所以  上式表明矩阵F 范数的平方等于矩阵的所有奇异值的平方和。对于一幅图像，通常用图像矩阵的 F范数来衡量图像的能量，所以图像的主要能量集中在矩阵哪些数值较大的奇异值上。

定理（Weyl定理）

设 A为一个大小为的矩阵，，是矩阵A 的一个扰动，假设矩阵 A，B 的奇异值分别为 和 ，，是矩阵 的最大奇异值，则有，其中 表示2—范数。

由Weyl定理可知，当图像被施加小的扰动时图像矩阵的奇异值的变化不会超过扰动矩阵的最大奇异值，因此基于矩阵奇异值分解的数字水印算法具有很好的稳定性，能够有效地抵御噪声对水印信息的干扰。由上述矩阵奇异值分解的性质可知，图像矩阵奇异值分解的稳定性非常好。当图像加入小的扰动时，其矩阵奇异值的变换不超过扰动矩阵的最大奇异值。基于矩阵奇异值分解的数字水印算法正是将想要嵌入的水印信息嵌入到图像矩阵的奇异值中，如果在嵌入水印的过程中选择一个嵌入强度因子来控制水印信息嵌入的程度，那么当嵌入强度因子足够小时，图像在视觉上不会产生明显的变化。

2  水印嵌入

设一副图像对应的矩阵A 大小为，需要嵌入的水印对应的矩阵 W大小为，自然地有 ，。在矩阵 A的左上角取一个大小为 的子块。首先对 进行奇异值分解，得到，其中 是的奇异值矩阵。我们的目标就是将水印 嵌入到矩阵中，在这里定义一个描述水印嵌入过程的参数，称为嵌入强度因子，则水印嵌入的过程表示为。可以看出，矩阵 包含了所有的水印信息，水印信息的能量反映在 的奇异值当中。对进行奇异值分解，得到 ，则反映了嵌入水印的图像的全部信息，由此得到子块 嵌入水印后的图像子块，这样就完成了水印的嵌入。

上述水印嵌入算法的基本过程可以表示为  嵌入强度因子 的意义在上式中一目了然，它衡量了水印对原图像的扰动情况。在水印嵌入时，选择合适的嵌入强度因子是十分重要的。小的嵌入强度因子有利于水印的透明性，但嵌入的水印信息容易受到外界噪声的干扰，如果噪声强度足够大，则可能使水印信息被噪声淹没而完全丢失，导致提取水印时无法得到水印的全部信息。大的嵌入强度因子有利于增强算法的鲁棒性，即使在噪声较强的情况下水印信息也不会受到很大的影响，但是过大的嵌入强度因子可能对原矩阵的奇异值产生较大的影响，有可能破坏水印的透明性，影响图像的质量。因此，在水印嵌入时要选择适当的嵌入强度因子使得水印图像的不可觉察性与鲁棒性达到最佳。

3  水印提取  

水印提取是上述水印嵌入过程的逆过程。在水印提取时，假设我们得到的是受扰动的图像矩阵，首先对 进行奇异值分解由此得到包含有全部水印信息的奇异值矩阵 ，然后利用水印嵌入时的矩阵 ，得到由上述水印嵌入的算法可得这样就得到了水印的信息，其中 为水印嵌入时所用的嵌入强度因子，为原图像 的奇异值矩阵。

四.实验步骤

1.打开matlab软件，新建脚本文件并保存

2.输入程序

3.编译通过后点击运行

五.实验代码

clc, clear
A=imread('tu5.bmp'); %读入载体文件
W=imread('tu6.bmp');  %读入水印文件
[m1,m2,m3]=size(W);  %给出矩阵W的维数
A0=A([1:m1],[1:m2],:); %在矩阵A的左上角选取与W同样大小的子块
A0=double(A0); W=double(W); %进行数据类型转换
a=0.05; %嵌入强度因子为0.05
for i=1:3
    [U1{i},S1{i},V1{i}]=svd(A0(:,:,i)); %对载体R，G，B层分别进行奇异值分解
    A1(:,:,i)=S1{i}+a*W(:,:,i);  %计算A1矩阵
    [U2{i},S2{i},V2{i}]=svd(A1(:,:,i)); %对A1的各层进行奇异值分解
    A2(:,:,i)=U1{i}*S2{i}*V1{i}'; %计算A2矩阵
end
AW=A;  %整体水印合成图片初始化
AW([1:m1],[1:m2],:)=A2; %左上角替换成水印合成子块，水印嵌入完成
AW=uint8(AW); W=uint8(W); %变换回原来的数据类型
subplot(1,3,1), imshow(A) %显示载体图片
subplot(1,3,2), imshow(W) %显示水印图片
subplot(1,3,3), imshow(AW) %显示嵌入水印的合成图片
%以下是水印的提出
AWstar=imnoise(AW,'gaussian',0,0.01); %加入高斯噪声
A2star=AWstar([1:m1],[1:m2],:);  %提出子块
A2star=double(A2star); %进行数据类型转换
for i=1:3
    [U3{i},S2star{i},V3{i}]=svd(A2star(:,:,i)); %奇异值分解
    A1star(:,:,i)=U2{i}*S2star{i}*V2{i}';  %计算A1*
    Wstar(:,:,i)=(A1star(:,:,i)-S1{i})/a;  %计算W*
end
for i=1:3
Wstar(:,:,i)=medfilt2(Wstar(:,:,i));  %对提取水印的R,G,B层分别进行中值滤波
end
Wstar=uint8(Wstar); %进行类型转换
figure, subplot(1,2,1), imshow(AWstar) %显示被噪声污染的合成图片
subplot(1,2,2), imshow(Wstar) %显示提出的水印图片