CN108596823B - 一种基于稀疏变换的数字盲水印的嵌入和提取方法 - Google Patents

一种基于稀疏变换的数字盲水印的嵌入和提取方法 Download PDF

Info

Publication number
CN108596823B
CN108596823B CN201810403777.0A CN201810403777A CN108596823B CN 108596823 B CN108596823 B CN 108596823B CN 201810403777 A CN201810403777 A CN 201810403777A CN 108596823 B CN108596823 B CN 108596823B
Authority
CN
China
Prior art keywords
sparse
watermark
image
transform
embedded
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810403777.0A
Other languages
English (en)
Other versions
CN108596823A (zh
Inventor
顾宇鑫
马小虎
边丽娜
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Suzhou University
Original Assignee
Suzhou University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Suzhou University filed Critical Suzhou University
Priority to CN201810403777.0A priority Critical patent/CN108596823B/zh
Publication of CN108596823A publication Critical patent/CN108596823A/zh
Application granted granted Critical
Publication of CN108596823B publication Critical patent/CN108596823B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T1/00—General purpose image data processing
    • G06T1/0021—Image watermarking
    • G06T1/005—Robust watermarking, e.g. average attack or collusion attack resistant
    • G06T1/0057—Compression invariant watermarking
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/214—Generating training patterns; Bootstrap methods, e.g. bagging or boosting
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F21/00—Security arrangements for protecting computers, components thereof, programs or data against unauthorised activity
    • G06F21/10—Protecting distributed programs or content, e.g. vending or licensing of copyrighted material ; Digital rights management [DRM]
    • G06F21/16—Program or content traceability, e.g. by watermarking

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Computation (AREA)
  • Evolutionary Biology (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Artificial Intelligence (AREA)
  • Multimedia (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Technology Law (AREA)
  • Computer Hardware Design (AREA)
  • Computer Security & Cryptography (AREA)
  • Editing Of Facsimile Originals (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种基于稀疏变换的数字盲水印的嵌入和提取方法，对原始图像进行8×8不重叠分块，用分割后的图像块结合稀疏变换模型，用稀疏的频域空间压缩嵌入范围，之后基于Alpha分布对各个频域分量系数进行估计，选择分布相近的频域分量进行水印嵌入。水印提取时，不需要原始图像，实现了盲提取。实验结果表明，本发明的方法能够有效对抗JPEG压缩、噪声、剪切等攻击，具有较好的不可见性和鲁棒性。

Description

一种基于稀疏变换的数字盲水印的嵌入和提取方法
技术领域
本发明涉及一种数字水印处理方法，具体涉及一种基于稀疏变换和Alpha稳定分布的数字水印嵌入和提取方法。
背景技术
随着信息技术和网络技术的飞速发展，数字多媒体信息的存储、复制与传播变得十分方便，对数字音乐和数字图像的编辑、修改、复制和散布涉及到数字媒体原创者的版权保护以及数字媒体的信息安全问题。盗版问题和版权纠纷问题已成为日益严重的社会问题。
数字水印技术是近年来兴起的前沿研究领域，在多媒体信息的版权保护和完整性认证方面得到迅猛发展。将含有特定版权信息的数字水印嵌入到音乐、图像或视频中，可以用于识别多媒体信息的发布者是否获得授权，因此，数字水印技术已成为版权保护的重要工具。
以数字图像水印为例，早期的水印算法大多集中在空域水印算法，如LSB、Patchwork等，算法的透明性较好，但对于滤波、噪声、JPEG压缩等常见攻击的鲁棒性较低且嵌入容量有限，不足以达到版权保护的要求。目前主流的数字水印算法主要集中在变换域水印算法，如将原始图像进行离散傅里叶变换（discrete Fourier transform, DFT）、离散余弦变换（discrete cosine transform, DCT）、离散小波变换等，在变换后的频域空间进行水印嵌入。近年来，随着信号和图像处理技术不断发展，常规变换域的拓展域也被应用在数字水印中，如同时包含信号时域和频域特征的分数阶傅里叶变换、能对信号进行多方向多尺度分解的轮廓波（Coutourlet）。
稀疏变换模型是近年来新出现的一个稀疏表达的模型，由于其在求解稀疏编码上较低的复杂度和较好的收敛性，在诸如图像去噪和磁共振成像（MRI）领域受到了研究者的广泛关注。与其他常见的变换域如DCT和小波变换等相比，稀疏变换的稀疏程度更高，对信号主要成分的表达更好，依托稀疏变换，能寻找更鲁棒的嵌入位置，因而可以考虑应用于数字水印。然而，稀疏域下的系数由于缺失了绝大部分频域分量，导致稀疏域系数中零项较多，且出现位置不存在明显规律，因此，通常的水印嵌入方式如加性嵌入、乘性嵌入以及相关性调制等方式均不适用。如何根据稀疏的特征构建新的水印嵌入方式，是水印研究者需要解决的问题。
发明内容
本发明的发明目的是提供一种基于稀疏变换的数字盲水印的嵌入和提取方法，以解决现有技术中难以实现基于稀疏变换的水印构建的问题，同时实现水印的盲提取，提供一种更好的数字水印的解决方案。
为达到上述发明目的，本发明采用的技术方案是：一种基于稀疏变换的数字盲水印的嵌入方法，包括以下步骤：
(1) 对水印图像
Figure 100002_DEST_PATH_IMAGE001
进行Arnold置乱，得到置乱后的水印图像并转换为一维序列
Figure 158458DEST_PATH_IMAGE002
，Arnold迭代次数为密钥key1；其中，M为p×q大小的水印图像，R表示实数空间；
(2) 将原始图像分成大小为8×8的N个不重叠分块，对每个分块进行白化处理后拓展为一维向量，按列组合成样本矩阵Y∈R64×N；
(3) 按稀疏变换算法对样本矩阵进行训练，获得变换矩阵W∈R64×64和稀疏域X∈R64×N，计算变换后的残差矩阵E=WY-X，E∈R64×N；
(4) 对稀疏域X的每一分量，用S
Figure 100002_DEST_PATH_IMAGE003
S分布拟合数据，得到各分量的分布参数，选择X中部分分量构建嵌入向量u i ，选择的分量坐标存储为密钥向量key2；i为对应图像块X稀疏分量的下标，i∈{1,2,…N}；
(5) 使用密钥key3生成扩展变换抖动调制(spread transform dithermodulation，STDM)中所需要的投影向量v，同时确定量化步长
Figure 612573DEST_PATH_IMAGE004
和随机抖动信号d m ；
(6) 利用STDM嵌入器在每个8×8分块的向量u i 中嵌入水印位m i ，方法为：
Figure 100002_DEST_PATH_IMAGE005
式中，s i 为携带水印信息的向量，抖动调制器QDM中，
Figure 643371DEST_PATH_IMAGE004
为量化步长，d m 是对应于欲嵌入水印信息m i 的随机抖动信号；
(7) s i 和稀疏域中未修改的部分共同构成含水印的稀疏域
Figure 182936DEST_PATH_IMAGE006
，经
Figure 100002_DEST_PATH_IMAGE007
逆变换后，重采样成原始分辨率大小，得到嵌入水印后的图像；
Figure 705185DEST_PATH_IMAGE008
为水印嵌入后的样本矩阵。
上述技术方案中，基于稀疏变换和Alpha稳定分布实现，改进了传统变换域水印算法在嵌入分量上的选取问题，对原始图像进行8×8不重叠分块，用分割后的图像块结合稀疏变换模型，用稀疏的频域空间压缩嵌入范围，之后基于Alpha分布对各个频域分量系数进行估计，选择分布相近的频域分量进行水印嵌入。由此，水印提取时，不需要原始图像，实现了盲提取。
上述技术方案中，步骤(3)中，所述稀疏变换算法为，
输入：采样向量矩阵
Figure 100002_DEST_PATH_IMAGE009
，其中，n=64，稀疏程度s，训练参数
Figure 13806DEST_PATH_IMAGE010
、
Figure DEST_PATH_IMAGE011
，梯度下降步长
Figure 212575DEST_PATH_IMAGE012
，算法整体迭代次数t，共轭梯度下降的迭代次数g；
输出：变换矩阵
Figure 100002_DEST_PATH_IMAGE013
，稀疏域
Figure 973858DEST_PATH_IMAGE014
；
用二维离散余弦变换矩阵初始化W；初始化算法整体迭代次数k=0；
循环执行以下步骤：
Step1，判断是否满足k=t，若满足则循环终止并输出W和X，否则执行Step2；
Step2，根据X=WY计算稀疏域X, 保留X中每列系数前s个最大的值, 其余系数设为0；
Step3，初始化当前梯度下降迭代次数i=0；
Step4，判断梯度下降迭代次数是否满足i=g，若满足则令k=k+1并执行Step1，否则执行Step5；
Step5，共轭梯度下降计算梯度：
Figure 100002_DEST_PATH_IMAGE015
，其中，G为损失函数的梯度矩阵；
计算搜索方向：
Figure 299797DEST_PATH_IMAGE016
，其中，
Figure 100002_DEST_PATH_IMAGE017
为共轭梯度矩阵，F表示矩阵的F范数；
Step6，稀疏变换矩阵更新：
Figure 462925DEST_PATH_IMAGE018
，返回执行Step4。
优选地，λ=μ=4×105，
Figure 848907DEST_PATH_IMAGE012
=1×10-8，t=300，g=128，s=16。
上述技术方案中，步骤(4)中，嵌入向量的构建方法为：
Figure 877912DEST_PATH_IMAGE020
其中
Figure 7542DEST_PATH_IMAGE022
为第j个图像块样本第i个稀疏分量对应系数，C是稀疏域中非零系数个数大于样本数一半的分量集合。
Figure 290756DEST_PATH_IMAGE024
是可选参数集，包含
Figure 113218DEST_PATH_IMAGE026
分布下估计的特性指数
Figure 114672DEST_PATH_IMAGE028
和尺度参数
Figure 782414DEST_PATH_IMAGE030
。统计C中各稀疏分量在
Figure 185713DEST_PATH_IMAGE026
分布下用极大似然估计的
Figure 179077DEST_PATH_IMAGE024
。对这些参数划分成10个区间统计直方图信息，L为第k个分量中根据参数最小值和最大值均匀划分的10个区间，N为落在该区间中样本数量，选择其中众数所在的区间
Figure 651516DEST_PATH_IMAGE032
，计算均值
Figure 388528DEST_PATH_IMAGE034
，作为参数标准量。最后选择C中与参数标准量差异小于阈值e的全部分量构成嵌入分量集合G。对每个图像块选择G中分量对应的稀疏域系数构建嵌入向量，之后通过扩展变换抖动调制进行水印嵌入。
上述技术方案中，统计C中各稀疏分量在
Figure 646334DEST_PATH_IMAGE026
分布下估计的
Figure 545020DEST_PATH_IMAGE024
，对这些参数划分成10个区间统计直方图信息，选择其中众数所在的区间
Figure 521066DEST_PATH_IMAGE032
，计算均值
Figure 796189DEST_PATH_IMAGE034
，作为参数标准量。最后选择C中与参数标准量差异小于阈值e的分量构建嵌入向量。
本发明同时提供了一种基于稀疏变换的数字盲水印的提取方法，用于提取按上述方法嵌入的水印，包括以下步骤：
(a) 将含水印图像按大小为8×8进行不重叠分块；
(b) 按照嵌入水印时相同的方法训练稀疏变换矩阵W’和稀疏域X’；
(c) 根据key2对X’的每个分块重新构建嵌入主向量u i ’；
(d) 利用密钥key3计算每个主向量u i ’在STDM调制中所需要的投影向量v；
(e) 根据v及已知的
Figure 908502DEST_PATH_IMAGE004
、d m ，通过下式按序提取每个分块图像对应的嵌入向量u i ’中所携带的水印位：
Figure DEST_PATH_IMAGE035
其中，
Figure 224427DEST_PATH_IMAGE036
；
由此获得水印序列
Figure DEST_PATH_IMAGE037
；
(f) 将提取出的水印序列
Figure 156611DEST_PATH_IMAGE037
转换为
Figure 235425DEST_PATH_IMAGE038
大小的二维矩阵，再通过key1进行Arnold逆置乱获得最终提取的水印图像M’。
由于上述技术方案运用，本发明与现有技术相比具有下列优点：
1、本发明基于稀疏变换和Alpha稳定分布，实现了数字水印的嵌入和提取，能够有效对抗JPEG压缩、噪声、剪切等攻击，具有较好的不可见性和鲁棒性。
2、本发明的水印提取时，不需要原始图像，实现了盲提取。
附图说明
图1是DCT域系数和稀疏域系数对比；
图2 是图像稀疏域系数在各分布下的log-scale PDF；
图3是实施例中水印嵌入流程示意图；
图4是实施例中水印提取流程示意图；
图5是实施例中的原始载体图像；
图6是实施例中的水印图像。
具体实施方式
下面结合附图及实施例对本发明作进一步描述：
实施例一：一种基于稀疏变换的数字盲水印的嵌入和提取方法，通过对原始图像进行8×8不重叠分块，用分割后的图像块结合稀疏变换模型，用稀疏的频域空间压缩嵌入范围，之后基于Alpha分布对各个频域分量系数进行估计，选择分布相近的频域分量进行水印嵌入。水印提取时，不需要原始图像，实现了盲提取。
为便于理解，首先对稀疏变换和Alpha稳定分布模型介绍如下：
1、稀疏变换
稀疏变换模型，假定信号
Figure 467823DEST_PATH_IMAGE040
在经过变换
Figure 708312DEST_PATH_IMAGE042
后可以被稀疏表示，
Figure 658950DEST_PATH_IMAGE044
，其中
Figure 275877DEST_PATH_IMAGE046
是稀疏编码，e是表示误差。和通常的稀疏字典学习不同的是，W是一个变换矩阵，且不是正交投影变换，每行向量没有正交和归一化的约束，意味着x不会受限于W的范围，这个性质允许稀疏变换模型有更宽的信号类别，式(1)为稀疏变换模型，其中s为稀疏程度。
Figure 612049DEST_PATH_IMAGE048
(1)
加上对W矩阵的非奇异约束，和尺度约束，目标函数变为：
Figure 23439DEST_PATH_IMAGE050
(2)
与其他常见的变换域如DCT和小波变换等相比，稀疏变换的稀疏程度更高，对信号主要成分的表达更好，依托稀疏变换，能寻找更鲁棒的嵌入位置。然而稀疏域下的系数由于缺失了绝大部分频域分量，通常的水印嵌入方式如加性嵌入、乘性嵌入以及相关性调制等方式均不适用。附图1所示为Lena图左上方8×8图像块的DCT域系数及稀疏域系数，其中，左侧为DCT域系数，右侧为稀疏域系数。显然，稀疏域系数中零项较多，且出现位置不存在明显规律。
2、
Figure 195794DEST_PATH_IMAGE026
分布
学术界对自然图像DCT交流系数（AC）的统计分布进行了长久的研究。Barni 等人采用广义高斯分布对170 张图像的整幅DCT 变换交流系数的分布进行拟合，得到其统计模型非常接近拉普拉斯分布。Reininger等用KS检验验证了图像DCT交流系数服从拉普拉斯（Laplacian）分布。Sadreazami基于Alpha稳定分布对轮廓波（contourlet）的各个频域子带进行分析进而设计水印检测算法。正因为通常的变换域系数都服从一个非随机的分布，假设以DCT矩阵作为初始点求解得到的稀疏变换，其各个频域分量也服从非随机分布，就能比较便捷地对稀疏域进行水印嵌入。
经实验发现，图像块的稀疏系数具有非高斯性质和重尾分布的情况。稀疏域系数的经验分布具有较大峰值，并且尾部比高斯PDF更重，即重尾PDF。有鉴于此，本发明用Alpha稳定分布模型(
Figure 616411DEST_PATH_IMAGE026
)去估计稀疏域各个分量的系数分布。该模型适用于描述具有非高斯估计和重尾的信号。Alpha稳定分布模型需要四个参数来描述其完整性：特征指数
Figure 557822DEST_PATH_IMAGE052
(
Figure 405692DEST_PATH_IMAGE054
)；偏斜参数
Figure 65344DEST_PATH_IMAGE056
(
Figure 24073DEST_PATH_IMAGE058
)；尺度参数
Figure 334837DEST_PATH_IMAGE060
，其意义类似于高斯分布时的方差；位置参数
Figure 353609DEST_PATH_IMAGE062
，表示概率密度函数(probability density function, PDF)在X轴的偏移。其中当
Figure 234977DEST_PATH_IMAGE064
为0时的Alpha稳定分布为对称
Figure 997397DEST_PATH_IMAGE066
稳定分布(symmetric alpha-stable，S
Figure 913400DEST_PATH_IMAGE066
S)，其概率密度函数可以表示为
Figure 103073DEST_PATH_IMAGE068
(3)
特征指数
Figure 737317DEST_PATH_IMAGE066
是决定分布形状的最重要的参数，
Figure 37848DEST_PATH_IMAGE066
的值越小，分布的尾部越重，这意味着随着S
Figure 57625DEST_PATH_IMAGE066
S分布具有小特征指数的随机变量是高度脉冲的。S
Figure 418200DEST_PATH_IMAGE066
S分布除了
Figure DEST_PATH_IMAGE069
和
Figure 743002DEST_PATH_IMAGE070
分别定义柯西和高斯分布外，没有闭合表达式。 尽管S
Figure 112803DEST_PATH_IMAGE066
S分布在近原点附近表现为近似高斯密度，但其尾部衰减速率较低。
3、稀疏域分量在S
Figure 3399DEST_PATH_IMAGE066
S分布下的估计
对称的
Figure 534874DEST_PATH_IMAGE066
稳定分布已经在重尾数据的建模(例如变换域图像系数)中引起了关注。为了验证图像的稀疏变换域系数的分布情况，使用了S
Figure 878131DEST_PATH_IMAGE066
S分布。为此，当给定的测试图像按照8×8不重叠分块，训练得到
Figure 772662DEST_PATH_IMAGE072
，其中
Figure 252185DEST_PATH_IMAGE074
是白化处理后的8×8不重叠分块向量。如图1，稀疏域分量以类似DCT系数的方式排列，共64分量，将出现非零项出现次数占总体样本1/3的分量，作为待估计的分量，去估计S
Figure 954562DEST_PATH_IMAGE066
S分布下的特征指数
Figure 785115DEST_PATH_IMAGE066
。表1是用最大似然方法估计
Figure 496719DEST_PATH_IMAGE066
获得的结果，对于少数测试图像。从表中可以看出，
Figure 96327DEST_PATH_IMAGE066
的值在0.6到1.9之间变化，表明稀疏域系数的重尾特性，并且分布不是高斯分布。因此，图像的稀疏域系数的分布可以用S
Figure 235185DEST_PATH_IMAGE066
S拟合。
表1 在S
Figure 287454DEST_PATH_IMAGE066
S估计下稀疏域不同分量的特征指数
Figure 786438DEST_PATH_IMAGE066
稀疏域系数 Lena Barbara Peppers Baboon Airplane
C<sub>0,1</sub> 0.985 1.322 0.968 1.825 0.688
C<sub>1,0</sub> 1.079 1.254 1.064 1.795 0.852
C<sub>1,1</sub> 1.100 1.375 1.187 1.927 0.864
C<sub>0,2</sub> 1.058 1.329 1.074 1.909 0.812
C<sub>2,0</sub> 1.247 1.446 1.102 1.746 0.913
C<sub>1,2</sub> 1.236 1.468 1.329 1.957 0.966
C<sub>2,1</sub> 1.189 1.518 1.307 1.917 0.900
C<sub>2,2</sub> 1.218 1.569 1.547 1.937 1.113
主要对数据的经验分布（empirical distribution）以及S
Figure 506132DEST_PATH_IMAGE066
S，广义高斯分布（generalized Gaussian distribution,GGD）和拉普拉斯分布（Laplacian distribution）这几个常见分布进行估计。由于S
Figure 550311DEST_PATH_IMAGE066
S稳定PDF不具有封闭形式表达，
Figure 89877DEST_PATH_IMAGE069
时对应柯西分布(Cauchy distribution)，所以也探究柯西分布下对稀疏域系数的估计情况。图2显示了Lena和Barbara两个图像的某一个稀疏域分量下的系数的估计情况。从图中可以看出，S
Figure 877704DEST_PATH_IMAGE003
S分布及其柯西成员对比经验分布比GGD和Laplacian分布拟合的更好。其他测试图像也能得到类似的结果。此外，为了量化PDF的拟合情况，采用Kolmogorov-Smirnov距离（KSD）描述分布的拟合度。
Figure 451905DEST_PATH_IMAGE076
(4)
其中
Figure 666986DEST_PATH_IMAGE078
是数据拟合的分布，
Figure 959427DEST_PATH_IMAGE080
是样本的经验分布。表2给出了20张分辨率为512×512的图像，每张图像按8×8划分成64×64个图像块样本，共20×4096个样本。计算这些图像小块的稀疏域系数在S
Figure 534634DEST_PATH_IMAGE066
S分布、Cauchy分布及GGD下的平均KSD距离。实验验证了S
Figure 228920DEST_PATH_IMAGE066
S分布对稀疏域系数拟合更准确。
表2 图像稀疏域系数在各分布下KSD均值
Figure 614902DEST_PATH_IMAGE082
4、根据S
Figure 129060DEST_PATH_IMAGE003
S分布参数构建水印嵌入向量
根据公式(5)选择稀疏域系数构建嵌入向量
Figure 258690DEST_PATH_IMAGE084
(5)
其中
Figure 807483DEST_PATH_IMAGE022
为第j个图像块样本第i个稀疏分量对应系数，C是稀疏域中非零系数个数大于样本数一半的分量集合。
Figure 364366DEST_PATH_IMAGE024
是可选参数集，包含
Figure 365820DEST_PATH_IMAGE026
分布下估计的特性指数
Figure 548409DEST_PATH_IMAGE028
和尺度参数
Figure 951708DEST_PATH_IMAGE030
。统计C中各稀疏分量在
Figure 679493DEST_PATH_IMAGE026
分布下用极大似然估计的
Figure 168243DEST_PATH_IMAGE024
。对这些参数划分成10个区间统计直方图信息，L为第k个分量中根据参数最小值和最大值均匀划分的10个区间，N为落在该区间中样本数量，选择其中众数所在的区间
Figure 905255DEST_PATH_IMAGE032
，计算均值
Figure 163061DEST_PATH_IMAGE034
，作为参数标准量。最后选择C中与参数标准量差异小于阈值e的全部分量构成嵌入分量集合G。对每个图像块选择G中分量对应的稀疏域系数构建嵌入向量，之后通过扩展变换抖动调制进行水印嵌入。
以下介绍本发明的一种具体实施方法：
1、一种基于稀疏变换的数字盲水印的嵌入方法，参见附图3所示，包括以下步骤：
Step1. 对水印图像
Figure 61747DEST_PATH_IMAGE086
进行Arnold置乱得到置乱后的水印图像 ,将M转换为一维序列
Figure 37793DEST_PATH_IMAGE088
， Arnold迭代次数则作为密钥
Figure 570973DEST_PATH_IMAGE090
保存.
Step2. 对原始载体图像I按8×8不重叠分块，每个图像块白化处理后拓展为一维向量，按列组合成样本矩阵
Figure 683286DEST_PATH_IMAGE092
。
Step3.按照上述稀疏变换算法训练变换
Figure 18452DEST_PATH_IMAGE094
和稀疏域
Figure 216215DEST_PATH_IMAGE096
, 计算变换后的残差矩阵
Figure 295030DEST_PATH_IMAGE098
，
Figure 261849DEST_PATH_IMAGE100
。
Step4. 对稀疏域X的每一分量，用
Figure 767917DEST_PATH_IMAGE102
分布拟合数据，得到各分量的分布参数，按稀疏编码、训练稀疏变换方法选择X中部分分量构建嵌入向量
Figure 452976DEST_PATH_IMAGE104
，选择的分量坐标存储为密钥向量
Figure 319170DEST_PATH_IMAGE106
。
Step5. 使用密钥
Figure 406074DEST_PATH_IMAGE108
生成STDM调制中所需要的投影向量v，同时确定量化步长
Figure 83043DEST_PATH_IMAGE004
和
Figure 255399DEST_PATH_IMAGE110
，且与检测端保持一致，以保证水印检测效果。
Step6. 利用STDM嵌入器向每个8×8分块的向量
Figure 676016DEST_PATH_IMAGE104
中嵌入水印位
Figure 617427DEST_PATH_IMAGE112
,嵌入后携带水印信息的向量
Figure 465297DEST_PATH_IMAGE114
可用下式导出：
Figure 124949DEST_PATH_IMAGE116
(8)
抖动调制器（QDM）中，
Figure 332945DEST_PATH_IMAGE004
为量化步长，
Figure 394442DEST_PATH_IMAGE110
是对应于欲嵌入水印信息
Figure 147634DEST_PATH_IMAGE112
的随机抖动信号。
Step7.
Figure 294582DEST_PATH_IMAGE114
和未修改部分构成含水印的稀疏域
Figure 57001DEST_PATH_IMAGE118
，经
Figure 973005DEST_PATH_IMAGE120
逆变换后，
Figure 162678DEST_PATH_IMAGE122
重采样成512×512的原始分辨率大小，得到水印嵌入后的图像
Figure 531342DEST_PATH_IMAGE124
。
其中，稀疏编码、训练稀疏变换方法如下：
输入. 采样向量矩阵
Figure DEST_PATH_IMAGE126
，稀疏程度s，训练参数
Figure 815562DEST_PATH_IMAGE127
，
Figure DEST_PATH_IMAGE128
，梯度下降步长
Figure 586072DEST_PATH_IMAGE012
，算法整体迭代次数t，共轭梯度下降的迭代次数g。
输出. 变换矩阵
Figure DEST_PATH_IMAGE130
，稀疏域
Figure DEST_PATH_IMAGE132
。
初始化. 用二维离散余弦变换矩阵初始化W；对载体图像进行
Figure 415487DEST_PATH_IMAGE133
不重叠分块，白化后拓展为n维向量，按列组合成样本矩阵Y；初始化算法整体迭代次数
Figure 271448DEST_PATH_IMAGE135
。
循环执行以下步骤：
Step1. 判断k是否满足
Figure 627867DEST_PATH_IMAGE137
，若满足则循环终止并输出W和X, 否则执行Step2；
Step2. 根据X=WY计算稀疏域X，保留X中每列系数前s个最大的值，其余系数设为0；
Step3. 初始化当前梯度下降迭代次数
Figure 518463DEST_PATH_IMAGE139
；
Step4. 判断梯度下降迭代次数是否满足
Figure 49939DEST_PATH_IMAGE141
，若满足则
Figure 393195DEST_PATH_IMAGE143
并执行Step1，否则执行Step5, Step6；
Step5. 共轭梯度下降计算梯度：
Figure 35529DEST_PATH_IMAGE145
(6)
计算搜索方向:
Figure 46210DEST_PATH_IMAGE147
(7)
Step6. 稀疏变换矩阵更新：
Figure 483008DEST_PATH_IMAGE149
，并执行Step4。
本发明实验参数设置为:
Figure 313561DEST_PATH_IMAGE151
。
2、水印提取算法
本发明的水印算法是盲数字水印, 即水印提取时不需要原始载体图像。
水印提取的具体过程如下：
Step1. 将含水印图像
Figure 25165DEST_PATH_IMAGE153
按8×8不重叠分块，按照嵌入时的方式训练稀疏变换矩阵
Figure 874041DEST_PATH_IMAGE155
和稀疏域
Figure 747319DEST_PATH_IMAGE157
。
Step2. 按保存下来的
Figure 65168DEST_PATH_IMAGE159
对
Figure 314884DEST_PATH_IMAGE157
的每个分块重新构建嵌入主向量
Figure 768999DEST_PATH_IMAGE161
。
Step3. 利用密钥
Figure 813178DEST_PATH_IMAGE163
生成每个主向量
Figure 618323DEST_PATH_IMAGE161
在STDM调制中所需要的投影向量v。
Step4. 根据v及确定的
Figure 406151DEST_PATH_IMAGE004
、
Figure 964040DEST_PATH_IMAGE165
，通过式(9)按序提取每个分块图像对应的嵌入向量
Figure 179120DEST_PATH_IMAGE161
，估计出中每个系数中所携带的水印位。
Figure 471561DEST_PATH_IMAGE167
(9)
其中，
Figure 63080DEST_PATH_IMAGE169
。
Step5. 将提取出的水印序列
Figure 491787DEST_PATH_IMAGE171
转换为
Figure 877769DEST_PATH_IMAGE173
大小的二维矩阵, 再通过
Figure 391927DEST_PATH_IMAGE175
进行Arnold逆置乱获得最终提取的水印图像
Figure 787136DEST_PATH_IMAGE177
。
3、对本实施例的方法进行实验验证，实验结果分析如下：
（1）实验载体
本发明实验平台为Matlab2015a，原始载体图像为512×512的Lena灰度图，如图5所示；大小32×32、含有“苏州大学”logo的二值图像作为水印图像，如图6所示。
（2）图像评价指标
Figure 585197DEST_PATH_IMAGE179
(10)
Figure 142080DEST_PATH_IMAGE181
(11)
Figure 143534DEST_PATH_IMAGE183
(12)
其中I和
Figure 76855DEST_PATH_IMAGE185
分别为原始图像和含水印图像，W和
Figure 480155DEST_PATH_IMAGE187
分别为原始水印图像和经过攻击后提取出的水印图像. 峰值信噪比(peak signal to noise ratio, PSNR, 单位dB) 用于检测含水印图像与原始载体图像的差异度，评测水印算法的透明性. 相关系数(normalcorrelation, NC) 和误码率(bit error ratio, BER) 用于衡量水印算法对各种信号攻击的鲁棒性. NC越大, BER越低, 则说明提取出的水印与原始水印相似度越高, 算法鲁棒性越好。
（3）图像攻击类型
实验时将含水印载体图像归一化,以便于攻击处理.使用的常规信号攻击参数说明:
(a)JPEG10：JPEG压缩质量为10；JPEG30, 即JPEG压缩质量为30；
(b) 3×3高斯滤波：标准差
Figure 207939DEST_PATH_IMAGE189
, 模板尺寸为3×3；5×5高斯滤波的标准差
Figure 696689DEST_PATH_IMAGE189
, 模板尺寸5×5；
(c) 3×3中值滤波：模板尺寸3×3；
(d)0.3%高斯噪声：均值
Figure 168122DEST_PATH_IMAGE191
，方差
Figure 672266DEST_PATH_IMAGE193
；0.5%高斯噪声即
Figure 836531DEST_PATH_IMAGE191
，
Figure 812577DEST_PATH_IMAGE195
；
(e)1%椒盐噪声：1%图像像素随机受椒盐噪声影响；
(f)1%斑点噪声(乘积性噪声)：
Figure 822122DEST_PATH_IMAGE197
，其中
Figure DEST_PATH_IMAGE198
为原图像素，
Figure 934434DEST_PATH_IMAGE199
为受攻击后的像素，
Figure DEST_PATH_IMAGE200
为均值为0, 方差为0.01的随机分布；
(g)剪切1/4: 将原图右下角1/4以像素值0代替；
(h)缩放1/2: 通过将512×512分辨率的图像采样成半分辨率256×256图像后,通过线性插值的方式恢复为原分辨率512×512。
（4）透明性分析
数字水印的透明性，决定了水印嵌入不会对载体图像进行大幅修改。因此，从视觉上几乎无法感知图像中嵌入的水印信息。为了量化嵌入水印前与嵌入水印后载体图像的视觉差异，实验以PSNR值作为衡量指标。表3所示为不同载体图像的嵌入水印后的PSNR值，可以看到本文算法在不同载体图像嵌入水印后的PSNR值均达到39dB以上，算法具有良好的透明性。
表3 不同载体图像嵌入水印后的PSNR值
载体图像 Lena Barbara Peppers Airplane
PSNR/dB 40.30 39.56 40.20 40.28
（5）鲁棒性分析
为了检测算法的鲁棒性，分别对含水印的图像分别进行了以下四类图像处理攻击：(1)格式压缩攻击：JPEG压缩；(2)滤波攻击：高斯低通滤波和中值滤波；(3)噪声攻击：高斯噪声、椒盐噪声和斑点噪声；(4)几何攻击：剪切和缩放。实验中，各含水印载体图像在图像攻击后，提取出水印的NC和BER如图表3所示。
由表4中鲁棒性实验结果可知：本实施例中的算法对JPEG压缩攻击、剪切攻击、缩放攻击的鲁棒性较好，NC值均达到了0.95以上，对应误码率BER均小于5%；对高斯滤波攻击、椒盐噪声、斑点噪声攻击的NC值也达到了0.89以上；但算法在遭受较高强度图像攻击，如5×5以上的中值滤波和0.5%高斯噪声攻击时鲁棒性稍差，无法提供较好的水印提取效果。
表4 不同载体图像嵌入水印后的鲁棒性结果
Figure DEST_PATH_IMAGE202
。
表5列出了以Lena为载体图像，“苏州大学”为水印图像，在相同图像攻击下本文算法与其他算法提取水印的鲁棒性对比，算法A为Jayalakshmi基于contourlet嵌入的水印算法；算法B为Lang J等人提出的基于分数阶傅里叶变化的水印算法；算法C为Duman O等人提出的基于分数阶傅里叶变换和小波变换相结合的数字水印算法。本文算法与其他算法对应的PSNR值分别为40.30dB、36.57dB、38.24dB、36.68dB。从表中可以看出本文算法在PSNR值高于这三个算法的前提下，对JPEG压缩、滤波攻击、剪切、缩放等常规信号攻击表现出优异的鲁棒性。JPEG压缩质量30攻击下的NC值也达到了0.99以上，提取出的水印图像与原始水印基本没有差异。对大尺度高斯滤波有着不错的鲁棒性，其NC值达到0.93以上。对去同步攻击如剪切和缩放的NC值也有0.96以上。虽然算法对噪声等随机性信号攻击的鲁棒性有所不足，但对强度一般的噪声攻击也存有0.90的NC值，仍能肉眼分辨出提取的水印信息。
表5 四种算法在相同攻击下的鲁棒性(NC)比较
Figure DEST_PATH_IMAGE204
。
本发明针对稀疏域下水印嵌入系数选择的问题，从数据统计的角度对稀疏分量的系数进行分析。用S
Figure 456551DEST_PATH_IMAGE205
S分布对各稀疏分量系数进行拟合，筛选分布参数相近的分量构建嵌入向量，每个图像分块的嵌入向量对应一个水印信息位。即使图像块系数由于稀疏的关系，导致嵌入向量中某些分量为零，但只要构建的嵌入向量不为零向量，就依旧能实现水印嵌入。在水印提取时，只需要变换到稀疏域就可以按序检测水印信息，实现了基于稀疏变换的盲数字水印算法。同时提出的算法较其他水印算法，在JPEG压缩、滤波、噪声、剪切、缩放等常规图像攻击中表现出较好的鲁棒性。

Claims (5)

1.一种基于稀疏变换的数字盲水印的嵌入方法，其特征在于，包括以下步骤：
(1) 对水印图像
Figure DEST_PATH_IMAGE001
进行Arnold置乱，得到置乱后的水印图像并转换为一维序列
Figure DEST_PATH_IMAGE002
，Arnold迭代次数为密钥key1；
(2) 将原始图像分成大小为8×8的N个不重叠分块，对每个分块进行白化处理后拓展为一维向量，按列组合成样本矩阵Y∈R64×N；
(3) 按稀疏变换算法对样本矩阵进行训练，获得变换矩阵W∈R64×64和稀疏域X∈R64×N，计算变换后的残差矩阵E=WY-X，E∈R64×N；
(4) 对稀疏域X的每一分量，用S
Figure DEST_PATH_IMAGE003
S分布拟合数据，得到各分量的分布参数，选择X中部分分量构建嵌入向量u i ，选择的分量坐标存储为密钥向量key2；
(5) 使用密钥key3生成扩展变换抖动调制(STDM)中所需要的投影向量v，同时确定量化步长
Figure DEST_PATH_IMAGE004
和随机抖动信号d m ；
(6) 利用扩展变换抖动调制(STDM)嵌入器在每个8×8分块的向量u i 中嵌入水印位m i ，方法为：
Figure DEST_PATH_IMAGE005
式中，s i 为携带水印信息的向量，抖动调制器QDM中，
Figure 879659DEST_PATH_IMAGE004
为量化步长，d m 是对应于欲嵌入水印信息m i 的随机抖动信号；
(7) s i 和稀疏域中未修改的部分共同构成含水印的稀疏域
Figure DEST_PATH_IMAGE006
，经
Figure DEST_PATH_IMAGE007
逆变换后，重采样成原始分辨率大小，得到嵌入水印后的图像。
2.根据权利要求1所述的基于稀疏变换的数字盲水印的嵌入方法，其步骤(3)中，所述稀疏变换算法为，
输入：采样向量矩阵
Figure DEST_PATH_IMAGE008
，其中，n=64，稀疏程度s，训练参数
Figure DEST_PATH_IMAGE009
、
Figure DEST_PATH_IMAGE010
，梯度下降步长
Figure DEST_PATH_IMAGE012
，算法整体迭代次数t，共轭梯度下降的迭代次数g；
输出：变换矩阵
Figure DEST_PATH_IMAGE013
，稀疏域
Figure DEST_PATH_IMAGE014
；
用二维离散余弦变换矩阵初始化W；初始化算法整体迭代次数k=0；
循环执行以下步骤：
Step1，判断是否满足k=t，若满足则循环终止并输出W和X，否则执行Step2；
Step2，根据X=WY计算稀疏域X, 保留X中每列系数前s个最大的值, 其余系数设为0；
Step3，初始化当前梯度下降迭代次数i=0；
Step4，判断梯度下降迭代次数是否满足i=g，若满足则令k=k+1并执行Step1，否则执行Step5；
Step5，共轭梯度下降计算梯度：
Figure DEST_PATH_IMAGE015
计算搜索方向：
Figure DEST_PATH_IMAGE016
Step6，稀疏变换矩阵更新：
Figure DEST_PATH_IMAGE017
，返回执行Step4。
3.根据权利要求2所述的基于稀疏变换的数字盲水印的嵌入方法，其特征在于：λ=μ=4×105，
Figure 228863DEST_PATH_IMAGE012
=1×10-8，t=300，g=128，s=16。
4.
Figure DEST_PATH_IMAGE019
根据权利要求1所述的基于稀疏变换的数字盲水印的嵌入方法，其特征在于：步骤(4)中，嵌入向量的构建方法为：
Figure DEST_PATH_IMAGE021
Figure DEST_PATH_IMAGE023
式中，i对应64个稀疏分量，j对应图像分块，C是稀疏域中非零系数个数大于样本数一半的分量集合，
Figure DEST_PATH_IMAGE024
是可选参数集，包含S
Figure 246498DEST_PATH_IMAGE003
S分布下估计的特征指数
Figure 334539DEST_PATH_IMAGE003
和尺度参数
Figure DEST_PATH_IMAGE026
，
Figure DEST_PATH_IMAGE028
为第j个图像块样本第i个稀疏分量对应系数，划分成10个区间统计直方图信息，L为第k个分量中根据参数最小值和最大值均匀划分的10个区间，N为落在该区间中样本数量，H为其中众数所在的区间，V为参数标准量，G为嵌入分量集合。
5.一种基于稀疏变换的数字盲水印的提取方法，用于提取按权利要求1-4中任一方法嵌入的水印，其特征在于，包括以下步骤：
(a) 将含水印图像按大小为8×8进行不重叠分块；
(b) 按照嵌入水印时相同的方法训练稀疏变换矩阵W’和稀疏域X’；
(c) 根据key2对X’的每个分块重新构建嵌入主向量u i ’；
(d) 利用密钥key3计算每个主向量u i ’在扩展变换抖动调制(STDM)中所需要的投影向量v；
(e) 根据v及已知的
Figure 411692DEST_PATH_IMAGE004
、d m ，通过下式按序提取每个分块图像对应的嵌入向量u i ’中所携带的水印位：
Figure DEST_PATH_IMAGE029
其中，
Figure DEST_PATH_IMAGE031
；
由此获得水印序列
Figure DEST_PATH_IMAGE032
；
(f) 将提取出的水印序列
Figure 876172DEST_PATH_IMAGE032
转换为
Figure DEST_PATH_IMAGE033
大小的二维矩阵，再通过key1进行Arnold逆置乱获得最终提取的水印图像M’。
CN201810403777.0A 2018-04-28 2018-04-28 一种基于稀疏变换的数字盲水印的嵌入和提取方法 Active CN108596823B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810403777.0A CN108596823B (zh) 2018-04-28 2018-04-28 一种基于稀疏变换的数字盲水印的嵌入和提取方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810403777.0A CN108596823B (zh) 2018-04-28 2018-04-28 一种基于稀疏变换的数字盲水印的嵌入和提取方法

Publications (2)

Publication Number Publication Date
CN108596823A CN108596823A (zh) 2018-09-28
CN108596823B true CN108596823B (zh) 2022-06-07

Family

ID=63619330

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810403777.0A Active CN108596823B (zh) 2018-04-28 2018-04-28 一种基于稀疏变换的数字盲水印的嵌入和提取方法

Country Status (1)

Country Link
CN (1) CN108596823B (zh)

Families Citing this family (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109493271B (zh) * 2018-11-16 2020-12-01 中国科学院自动化研究所 图像差值量化水印嵌入方法、提取方法、设备及存储介质
CN109714284B (zh) * 2018-11-27 2020-06-30 华中科技大学 一种基于k-s检测的射频水印检测方法
CN110009663B (zh) * 2019-04-10 2023-06-09 苏州大学 一种目标跟踪方法、装置、设备及计算机可读存储介质
CN110378832B (zh) * 2019-07-23 2023-01-24 上饶师范学院 水印嵌入方法、装置及电子设备
CN111311472B (zh) * 2020-01-15 2023-03-28 中国科学技术大学 一种图像处理模型与图像处理算法的产权保护方法
CN111340675B (zh) * 2020-02-12 2023-04-25 中南林业科技大学 一种基于稀疏表示的彩色图案水印嵌入和提取方法
CN111951148B (zh) * 2020-07-13 2022-11-11 清华大学 一种pdf文档水印生成方法和水印提取方法
CN113034335B (zh) * 2021-03-29 2024-01-05 杨帆 水印嵌入方法、系统及含水印图像的水印提取方法、系统
CN113139890B (zh) * 2021-05-18 2024-04-19 辽宁工程技术大学 一种基于子块均值的空域盲水印算法
CN114445256A (zh) * 2022-01-27 2022-05-06 腾讯科技（深圳）有限公司 一种数字水印的训练方法、装置、设备及存储介质

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102682417A (zh) * 2012-05-08 2012-09-19 常熟南师大发展研究院有限公司 一种基于压缩感知的矢量地理数据水印信息生成方法
CN103700061A (zh) * 2013-12-30 2014-04-02 东北大学 一种基于压缩感知的数字图像水印嵌入和提取方法
CN103971322A (zh) * 2014-05-13 2014-08-06 安徽师范大学 宿主图像的水印嵌入方法及含水印图像的水印提取方法
CN106709291A (zh) * 2016-12-02 2017-05-24 杭州师范大学 一种基于结构化贝叶斯压缩感知的数字水印方法
CN106772275A (zh) * 2015-12-22 2017-05-31 中国电子科技集团公司第二十研究所 一种基于Alpha稳定分布的低秩矩阵恢复检测方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102682417A (zh) * 2012-05-08 2012-09-19 常熟南师大发展研究院有限公司 一种基于压缩感知的矢量地理数据水印信息生成方法
CN103700061A (zh) * 2013-12-30 2014-04-02 东北大学 一种基于压缩感知的数字图像水印嵌入和提取方法
CN103971322A (zh) * 2014-05-13 2014-08-06 安徽师范大学 宿主图像的水印嵌入方法及含水印图像的水印提取方法
CN106772275A (zh) * 2015-12-22 2017-05-31 中国电子科技集团公司第二十研究所 一种基于Alpha稳定分布的低秩矩阵恢复检测方法
CN106709291A (zh) * 2016-12-02 2017-05-24 杭州师范大学 一种基于结构化贝叶斯压缩感知的数字水印方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
A Study of Multiplicative Watermark Detection in the Contourlet Domain Using Alpha-Stable Distributions;Hamidreza Sadreazami，etc;《IEEE Transactions on Image Processing》;20140716;第23卷(第10期);4348 - 4360 *
基于Hadamard变化和LU分解的盲水印算法;吴捷，马小虎;《计算机技术与发展》;20160229;第26卷(第2期);P82-86 *
基于压缩感知的数字图像水印技术;唐沛钦;《中国优秀博硕士学位论文全文数据库(硕士)信息科技辑 》;20180115(第01期);I138-1787 *

Also Published As

Publication number Publication date
CN108596823A (zh) 2018-09-28

Similar Documents

Publication Publication Date Title
CN108596823B (zh) 一种基于稀疏变换的数字盲水印的嵌入和提取方法
Andalibi et al. Digital image watermarking via adaptive logo texturization
Cotting et al. Robust watermarking of point-sampled geometry
Wang et al. Blind optimum detector for robust image watermarking in nonsubsampled shearlet domain
Bhatnagar et al. A new robust adjustable logo watermarking scheme
CN102903075B (zh) 基于图像特征点全局校正的鲁棒水印方法
Jia et al. A novel color image watermarking scheme based on DWT and QR decomposition
Liu et al. Robust blind image watermarking based on chaotic mixtures
Thajeel et al. A NEW COLOR IMAGE WATERMARKING TECHNIQUE USING MULTIPLE DECOMPOSITIONS.
Dharwadkar et al. Watermarking scheme for color images using wavelet transform based texture properties and secret sharing
Ou et al. A robust watermarking method for stereo-pair images based on unmatched block bitmap
Prajwalasimha et al. Digital Image Watermarking based on Successive Division
Saneie et al. Introducing a new method of robust digital image watermarking against cropping and salt & pepper noise using Sudoku
CN111340675B (zh) 一种基于稀疏表示的彩色图案水印嵌入和提取方法
Bhatnagar et al. Robust reference-watermarking scheme using wavelet packet transform and bidiagonal-singular value decomposition
Chongtham et al. A copyright protection scheme for videos based on the SIFT
Shahdoosti et al. A new digital image watermarking using optimal embedding weights
Dhar et al. Color image watermarking based on radon transform and jordan decomposition
Chu et al. Region of interest fragile watermarking for image authentication
Bas et al. Vulnerability of dm watermarking of non-iid host signals to attacks utilising the statistics of independent components
Gaata Robust watermarking scheme for GIS vector maps
Thajeel Robust Video Watermarking of Hybrid Based Techniques
Maharjan et al. Video watermarking algorithm: reducing vulnerability to geometric attacks
Mehta et al. Sub-band discrete cosine transform-based greyscale image watermarking using general regression neural network
Aherrahrou et al. The efficiency of PDE decomposition in images watermarking

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant