CN104732474A - 基于多级小波系数加权和量化的数字水印嵌入及提取方法 - Google Patents

Abstract

本发明提供一种基于多级小波系数加权和量化的数字水印嵌入及提取方法。所提方法通过掩体图像的MD5值和用户给定的密钥产生水印加密环节和嵌入环节的初始值和参数，用于驱动混沌映射生成2值伪随机矩阵来加密2值水印，并将混沌映射产生的中间参数应用到水印嵌入环节的系数选择中，从而将水印加密环节和嵌入环节与掩体图像和用户给定的密钥紧密耦合，提高了水印和水印嵌入的安全性。所提方法采用多级小波系数加权的嵌入策略，权重根据对JPEG压缩的小波系数绝对变化量的算术平均数进行设置，从而不同级的小波系数权重各不相同，减少了单个小波系数对水印的影响，进一步提高了图像视觉质量和嵌入水印的安全性。在水印提取时，通过滤除孤立黑点来提高水印提取质量，从而有效地减小水印的恢复误差。

Description

基于多级小波系数加权和量化的数字水印嵌入及提取方法

技术领域

本发明属于信息安全和数字图像信号处理交叉研究领域，涉及一种数字水印嵌入及提取方法，特别涉及一种基于多级小波系数加权和量化的数字水印嵌入及提取方法。

背景技术

多媒体和网络技术的发展，使人们能够更加方便地获得各种数字图像信息。数字图像具有易存储性、易传输性以及易操作性等优点，但同时也使得盗版者可以对其进行未授权的复制和传播；恶意攻击者可容易地对其进行篡改或伪造，从而易造成严重后果和导致重大经济损失。

针对图像信息安全，人们已提出多种方法，如：①将明文图像转换为密文图像的数字图像加密技术；②将机密图像隐藏于非相关载体中的数字图像隐写技术；③结合现有技术特点对图像的真实性和完整性进行甄别的数字图像盲取证技术；④将机密图像拆分为影子图像，利用部分分发影子图像重构秘密图像的数字图像分存技术以及⑤给图像添加版权认证标识，对图像的真实性和完整性，来源者和使用者进行鉴别认证的数字图像水印技术。其中数字图像水印技术，已成为图像信息安全的研究热点。

Van Schyndel等较早提出数字水印(A digital watermark.Proceedings of theIEEE International Conference on Image Processing,1994,2:86-90.)。

传统的水印技术通常关注的是水印如何嵌入，以减少对掩体视觉质量或听觉质量的影响，而对水印嵌入位置的安全性却较少涉及，通常采用将水印嵌入到频域变换的固定位置。如：Lai等将水印分成两部分，分别通过修改掩体图像1级Haar小波分解的LH和HL子带的奇异值来嵌入水印(Digital imagewatermarking using discrete wavelet transform and singular value decomposition.IEEE Transactions on Instrumentation and Measurement,2010,59(11):3060-3063.)；Tao等将经过仿射变换的水印嵌入到掩体图像3级小波分解的HL₃子带系数的各个分块中，并利用粒子群优化算法来优化水印嵌入强度(Awavelet-based particle swarm optimization algorithm for digital image watermarking.Integrated Computer-Aided Engineering,2012,19(1):81-91.)；Zope等将Arnold变换加密的水印嵌入到掩体图像3级小波分解的LL₃和HH₃频带的小波系数中(Robust copyright protection of raster images using wavelet based digitalwatermarking.2014IEEE International Geoscience and Remote Sensing Symposium(IGARSS),2014:3129-3132.)；Preda等用2级整数小波分解的LL₂频带的系数生成水印，对其进行随机置换后嵌入HL₁、LH₁和HH₁频带的小波系数中(Self-recovery of unauthentic images using a new digital watermarking approach inthe wavelet domain.2014 10th IEEE International Conference on Communications,2014:1-4.)。

以上方法关注的是将水印嵌入到载体中，使其满足视觉上的不可见性，而对水印嵌入策略的安全性没有涉及，采用固定的频域变换导致水印在掩体上的嵌入位置固定。一些方法尽管试图通过对水印进行置乱来增强对嵌入的水印进行保护,但依然使用固定的频域变换来进行水印嵌入，由此在掩体上的嵌入位置也依然固定。固定的嵌入位置导致水印嵌入策略一旦公开，所嵌入的水印就可被轻易地剔除或篡改掉。

为克服传统的水印技术使用固定的小波基来实现小波变换这一缺陷。Meerwald等将参数化小波变换应用到水印系统，提出了基于参数化离散小波变换的安全水印方法(Watermark security via wavelet filter parametrization.Proceedings of the IEEE International Conference on Image Processing,2001,3:1027-1030.)，其基本思想是通过构造参数化滤波器来试图实现不同的小波参数对应不同的变换域，并通过对参数进行保密而保护水印的嵌入位置，以提高数字水印的安全性。在此基础上，Wu等提出一种基于参数化整数小波的半脆弱水印图像认证方法(A secure semi-fragile watermarking for image authentication basedon integer wavelet transform with parameters.Proceedings of the Australasianworkshop on Grid computing and e-research,2005,44:75-80.)，为同时获得安全性和较低的计算复杂性，使用提升算法构造参数化整数小波变换。该方法将掩体图像进行3级参数化整数小波分解，然后通过修改低频子带LL₃中系数的最低五位嵌入水印。Chamlawi等在Wu的基础上添加了Piva等的方法(Self recoveryauthentication of images in the DWT domain.International Journal of Image andGraphics,2005,5(1):149-165.)提出了基于参数化小波的图像认证和恢复方法(Wavelet based image authentication and recovery.Journal of Computer Science andTechnology,2007,22(6):795-804.)，通过嵌入双重水印达到图像认证和恢复的功能。

但将参数化小波用于决定和影响水印的嵌入位置以提高水印嵌入策略的安全性，以上方法都没有进行理论上的分析，其所依据的只是实验个例。经实验检验，使用参数化整数小波对于参数值改变的敏感性难以满足实际需求，当参数化小波的参数发生改变时，依然能够近似提取嵌入的水印，因此嵌入的水印的安全性依然难以保证。

发明内容

本发明目的在于克服现有技术缺陷，提供一种基于多级小波系数加权和量化的数字水印嵌入及提取方法，通过掩体图像MD5值，用户密钥生成初始值和参数驱动混沌映射生成2值伪随机矩阵来加密2值水印，并将混沌映射产生的中间参数应用到水印嵌入环节的系数选择中，再利用选取的多级小波系数加权、量化嵌入水印，提高水印和水印嵌入的安全性。

为实现上述目的，本发明采用以下技术方案：

基于多级小波系数加权和量化的数字水印嵌入方法，包括以下步骤：

第1步：记掩体图像A＝(P_i,j)_m×n对应的128位MD5值为A_md5，从A_md5的左边最高位开始依次将其转换为32个16进制数构成16进制序列，记为S_a＝＜a₀,…a_k,…,a₃₁＞,a_k∈{0,1,…,15}；

第2步：由用户给定x_init∈(0,1)和μ_init∈[3.57,4]作为式(1)中的初始值x和参数μ并驱动式(1)产生长度为32的伪随机序列S_b＝＜b₀,…b_k,…,b₃₁＞,b_k∈(0,1)，在迭代产生S_b时，为消除暂态效应，将迭代前num≥0次结果抛弃，从num+1起产生S_b，num由用户给定；

x＝μx(1-x)            (1)

第3步：利用S_a和S_b计算用于对水印和嵌入位置进行保护的初始值x₀和参数μ₀；

第4步：记2值水印图为D＝(d_i,j)_mm×nn，将x₀和μ₀代入式(1)产生长度为mm×nn的伪随机序列S_c＝＜c₀,…,c_k,…,c_mm×nn-1＞,c_k∈(0,1)，同样为消除暂态效应，将迭代前num≥0次结果抛弃；

第5步：将伪随机序列S_c转化为2值矩阵E＝(e_i,j)_mm×nn，对D＝(d_i,j)_mm×nn按式(7)得到异或加密后的水印图D′＝(d′_i,j)_mm×nn；

$d_{i,j}^{\′}=d_{i,j}\⊕e_{i,j}---\left(7\right)$

第6步：对掩体图像A＝(P_i,j)_m×n进行l级小波金字塔变换形成3l+1个子带，包括1个低频子带LL_l和3l个高频子带HL_k,LH_k,HH_k,k＝1,2,…，l，l的取值由用户给定且满足HL_l和LH_l频带所有系数的总个数≥mm×nn；

第7步：利用x₀和μ₀产生水印嵌入环节的中间参数x_k,μ_k,k＝1,2,…l；

第8步：对第k级小波变换的HL_k和LH_k,k＝1,2,…,l频带所有系数，按行优先顺序扫描完HL_k再扫描LH_k，构成长度为len_k的小波系数序列，记为 $S_h^k=<h_0^k,...,h_i^k,...,h_{{\mathrm{len}}_k-1}^k>,{\mathrm{len}}_k=m\&times;n/2^{2k-1};$

第9步：对第k级小波变换，将x_k和μ_k代入式(1)产生长度为len_k的伪随机序列，记为 $S_r^k=<r_0^k,...,{r_i}^k,...,r_{{\mathrm{len}}_k-1}^k>,{r_i}^k\&Element;\left(\mathrm{0,1}\right),k=\mathrm{1,2},...,l,$ 同样为消除暂态效应，将迭代前num≥0次结果抛弃；

第10步：对k＝1,2,…，l按大小进行排序，记元素下标索引对应的置换为 $<0,...,i,...,{\mathrm{len}}_k-1>\&RightArrow;<s_0,...,s_i,...,s_{{\mathrm{len}}_k-1}>,$ 将其用于对置换，选取置换后的前mm×nn个小波系数，记为 表示选出的第k级小波系数的第ii个；

第11步：利用JPEG压缩估算筛选出的每级小波系数k＝1,2,…,l对应的权重w_k,k＝1,2,…,l；

第12步：按式(14)计算小波系数加权均值weight_ii,ii＝0,1,…,mm×nn-1；

${\mathrm{weight}}_{\mathrm{ii}}=\underset{k=1}{\overset{l}{\mathrm{\&Sigma;}}}{w}_k{v}_{\mathrm{ii}}^k---\left(14\right)$

第13步：通过修改小波系数加权均值weight_ii,ii＝0,1,…,mm×nn-1，嵌入水印；

第14步：将l级小波金字塔变换和嵌入水印后对应的矩阵记为A′＝(P′_i,j)_m×n，对A′＝(P′_i,j)_m×n其进行l级小波逆金字塔变换重构含水印图像A″＝(P″_i,j)_m×n。

作为进一步的优选方案，所述第3步中通过式(2)和式(3)计算用于水印和嵌入位置进行保护的初始值x₀和参数μ₀，mod1表示取小数运算；

$x_0=\left(\underset{k=0}{\overset{15}{\mathrm{\&Sigma;}}}{b}_k{\log }_{(15+\mathrm{\&alpha;}+\mathrm{\&beta;})}(a_k+\mathrm{\&alpha;})\right)\mathrm{mod}1---\left(2\right)$

${\mathrm{\&mu;}}_0=4-\left(\underset{k=16}{\overset{31}{\mathrm{\&Sigma;}}}{b}_k{\log }_{(15+\mathrm{\&alpha;}+\mathrm{\&beta;})}(a_k+\mathrm{\&alpha;})\right)\mathrm{mod}1\&times;0.1---\left(3\right)$

式(2)，式(3)中，a_k∈S_a,b_k∈S_b，α∈(1,15),β∈(0,15)由式(4)和式(5)产生，式(2)和式(3)的运算性质保证了x₀∈(0,1)，μ₀∈(3.9,4]，即式(1)处于混沌状态；

$\mathrm{\&alpha;}=(\left(\underset{k=0}{\overset{15}{\mathrm{\&Sigma;}}}{a}_k{b}_k\right)\&times;{10}^2)\mathrm{mod}13+b_0+1---\left(4\right)$

$\mathrm{\&beta;}=(\left(\underset{k=16}{\overset{31}{\mathrm{\&Sigma;}}}{a}_k{b}_k\right)\&times;{10}^2)\mathrm{mod}14+b_{31}---\left(5\right).$

作为进一步的优选方案，所述第5步中按式(6)将S_c转化为2值矩阵E＝(e_i,j)_mm×nn,对D＝(d_i,j)_mm×nn进行异或加密得到加密后的水印图D′＝(d′_i,j)_mm×nn；

$e_{i,j}=\left\{\begin{array}{cc}1,& c_k\&GreaterEqual;0.5\\ 0,& c_k<0.5\end{array}\right.,k=\mathrm{0,1},...,\mathrm{mm}\&times;\mathrm{nn}-1,i=k/\mathrm{mm},j=k\mathrm{mod}\mathrm{mm}---\left(6\right).$

作为进一步的优选方案，所述第6步对A＝(P_i,j)_m×n进行l级Haar小波金字塔变换形成3l+1个子带；

Haar小波变换公式如式(8)所示：

$\left[\begin{array}{cc}{f}_0^{\&prime;}& {f_1}^{\&prime;}\\ f_2^{\&prime;}& {f_3}^{\&prime;}\end{array}\right]=\frac{1}{2}\left[\begin{array}{cc}1& 1\\ 1& -1\end{array}\right]\left[\begin{array}{cc}{f}_0& f_1\\ f_2& f_3\end{array}\right]\left[\begin{array}{cc}1& 1\\ 1& -1\end{array}\right]---\left(8\right)$

式(8)中，f₀,f₁,f₂,f₃是变换前矩阵2×2小块上下左右相邻的4个元素，f′₀,f′₁,f′₂,f′₃为变换后的系数，将变换后的所有系数按所处的频带放置在一起，即构成对原矩阵元素的一级金字塔分解，其中f′₀,f′₁,f′₂,f′₃分别对应为LL,HL,LH,HH子带；

所述第14步中将l级小波金字塔变换和嵌入水印后对应的矩阵记为A′＝(P′_i,j)_m×n，对A′＝(P′_i,j)_m×n其进行l级小波逆金字塔变换重构含水印图像A″＝(P″_i,j)_m×n，其中l级Haar小波逆金字塔变换对应的2D Haar小波逆变换如式(17)所示:

$\left[\begin{array}{cc}{f}_0& f_1\\ f_2& f_3\end{array}\right]=\frac{1}{2}\left[\begin{array}{cc}1& 1\\ 1& -1\end{array}\right]\left[\begin{array}{cc}{f}_0^{\&prime;}& {f_1}^{\&prime;}\\ f_2^{\&prime;}& {f_3}^{\&prime;}\end{array}\right]\left[\begin{array}{cc}1& 1\\ 1& -1\end{array}\right]---\left(17\right).$

作为进一步的优选方案，所述第7步中将x₀和μ₀代入式(1)产生长度为2l的伪随机序列S_g＝＜g₀,…,g_i,…,g_2l-1＞,g_i∈(0,1)，同样为消除暂态效应，将迭代前num≥0次结果抛弃，并按式(9)，式(10)产生水印嵌入环节的中间参数x_k∈(0,1),μ_k∈(3.9,4],k＝1,2,…l；

x_k＝g_k-1               (9)

μ_k＝4-(g_k+l-1×10⁸)mod1×0.1            (10)。

作为进一步的优选方案，所述第11步中对A＝(P_i,j)_m×n进行l级Haar小波金字塔变换，记为对A＝(P_i,j)_m×n进行质量因子Q＝80的JPEG压缩之后再进行l级Haar小波金字塔变换，记为按式(11)求出第k级小波变换的HL_k和LH_k频带中每个系数的绝对变化量按式(12)求出第k级小波变换的HL_k和LH_k频带所有系数的绝对变化量的算术平均数按式(13)计算w_k，w_k按四舍五入法保留小数点后两位；

${\mathrm{\&sigma;}}_{i,j}^k=|P_{i,j}^1-P_{i,j}^2|,P_{i,j}^1,P_{i,j}^2\&Element;{\mathrm{HL}}_k|{\mathrm{LH}}_k---\left(11\right)$

${\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}^k=\frac{1}{{\mathrm{len}}_k}\underset{{\mathrm{\&sigma;}}_{i,j}^k\&Element;{\mathrm{HL}}_k|{\mathrm{LH}}_k}{\mathrm{\&Sigma;}}{\mathrm{\&sigma;}}_{i,j}^k---\left(12\right)$

$w_k={\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}^1/{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}^k---\left(13\right)$

式(12)中，表示落入HL_k和LH_k子带系数的变化量，len_k为HL_k和LH_k子带所有系数的数量，式(13)中，表示第1级小波变换的HL₁和LH₁频带所有系数的绝对变化量的算术平均数。

作为进一步的优选方案，所述第13步中通过对小波系数加权均值weight_ii,ii＝0,1,…,mm×nn-1按式(15)进行修改，按式(16)计算变化量Δ_ii，嵌入水印，q为量化步长：

i＝ii/mm,j＝iimodmm

$v_{\mathrm{ii}}^{\&prime;k}=v_{\mathrm{ii}}^k+2^{k-1}{\mathrm{\&Delta;}}_{\mathrm{ii}}/\underset{k=1}{\overset{l}{\mathrm{\&Sigma;}}}{2}^{k-1}{w}_k,{\mathrm{\&Delta;}}_{\mathrm{ii}}={\mathrm{weight}}_{\mathrm{ii}}^{\&prime;}-{\mathrm{weight}}_{\mathrm{ii}}---\left(16\right).$

基于多级小波系数加权和量化的数字水印提取方法，包括以下步骤：

第1步：输入掩体图像对应的128位MD5值A_md5，从左边最高位开始依次将其转换为32个16进制数构成16进制序列，记为S_a＝＜a₀,…a_k,…,a₃₁＞,a_k∈{0,1,…,15}；

第2步：输入x_init∈(0,1)和μ_init∈[3.57,4]作为式(1)的初始值x和参数μ并驱动式(1)产生长度为32的伪随机序列S_b＝＜b₀,…b_k,…,b₃₁＞,b_k∈(0,1)，在迭代产生S_b时，为消除暂态效应，将迭代前num≥0次结果抛弃，从num+1起产生S_b，num由用户输入；

x＝μx(1-x)                (1)

第3步：利用S_a和S_b计算用于对水印和嵌入位置进行保护的初始值x₀和参数μ₀；

第4步：利用x₀和μ₀产生水印嵌入环节的中间参数x_k,μ_k,k＝1,2,…l；

第5步：假设待提取水印的图像为A″＝(P″_i,j)_m×n，对其进行l级小波金字塔变换形成3l+1个子带，将变换后的所有系数按所处的频带放置在一起，即构成对原矩阵元素的一级金字塔分解，l的取值由用户给定；

第6步：对第k级小波变换的HL_k和LH_k,k＝1,2,…，l频带所有系数，按行优先顺序扫描完HL_k再扫描LH_k，构成长度为len_k的小波系数序列 $S_h^k=<h_0^k,...,h_i^k,...,h_{{\mathrm{len}}_k-1}^k>,{\mathrm{len}}_k=m\&times;n/2^{2k-1};$

第7步：对第k级小波变换，将x_k和μ_k代入式(1)产生长度为len_k的伪随机序列 $S_r^k=<r_0^k,...,{r_i}^k,...,r_{{\mathrm{len}}_k-1}^k>,{r_i}^k\&Element;\left(\mathrm{0,1}\right),k=\mathrm{1,2},...,l,$ 同样为消除暂态效应，将迭代前num≥0次结果抛弃；

第8步：对k＝1,2,…,l其按大小进行排序，记元素下标索引对应的置换为 $<0,...,i,...,{\mathrm{len}}_k-1>\&RightArrow;<s_0,...,s_i,...,s_{{\mathrm{len}}_k-1}>,$ 将其用于对置换，选取置换后的前mm×nn个小波系数，记为

第9步：将w_k和按式(14)计算小波系数加权均值weight_ii,ii＝0,1,…,mm×nn-1；

${\mathrm{weight}}_{\mathrm{ii}}=\underset{k=1}{\overset{l}{\mathrm{\&Sigma;}}}{w}_k{v}_{\mathrm{ii}}^k---\left(14\right)$

第10步：提取出加密的水印图D′＝(d′_i,j)_mm×nn；

第11步：将x₀和μ₀代入式(1)产生长度为mm×nn的伪随机序列S_c＝＜c₀,…,c_k,…,c_mm×nn-1＞,c_k∈(0,1)，同样为消除暂态效应，将迭代前num≥0次结果抛弃；

第12步：将S_c转化为2值矩阵E＝(e_i,j)_mm×nn,对加密的水印图D′＝(d′_i,j)_mm×nn根据式(19)求得解密后的水印图D＝(d_i,j)_mm×nn；

$d_{i,j}^{\&prime;}=d_{i,j}^{\&prime;}\&CirclePlus;e_{i,j}---\left(19\right)$

第13步：滤除水印图D＝(d_i,j)_mm×nn的孤立黑点。

作为进一步的优选方案，所述第10步中按式(18)提取出加密的水印图D′＝(d′_i,j)_mm×nn；

作为进一步的优选方案，所述第13步中按式(20)滤除水印图D＝(d_i,j)_mm×nn的孤立黑点；

${\hat{d}}_{i,j}=\left\{\begin{array}{cc}1,& d_{i,j}=0,\underset{(x,y)\&Element;{\mathrm{\&psi;}}_{i,j}}{\mathrm{\&Sigma;}}{d}_{x,y}=8\\ d_{i,j}& \end{array}\right.---\left(20\right)$

式(20)中，ψ_i,j＝{(i-1,j-1),(i-1,j),(i-1,j+1),(i,j-1),(i,j+1),(i+1,j-1),(i+1,j),(i+1,j+1)}。

本发明同现有技术相比优点如下：

①本发明给出的是一种嵌入和提取策略可公开的水印构造方法，在该方法中，将掩体图像的MD5值和用户给定的密钥产生水印加密环节和嵌入环节的初始值和参数，从而将水印加密环节和嵌入环节与掩体图像和用户给定的密钥紧密耦合。同传统水印不同，即使嵌入和提取算法完全公开，采用不同的密钥，不同的掩体图像，都将导致嵌入的水印经历不同的水印加密和嵌入环节，从而使得攻击者试图通过选择特殊的掩体图像猜测出水印的加密参数和找到水印对应的嵌入位置进而对水印进行篡改和修改的策略彻底失效。

②本发明给出的是一种基于小波加权均值的水印嵌入和提取方法，同现有方法嵌入不同，采用了多级小波系数加权的嵌入策略，权重根据对JPEG压缩的小波系数绝对变化量的算术平均数进行设置，从而不同级的小波系数权重各不相同，减少了单个小波系数对水印的影响，提高了图像的视觉质量和安全性。在水印提取时，通过滤除孤立黑点来提高水印的提取质量，从而有效地减小水印的恢复误差。

③本发明给出的是一种真正意义的敏感性水印构造方法，对用户提供的密钥和掩体图像特征极度敏感，并且这种敏感性是理论上可以证实和保证的。而使用参数化整数小波对于参数值改变的敏感性难以满足实际需求，经实验检验，当参数化小波的参数发生改变时，依然能够近似提取嵌入的水印，因此嵌入的水印的安全性依然难以保证。

④本发明给出的尽管是一种敏感性水印构造方法，但依然具备传统鲁棒水印的特征，对JPEG压缩、剪裁、噪声、覆盖、涂鸦具有一定的鲁棒性，因此相对于传统鲁棒水印方法，具备更高的安全性和更好的实际应用价值。

附图说明

图1是水印嵌入方法流程图；

图2是水印提取方法流程图；

图3是实施例：掩体图像，为512×512分辨率的8位标准灰度图像truck；

图4是实施例：嵌入水印后的含水印图像，为512×512分辨率的8位灰度图像，相对图3的PSNR＝46.77；

图5是实施例：水印图像，为64×64分辨率的8位2值图像；

图6是实施例：从图4中提取的水印图像，为64×64分辨率的8位2值图像；

图7是实施例：提供错误的MD5值(以16进制数表示为：0xb7fd46987317ca3c1141f7750ed0ebbc)提取的水印图像，为64×64分辨率的8位2值图像，相对图5的NC＝0.49；

图8是实施例：修改MD5值中1bit(以16进制数表示为：0x18b8b036c97091b4e8a09b5066e085f7)提取的水印图像，分辨率为64×64的8位2值图像，相对图5的NC＝0.49；

图9是实施例：提供错误的x_init值(将0.5836541938742647修改为0.5836541938742646)提取的水印图像，为64×64分辨率的8位2值图像，相对图5的NC＝0.51；

图10是实施例：提供错误的x_init值(将0.5836541938742647修改为0.5836541938742648)提取的水印图像，为64×64分辨率的8位2值图像，相对图5的NC＝0.50；

图11是实施例：提供错误的μ_init值(将3.9159734698471825修改为3.9159734698471824)提取的水印图像，为64×64分辨率的8位2值图像，相对图5的NC＝0.51；

图12是实施例：提供错误的μ_init值(将3.9159734698471825修改为3.9159734698471826)提取的水印图像，为64×64分辨率的8位2值图像，相对图5的NC＝0.52；

图13是实施例：提供错误的num值(将2000修改为1999)提取的水印图像，为64×64分辨率的8位2值图像，相对图5的NC＝0.51；

图14是实施例：提供错误的l值(将3修改为2)提取的水印图像，为64×64分辨率的8位2值图像，相对图5的NC＝0.50；

图15是实施例：按质量因子Q＝80的JPEG压缩后的含水印图像，为512×512分辨率的8位灰度图像，相对图3的PSNR＝36.54；

图16是实施例：剪裁后的含水印图像，为512×512分辨率的8位灰度图像，相对图3的PSNR＝11.94；

图17是实施例：添加0.4％白噪声后的含水印图像，为512×512分辨率的8位灰度图像，相对图3的PSNR＝28.54；

图18是实施例：覆盖后的含水印图像，为512×512分辨率的8位灰度图像，相对图3的PSNR＝21.64；

图19是实施例：涂鸦后的含水印图像，为512×512分辨率的8位灰度图像，相对图3的PSNR＝20.24；

图20是实施例：图4完全篡改水印后的图像，为512×512分辨率的8位灰度图像，相对图3的PSNR＝26.45；

图21是实施例：按质量因子Q＝80的JPEG压缩后的含水印图像中提取的水印图像，为64×64分辨率的8位2值图像，相对图5的NC＝0.97；

图22是实施例：剪裁后的含水印图像中提取的水印图像，为64×64分辨率的8位2值图像，相对图5的NC＝0.91；

图23是实施例：添加0.4％白噪声后的含水印图像中提取的水印图像，为64×64分辨率的8位2值图像，相对图5的NC＝0.87；

图24是实施例：覆盖后的含水印图像中提取的水印图像，为64×64分辨率的8位2值图像，相对图5的NC＝0.89；

图25是实施例：涂鸦后的含水印图像中提取的水印图像，为64×64分辨率的8位2值图像，相对图5的NC＝0.94；

具体实施方式

以下结合具体实施例对本发明进行说明：

以JAVA jdk1.7.0_09作为案例实施环境，结合附图对本发明实施方式进行详细说明，但不局限于本实施案例，其中图1是水印嵌入方法流程图，图2是水印提取方法流程图。

嵌入过程：

第1步：选取图像A如图3所示，为512×512的8位标准灰度图像truck.bmp，记掩体图像A＝(P_i,j)_512×512对应的128位MD5值为A_md5＝00011000101110001011000000110110110010010111000010010001101101001110100010100000100110110101000001100110111000001000010111111000，从A_md5的左边最高位开始依次将其转换为32个16进制数构成16进制序列，记为S_a＝＜1,8,11,8,11,0,3,6,12,9,7,0,9,1,11,4,14,8,10,0,9,11,5,0,6,6,14,0,8,5,15,8＞；

第2步：由用户给定x_init∈(0,1)和μ_init∈[3.57,4]作为式(1)的初始值x和参数μ并驱动式(1)产生长度为32的伪随机序列S_b＝＜b₀,…b_k,…,b₃₁＞,b_k∈(0,1)，在迭代产生S_b时，为消除暂态效应，将迭代前num≥0次结果抛弃，从num+1起产生S_b，num由用户给定；

例如：x_init＝0.5836541938742647，μ_init＝3.9159734698471825，num＝2000，产生长度为32的伪随机序列S_b＝＜b₀,…b_k,…,b₃₁＞,b_k∈(0,1)，b₀＝0.4454088439071509,b₁＝0.9673230055557697,b₂＝0.1237808218039777,b₃＝0.4247230754763154,b₄＝0.9568030520257117,b₅＝0.1618509885035459,b₆＝0.5312223444754627,b₇＝0.9751759402688756,b₈＝0.09479720355471674,b₉＝0.3360324001656252,b₁₀＝0.8737109569519057,b₁₁＝0.4320889851411571,b₁₂＝0.9609332661586223,b₁₃＝0.1470076966090752,b₁₄＝0.4910491077657985,b₁₅＝0.9786796256518280,b₁₆＝0.08170998182720710,b₁₇＝0.2938290414402912,b₁₈＝0.8125391815399567,b₁₉＝0.5964781811158950,b₂₀＝0.9425433319920187,b₂₁＝0.2120711069445983,b₂₂＝0.6543472330541910,b₂₃＝0.8857028638270253,b₂₄＝0.3964269203373903,b₂₅＝0.9369852208954993,b₂₆＝0.2312144114271247,b₂₇＝0.6960811518336486,b₂₈＝0.8284327317914145,b₂₉＝0.5565849089522140,b₃₀＝0.9664550002842699,b₃₁＝0.1269548131911206；

第3步：按式(4)，式(5)计算α,β，再按式(2)，式(3)计算x₀和μ₀，mod1表示取小数运算；

例如： $\mathrm{\&alpha;}=(\left(\underset{k=0}{\overset{15}{\mathrm{\&Sigma;}}}{a}_k{b}_k\right)\&times;{10}^2)\mathrm{mod}13+b_0+1=8.747832146115547,$ $\mathrm{\&beta;}=(\left(\underset{k=16}{\overset{31}{\mathrm{\&Sigma;}}}{a}_k{b}_k\right)\&times;{10}^2)\mathrm{mod}14+b_{31}=5.994651551546104,$ b₀log_(15+α+β)(a₀+α)＝0.4454088439071509×log_{29.742483697661651}(1+8.747832146115547),b₁log_(15+α+β)(a₁+α)＝0.9673230055557697×log_{29.742483697661651}(8+8.747832146115547),…,b₃₁log_(15+α+β)(a₃₁+α)＝0.1269548131911206×log_{29.742483697661651}(8+8.747832146115547)，按式(2)，式(3)计算可得到： $x_0=\left(\underset{k=0}{\overset{15}{\mathrm{\&Sigma;}}}{b}_k{\log }_{(15+\mathrm{\&alpha;}+\mathrm{\&beta;})}(a_k+\mathrm{\&alpha;})\right)\mathrm{mod}1=0.103478535807934540367326244225914,$ ${\mathrm{\&mu;}}_0=4-\left(\underset{k=16}{\overset{31}{\mathrm{\&Sigma;}}}{b}_k{\log }_{(15+\mathrm{\&alpha;}+\mathrm{\&beta;})}(a_k+\mathrm{\&alpha;})\right)\mathrm{mod}1\&times;0.1=3.9698439504568444257327859293425170;$

第4步：选取图像D如图5所示，为64×64分辨率的8位2值图像，将x₀和μ₀代入式(1)产生长度为64×64的伪随机序列S_c＝＜c₀,…,c_k,…,c₄₀₉₅＞,c_k∈(0,1)，将迭代前num＝2000次结果抛弃，其中，c₀＝0.4104962598890171,c₁＝0.9606588873239365,…,c₄₀₉₅＝0.3082923173578490；

第5步：按式(6)将S_c转化为2值矩阵E＝(e_i,j)_64×64，对D＝(d_i,j)_64×64按式(7)得到加密后的水印图D′＝(d′_i,j)_mm×nn；

例如：d_0,0＝1,d_0,1＝1,…,d_63,63＝1,e_0,0＝0,e_0,1＝1,…,e_63,63＝0, $d_{\mathrm{0,1}}^{\&prime;}=1\&CirclePlus;1=0,d_{\mathrm{63,63}}^{\&prime;}=1\&CirclePlus;0=1;$

第6步：若l＝3，对A＝(P_i,j)_m×n进行3级Haar小波金字塔变换形成10个子带，形成1个低频子带LL₃和9个高频子带HL_l,LH_l和HH_l,l＝1,2,3；

若A左上角的四个系数为 $\left[\begin{array}{cc}{f}_0& f_1\\ f_2& f_3\end{array}\right]=\left[\begin{array}{cc}77& 89\\ 93& 92\end{array}\right],$ 按式(8)进行1级Haar小波变换为 $\left[\begin{array}{cc}175.5& -5.5\\ -9.5& -6.5\end{array}\right]=\frac{1}{2}\left[\begin{array}{cc}1& 1\\ 1& -1\end{array}\right]\left[\begin{array}{cc}77& 89\\ 93& 92\end{array}\right]\left[\begin{array}{cc}1& 1\\ 1& -1\end{array}\right],$ 同理对所有系数块按式(8)进行运算，将变换后的所有系数按所处的频带放置在一起，即构成对原矩阵元素的1级金字塔分解，得到LL₁,HL₁,LH₁,HH₁四个频带的系数，再对左上角LL₁频带系数进行小波变换，得到LL₂,HL₂,LH₂,HH₂四个频带的系数，然后，对变换后的左上角LL₂频带系数进行小波变换，得到LL₃,HL₃,LH₃,HH₃四个频带的系数，完成3级Haar小波金字塔变换；

第7步：将x₀和μ₀代入式(1)产生长度为2l的伪随机序列S_g＝＜g₀,…,g_i,…,g_2l-1＞,g_i∈(0,1)，同样为消除暂态效应，将迭代前num≥0次结果抛弃，并按式(9)，式(10)产生水印嵌入环节的中间参数x_k∈(0,1),μ_k∈(3.9,4],k＝1,2,…l；

例如：产生长度为2l的伪随机序列S_g＝＜g₀,…,g_i,…,g_2l-1＞,g_i∈(0,1)，g₀＝0.4104962598890171,g₁＝0.9606588873239365,…,g₄₀₉₅＝0.03030914664336248，x₁＝g₀＝0.4104962598890171,μ₁＝4-(g₃×10⁸)mod1×0.1＝3.909866713，x₂＝g₁＝0.9606588873239365,μ₂＝4-(g₄×10⁸)mod1×0.1＝3.941740829，x₃＝g₂＝0.1500338587908365,μ₃＝4-(g₅×10⁸)mod1×0.1＝3.9335663752；

第8步：对第1级小波变换的HL₁和LH₁频带所有系数，若 ${\mathrm{HL}}_1=\left[\begin{array}{cccc}-5& 8.5& ...& 4.5\\ 0& -1.5& ...& -7.0\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ -2.5& -3.0& ...& -2.0\end{array}\right],{\mathrm{LH}}_1=\left[\begin{array}{cccc}-9.5& 6.5& ...& -5.5\\ -12.0& 8.5& ...& -15.0\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ 4.5& 9.0& ...& 15.0\end{array}\right],$ 按行优先顺序扫描完HL₁再扫描LH₁，构成长度为len₁＝131072的小波系数序列 $S_h^1=<-\mathrm{5,8.5},...,-2.0,-\mathrm{9.5,6.5},...,15.0>,$ 同理可得到

第9步：对第1级小波变换，将x₁和μ₁代入式(1)产生长度为len₁的伪随机序列将迭代前num次结果抛弃,其中， $r_0^1=0.5529855640821721,{r_1}^1=0.9664898447449297,...,r_{{\mathrm{len}}_1-1}^1=0.1911986945616335,$ 同理可产生

第10步：对按大小进行排序，元素下标索引对应的置换为＜0,1,…,131071＞→＜62130,118422,…,18452＞，＜0,1,…,32767＞→＜10957,29217,…,23957＞，＜0,1,…,8191＞→＜6377,2429,…,1747＞，将其用于对置换，选取置换后的前64×64个小波系数，将其记为 $S_v^1=<v_0^1,...,v_k^1,...,v_{4095}^1>,S_v^2=<v_0^2,...,v_k^2,...,v_{4095}^2>,S_v^3=<v_0^3,...,v_k^3,...,v_{4095}^3>,$ 其中， $v_0^1=-3.5,v_1^1=-7.5,...,v_{4095}^1=18.0,v_0^2=-2.5,v_1^2=-11.0,...,v_{4095}^2=-2.0,v_0^3=27.875,$ $v_1^3=-0.125,...,v_{4095}^3=20.125;$

第11步：对A＝(P_i,j)_512×512进行3级Haar小波金字塔变换，记为对A＝(P_i,j)_512×512进行质量因子Q＝80的JPEG压缩之后再进行3级Haar小波金字塔变换，记为按式(11)计算得 ${\mathrm{\&sigma;}}_{\mathrm{0,256}}^1=|P_{\mathrm{0,256}}^1-P_{\mathrm{0,256}}^2|=|-5.5-(-8.5)|=3,$ 类似求出第k级小波变换的HL_k和LH_k频带中每个系数的绝对变化量对第k级小波变换的HL_k和LH_k频带所有系数的绝对变化量求其算术平均数若算术平均数 ${\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}^1=3.0601730346679688,{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}^2=2.0460662841796875,{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}^3=1.416412353515625,$ 按式(13)计算w_k，其中 $w_1={\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}^1/{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}^1=1.00,w_2={\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}^1/{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}^2=1.50,w_3={\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}^1/{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}^3=2.16,$ w_k按四舍五入法保留小数点后两位；

第12步：按式(14)计算小波系数加权均值weight_ii,ii＝0,1,…,4095，其中，w₁,w₂,w₃为各级小波系数的权重；

例如： ${\mathrm{weight}}_0=w_1{v}_0^1+w_2{v}_0^2+w_3{v}_0^3=1.00\&times;-3.5+1.50\&times;-2.5+2.16\&times;27.875=52.96;$

第13步：通过对小波系数加权均值weight_ii,ii＝0,1,…,mm×nn-1按式(15)进行修改，按式(16)计算变化量Δ_ii，通过修改小波系数嵌入水印；

例如：若按式(15)进行计算按式(16)计算变化量Δ₀＝weight′₀-weight₀＝45.0-52.96＝-7.96，修改小波系数

$\begin{array}{c}{v}_0^{\&prime;1}=v_0^1+2^0{\mathrm{\&Delta;}}_0/\underset{l=1}{\overset{3}{\mathrm{\&Sigma;}}}{2}^{l-1}{w}_l=-3.5+(-7.96)/(2^0\&times;1.00+2^1\&times;1.50+2^2\&times;2.16)=\\ -4.129746835443038,v_0^{\&prime;2}=v_0^2+2^1{\mathrm{\&Delta;}}_0/\underset{l=1}{\overset{3}{\mathrm{\&Sigma;}}}{2}^{l-1}{w}_l=-2.5+2\&times;(-7.96)/(2^0\&times;1.00+\\ 2^1\&times;1.50+2^2\&times;2.16)=-3.759493670886076,v_0^{\&prime;3}=v_0^3+2^2{\mathrm{\&Delta;}}_0/\underset{l=1}{\overset{3}{\mathrm{\&Sigma;}}}{2}^{l-1}{w}_l=27.875+\\ 4\&times;(-7.96)/(2^0\&times;1.00+2^1\&times;1.50+2^2\&times;2.16)=25.356012658227849;\end{array}$

第14步：将3级Haar小波金字塔变换和按式(16)嵌入水印后对应的矩阵记为A′＝(P′_i,j)_512×512，对A′＝(P′_i,j)_512×512其进行3级Haar小波逆金字塔变换可重构含水印图像A″＝(P″_i,j)_512×512；

若A′的LL₃,HL₃,LH₃,HH₃频带左上角的系数分别为640.125，2.375，30.7753164556962，-18.625，将其构成矩阵 $\left[\begin{array}{cc}640.125& 2.375\\ 30.7753164556962& -18.625\end{array}\right],$ 按式(12)进行1级Haar小波逆金字塔变换为 $\left[\begin{array}{cc}327.3251582278481& 343.5751582278481\\ 315.1748417721519& 294.1748417721519\end{array}\right]=\frac{1}{2}\left[\begin{array}{cc}1& 1\\ 1& -1\end{array}\right]\left[\begin{array}{cc}640.125& 2.375\\ 30.7753164556962& -18.625\end{array}\right]\left[\begin{array}{cc}1& 1\\ 1& -1\end{array}\right],$ 将其作为LL₂频带左上角的四个系数，类似的运算可以得到所有LL₂频带的所有系数，对LL₂,HL₂,LH₂,HH₂频带系数重复上述步骤可得到LL₁频带的所有系数，再对LL₁,HL₁,LH₁,HH₁频带系数重复上述步骤可得到进行3级Haar小波逆金字塔变换的含水印图像A″＝(P″_i,j)_512×512，如图4所示。

基于多级小波系数加权和量化的数字水印提取方法，包括以下步骤：

第1步：输入掩体图像A＝(P_i,j)_m×n对应的128位MD5值为A_md5＝00011000101110001011000000110110110010010111000010010001101101001110100010100000100110110101000001100110111000001000010111111000，从A_md5的左边最高位开始依次将其转换为32个16进制数构成16进制序列，记为S_a＝＜1,8,11,8,11,0,3,6,12,9,7,0,9,1,11,4,14,8,10,0,9,11,5,0,6,6,14,0,8,5,15,8＞；

第2步：输入x_init∈(0,1)和μ_init∈[3.57,4]作为式(1)的初始值和参数并驱动式(1)产生长度为32的伪随机序列S_b＝＜b₀,…b_k,…,b₃₁＞,b_k∈(0,1)，在迭代产生S_b时，为消除暂态效应，将迭代前num≥0次结果抛弃，从num+1起产生S_b，num由用户输入；

例如：x_init＝0.5836541938742647，μ_init＝3.9159734698471825，num＝2000，产生长度为32的伪随机序列S_b＝＜b₀,…b_k,…,b₃₁＞,b_k∈(0,1)，b₀＝0.4454088439071509,b₁＝0.9673230055557697,b₂＝0.1237808218039777,b₃＝0.4247230754763154,b₄＝0.9568030520257117,b₅＝0.1618509885035459,b₆＝0.5312223444754627,b₇＝0.9751759402688756,b₈＝0.09479720355471674,b₉＝0.3360324001656252,b₁₀＝0.8737109569519057,b₁₁＝0.4320889851411571,b₁₂＝0.9609332661586223,b₁₃＝0.1470076966090752,b₁₄＝0.4910491077657985,b₁₅＝0.9786796256518280,b₁₆＝0.08170998182720710,b₁₇＝0.2938290414402912,b₁₈＝0.8125391815399567,b₁₉＝0.5964781811158950,b₂₀＝0.9425433319920187,b₂₁＝0.2120711069445983,b₂₂＝0.6543472330541910,b₂₃＝0.8857028638270253,b₂₄＝0.3964269203373903,b₂₅＝0.9369852208954993,b₂₆＝0.2312144114271247,b₂₇＝0.6960811518336486,b₂₈＝0.8284327317914145,b₂₉＝0.5565849089522140,b₃₀＝0.9664550002842699,b₃₁＝0.1269548131911206；

第3步：按式(4)，式(5)计算α,β，再按式(2)，式(3)计算x₀和μ₀，mod1表示取小数运算；

例如： $\mathrm{\&alpha;}=(\left(\underset{k=0}{\overset{15}{\mathrm{\&Sigma;}}}{a}_k{b}_k\right)\&times;{10}^2)\mathrm{mod}13+b_0+1=8.747832146115547,$ $\mathrm{\&beta;}=(\left(\underset{k=16}{\overset{31}{\mathrm{\&Sigma;}}}{a}_k{b}_k\right)\&times;{10}^2)\mathrm{mod}14+b_{31}=5.994651551546104,$ b₀log_(15+α+β)(a₀+α)＝0.4454088439071509×log_{29.742483697661651}(1+8.747832146115547),b₁log_(15+α+β)(a₁+α)＝0.9673230055557697×log_{29.742483697661651}(8+8.747832146115547),…,b₃₁log_(15+α+β)(a₃₁+α)＝0.1269548131911206×log_{29.742483697661651}(8+8.747832146115547)，按式(2)，式(3)计算得到： $x_0=\left(\underset{k=0}{\overset{15}{\mathrm{\&Sigma;}}}{b}_k{\log }_{(15+\mathrm{\&alpha;}+\mathrm{\&beta;})}(a_k+\mathrm{\&alpha;})\right)\mathrm{mod}1=0.103478535807934540367326244225914,$ ${\mathrm{\&mu;}}_0=4-\left(\underset{k=16}{\overset{31}{\mathrm{\&Sigma;}}}{b}_k{\log }_{(15+\mathrm{\&alpha;}+\mathrm{\&beta;})}(a_k+\mathrm{\&alpha;})\right)\mathrm{mod}1\&times;0.1=3.9698439504568444257327859293425170;$

第4步：将x₀和μ₀代入式(1)产生长度为2l的伪随机序列S_g＝＜g₀,…,g_i,…,g_2l-1＞,g_i∈(0,1)，同样为消除暂态效应，将迭代前num≥0次结果抛弃，并按式(9)，式(10)产生水印嵌入环节的中间参数x_k∈(0,1),μ_k∈(3.9,4],k＝1,2,…l；

例如：产生长度为2l的伪随机序列S_g＝＜g₀,…,g_i,…,g_2l-1＞,g_i∈(0,1)，g₀＝0.4104962598890171,g₁＝0.9606588873239365,…,g₄₀₉₅＝0.03030914664336248，x₁＝g₀＝0.4104962598890171,μ₁＝4-(g₃×10⁸)mod1×0.1＝3.909866713，x₂＝g₁＝0.9606588873239365,μ₂＝4-(g₄×10⁸)mod1×0.1＝3.941740829，x₃＝g₂＝0.1500338587908365,μ₃＝4-(g₅×10⁸)mod1×0.1＝3.9335663752；

第5步：假设待提取水印的图像A″＝(P″_i,j)_512×512如图4所示，对其进行3级Haar小波金字塔变换形成10个子带，形成1个低频子带LL₃和9个高频子带HL_l,LH_l和HH_l,l＝1,2,3；

若A″左上角的四个系数为 $\left[\begin{array}{cc}{f}_0& f_1\\ f_2& f_3\end{array}\right]=\left[\begin{array}{cc}78& 90\\ 94& 93\end{array}\right],$ 按式(8)进行1级Haar小波变换为 $\left[\begin{array}{cc}177.5& -5.5\\ -9.5& -6.5\end{array}\right]=\frac{1}{2}\left[\begin{array}{cc}1& 1\\ 1& -1\end{array}\right]\left[\begin{array}{cc}78& 90\\ 94& 93\end{array}\right]\left[\begin{array}{cc}1& 1\\ 1& -1\end{array}\right],$ 同理对所有系数块按式(8)进行运算，将变换后的所有系数按所处的频带放置在一起，即构成对原矩阵元素的1级金字塔分解，得到LL₁,HL₁,LH₁,HH₁四个频带的系数，再对左上角LL₁频带系数进行小波变换，得到LL₂,HL₂,LH₂,HH₂四个频带的系数，然后，对变换后的左上角LL₂频带系数进行小波变换，得到LL₃,HL₃,LH₃,HH₃四个频带的系数，完成3级Haar小波金字塔变换；

第5步：对第1级小波变换的HL₁和LH₁频带所有系数，若 ${\mathrm{HL}}_1=\left[\begin{array}{cccc}-5.5& 8.5& ...& 4.5\\ 0& -1.5& ...& -7\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ -2.5& -3.0& ...& -2.0\end{array}\right],{\mathrm{LH}}_1=\left[\begin{array}{cccc}-9.5& 6.5& ...& -5.5\\ -12.0& 8.5& ...& -15.0\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ 4.5& 9.0& ...& 15.0\end{array}\right],$ 按行优先顺序扫描完HL₁再扫描LH₁，构成长度为len₁＝131072的小波系数序列 $S_h^1=<-\mathrm{5.5,8.5},...,-2.0,-\mathrm{9.5,6.5},...,15.0>,$ 同理可得到

第6步：对第1级小波变换，将x₁和μ₁代入式(1)产生长度为len₁的伪随机序列将迭代前num次结果抛弃,其中， $r_0^1=0.5529855640821721,{r_1}^1=0.9664898447449297,...,r_{{\mathrm{len}}_1-1}^1=0.1911986945616335,$ 同理可产生

第7步：对按大小进行排序，元素下标索引对应的置换为＜0,1,…,131071＞→＜62130,118422,…,18452＞，＜0,1,…,32767＞→＜10957,29217,…,23957＞，＜0,1,…,8191＞→＜6377,2429,…,1747＞，将其用于对置换，选取置换后的前64×64个小波系数，将其记为 $S_v^1=<v_0^1,...,v_k^1,...,v_{4095}^1>,S_v^2=<v_0^2,...,v_k^2,...,v_{4095}^2>,S_v^3=<v_0^3,...,v_k^3,...,v_{4095}^3>,$ 其中， $v_0^1=-4.5,v_1^1=-4.5,...,v_{4095}^1=21.0,v_0^2=-4.5,v_1^2=-5.0,...,v_{4095}^2=10.0,v_0^3=26.625,$ $v_1^3=11.875,...,v_{4095}^3=33.125;$

第8步：将w_k和按式(14)计算小波系数加权均值weight_ii,ii＝0,1,…,4095；

例如：w₁＝1.00,w₂＝1.50,w₃＝2.16，按式(14)计算， ${\mathrm{weight}}_0=w_1{v}_0^1+w_2{v}_0^2+w_3{v}_0^3=1.00\&times;-4.5+1.50\&times;-4.5+2.16\&times;26.625=46.26;$

第9步：按式(18)提取出加密的水印图D′＝(d′_i,j)_64×64；

例如：

第10步：将x₀和μ₀代入式(1)产生长度为64×64的伪随机序列S_c＝＜c₀,…,c_k,…,c₄₀₉₅＞,c_k∈(0,1)，将迭代前num次结果抛弃，其中，c₀＝0.4104962598890171,c₁＝0.9606588873239365,…,c₄₀₉₅＝0.3082923173578490；

第11步：按式(6)将S_c转化为2值矩阵E＝(e_i,j)_64×64,对加密的水印图D′＝(d′_i,j)_64×64，根据式(19)求得解密后的水印图D＝(d_i,j)_64×64；

例如：d′_0，0＝1，e_0，0＝0，

第12步：按式(20)滤除水印图D＝(d_i,j)_64×64的孤立黑点，如图6所示。

例如：d_0,0＝1，按式(20)，则

为验证所提出的方法对掩体图像MD5值，用户密钥的敏感性，图7～14分别为修改其中之一相对应的提取的水印；从提取的水印可以看出本发明方法，对用户提供的密钥和掩体图像特征极度敏感。

为验证所提出的方法对一些图像攻击操作的鲁棒性，图15～19分别为经过JPEG压缩、剪裁、噪声、覆盖、涂鸦的含水印图，图21～25分别为相对应的提取的水印；从提取的水印可以看出本发明方法，依然具备传统鲁棒水印的特征，因此相对于传统鲁棒水印方法，具备更高的安全性和更好的实际应用价值。

为验证所提出的方法当用户试图完全篡改水印时，嵌入水印后图像视觉质量严重下降。图20为用户试图完全剔除嵌入的水印，即将图4各级Haar小波分解HL子带和LH子带系数均置为0得到的图像；从得到的图像可以看出本发明方法，若用户试图完全篡改水印，只能得到严重降质图像。

Claims (10)

1.基于多级小波系数加权和量化的数字水印嵌入方法，其特征在于包括以下步骤：

第1步：记掩体图像A＝(P_i,j)_m×n对应的128位MD5值为A_md5，从A_md5的左边最高位开始依次将其转换为32个16进制数构成16进制序列，记为S_a＝＜a₀,…a_k,…,a₃₁＞,a_k∈{0,1,…,15}；

第2步：由用户给定x_init∈(0,1)和μ_init∈[3.57,4]作为式(1)中的初始值x和参数μ并驱动式(1)产生长度为32的伪随机序列S_b＝＜b₀,…b_k,…,b₃₁＞,b_k∈(0,1)，在迭代产生S_b时，为消除暂态效应，将迭代前num≥0次结果抛弃，从num+1起产生S_b，num由用户给定；

x＝μx(1-x)                   (1)

第3步：利用S_a和S_b计算用于对水印和嵌入位置进行保护的初始值x₀和参数μ₀；

第4步：记2值水印图为D＝(d_i,j)_mm×nn，将x₀和μ₀代入式(1)产生长度为mm×nn的伪随机序列S_c＝＜c₀,…,c_k,…,c_mm×nn-1＞,c_k∈(0,1)，同样为消除暂态效应，将迭代前num≥0次结果抛弃；

第5步：将伪随机序列S_c转化为2值矩阵E＝(e_i,j)_mm×nn，对D＝(d_i,j)_mm×nn按式(7)得到异或加密后的水印图D′＝(d′_i,j)_mm×nn；

$d_{i,j}^{\&prime;}=d_{i,j}\&CirclePlus;e_{i,j}---\left(7\right)$

第6步：对掩体图像A＝(P_i,j)_m×n进行l级小波金字塔变换形成3l+1个子带，包括1个低频子带LL_l和3l个高频子带HL_k,LH_k,HH_k,k＝1,2,…,l，l的取值由用户给定且满足HL_l和LH_l频带所有系数的总个数≥mm×nn；

第7步：利用x₀和μ₀产生水印嵌入环节的中间参数x_k,μ_k,k＝1,2,…l；

第8步：对第k级小波变换的HL_k和LH_k,k＝1,2,…,l频带所有系数，按行优先顺序扫描完HL_k再扫描LH_k，构成长度为len_k的小波系数序列，记为 $S_h^k=<h_0^k,...,h_i^k,...,h_{{\mathrm{len}}_k-1}^k>,$ len_k＝m×n/2^2k-1；

第9步：对第k级小波变换，将x_k和μ_k代入式(1)产生长度为len_k的伪随机序列，记为 $S_r^k=<r_0^k,...,r_i^k,...,r_{{\mathrm{len}}_k-1}^k>,$ $r_i^k\&Element;\left(\mathrm{0,1}\right),k=\mathrm{1,2},...,l,$ 同样为消除暂态效应，将迭代前num≥0次结果抛弃；

第10步：对按大小进行排序，记元素下标索引对应的置换为 $<0,...,i,...,{\mathrm{len}}_k-1>\&RightArrow;<s_0,...,s_i,...,s_{{\mathrm{len}}_k-1}>,$ 将其用于对置换，选取置换后的前mm×nn个小波系数，记为 表示选出的第k级小波系数的第ii个；

第11步：利用JPEG压缩估算筛选出的每级小波系数对应的权重w_k,k＝1,2,…,l；

第12步：按式(14)计算小波系数加权均值weight_ii,ii＝0,1,…,mm×nn-1；

${\mathrm{weight}}_{\mathrm{ii}}=\underset{k=1}{\overset{l}{\mathrm{\&Sigma;}}}{w}_k{v}_{\mathrm{ii}}^k---\left(14\right)$

第13步：通过修改小波系数加权均值weight_ii,ii＝0,1,…,mm×nn-1，嵌入水印；

第14步：将l级小波金字塔变换和嵌入水印后对应的矩阵记为A′＝(P′_i,j)_m×n，对A′＝(P′_i,j)_m×n其进行l级小波逆金字塔变换重构含水印图像A″＝(P″_i,j)_m×n。

2.如权利要求1所述的基于多级小波系数加权和量化的数字水印嵌入方法，其特征在于：所述第3步中通过式(2)和式(3)计算用于水印和嵌入位置进行保护的初始值x₀和参数μ₀，mod1表示取小数运算；

$x_0=\left(\underset{k=0}{\overset{15}{\mathrm{\&Sigma;}}}{b}_k{\log }_{(15+\mathrm{\&alpha;}+\mathrm{\&beta;})}(a_k+\mathrm{\&alpha;})\right)\mathrm{mod}1---\left(2\right)$

${\mathrm{\&mu;}}_0=4-\left(\underset{k=16}{\overset{31}{\mathrm{\&Sigma;}}}{b}_k{\log }_{(15+\mathrm{\&alpha;}+\mathrm{\&beta;})}(a_k+\mathrm{\&alpha;})\right)\mathrm{mod}1\&times;0.1---\left(3\right)$

式(2)，式(3)中，a_k∈S_a,b_k∈S_b，α∈(1,15),β∈(0,15)由式(4)和式(5)产生，式(2)和式(3)的运算性质保证了x₀∈(0,1)，μ₀∈(3.9,4]，即式(1)处于混沌状态；

$\mathrm{\&alpha;}=(\left(\underset{k=0}{\overset{15}{\mathrm{\&Sigma;}}}{a}_k{b}_k\right)\&times;{10}^2)\mathrm{mod}13+b_0+1---\left(4\right)$

$\mathrm{\&beta;}=(\left(\underset{k=0}{\overset{31}{\mathrm{\&Sigma;}}}{a}_k{b}_k\right)\&times;{10}^2)\mathrm{mod}14+b_{31}---\left(5\right).$

3.如权利要求1所述的基于多级小波系数加权和量化的数字水印嵌入方法，其特征在于：所述第5步中按式(6)将S_c转化为2值矩阵E＝(e_i,j)_mm×nn,对D＝(d_i,j)_mm×nn进行异或加密得到加密后的水印图D′＝(d′_i,j)_mm×nn；

$e_{i,j}=\left\{\begin{array}{cc}1,& c_k\&GreaterEqual;0.5\\ 0,& c_k<0.5\end{array}\right.,k=\mathrm{0,1},...,\mathrm{mm}\&times;\mathrm{nn}-1,i=k/\mathrm{mm},j=k\mathrm{mod}\mathrm{mm}---\left(6\right).$

4.如权利要求1所述的基于多级小波系数加权和量化的数字水印嵌入方法，其特征在于：所述第6步对A＝(P_i,j)_m×n进行l级Haar小波金字塔变换形成3l+1个子带；

Haar小波变换公式如式(8)所示：

$\left[\begin{array}{cc}{f}_0^{\&prime;}& f_1^{\&prime;}\\ f_2^{\&prime;}& f_3^{\&prime;}\end{array}\right]=\frac{1}{2}\left[\begin{array}{cc}1& 1\\ 1& -1\end{array}\right]\left[\begin{array}{cc}{f}_0& f_1\\ f_2& f_3\end{array}\right]\left[\begin{array}{cc}1& 1\\ 1& -1\end{array}\right]---\left(8\right)$

式(8)中，f₀,f₁,f₂,f₃是变换前矩阵2×2小块上下左右相邻的4个元素，f₀′,f₁′,f₂′,f₃′为变换后的系数，将变换后的所有系数按所处的频带放置在一起，即构成对原矩阵元素的一级金字塔分解，其中f₀′,f₁′,f₂′,f₃′分别对应为LL,HL,LH,HH子带；

所述第14步中将l级小波金字塔变换和嵌入水印后对应的矩阵记为A′＝(P′_i,j)_m×n，对A′＝(P′_i,j)_m×n其进行l级小波逆金字塔变换重构含水印图像A″＝(P″_i,j)_m×n，其中l级Haar小波逆金字塔变换对应的2D Haar小波逆变换如式(17)所示:

$\left[\begin{array}{cc}{f}_0& f_1\\ f_2& f_3\end{array}\right]=\frac{1}{2}\left[\begin{array}{cc}1& 1\\ 1& -1\end{array}\right]\left[\begin{array}{cc}{f}_0^{\&prime;}& f_1^{\&prime;}\\ f_2^{\&prime;}& f_3^{\&prime;}\end{array}\right]\left[\begin{array}{cc}1& 1\\ 1& -1\end{array}\right]---\left(17\right).$

5.如权利要求1所述的基于多级小波系数加权和量化的数字水印嵌入方法，其特征在于：所述第7步中将x₀和μ₀代入式(1)产生长度为2l的伪随机序列S_g＝＜g₀,…,g_i,…,g_2l-1＞,g_i∈(0,1)，同样为消除暂态效应，将迭代前num≥0次结果抛弃，并按式(9)，式(10)产生水印嵌入环节的中间参数x_k∈(0,1),μ_k∈(3.9,4],k＝1,2,…l；

x_k＝g_k-1                     (9)

μ_k＝4-(g_k+l-1×10⁸)mod1×0.1           (10)。

6.如权利要求1所述的基于多级小波系数加权和量化的数字水印嵌入方法，其特征在于：所述第11步中对A＝(P_i,j)_m×n进行l级Haar小波金字塔变换，记为对A＝(P_i,j)_m×n进行质量因子Q＝80的JPEG压缩之后再进行l级Haar小波金字塔变换，记为按式(11)求出第k级小波变换的HL_k和LH_k频带中每个系数的绝对变化量按式(12)求出第k级小波变换的HL_k和LH_k频带所有系数的绝对变化量的算术平均数按式(13)计算w_k，w_k按四舍五入法保留小数点后两位；

${\mathrm{\&sigma;}}_{i,j}^k=|P_{i,j}^1-P_{i,j}^2|,P_{i,j}^1,P_{i,j}^2\&Element;{\mathrm{HL}}_k|{\mathrm{LH}}_k---\left(11\right)$

${\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}^k=\frac{1}{{\mathrm{len}}_k}\underset{{\mathrm{\&sigma;}}_{i,j}^k\&Element;{\mathrm{HL}}_k|{\mathrm{LH}}_k}{\mathrm{\&Sigma;}}{\mathrm{\&sigma;}}_{i,j}^k---\left(12\right)$

$w_k={\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}^1/{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}^k---\left(13\right)$

式(12)中，表示落入HL_k和LH_k子带系数的变化量，len_k为HL_k和LH_k子带所有系数的数量，式(13)中，表示第1级小波变换的HL₁和LH₁频带所有系数的绝对变化量的算术平均数。

7.如权利要求1所述的基于多级小波系数加权和量化的数字水印嵌入方法，其特征在于：所述第13步中通过对小波系数加权均值weight_ii,ii＝0,1,…,mm×nn-1按式(15)进行修改，按式(16)计算变化量Δ_ii，嵌入水印，q为量化步长：

i＝ii/mm,j＝iimodmm

$v_{\mathrm{ii}}^{\&prime;k}=v_{\mathrm{ii}}^k+2^{k-1}{\mathrm{\&Delta;}}_{\mathrm{ii}}/\underset{k=1}{\overset{l}{\mathrm{\&Sigma;}}}{2}^{k-1}{w}_k,{\mathrm{\&Delta;}}_{\mathrm{ii}}={\mathrm{weight}}_{\mathrm{ii}}^{\&prime;}-{\mathrm{weight}}_{\mathrm{ii}}---\left(16\right).$

8.与权利要求1对应的基于多级小波系数加权和量化的数字水印提取方法，其特征在于包括以下步骤：

第1步：输入掩体图像对应的128位MD5值A_md5，从左边最高位开始依次将其转换为32个16进制数构成16进制序列，记为S_a＝＜a₀,…a_k,…,a₃₁＞,a_k∈{0,1,…,15}；

第2步：输入x_init∈(0,1)和μ_init∈[3.57,4]作为式(1)的初始值x和参数μ并驱动式(1)产生长度为32的伪随机序列S_b＝＜b₀,…b_k,…,b₃₁＞,b_k∈(0,1)，在迭代产生S_b时，为消除暂态效应，将迭代前num≥0次结果抛弃，从num+1起产生S_b，num由用户输入；

x＝μx(1-x)                  (1)

第3步：利用S_a和S_b计算用于对水印和嵌入位置进行保护的初始值x₀和参数μ₀；

第4步：利用x₀和μ₀产生水印嵌入环节的中间参数x_k,μ_k,k＝1,2,…l；

第5步：假设待提取水印的图像为A″＝(P″_i,j)_m×n，对其进行l级小波金字塔变换形成3l+1个子带，将变换后的所有系数按所处的频带放置在一起，即构成对原矩阵元素的一级金字塔分解，l的取值由用户输入；

第6步：对第k级小波变换的HL_k和LH_k,k＝1,2,…,l频带所有系数，按行优先顺序扫描完HL_k再扫描LH_k，构成长度为len_k的小波系数序列 $S_h^k=<h_0^k,...,h_i^k,...,h_{{\mathrm{len}}_k-1}^k>,$ len_k＝m×n/2^2k-1；

第7步：对第k级小波变换，将x_k和μ_k代入式(1)产生长度为len_k的伪随机序列 $S_r^k=<r_0^k,...,r_i^k,...,r_{{\mathrm{len}}_k-1}^k>,$ $r_i^k\&Element;\left(\mathrm{0,1}\right),k=\mathrm{1,2},...,l,$ 同样为消除暂态效应，将迭代前num≥0次结果抛弃；

第8步：对其按大小进行排序，记元素下标索引对应的置换为 $<0,...,i,...,{\mathrm{len}}_k-1>\&RightArrow;<s_0,...,s_i,...,s_{{\mathrm{len}}_k-1}>,$ 将其用于对置换，选取置换后的前mm×nn个小波系数，记为

第9步：将w_k和按式(14)计算小波系数加权均值weight_ii,ii＝0,1,…,mm×nn-1；

${\mathrm{weight}}_{\mathrm{ii}}=\underset{k=1}{\overset{l}{\mathrm{\&Sigma;}}}{w}_k{v}_{\mathrm{ii}}^k---\left(14\right)$

第10步：提取出加密的水印图D′＝(d′_i,j)_mm×nn；

第11步：将x₀和μ₀代入式(1)产生长度为mm×nn的伪随机序列S_c＝＜c₀,…,c_k,…,c_mm×nn-1＞,c_k∈(0,1)，同样为消除暂态效应，将迭代前num≥0次结果抛弃；

第12步：将S_c转化为2值矩阵E＝(e_i,j)_mm×nn,对加密的水印图D′＝(d′_i,j)_mm×nn根据式(19)求得解密后的水印图D＝(d_i,j)_mm×nn；

$d_{i,j}=d_{i,j}^{\&prime;}\&CirclePlus;e_{i,j}---\left(19\right)$

第13步：滤除水印图D＝(d_i,j)_mm×nn的孤立黑点。

9.如权利要求8所述的基于多级小波系数加权和量化的数字水印提取方法，其特征在于：所述第10步中按式(18)提取出加密的水印图D′＝(d′_i,j)_mm×nn；

                                 (18)。

i＝ii/mm,j＝iimodmm

10.如权利要求8所述的基于多级小波系数加权和量化的数字水印提取方法，其特征在于：所述第13步中按式(20)滤除水印图D＝(d_i,j)_mm×nn的孤立黑点；

${\hat{d}}_{i,j}=\left\{\begin{array}{cc}1,& d_{i,j}=0,\underset{(x,y)\&Element;{\mathrm{\&psi;}}_{i,j}}{\mathrm{\&Sigma;}}{d}_{x,y}=0\\ d_{i,j}& \end{array}\right.---\left(20\right)$

式(20)中，ψ_i,j＝{(i-1,j-1),(i-1,j),(i-1,j+1),(i,j-1),(i,j+1),(i+1,j-1),(i+1,j),(i+1,j+1)}。