基于二元Weibull分布的非下采样Shearlet域图像水印算法

作为图像作品版权保护的有效手段，数字水印技术已成为国际学术界研究的热点之一.不可感知性、鲁棒性、水印容量是衡量一个图像水印算法优劣的最重要指标[1]，而这3方面却又存在着固有的相互矛盾关系，三者之间的最佳平衡也成为图像水印算法所共同追求的目标.可有效保持不可感知性、鲁棒性、水印容量之间良好平衡的高性能数字图像水印算法研究是一项富有挑战性的工作[2].

数字图像水印技术发展到今天，已有大量不同的图像水印算法.而依据水印嵌入方式，人们习惯将数字图像水印算法划分为3类，即加性水印、乘性水印、量化水印[3-6].相比之下，乘性嵌入方式在保证含水印图像质量在可接受水平内，能够接受更大的水印嵌入强度和水印容量，即：乘性嵌入方法在保证不可感知性的同时，能够使得水印的嵌入强度随着载体信号强度成一定比例变化，让人更不容易察觉，获得更高的鲁棒性和更大的水印容量，从而最接近于三者之间的良好平衡.一般说来，为全面提高图像水印的不可感知性、鲁棒性和水印容量，应该充分结合人眼视觉掩蔽特性与图像自身统计特性而进行数字水印信息的嵌入与检测.基于统计模型的变换域乘性水印较好地体现了上述思想，为有效解决不可感知性、鲁棒性、水印容量之间良好平衡问题提供了可能的解决方向[7].基于统计模型的变换域乘性水印工作原理为[7]：在水印嵌入环节，利用简单的乘性策略调制原始载体图像信号，以保证水印嵌入强度与载体图像信号强度成一定比例变化，让人更不容易察觉，从而最大程度地平衡鲁棒性和不可感知性；在水印检测环节，结合能够充分体现图像自身特性的多尺度变换特性(多分辨率性、能量聚集性等)，有效利用变换系数的统计特性而构造数字水印检测器，以检测和提取数字水印信息.其中，基于统计模型的水印检测器构造是整个变换域乘性图像水印方案的核心.

截止到目前，人们主要采用2类统计模型设计变换域乘性图像水印方案，分别为多尺度联合统计模型[8](multiscale joint statistical model)和多尺度边缘分布模型[9](multiscale marginal distribution model).

多尺度联合统计模型也称多尺度统计相关模型，其着重考察多尺度变换系数之间的相互关系，利用基于多尺度变换的各种数据结构，如隐Markov模型(hidden Markov models, HMM)描述并确定多尺度变换系数的隐含状态(隐含未知参数)，进而通过隐含状态联系建立多尺度变换系数关联，并进一步利用隐含未知参数来分析处理图像(如图像水印检测器构造).Wang等人[10]提出了基于小波域先验HMM模型参数的数字图像水印算法.该算法首先对原始载体图像进行3级小波分解并利用高频子带构造系数向量树，然后结合小波Watson视觉失真测度与先验HMM模型参数自适应确定水印嵌入强度，进而将水印信息嵌入到小波系数向量树中.在水印检测阶段，该算法结合泰勒级数近似局部最优检验理论，利用先验HMM模型参数构造了数字水印检测器.Amini等人[11]结合小波域向量HMM联合统计建模理论，提出了一种局部最优盲图像水印检测算法.该算法选取具有最大系数方差的高频子带作为水印嵌入区，并采用线性乘性方法嵌入水印信息.在水印提取阶段，算法首先利用向量HMM对小波系数进行了多尺度联合统计建模，然后结合对数似然比检验理论构造了数字图像水印检测器.总体说来，现有多尺度联合统计模型普遍采用高斯混合模型(Gaussian mixture model, GMM)描述变换系数的边缘分布，并利用隐Markov树捕捉变换系数间的相关特性.然而仿真实验表明，对具有明显高尖峰重拖尾分布的变换系数而言，高斯混合模型难以精确描述其非高斯分布特性.也就是说，现有基于GMM的多尺度联合统计模型尚无法高效捕获变换系数的边缘分布与变换系数间的相关特性，而这必然严重影响图像水印的检测性能.

多尺度边缘分布模型也称多尺度随机过程模型，其是将各多尺度变换系数看作是一个随机变量，进而考察各单个多尺度变换系数的边缘分布特点，即利用各种概率密度函数(probability density function, PDF)反映单个多尺度变换系数分布，描述多尺度变换系数统计特性，并进一步利用概率密度函数参数分析处理图像(如图像水印检测器构造).Akhaee等人[12]选择将水印信息嵌入到Contourlet域高频子带内，同时利用广义高斯分布(generalized Gaussian distribution, GGD)对Contourlet域高频系数进行统计建模，并进一步结合极大似然方法估计GGD参数并构造数字水印检测器.然而，该算法估计GGD模型参数时需要原始载体信号，故不利于实际应用.Bian等人[13]利用贝塞尔K(Bessel K form, BKF)概率密度函数对图像小波系数进行边缘分布建模，并根据极大似然理论构造了局部最优水印检测器.同时，结合“渐进相对效率”和“最大渐进相对效率”等客观评价指标证明了该检测器的有效性.然而，该局部最优水印检测器的抗攻击能力较弱.Etemad等人[14]首先依据Kolmogorov-Smirnov检验理论，证明t location-scale概率密度函数可高效描述Contourlet系数的边缘分布特点，进而结合似然比检验与t location-scale分布理论，构造了最优乘性数字水印检测器.Sadreazami等人[15]利用正态逆高斯(normal inverse Gaussian, NIG)分布刻画图像Contourlet域系数统计特性，并使用局部窗口自适应估计每个子带的NIG模型参数，同时结合Bayesian最大后验概率准则构造了水印检测器.但NIG分布无封闭形式解，且参数估计较复杂.Wang等人[16]结合非下采样剪切波变换与BKF统计建模理论，首先通过乘性方法将水印信息嵌入到非下采样剪切波变换重要系数中；然后，利用BKF边缘分布对非下采样剪切波变换域系数进行建模，并结合子带内系数相关性估计出BKF模型参数；最后，结合BKF统计模型与局部最大势能检验理论构造出局部最优检测器并进行水印提取.张悦等人[17]提出了一种基于共轭对称列率复数哈达码变换(conju-gate symmetric sequency-ordered complex Hadamard transform, CS-SCHT)的数字图像水印算法.该算法首先选择CS-SCHT系数幅值作为载体并进行乘性水印嵌入，然后结合实验对CS-SCHT系数幅值分布进行分析并确定出其近似分布——韦伯(Weibull)分布，最后设计出基于Weibull分布的乘性水印局部最优检测器.Dong等人[18]提出了一种基于Weibull边缘分布的DCT域图像水印算法.该算法利用Weibull分布对DCT直流系数和交流系数进行了统计建模，并结合极大似然方法构造出了数字水印检测器.由于上述多尺度边缘分布模型图像水印仅仅结合各种概率密度函数，考察了各单个多尺度变换系数的边缘分布特点，而忽略了多尺度变换系数间的重要相关特性(如尺度间相关性)，故数字水印检测性能并不理想.尽管Sadreazami等人[19-20]尝试引入了二元α-稳定分布、二元Cauchy分布等多元Cauchy统计建模理论，通过初步考虑变换系数尺度间相关性，一定程度上提升了多尺度边缘分布模型的建模能力，但数字水印的检测性能并未得到明显改善.

本文以非下采样Shearlet变换(nonsubsampled Shearlet transform, NSST)与二元Weibull分布理论为基础，提出了一种基于二元Weibull统计建模的非下采样Shearlet域数字图像水印算法.由于非下采样Shearlet变换拥有鲁棒而接近最优的图像几何结构特征捕获能力，同时二元Weibull分布可精确描述NSST域图像信号的边缘分布与相关特性，故利用二元Weibull边缘分布对NSST域高熵块奇异值进行统计建模并估计出其模型参数，同时利用最大似然决策构造出二元数字水印检测器并盲提取水印信息，可获得较好的水印检测性能.

1 多元Weibull分布的NSST域统计建模

1.1 NSST系数统计特性分析

2008年Easley等人[21]提出了接近最优(即“最稀疏”)图像表示方法——剪切波(Shearlet)变换理论，其不仅可以鲁棒捕获图像几何结构特征，更具有方向敏感、重构准确、结构简单、表示稀疏等诸多优良特性.NSST是Shearlet变换的拓展表示形式[22]，其消除了采样步骤，将非下采样拉普拉斯金字塔(nonsubsampled Laplacian pyramid, NSLP)变换与一系列剪切滤波器进行结合，通过迭代操作实现.

NSST继承了Shearlet变换优点，不仅有效解决了Shearlet变换所存在的“混沌”现象，使同一方向信息不会在不同方向子带中重复出现，而且为二维图像提供了最优的渐近表示，对图像重要特征具有极强的捕获能力.图1给出了标准灰度图像Barbara的2级NSST示意图.这里，尺度2和尺度1均采用了4个分解方向.

为全面客观了解NSST系数统计特性，图2给出了标准灰度图像Barbara的部分NSST域高频子带系数分布直方图.从图2中可以看出，绝大多数NSST高频子带系数接近于0，只有少量NSST高频子带系数幅值较大，即体现了明显的高尖峰重拖尾分布特点.同时，4个NSST域高频子带的系数峰度值分别为23.754 4，22.901 7，22.931 0，12.126 7，远大于标准高斯分布的峰度值3，其表明NSST域高频子带系数分布呈现出了明显的非高斯性.

近年来，研究者普遍关注并利用了变换域系数的高尖峰重拖尾非高斯性，但不同程度忽略了变换系数之间的强依赖性(多种相关性).事实上，同一子带内的系数相关性、不同方向间的系数相关性及不同尺度间的系数相关性是NSST域高频子带客观存在的3种重要相关性.这里，我们采用chi-plot工具[23]分析了NSST域高频子带的系数相关特性(包括子带内、方向间、尺度间相关性).chi-plot是Scatterplot的扩展，它能更好地并且以可视化方式评估二维变量之间的相关程度.在chi-plot中，n对点(Xi,Yi)被变换成n对(λn i,χn i)(i=1,2,…,n).χn i表示在采样点(Xi,Yi)处二维分布函数能被分解为2个边缘分布的失败情况.λn i表示(Xi,Yi)到二维中值的距离，如果X和Y相关，则λn i的值往往比较集中.chi-plot图在二维坐标中的范围是[-1,1]×[-1,1].(λn i,χn i)偏离水平线χ=0的程度反映了X和Y变量的相关程度.图3给出了NSST域高频子带的子带内、方向间、尺度间chi-plot图.

从图3可以看出，3幅chi-plot图中的大多数点比较集中且较大地偏离水平线χ=0，其说明同一子带内、不同方向间、不同尺度间的NSST域高频子带系数均存在明显的相关特性.需要指出的是，计算χn i值时会因为采样的不稳定性而产生偏差，在chi-plot图中用灰色部分予以表示.

1.2 基于二元Weibull分布的NSST域统计建模

为充分利用NSST域的尺度间相关性，本文利用二元Weibull分布对NSST域尺度间父子系数进行建模，二元Weibull分布的PDF具体表示参见附录A.图4(b)给出了二元Weibull的父子联合概率分布P(X1,X2)拟合图.为了能够说明二元Weibull统计模型的精确性，这里又利用二元Cauchy分布、二元BKF分布、二元标准正态分布对NSST域尺度间父子系数进行了建模，如图4(c)(d)(e)所示.

图4(a)给出了NSST域尺度间父子系数联合分布直方图.不难看出，相比其他二元分布，利用二元Weibull分布进行统计建模，可以更加精确地拟合出NSST域尺度间父子系数的分布.

为进一步比较4种二元统计模型的拟合精度，本文选择KL距离(Kullback-Leibler divergence)来评判各种分布模型对NSST系数建模的拟合效果，KL距离计算:

其中，Pi表示联合分布直方图的bin值，Qi表示在每个bin中由二元分布PDF求出的概率值.在计算中，要先将Pi和Qi进行归一化处理.以标准灰度图像Lena，Barbara，Mandrill为例，首先获得NSST域第1尺度4个方向和第2尺度4个方向的高频子带父子系数，然后利用不同二元分布对每组父子系数进行建模，并计算出每种模型的平均KL距离，表1给出了测试结果.

Note: The black bold font provides a better performance than other commonly used distributions.

从表1可以看出，相比其他3种二元分布，二元Weibull分布可以更好地对图像NSST域尺度间父子系数相关性进行建模，因此本文利用二元Weibull统计模型构造水印检测器.

2 数字水印嵌入

本文结合人眼视觉感知特性选取非线性嵌入强度函数，采用乘性策略将水印信息嵌入到鲁棒的NSST域高频系数块奇异值内.假设I={f(x,y),0≤x<P,0≤y<Q)}表示原始载体图像,W={w(i)∈{0,1},0≤i≤N}表示随机生成的二值水印序列.则整个图像水印的嵌入过程可描述如下：

1) 原始载体的非下采样Shearlet变换

首先对原始载体图像进行2级非下采样Shearlet变换，以获得1个低频子带及第1尺度4个方向、第2尺度4个方向的高频子带.选取第1尺度和第2尺度中具有最高能量的方向子带(以下简称“重要子带”)用于水印嵌入，子带能量的计算方法:

其中，Dk,j表示位于第k尺度、第j个方向的子带.

2) 重要Shearlet系数块选取及其奇异值分解

首先，将所选取出的“重要子带”划分成不重叠的Shearlet系数块(每块大小为L×L)，并计算出每个Shearlet系数块的熵值.然后，根据熵值，对第1与第2尺度“重要子带”的Shearlet系数块进行各自排序处理，进而分别选取前N个高熵Shearlet系数块(以下简称“重要Shearlet系数块”)用于嵌入水印.由于高熵系数块集中了大量信息，人眼对其不敏感，故将水印信息嵌入到高熵Shearlet系数块中，可以增强水印不可感知性.最后，对所选取的N个重要Shearlet系数块分别进行奇异值分解(singular value decomposition, SVD)，并选取每个重要Shearlet系数块的第2个奇异值作为水印嵌入位置.由于奇异值的稳定性较好，故通过对奇异值的微小改变来嵌入水印信息，可以提高数字水印的不可感知性及鲁棒性.

3) 数字水印嵌入

现有统计模型图像水印方案普遍采用线性单调函数刻画水印嵌入强度，并用于自适应调节载体信号[11,13-15,18-19].然而，水印嵌入强度增量与载体信号改变量间并非具有简单的线性单调关系，而是呈现出明显的非线性单调关系.为此，本文以刻画能力更强的非线性单调函数为基础，构造了自适应高阶水印嵌入强度(函数)，并进一步给出了自适应图像内容的高性能乘性水印嵌入方法.

设Bi(i=0,1,…,N-1)表示原始重要Shearlet系数块，

表示含水印重要Shearlet系数块，式(3)给出了自适应图像内容的数字水印嵌入方法.

其中，xi∈Bi和

分别表示原始重要Shearlet系数块奇异值和含水印重要Shearlet系数块奇异值.f1(xi)和f0(xi)分别表示待嵌入水印位为“1”时的嵌入强度函数和待嵌入水印位为“0”时的嵌入强度函数.

这里，a1=100,b1=1.2,a2=150,b2=1.5.图5给出了水印位为“1”和“0”时的嵌入强度函数.

与现有统计模型水印方案普遍采用线性水印嵌入强度不同，本文采用刻画能力更强的非线性单调函数——反三角函数，构造了自适应水印嵌入强度函数.从图5可以看出，本文所选取的非线性单调函数具有使嵌入强度增量与载体信号改变量保持非线性关系的性质，即：当嵌入强度增加时，使得值大的重要Shearlet系数块奇异值发生大的改变，而值小的重要Shearlet系数块奇异值发生小的改变，反之亦然.

重复步骤3)，直到所有重要Shearlet系数块处理完毕为止.最后，对每个含水印奇异值进行逆奇异值分解即可得到含水印重要Shearlet系数块，进一步对含水印Shearlet域重要子带与其他子带一起进行逆非下采样Shearlet变换，即可得到含水印图像.

3 数字水印检测

3.1 二元水印检测器构造

假设水印检测器端接收到的含水印重要Shearlet系数块奇异值yi受到零均值的高斯白噪声

其中，f表示2种水印嵌入强度函数.由于重要Shearlet系数块奇异值与噪声项是相互独立的，因此可获得含噪声重要Shearlet系数块奇异值的概率密度函数公式：

为简化式(6)，本文利用three-sigma准则估计高斯分布的概率密度函数fn(τ)：

由式(7)可见，存在2个非0的闭合区间：[-3σn,0]和[0,3σn]，因此式(6)可表示成式(7)的形式：

根据Simpson’s rule，将式(7)带入式(8)中，化简后可以得到：

其中，fBWD表示二元Weibull分布概率密度函数[24].

本文中，假设二值水印序列与重要Shearlet系数块奇异值均为服从独立同分布的随机变量，于是有重要Shearlet系数块奇异值所组成的y11,y12,…,y1N(来自第1尺度重要子带)和y21,y22,…,y2N(来自第2尺度重要子带)，在嵌入水印位为“1”和嵌入水印位为“0”两种假设下的检测概率.

其中，yi=[y1i,y2i],i=1,2,…,N,fBWD(yi|1)和fBWD(yi|0)分别是在使用嵌入强度函数的情况下被计算出来的：

其中，g1(·)和g0(·)分别代表嵌入强度函数f1(xi)和f0(xi)的反函数，ρ表示y1i和y2i之间的相关系数(i=1,2,…,N),C1和k1表示y1i的形状参数和尺度参数；C2和k2表示y2i的形状参数和尺度参数.这里我们采用最大似然方法对4个参数C1,k1,C2,k2进行估计，并进一步可通过最大似然决策规则获得最优水印检测器.

代入式(10)和式(11)，经过代数操作后，Zi(y)可以表示为

于是，可以从第i个重要Shearlet系数块中提取出水印信息位

其中，

和fBWD(yi|0)由式(12)(13)求得.

3.2 数字水印提取

首先对含水印图像进行2级非下采样Shearlet变换，并按照第2节步骤1)和步骤2)方法确定出“重要Shearlet系数块”，同时对其进行奇异值分解.然后，利用3.1节所构造的数字水印检测器，从每个重要Shearlet系数块奇异值所组成的系数对yi=[y1i,y2i]中提取出水印信息位.

3.3 误码概率分析

误码概率(bit error probability, BEP)用于分析提出的水印检测器的工作性能.在无攻击检测条件下BEP的计算：

在无噪检测条件下，含水印重要Shearlet系数块奇异值都嵌入“0”的H0假设下的检测概率Zi(y):

由于Zi(y|H0)是大量的独立随机变量之和，因此，Zi(y|H0)近似地服从正态分布，且在2种假设下，具有有限的数学期望和方差：(μ0,σ0)和(μ1,σ1).具体地，在H0假设下，其期望μ0的计算:

在H0假设下，其期望

的计算：

由于Zi(y|H0)=-Zi(y|H1)，可以得到μ1=-μ0和

故可求得检测1位水印信息位的误码概率

其中，

因此，如果嵌入到原始载体图像的二进制序列“0”或“1”具有相同的概率，总的BEP计算：

应用式(23)，即可对所提出的数字水印检测器进行性能评估.

4 实验结果

为验证本文图像水印算法的高效性，分别给出了本文算法的数字水印检测性能、不可感知性和鲁棒性的测试结果，并与文献[12,15-16,19]算法进行了对比.文献[12,15-16,19]均为基于统计模型的变换域乘性水印算法，且与本文算法密切相关.其中，文献[15,19]采用了线性单调函数刻画水印嵌入强度，文献[12,16]采用了非线性单调函数(指数函数)刻画水印嵌入强度；此外，文献[12,15-16]利用了单元概率密度函数建立多尺度边缘分布模型，文献[19]则是利用了多元概率密度函数建立多尺度边缘分布模型.实验中，原始载体图像选用了512×512×8 b的标准灰度图像Lena，Mandrill，Barbara，Peppers，数字水印则采用了128 b，256 b，512 b，1 024 b的伪随机序列.此外，我们利用峰值信噪比(peak signal to noise ratio, PSNR)和失真率(bit error ratio, BER)评价了图像水印算法的不可感知性和鲁棒性.同时，非下采样Shearlet变换的分解级数K=2，第1尺度和第2尺度均采用了4个分解方向，Shearlet系数块大小为8×8.

4.1 水印检测器的性能测试

本文结合非下采样Shearlet域二元Weibull边缘分布模型与最大似然决策理论，构造了二元数字水印检测器并盲提取了数字水印信息.本部分以标准灰度图像Barbara为例，对所构造的数字水印检测器的工作性能予以测试.图6分别给出了在嵌入水印位为“1”和“0”2种情况下数字水印检测器的检测响应.其中，横坐标为接收端含水印图像信号序列，纵坐标表示水印检测器的响应值，虚线表示检测器关于含水印图像信号的检测响应值，实线表示关于含水印图像信号的检测阈值.

通过图6可以看出，本文所提出的二元数字图像水印检测器能够准确地提取出数字水印信息(包括各种攻击下)，即当嵌入水印位为“1”时，有SML>ST；当嵌入水印位为“0”时，有ST>SML.这里，S表示曲线(直线)与x轴和y轴围成的面积.根据数字水印检测原理，当嵌入水印位为“1”时，检测器(ML detector， ML)的曲线面积SML越大于检测阈值(threshold， T)的曲线面积ST，则说明检测响应效果越好，即检测器工作能力越强.同理，当嵌入水印位为“0”时，检测器的曲线面积越小于检测阈值的曲线面积，则说明检测响应效果越好，即检测器工作能力越强.

表2给出了不同水印容量(包括128 b，256 b，512 b，1 024 b)情况下，整个图像水印系统的PSNR、BER、水印嵌入时间以及水印提取时间等工作性能.

4.2 不可感知性测试

理论分析与实验结果均表明:水印嵌入强度对算法工作性能有一定影响，为了取得不可感知性和鲁棒性的良好平衡，水印嵌入强度应自适应于载体信号内容，即水印嵌入强度应该随着载体信号强度(值)的增大而增大.现有算法普遍采用线性单调函数刻画水印嵌入强度，并用于自适应调节载体信号(即嵌入水印信息)[15,19].然而，水印嵌入强度增量与载体信号改变量间并非具有简单的线性单调关系，而是呈现出明显的非线性单调关系，为此本文算法采用刻画能力更强的非线性单调函数——反三角函数，构造了自适应水印嵌入强度函数.与指数函数相比[12,16]，反三角函数具有更强刻画水印嵌入强度增量与载体信号改变量间非线性关系的性质，更好地获得了不可感知性和鲁棒性的良好平衡.

图7给出了部分标准灰度图像的不可感知性主观测试结果，包括原始载体图像、含水印图像及放大10倍的差值图像.不难看出，利用本文所提出算法嵌入数字水印后，人眼无法察觉到含水印图像和原图像之间的差别.

表3给出了不同统计模型图像水印算法的感知透明性客观评价结果.不难看出，本文所提出的统计模型图像水印算法具有较好的不可感知性能.

4.3 鲁棒性测试

为验证本文算法的鲁棒性，我们以标准灰度图像Barbara和Mandrill为原始载体，128 b伪随机序列作为数字水印，对含水印图像进行了一系列攻击实验.表4给出了本文算法与典型统计模型图像水印算法(包括文献[12,15-16,19]算法)的鲁棒性测试对比结果.不难看出，相对于3种典型统计模型图像水印算法，本文所提出的图像水印方案具有更好的鲁棒性能，可以有效抵抗常规攻击和几何攻击.

5 结束语

本文以非下采样Shearlet变换与二元Weibull分布理论为基础，提出了一种基于二元Weibull统计建模的非下采样Shearlet域数字图像水印算法.由于非下采样Shearlet变换拥有鲁棒而接近最优的图像几何结构特征捕获能力，同时二元Weibull分布可精确描述非下采样Shearlet域图像信号的边缘分布与相关特性，故使得本文图像水印算法可较好获得不可感知性、鲁棒性、水印容量之间的良好平衡.同时，仿真实验结果也证明了本文算法的有效性.

本文算法虽然具有较强的创新性和重要的理论价值，且亦可以应用到彩色图像处理等其他领域，但仍存在些许不足之处.例如本算法仅初步考虑了变换系数尺度间相关性，故多尺度边缘分布模型的建模能力未得到明显改善.在未来的工作中，我们将充分利用多尺度变换系数间的多种重要相关特性(子带内、方向间、尺度间相关性)，进一步提升多尺度边缘分布模型的建模能力.

[1]Shih F Y. Digital Watermarking and Steganography: Fundamentals and Techniques[M]. 2nd ed. Boca Raton, USA: CRC, 2017

[2]Qasim A F, Meziane F, Aspin R. Digital watermarking: Applicability for developing trust in medical imaging workflows state of the art review[J]. Computer Science Review, 2018, 27: 45-60

[3]Wang Chunpeng, Wang Xingyuan, Zhang Chuan, et al. Stereo image zero-watermarking algorithm based on ternary polar harmonic Fourier moments and chaotic mapping[J]. Scientia Sinica Informationis, 2018, 48(1): 79-99 (in Chinese)(王春鹏, 王兴元, 张川, 等. 基于三元数极谐-Fourier矩和混沌映射的立体图像零水印算法[J]. 中国科学: 信息科学, 2018, 48(1): 79-99)

[4]Panah A S, Schyndel R V, Sellis T. On the properties of non-media digital watermarking: A review of state of the art techniques[J]. IEEE Access, 2016, 4: 2670-2704

[5]Wang Xiangyang, Yang Hongying, Niu Panpan, et al. Quaternion exponent moments based robust color image watermarking[J]. Journal of Computer Research and Development, 2016, 53(3): 651-665 (in Chinese)(王向阳, 杨红颖, 牛盼盼, 等. 基于四元数指数矩的抗几何攻击彩色图像水印算法[J]. 计算机研究与发展, 2016, 53(3): 651-665)

[6]Shi Hui, Li Mingchu. Safe variable-capacity self-recovery watermarking scheme[J]. Journal of Computer Research and Development, 2014, 51(8): 1715-1726 (in Chinese)(石慧, 李明楚. 安全的变容量恢复水印方案[J]. 计算机研究与发展, 2014, 51(8): 1715-1726)

[7]Barni M, Bartolini F, DeRosa A, et al. Optimum decoding and detection of multiplicative watermarks[J]. IEEE Transactions on Signal Processing, 2003, 51(4): 1118-1123

[8]Do M N, Vetterli M. Rotation invariant texture characterization and retrieval using steerable wavelet-domain hidden Markov models[J]. IEEE Transactions on Multimedia, 2012, 4(4): 517-527

[9]Do M N, Vetterli M. Wavelet-based texture retrieval using generalized Gaussian density and Kullback-Leibler distance[J]. IEEE Transactions on Image Processing, 2002, 11(2): 146-158

[10]Wang Chuntao, Ni Jiangqun, Huang Jiwu. An informed watermarking scheme using hidden Markov model in the wavelet domain[J]. IEEE Transactions on Information Forensics and Security, 2012, 7(3): 853-867

[11]Amini M, Ahmad M O, Swamy M N S. A new locally optimum watermark detection using vector-based hidden Markov model in wavelet domain[J]. Signal Processing, 2017, 137: 213-222

[12]Akhaee M A, Sahraeian S M E, Marvasti F. Contourlet-based image watermarking using optimum detector in a noisy environment[J]. IEEE Transactions on Image Processing, 2010, 19(4): 967-980

[13]Bian Yong, Liang S. Locally optimal detection of image watermarks in the wavelet domain using Bessel K form distribution[J]. IEEE Transactions on Image Processing, 2013, 22(6): 2372-2384

[14]Etemad S, Amirmazlaghani M. A new multiplicative watermark detector in the contourlet domain using t location-scale distribution[J]. Pattern Recognition, 2018, 77: 99-112

[15]Sadreazami H, Ahmad M O, Swamy M N S. Multiplicative watermark decoder in contourlet domain using the normal inverse Gaussian distribution[J]. IEEE Transactions on Multimedia, 2016, 18(2): 196-207

[16]Wang Xiangyang, Liu Yunan, Xu Huan, et al. Blind optimum detector for robust image watermarking in nonsubsampled shearlet domain[J]. Information Sciences, 2016, 372: 634-654

[17]Zhang Yue, Shu Huazhong, Wu Jiasong, et al. A robust watermarking algorithm based on a new orthogonal complex transform[J]. Acta Electronica Sinica, 2013, 41(8): 1574-1579 (in Chinese)(张悦, 舒华忠, 伍家松, 等. 一种基于新的正交复数变换的鲁棒水印算法[J]. 电子学报, 2013, 41(8): 1574-1579)

[18]Dong Luan, Yan Qin, Lü Yong, et al. Full band watermarking in DCT domain with Weibull model[J]. Multmedia Tools and Applications, 2017, 76(2): 1983-2000

[19]Sadreazami H, Ahmad M O, Swamy M N S. A study of multiplicative watermark detection in the contourlet domain using alpha-stable distributions[J]. IEEE Transactions on Image Processing, 2014, 23(10): 4348-4360

[20]Sadreazami H, Ahmad M O, Swamy M N S. A robust multiplicative watermark detector for color images in sparse domain[J]. IEEE Transactions on Circuits and System, 2015, 62(12): 1159-1163

[21]Easley G, Labate D, Lim W Q. Sparse directional image representations using the discrete shearlet transform[J]. Applied and Computational Harmonic Analysis, 2008, 25(1): 25-46

[22]Lim W Q. A new directional transform and compactly supported shearlet frames[J]. IEEE Transactions on Image Processing, 2010, 19(5): 1166-1180

[23]Fisher N I, Switzer P. Chi-plots for assessing dependence[J]. Biometrica, 1985, 72(2): 253-265

[24]Villanueva D, Feijóo A, Pazos J L. Multivariate Weibull distribution for wind speed and wind power behavior assessment[J]. Resources, 2013, 2(3): 370-384