首页> 中国专利> 变方向差值扩展和同步嵌入的可逆水印嵌入和提取方法

变方向差值扩展和同步嵌入的可逆水印嵌入和提取方法

摘要

本发明提供一种变方向差值扩展和同步嵌入的可逆水印嵌入和提取方法,采用变方向差值扩展使得嵌入容量与嵌入数据无关且阈值选择可最大嵌入容量;为避免规范化梯度和降低分类精度,采用梯度和直接分类以提高分类精度和最大嵌入容量。使用压缩位置图来避免泛差值扩展与局部复杂性可逆水印方法存在的记录溢出像素位置大量消耗嵌入容量和附加数据过大。为避免嵌入负载数据后直接嵌入附加数据导致的不可逆,给出了同步附加备份数据的嵌入策略来保证嵌入数据可完全可逆并给出了基于排序与枚举的嵌入参数选择方法以降低计算复杂性。同基于泛差值扩展和局部复杂度可逆水印方法相比,所提方法可完全可逆且参数选择时间大为减少,具有更大的最大嵌入容量。

著录项

  • 公开/公告号CN106067157A

    专利类型发明专利

  • 公开/公告日2016-11-02

    原文格式PDF

  • 申请/专利权人 陕西师范大学;

    申请/专利号CN201610364977.0

  • 发明设计人 邵利平;陈文鑫;师军;

    申请日2016-05-27

  • 分类号G06T1/00;

  • 代理机构西安通大专利代理有限责任公司;

  • 代理人陆万寿

  • 地址 710062 陕西省西安市长安南路199号

  • 入库时间 2023-06-19 00:43:59

法律信息

  • 法律状态公告日

    法律状态信息

    法律状态

  • 2019-06-11

    授权

    授权

  • 2016-11-30

    实质审查的生效 IPC(主分类):G06T1/00 申请日:20160527

    实质审查的生效

  • 2016-11-02

    公开

    公开

说明书

技术领域

本发明属于图像信息安全和数字图像信号处理交叉领域,涉及一种可逆水印嵌入和提取方法,特别涉及一种变方向差值扩展和同步嵌入的可逆水印嵌入和提取方法。

背景技术

可逆水印是指水印提取后,嵌入载体可完整恢复的一类特殊水印。相对于传统水印,可逆水印对嵌入载体无损恢复有着严格的要求,一般用于重要图像的无失真保护,在军事图像、医学图像和遥感图像上有着重要的应用价值。

Tian等提出的差值扩展可逆水印方法是图像可逆水印方法的典型方法,该方法将邻近像素对进行Haar整数小波变换,将变换后的差值进行扩展以实现水印嵌入,其理论嵌入容量上限是0.5bpp,即2个像素最多能嵌入1比特(TIAN J.Reversible data embedding using a difference expansion[J].IEEE Transactions on Circuits and Systems for Video Technology,2003,13(8):890-896)。Alattar将Tian等的差值扩展推广至多个像素,使理论嵌入容量达到1bpp,但该方法在提高嵌入容量的同时也降低了嵌入水印后载体视觉质量(ALATTAR A M.Reversible watermark using the difference expansion of a generalized integer transform[J].IEEE Transactions on Image Processing,2004,13(8):1147-1156)。Thodi等使用MED(Median Edge Detection)预测器(Weinberger M J,Seroussi G,Sapiro G.The LOCO-I lossless image compression algorithm:principles and standardization into JPEG-LS[J].Image Processing,IEEE Transactions on,2000,9(8):1309-132)预测差值代替像素间差值来对差值扩展水印方法进行改进,在提高嵌入容量同时也有较好的嵌入视觉质量(THODI D M,RODRIGUEZ J J.Prediction-error based reversible watermarking[C]//Proceedings of IEEE International Conference on Image Processing.Singapore:IEEE,2004:1549-1552)。Thodi等给出的传统基于预测差值的可逆水印方法是:

对于任意像素x和其周围邻近的像素v1,v2,v3按MED计算其预测值如式(19)所示:

x^=max(v1,v2)v3min(v1,v2)min(v1,v2)v3max(v1,v2)v1+v2-v3otherwise---(19)

在式(19)中,v1,v2,v3和x的关系由式(20)所示:

xv1v2v3---(20)

得到后,可按式(3)计算其预测差值e,然后按式(21)进行差值扩展以嵌入1位数据b∈{0,1}或平移得到修改后预测差值e′,并按式(22)得到嵌入数据后像素x′。

e=2e+b-Te<Te+TeTe-Te<-T---(21)

式(21)中,T为预设阈值。

x=e+x^---(22)

按式(21)与式(22)进行嵌入,可能会导致像素溢出。对溢出像素,基于预测差值的可逆水印算法不做任何操作,通过在位置图上的标记来对嵌入的像素和未被改变的像素加以区分。

基于预测差值的可逆水印算法的提取过程是嵌入过程的逆过程,若对x′进行恢复,则需保证x′的邻近像素v1,v2,v3已恢复,此时可由v1,v2,v3按(19)计算得到由式(3)计算得到e′,并按式(23)提取出嵌入的数据b并按式(24)还原出预测差值e,然后按式(22)即可恢复得到原始像素x′。

b=e′&1,-2T≤e′<2T(23)

经典的基于预测差值的可逆水印算法可进行多次嵌入来提高最大嵌入容量,其最大嵌入容量可达到1bpp(bits per pixel),但每次嵌入时都存在位置图所带来的嵌入容量消耗,且在反复嵌入后的载体图像中按预测器计算预测差值的精度会逐渐偏低,在影响载体图像的视觉质量的同时,所能提供的最大嵌入容量依然十分有限。

为进一步提高基于预测差值的可逆水印算法的最大嵌入容量和嵌入后载体视觉质量,Li等(Li X,Yang B,Zeng T.Efficient reversible watermarking based on adaptive prediction-error expansion and pixel selection[J].Image Processing,IEEE Transactions on,2011,20(12):3524-3533.)对基于预测差值的可逆水印算法进行了分析,指出了对于预测差值较小的像素嵌入2位数据对视觉质量造成的影响仍可能小于在预测差值较大的像素中嵌入1位数据造成的影响,并将式(21)推广为式(25):

e=4e+b-Te<Te+3TeTe-3Te<-T---(25)

在式(25)中,式(21)中的1位嵌入数据b∈{0,1}被推广为b∈{0,1,2,3},并将图像中的像素划分为平滑像素与复杂像素两种类别,对于复杂像素,至多按式(21)嵌入1位数据,对于平滑像素则可按式(25)嵌入2位数据,同时采用式(25)直接对平滑像素嵌入2位数据也避免了对平滑数据多次嵌入所带来的预测器精度下降和多次记录位置图所带来的嵌入容量损失,因此相对于传统的基于预测差值的可逆水印算法具备更大的最大嵌入容量和在同等嵌入容量下避免了视觉质量下降。

进一步,Gui等(Gui X,Li X,Yang B.A high capacity reversible data hiding scheme based on generalized prediction-error expansion and adaptive embedding[J].Signal Processing,2014,98:370-380.)又将式(25)推广至一般形式,给出了泛差值扩展公式如式(26)所示:

在式(26)中,p为整数,表示预测差值e对应的类别,tp为类别p下的差值扩展参数,T为用于控制差值扩展区域的预设阈值,b∈{0,1,…,tp}为类别p预测差值可嵌入的数据,log2(tp+1)是对应的2进制嵌入位数。对于式(26),若tp=3,则式(26)退化为式(25);若tp=1,则退化为式(21)。

基于泛差值扩展公式,Gui等给出了结合局部复杂度和泛差值扩展的可逆水印算法。同传统的基于预测差值的可逆水印算法不同,Gui等通过每个可嵌入的像素x周围邻近的8个像素v1,v2,…,v8按式(4)构造GAP预测器来计算x的预测值并进一步增补v9,v10按式(1)进行局部复杂度C*(x)计算,像素x和周围10个像素v1,v2,…,v10之间的关系如式(2)所示,并通过C*(x)来对局部复杂度C(x)进行规范化,如式(27)所示:

在式(27)中,是所有局部复杂度C*(x)中的最大值,使用式(27)可将C*(x)规范为{0,1,…,255},C(x)越高,说明像素x的预测差值越大,越不易嵌入较多数据,因此Gui等按C(x)对像素进行分类,以决定像素可嵌入的数据大小,如式(28)所示:

p=LC(x)<sL...isi+1C(x)<si...0s1C(x)---(28)

式(28)中,L为最大类别值,si,i=0,1,…,L为规范化局部复杂度控制参数且满足0≤sL≤...≤s1≤256。

得到像素类别p后,即可对负载数据β进行嵌入,为实现盲检测,附加数据α也需嵌入到载体图像中,包括嵌入时用到的参数T、{s1,s2...,sL}、{t1,t2...,tL}、嵌入结束时位置和同时由于按式(26)与式(22)进行修改的像素可能会发生溢出,因此这些的像素的坐标和数量也被记录至α中,每个记录的坐标消耗的嵌入容量为位。

基于此,Gui等给出的基于局部复杂性与泛差值扩展的可逆水印嵌入方法:

第1步:自顶向下自左向右对p≠0的像素x按式(26)与式(22)进行修改,若发生溢出则记录溢出坐标位置,反之则按式(26)与式(22)嵌入log2(tp+1)比特数据或平移,直到嵌入完负载数据β;

第2步:根据所有溢出位置信息和嵌入时用到的参数生成附加数据α,并计算α的长度len(α),同时记录前len(α)个像素的最低有效位得到备份数据ξ;

第3步:按第1步相同的方式嵌入备份数据ξ到后续的像素中;

第4步:使用最低比特位替换的方式嵌入附加数据α到前len(α)个像素中,至此得到嵌入水印后载体图像。

除了嵌入方法以外,Gui等还给出了嵌入参数{s1,s2...,sL}和{t1,t2...,tL}的选择方法。记嵌入参数集合为Φ={T,<s1,t1>,...,<si,ti>,...,<sL,tL>},为确定Φ,Gui等首先给出了在参数集合Φ的容量估值公式和嵌入后视觉载体质量估值公式,如式(29)和式(30)所示:

EED(Φ)=Σi=1L(Σs=si+1si-1Σe=-255255gi(e)h(s,e))---(30)

在式(29)与式(30)中,h(s,e)表示载体中局部复杂度C(x)=s,预测差值为e的像素总数,式(30)中gi(e)表示的预测差值为e,像素类别为i的像素,修改后视觉质量影响估值,即像素修改后对MSE(Mean>

式(29)实际为忽略溢出情况下所有像素可能提供嵌入容量之和,式(30)为所有可能需要修改的像素对MSE影响的期望值之和。

对于式(29)和式(30),若参数集合Φ*满足对于任意参数集合Φ′(Φ′≠Φ*),EEC(Φ′)≥EEC*)与EED*)≥EED(Φ′)不同时成立,称Φ*称为理想参数集合,并通过枚举参数{<s1,t1>,...,<si,ti>,...,<sL,tL>}以确定理想参数集合后,选择容量大于len(β)的理想参数作为嵌入参数。

Gui等所提方法充分利用了各类别的像素以嵌入数据,达到了较大的最大嵌入容量,但在该方法中,存在下列问题:

①容量的期望值并不代表实际容量,由式(26),若tp+1不为2的幂值,对于任意e,仅有当时,可嵌入数据,否则只能嵌入数据,但容量确定前无法确定数据被嵌入的位置,因此对于任意可嵌入像素,其实际容量为而非log2(tp+1);

②式(26)中,预设阈值T相关于差值扩展区域和平移区域大小,由于T和最大嵌入容量不存在线性关系,因此为获取更大嵌入容量,T应该被参数选择方法确定而非预设,同时T的选择相关于所有类别像素,无法保证所有类别的像素提供的嵌入容量均达到最大;

③由式(27)和式(28),规范化局部复杂度决定影响像素类别,但按式(28)计算的局部复杂度为规范后的梯度和,不同的梯度和可能规范后相同,因此分类精度低于直接使用梯度和作为局部复杂度的情况,将影响最大容量与嵌入后载体视觉质量;

④由式(29)与式(30),确定理想参数集合的过程未考虑溢出像素的影响,难以保证选择的理想嵌入参数仍具有较大的真实容量;

⑤按记录位置的方式记录溢出像素则每个溢出的像素将消耗位容量,溢出像素较多时会对最大嵌入容量造成较大影响;

⑥为获得较大嵌入容量,L通常取3,由0≤sL≤...≤s1≤256,0<t1<...<tL<2L,为确定Φ则需枚举大量的参数集合,且需多次计算式(29)与式(30),直接暴力枚举的算法复杂度过高,需较长运行时间;

⑦负载数据β嵌入完成并得到备份数据ξ后,若嵌入备份数据ξ的过程中遇到了新的溢出像素,则这些溢出无法被记录,导致存在未备份的最低有效位,从而该方法不可逆;

⑧嵌入负载数据β后,通过提取前len(α)个像素的最低有效位生成备份数据ξ,若len(β)小于len(α),可能导致部分像素的最低有效位被备份到备份数据ξ后又用作嵌入ξ的数据,即进行了错误的最低有效位备份,使得水印方法不可逆。

发明内容

本发明的目的在于克服现有技术缺陷,提供一种变方向差值扩展和同步嵌入的可逆水印嵌入和提取方法,同基于泛差值扩展和局部复杂度可逆水印方法相比,所提方法可完全可逆且参数选择时间大为减少,具有更大的最大嵌入容量。

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

一种变方向差值扩展和同步嵌入的可逆水印嵌入方法,包括以下步骤:

第1步:计算所有存在预测差值像素的局部复杂度和预测差值;

第2步:确定变方向差值扩展和同步嵌入的可逆水印嵌入方法的嵌入参数;

第3步:对所有存在预测差值的像素进行分类,分为L+1个类别,其中L为大于0的整数;

第4步:按嵌入类别相关的嵌入参数生成记录溢出像素位置的位置图,将位置图进行压缩存储得到压缩后的位置图数据;

第5步:计算嵌入时结束位置并生成附加数据;

第6步:对附加数据与负载数据进行同步嵌入得到嵌入后载体图像。

进一步,在第1步中计算所有存在预测差值像素局部复杂度的具体方法是式(1):

C*(x)=|v1-v2|+Σi=35(|vi-vi+1|+|vi+4-vi+5|)+|v1-v5|+|v2-v6|+Σi=36|vi-vi+4|---(1)

式(1)中,C*(x)为像素x的局部复杂度,C*(x)包含13组邻近像素差值,邻近像素差值绝对值的最大值为255,故C*(x)的最大值为255×13=3155,像素x和周围像素v1,v2,…,v10的邻近关系如式(2)所示:

xv1v2v3v4v5v6v7v8v9v10---(2);

在第1步中计算所有存在预测差值像素x预测差值e的具体方法是式(3):

e=x-x^---(3)

式(3)中,是由式(2)像素x和周围像素v1,v2,…,v8按GAP预测器计算像素x预测值,具体的计算方法如式(4)所示:

式(4)中,δ为水平方向梯度与垂直方向梯度绝对值之差,为平滑上下文预测值,分别按式(5)与式(6)计算:

δ=|v1-v5|+|v3-v7|+|v4-v8|-|v1-v2|-|v3-v4|-|v4-v5|(5)

在第3步中对所有存在预测差值的像素分为L+1个类别的具体方法为式(7):

p=LC*(x)sL...isi+1<C*(x)si...0s1<C*(x)---(7)

式(7)中,局部复杂度控制参数s1,s2,…,sL,且满足0≤sL≤...≤s1≤3315,3315是C*(x)的最大值。

进一步,在第2步中变方向差值扩展的具体计算式如式(8)所示,式(8)中类别为i的像素可嵌入i位数据b=(b1b2…bq…bi)2,bq∈{0,1}表示待嵌入的第q位数据,ei,q表示按差值扩展嵌入q位数据或平移后的e,Ti且Ti≥0为对应类别下的差值扩展区域控制参数,Qi=(Qi,1Qi,2...Qi,i)2为对应类别i下i位二进制的调整方向控制参数,Qi,q则表明类别为i的像素进行第q次差值扩展时的调整方向;

ei,q=eq=02ei,q-1+(-1)Qi,q+1bqq{1,2,...,p},-TieTie+(2i-1)Ti+QiTi<ee-(2i-1)(Ti+1)+Qie<-Ti---(8)

在式(8)中,q=0,即没有数据被嵌入;当-Tp≤e≤Tp时,e被用作差值扩展,且嵌入第q位数据bq时应按Qi,q确定的调整方向对已嵌入q-1位的预测差值ei,q-1进行差值扩展;当Ti<e或e<-Ti时,对对应的预测差值进行平移。

进一步,在第2步中确定变方向差值扩展和同步嵌入的可逆水印嵌入方法嵌入参数的具体方法包括以下步骤:

第2.1步:初始化记录溢出像素的空间代价估值V=8,最大嵌入容量MC=0,嵌入参数集合Φ={};

第2.2步:确定调整方向控制参数集合{Q1,Q2,…,QL},其中Qi,i=1,…,L为第i个类别的调整方向控制参数;

第2.3步:确定局部复杂度控制参数集合{s1,s2,…,sL},其中si,i=1,…,L为第i个类别的局部复杂度控制参数;

第2.4步:确定差值扩展区域控制参数集合{T1,T2,…,TL},其中Ti,i=1,…,L为第i个类别的差值扩展区域控制参数;

第2.5步:在候选嵌入参数集合Φ′={<Q1,s1,T1>,...,<QL,sL,TL>}下,计算每个存在预测的像素能提供容量和,并生成位置图且进行压缩,然后按容量和减去压缩后位置图数据大小与附加信息的基本大小得到实际嵌入容量Cap,若Cap>MC则更新MC=Cap,Φ=Φ′;

第2.6步:令V=V-1,若V>0则转第2.3步,否则嵌入参数Φ与最大嵌入容量MC均已确定,并结束嵌入参数选择方法。

进一步,在第2.2步中确定调整方向控制参数集合{Q1,Q2,…,QL}的具体方法包括以下步骤:

第2.2.1步:初始化当前类别i=L;

第2.2.2步:初始化当前方向调整参数候选值TQ=2i-1,初始化类别i下变方向差值扩展且不发生溢出的最大像素数ME=0,方向调整参数初始化为Qi=TQ;

第2.2.3步:初始化按式(8)进行变方向差值扩展且不发生溢出的像素数EXP=0,对X中所有存在预测值的像素x计算其预测值是否同时满足式(9)对应的不溢出条件,若满足则更新EXP=EXP+1直至处理完所有存在预测值的像素,若EXP>ME则更新ME=EXP,Qi=TQ;

0x^+2ie+TQ255,0x^+2i(e-1)+TQ+1255---(9)

第2.2.4步:令TQ=TQ-1,若TQ<0则已对应类别为i的调整方向控制参数确定,并使i=i-1,否则转第2.2.3步;

第2.2.5步:若i=0,输出所有调整方向控制参数{Q1,Q2,…,QL},否则转第2.2.2步;

在第2.3步中确定局部复杂度控制参数集合{s1,s2,…,sL}的具体方法包括以下步骤:

第2.3.1步:对X中所有存在预测值的像素局部复杂度从小到大排序并去除重复值,得到序列其中A为不同局部复杂度值的数量,初始化当前类别i=L;

第2.3.2步:初始化sL+1=-1,初始化上一个已确定局部复杂度控制参数值在C*中的下标索引j=0,即

第2.3.3步:初始化si的候选值在C*中下标索引TJ=j并初始化约定enp(i,si+1,si+1,-255,255)=0,ovf(i,si+1,si+1,-255,255)=0,初始化所有当前像素类别能提供的最大嵌入容量MC=0,其中符号enp(i,sa,sb,ta,tb)和ovf(i,sa,sb,ta,tb)分别表示X中满足C*(x)∈(sa,sb],e∈[ta,tb]且按类别为i进行式(8)变方向差值扩展嵌入数据后不会溢出或可能发生溢出的像素数,并且满足当sa=sb或tb<0时,enp(i,sa,sb,ta,tb)=0和ovf(i,sa,sb,ta,tb)=0;

第2.3.4步:根据式(10)与式(11)按增量更新计算与更新TJ=TJ+1,

enp(i,si+1,CTJ+1*,-255,255)=enp(i,si+1,CTJ*-255,255)+enp(i,CTJ*,CTJ+1*,-255,255)---(10)

ovf(i,si+1,CTJ+1*,-255,255)=ovf(i,si+1,CTJ*-255,255)+ovf(i,CTJ*,CTJ+1*,-255,255)---(11)

第2.3.5步:计算当时,当前类别像素可提供的忽略Ti影响嵌入容量估值若Cap>MC则更新MC=Cap,j=TJ;

第2.3.6步:若TJ=A则设置否则转第2.3.4步;

第2.3.7步:更新i=i-1,若i>0则准备计算下一个局部复杂度控制参数,并转第2.3.3步;

第2.3.8步:若i=0,则输出所有局部复杂度控制参数集合{s1,s2,…,sL};

在第2.4步中确定差值扩展区域控制参数集合{T1,T2,…,TL}具体方法包括以下步骤:

第2.4.1步:初始化当前类别i=L;

第2.4.2步:对X中所有类别为i的像素的预测差值绝对值从小到大排序,记其中最大预测值绝对值为Me;

第2.4.3步:初始化Ti及对应的候选值TT均为-1,此时有enp(i,si+1,si,-TT,TT)=0,ovf(i,si+1,si,-TT,TT)=0,初始化所有当前类别像素能提供的最大嵌入容量MC=0;

第2.4.4步:按式(12)和式(13)增量更新enp(i,si+1,si,-TT-1,TT+1)与ovf(i,si+1,si,-TT-1,TT+1),更新TT=TT+1;

enp(i,si+1,si-TT-1,TT+1)=enp(i,si+1,si-TT-1,-TT-1)+enp(i,si+1,si-TT,TT)+enp(i,si+1,siTT+1,TT+1)---(12)

ovf(i,si+1,si-TT-1,TT+1)=ovf(i,si+1,si-TT-1,-TT-1)+ovf(i,si+1,si-TT,TT)+ovf(p,sa,sb,TT+1,TT+1)---(13)

第2.4.5步:当Ti=TT时,当前类别像素可提供的嵌入容量估值Cap=enp(i,si+1,si,-TT,TT)·i-ovf(i,si+1,si,-TT,TT)·V,若Cap>MC则更新MC=Cap,Ti=TT;

第2.4.6步:当时,结束当前类别Ti枚举,否则转第2.4.4步;

第2.4.7步:更新i=i-1,若i>0,则转第2.4.2步;

第2.4.8步:若i=0,则输出所有差值扩展区域控制参数集合{T1,T2,…,TL};

第2.5步在候选嵌入参数集合Φ′={<Q1,s1,T1>,...,<QL,sL,TL>}下,计算每个存在预测差值的像素能提供的容量和的具体方法是:记载体图像X=(Xi,j)M×N存在预测差值的像素对应的像素类别矩阵为P=(Pi,j)(M-2)×(N-3),则Cap的具体计算方法如式(14)所示:

Cap=∑Pi,j,Pi,j∈{Pi,j|Pi,j>0,Xi,j按式(8)将进行差值扩展且不发生溢出}(14)。

进一步,在第2.5步与第4步中生成记录溢出像素位置的位置图的具体方法在于包括以下步骤:

第4.1步:对载体图像X=(Xi,j)M×N存在预测差值的像素区域进行位置图L=(Li,j)(M-2)×(N-3)上的初始化设置,将所有元素初始化为Li,j=0;

第4.2步:记存在预测差值的像素区域像素Xi,j对应的像素类别为Pi,j,若Pi,j>0且Xi,j按式(8)进行变方向差值扩展发生溢出,即不同时满足式(9)所对应的条件,则置Li,j-1=1;

在第2.5步与第4步中将位置图进行压缩存储得到压缩后的位置图数据的具体方法是按JBIG压缩方法进行压缩。

进一步,在第5步中:计算嵌入时结束位置的具体方法在于包括以下步骤:

第5.1步:记嵌入时结束位置坐标为(epi,epj),记嵌入结束位置像素的位置为epk,初始化epi=0,epj=1,epk=0和待嵌入数据大小Remain=len(α)+len(β),其中α为附加数据,β为嵌入负载数据,len(α)和len(β)分别为α和β的长度;

第5.2步:按自顶向下自左向右的扫描序对存在预测差值的像素进行处理,若像素x可提供p位嵌入容量且p≥Remain,则设置(epi,epj)为当前坐标,epk=Remain并结束;

第5.3步:若像素x可提供p位嵌入容量且p<Remain,则Remain=Remain-p,并转到第5.2步处理后续像素;

在第5步中生成的附加数据具体应包含以下内容:①嵌入结束位置(epi,epj,epk),其中epi,epj需epk需位;②调整方向控制参数Q1,Q2,…,QL,其中Qp需要p位,存储Q1,Q2,…,QL总共需(1+L)·L/2位;③局部复杂度控制参数s1,s2…sL,其中sp<3315<4096=212,存储s1,s2…sL总共需12L位;④差值扩展区域控制参数T1,T2…TL,其中Tp<256=28,存储T1,T2…TL总共需8L位;⑤位置图压缩数据长度len(τ),需位和⑥位置图压缩数据,附加数据的总长度len(α)的具体计算方法如式(15)所示:

进一步,在第6步中对附加数据与负载数据进行同步嵌入得到嵌入后载体图像的具体方法在于包括以下步骤:

第6.1步:按自顶向下自左向右的扫描序,以1位为单位,按式(8)嵌入负载数据β到载体图像中;

第6.2步:将附加数据α使用最低L位有效位替换的方式嵌入到X0,0及存在预测差值的像素中,每次嵌入1位数据前,先提取原有数据并同步按式(8)嵌入到第1步完成后的预测差值中,直至α嵌入完成。

一种变方向差值扩展和同步嵌入的可逆水印提取方法,包括以下步骤:

第1步:按嵌入时相同的扫描序提取用于嵌入附加数据的像素最低L位有效位,得到不含压缩位置图数据的附加数据,并从中获得压缩位置图数据大小和继续提取出压缩位置图数据;

第2步:从附加数据中得到所有嵌入参数与嵌入结束时位置,并按位置图压缩方法对应的解压缩方法解压缩位置图得到解压缩位置图数据;

第3步:从嵌入结束时位置开始,按嵌入时逆扫描序进行备份数据提取与像素还原,每提取1位数据即按最低L位有效位替换对嵌入过附加数据的像素最低L位有效位进行恢复,直到嵌入过附加数据的像素最低L位有效位全部恢复;

第4步:按嵌入时逆扫描序对剩余像素进行数据提取与还原,得到负载数据并恢复原始载体图像。

进一步,在第2步中按位置图压缩方法对应的解压缩方法为JBIG解压缩方法;

在第2步中按位置图压缩方法对应的解压缩方法为JBIG解压缩方法;

在第3步中,按嵌入时逆扫描序进行备份数据提取与像素还原和第4步中按嵌入时逆扫描序对剩余像素进行还原的具体方法为式(16):

e=ep,0-2pTp+Qp-2p+1ep,q2pTp+Qpep,q-(2p-1)Tp-Qpep,q>2pTp+Qpep,q+(2p-1)(Tp+1)-Qpep,q<-2pTp+Qp-2p+1---(16)

式(16)中,ep,q表示类别为p的像素,按差值扩展嵌入q位数据或平移后的差值e,可按迭代的方式结合式(17)与式(18)进行嵌入数据的提取与计算:

bq=ep,q&1,p≥0(17)

式(17)和式(18)中,bq∈{0,1}表示待提取的第q位数据。

本发明同现有技术优点分析:

①现有的基于局部复杂性与泛差值扩展的可逆水印方法可充分地利用各类别像素进行数据嵌入,相比于之前的基于预测差值的可逆水印嵌入方法具有较大的最大嵌入容量。然而在该水印嵌入方法中,对于差值扩展参数不为2的幂次时,容量期望并不代表实际容量并且于具体嵌入的数据密切相关,因此存在预估容量足够而实际嵌入容量不足导致的嵌入失败。同时在泛差值扩展公式中,预设阈值和差值扩展区域的大小以及平移区域像素的平移量绑定,将影响所有类别像素,从而无法保证所有类别像素都能达到最大的嵌入容量。同以上方法不同,本发明对不同类别像素对应的扩展区域及平移区域分别确定,并且不同类别的像素对应不同的控制阈值,在此基础上给出了变方向差值扩展,通过找到每个类别像素对应的最适合调整方向使得所有类别的可嵌入像素尽可能多并找到每个类别像素最适合的调整方向来对数据进行逐位嵌入,在嵌入过程中始终按2的幂次进行扩展,因此嵌入容量与具体的负载数据无关,且在参数选择得当的情况下,使用变方向差值扩展可有效地减少溢出像素,提升嵌入容量并进一步提高嵌入后载体视觉质量,且在所给出的方法中,每个像素类别的差值扩展区域控制参数分别设置,从而避免了传统泛差值扩展方法对所有类别像素差值扩展区域统一设置导致的实际嵌入容量难以同时达到最大的设计缺陷。

②现有的基于局部复杂性与泛差值扩展的可逆水印方法为避免计算复杂度过大,将局部复杂度C*(x)规范到[0,255]之内的整数,从而降低了分类精度,将影响最大容量与嵌入后载体视觉质量。同以上方法不同,本发明对参数选择方法的枚举环节进行了改进,引入了调整方向控制参数、局部复杂度控制参数和差值扩展区域控制参数枚举方法,从而可以在更高的精度下对像素进行分类,进一步确定最大容量和提高嵌入后载体视觉质量。

③在现有的基于局部复杂性与泛差值扩展的可逆水印方法中,对于分辨率为M×N大小的载体图像,每个可能溢出的像素将消耗位容量,由于大容量嵌入时易存在大量可能发生溢出的像素,记录像素溢出位置的方式处理溢出像素将严重影响最大嵌入容量。同以上方法不同,本发明采用位置图策略来处理像素溢出,并通过JBIG压缩算法来节省存储空间,在大容量嵌入时,依然不会对最大容量造成太大影响。

④现有的基于局部复杂性与泛差值扩展的可逆水印方法,一方面在嵌入负载数据β时才进行附加数据α的溢出位置记录,但β嵌入完成后,还涉及到备份数据ξ和附加数据α的嵌入,其中附加数据α是采用最低有效位的替换方式进行嵌入,而备份数据ξ即被替换出的最低有效位,采用的是与β相同的嵌入方法,但此时α的容量已经确定,因此ξ的溢出位置标记将无法记录在α中,从而在嵌入ξ时,存在未被记录的溢出像素,导致错误的数据被提取。而另一方面,若len(β)小于len(α),则嵌入β后直接从前端像素中提取ξ的方式可能导致错误的最低有效位备份,在水印提取时,这些像素的预测差值将无法被正确地被恢复。以上2种情况都将导致所提方法不可逆。同以上方法不同,本发明通过使用位置图标记溢出位置,保证在嵌入任意数据前所有溢出位置均得以记录,同时同步嵌入附加数据与备份数据避免进行错误的最低有效位备份,所提嵌入策略完全可逆,在嵌入容量足够的情况下可保证数据成功嵌入。

⑤现有的基于局部复杂性与泛差值扩展的可逆水印方法为获得较大的嵌入容量,确定参数集合需进行大量枚举,而直接通过暴力枚举方法的复杂度较高,需较长运行时间,且在计算的过程中,未考虑溢出像素的影响。同以上方法不同,本发明引入排序方法,首先对调整方向控制参数进行优化,然后采用增量更新的机制,降低了局部复杂度参数与差值扩展区域控制参数的计算复杂度,同时通过引入溢出像素的容量代价估值使得选择的嵌入参数更为精确,因而能达到更大的嵌入容量。

附图说明

图1嵌入算法流程图

图2提取算法流程图

图3原始载体图像

图4原始载体图像对应的局部复杂度矩阵

图5原始载体图像对应的预测值矩阵

图6原始载体图像对应的预测差值矩阵

图7原始载体图像对应的分类矩阵

图8原始载体图像对应的位置图

图9负载数据嵌入后载体图像

图10完成嵌入后载体图像

图11附加数据与备份数据提取与还原后载体图像

图12所有数据提取与还原后载体图像

具体实施方式

以下结合附图具体实施例对本发明方法进行详细描述,其中图1是嵌入方法流程图,图2是提取方法流程图。

以L=3,负载数据为β=(0111)2为例,给出嵌入过程如下:

第1步:如图3所示分辨率为5×6的原始载体图像X,其存在预测差值的区域为图3中的实线所表示的区域,计算该区域中所有像素的局部复杂度与预测差值,这里以X0,1局部复杂度与预测差值为例:由式(1):C*(X0,1)=|X0,2-X0,3|+|X1,0-X1,1|+|X1,1-X1,2|+|X1,2-X1,3|+|X2,0-X2,1|+|X2,1-X2,2|+|X2,2-X2,3|+|X0,2-X1,2|+|X0,3-X1,3|+|X1,0-X2,0|+|X1,1-X2,1|+|X1,2-X2,2|+|X1,3-X2,3|=|203-223|+|182-194|+|194-222|+|222-244|+|186-228|+|228-245|+|245-241|+|203-222|+|223-244|+|182-186|+|194-228|+|222-245|+|244-241|=20+12+28+22+42+17+4+19+21+4+34+23+3=249,由式(5),X0,1对应的水平方向梯度与垂直方向梯度绝对值之差δ0,1=|X0,2-X1,2|+|X1,0-X2,0|+|X1,1-X2,1|-|X0,2-X0,3|-|X1,0-X1,1|-|X1,1-X1,2|=|203-222|+|182-186|+|194-228|-|203-223|-|182-194|-|194-222|=19+4+34-20-12-28=-3,由式(6),X0,1对应的平滑上下文预测值此时按式(4)计算X0,1的预测值由于δ0,1∈[-8,8],因此最后按式(3)计算X0,1对应的预测试值计算完成后,得到的对应的局部复杂度矩阵、预测值矩阵、预测差值矩阵分别如图4、图5、图6所示;

第2步:按排序与枚举确定变方向差值扩展和同步嵌入的可逆水印嵌入方法的嵌入参数,如步骤2.1~2.6所示:

第2.1步:初始化记录溢出像素的空间代价估值V=8,最大嵌入容量MC=0,嵌入参数集合Φ={};

第2.2步:按排序的方式确定调整方向控制参数集合{Q1,Q2,…,QL}即{Q1,Q2,Q3},以Q1的确定为例,对应如步骤2.2.1~2.2.3所示,例中当前类别i=1:

第2.2.1步:初始化当前方向调整参数候选值TQ=2i-1=1,初始化类别i下变方向差值扩展且不发生溢出的最大像素数ME=0,方向调整参数初始化为Q1=TQ=1;

第2.2.2步:初始化按式(8)进行变方向差值扩展且不发生溢出的像素数EXP=0,对X中所有存在预测值的像素x计算其预测值是否同时满足式(9)对应的不溢出条件,若满足则更新EXP=EXP+1直至处理完所有存在预测值的像素,以X0,1=185和X2,3=241为例:由结合式(9),分别计算因此X0,1满足式(9),并更新EXP=EXP+1=1,由计算得到因此X2,3不满足式(9),EXP保持不变;最后可得到EXP=8,由EXP>ME,因此更新ME=EXP=8,Q1=TQ=1;

第2.2.3步:更新TQ=TQ-1=0,并重新按2.2.2步同样的方式计算得到EXP=8,由EXP≤ME因此不更新Q1,此时已确定Q1=1;

假定完成2.2步后,所有的调整方向控制参数集合{Q1,Q2,Q3}={1,3,7}。

第2.3步:对图4所有局部复杂度排序后得到序列此时有A=9,得到C*后可按枚举的方式确定所有局部复杂度控制参数{s1,s2,s3},以s3为例,即对应的当前类别i=3,对应步骤如2.3.1~2.3.5所示:

第2.3.1步:初始化上一个已确定局部复杂度控制参数值在C*中的下标索引j=0,即

第2.3.2步:初始化si的候选值在C*中下标索引TJ=j=0并初始化按照约定有初始化所有当前像素类别能提供的最大嵌入容量MC=0;

第2.3.3步:根据式(10)与式(11)按增量更新计算与即计算局部复杂度处于(-1,187]的像素在对应类别为3,差值扩展区域为[-255,255]的情况下按式(8)变方向差值扩展嵌入数据后不会溢出或可能发生溢出的像素数并更新enp(3,-1,187,-255,255)与ovf(3,-1,187,-255,255),此时满足对应局部复杂度条件的像素仅有X0,3=223,对应的预测值与预测差值分别为由式(9)且Q3=7,计算得到即满足式(9)的不溢出条件,因此enp(3,-1,187,-255,255)=enp(3,-1,-1,-255,255)+1=1,ovf(3,-1,187,-255,255)=ovf(3,-1,-1,-255,255)=0,更新TJ=TJ+1=1;

第2.3.4步:计算当时,可提供的忽略T3影响嵌入容量估值由于Cap=3>MC=0,因此更新MC=Cap=3,j=TJ=1;

第2.3.5步:由于TJ=1≠A=9,因此重复2.3.3与2.3.4步,直至TJ=A=9,此时已确定j=7,

假定完成2.3步后,所有的局部复杂度控制参数{s1,s2,s3}={700,345,203}。

第2.4步:确定所有差值扩展区域控制参数{T1,T2,T3},以T3为例,即对应的当前类别i=3,对应步骤如2.4.1~2.4.5所示:

第2.4.1步:对X中所有类别为i=3的像素的预测差值绝对值从小到大排序,由式(7)且s3=203,仅有局部复杂度不大于s3的像素被分为类别3,即X0,2=203,X0,3=223为类别3的像素,对应的预测差值分别为e0,2=-10,e0,3=-15,按预测差值绝对值从小到大排序后得到序列{10,15},因此其中最大预测值绝对值为Me=15;

第2.4.2步:初始化T3及对应的候选值TT均为-1,此时有enp(i,si+1,si,-TT,TT)=enp(3,-1,203,1,-1)=0,ovf(i,si+1,si,-TT,TT)=ovf(3,-1,203,1,-1)=0,初始化所有当前类别像素能提供的最大嵌入容量MC=0;

第2.4.3步:按式(12)和式(13)增量更新enp(i,si+1,si,-TT-1,TT+1)=enp(3,-1,203,0,0)与ovf(i,si+1,si,-TT-1,TT+1)=ovf(3,-1,203,0,0),即计算类别为3的像素中预测差值绝对值为0的像素按式(8)变方向差值扩展嵌入数据后不会溢出或可能发生溢出的像素数并更新enp(3,-1,203,0,0)与ovf(3,-1,203,0,0),由于X0,2,X0,3对应的预测差值均不为0,因此enp(3,-1,203,0,0)=enp(3,-1,203,1,-1)=0,ovf(3,-1,203,0,0)=ovf(3,-1,203,-1,1)=0,更新TT=TT+1=0;

第2.4.4步:计算当T3=TT=0时,当前类别像素可提供的嵌入容量估值Cap=enp(i,si+1,si,-TT,TT)·i-ovf(i,si+1,si,-TT,TT)·V=enp(3,-1,203,0,0)·3-ovf(3,-1,203,0,0)·8=0,由于Cap=0≤MC=0,因此MC与T3保持不变;

第2.4.5步:当时,结束T3的枚举,否则转第2.4.3步;

假定完成2.4步后,所有的差值扩展区域控制参数{T1,T2,T3}={14,3,10}。

第2.5步:在候选嵌入参数集合Φ′={<Q1,s1,T1>,...,<QL,sL,TL>}={<1,700,14>,<3,345,3>,<7,203,10>}下,计算载体图像X存在预测差值的像素对应的像素类别矩阵为P=(Pi,j)3×3,以X0,1为例:由C*(X0,1)=249,结合式(7)有s3=203<C*(X0,1)≤s2=345,因此X0,1对应的类别P0,0=2;得到的P如图7所示。然后可按式(14)计算所有像素能提供的容量Cap,以X0,1为例:X0,1对应的类别为P0,0=2,对应的预测差值e0,1=-3,由式(8)且即e0,1∈[-T2,T2],因此X0,1将用于差值扩展,由式(9)计算可得X0,1在差值扩展嵌入数据后不会发生溢出,因此X0,1可提供的嵌入容量为P0,0=2bit;最后可得到在对应嵌入参数集合下Cap=8bit。此时可生成位置图,以X0,1为例:由于X0,1可按式(9)进行修改且不发生溢出,因此L0,0保持不变且L0,0=0;得到的位置图如图8所示,随后按JBIG压缩算法进行压缩得到压缩后位置图数据,并假设其值为(01)2,同时假设附加数据的基本大小为2bit,因此在对应嵌入参数下X的实际嵌入容量为Cap=Cap-2-2=4bit。由于实际容量Cap=4>MC=0,因此更新MC=Cap=4,Φ=Φ′={<1,700,14>,<3,345,3>,<7,203,10>};

第2.6步:令V=V-1=7,若V>0则转第2.3步,否则嵌入参数Φ与最大嵌入容量MC均已确定,并结束嵌入参数选择方法;

假定第2步完成时,最优嵌入参数集合为{<1,700,14>,<3,345,3>,<7,203,10>},对应的最大嵌入容量MC=4bit。

第3步:对所有存在预测差值的像素进行分类,由于第2步结束时的最优嵌入参数等于第2.5步例举的嵌入参数集合,因此对应的P仍如图7所示;

第4步:生成记录溢出像素位置的位置图,将位置图进行压缩存储得到压缩后的位置图数据,由于第2步结束时的最优嵌入参数等于第2.5步例举的嵌入参数集合,因此对应的L仍如图8所示,且压缩后位置图数据仍假设为(01)2

第5步:计算嵌入时结束位置并生成附加数据,分别如步骤5.1~5.3与步骤5.4所示:

第5.1步:初始化epi=0,epj=1,epk=0和待嵌入数据大小Remain=len(α)+len(β)=8bit;

第5.2步:按自顶向下自左向右的扫描序对存在预测差值的像素进行处理,由于X0,1可提供P0,0=2bit嵌入容量,但P0,1=2<Remain=8,因此直接执行下一步;

第5.3步:由于P0,0=2<Remain=8,因此更新Remain=Remain-P0,0=6,并转到5.2步对后续像素进行处理;

第5.1~-5.3步执行完成后可得到epi=2,epj=2,epk=1。

第5.4步:由附加数据具体应包含的内容①嵌入结束位置(epi,epj,epk)=(2,2,1),其中epi,epj需对应的二进制表示为epi·N+epj=2·6+2=14=(01110)2,epk需即对应的二进制表示为(01)2;②调整方向控制参数Q1,Q2,Q3=1,3,7,其中Qp需要p位,对应的二进制表示分别为(1)2,(11)2,(111)2;③局部复杂度控制参数s1,s2,s3,其中sp需要12bit,以s1=700为例,其二进制表示为(001010111100)2;④差值扩展区域控制参数T1,T2,T3,其中Tp需要8bit,以T1=14为例,对应的二进制表示为(00001110)2;⑤位置图压缩数据长度len(τ)=2,需位对应的二进制表示为(00010)2和⑥位置图压缩数据(01)2;可得到附加数据为:

附加数据的实际大小可按式(15)计算但为便于示例说明,仍采用步骤2.5步中的假定,即不包含位置图压缩数据的附加数据的基本大小为2bit,并假定该部分数据值为(10)2,即假定生成的完整的附加数据为:(1001)2

第6步:对附加数据与负载数据进行同步嵌入得到嵌入后载体图像,如步骤6.1~6.2所示:

第6.1步:按自顶向下自左向右的扫描序,以1位为单位,按式(8)嵌入负载数据β=(0111)2到载体图像中,如步骤6.1.1~6.1.2所示:

第6.1.1步:由L0,0=0,P0,0=2,-T2=-3≤e0,1=-3≤T2=3,即X0,1可按差值扩展提供P0,0=2bit嵌入容量,因此β的前2bit数据(01)2将嵌入到X0,1中,X0,1对应的预测值与预测差值分别为对应的调整方向控制参数因此嵌入第1bit数据0后,嵌入第2bit数据1后嵌入后X0,1可按式(22)计算,

第6.1.2步:由L0,1=0,P0,1=3,-T3=-10≤e0,2=-10≤T3=10,即X0,2可按差值扩展提供P0,1=3bit嵌入容量因此β中剩余的后2bit数据(11)2将嵌入到X0,2中,X0,2对应的预测值与预测差值分别为对应的调整方向控制参数因此嵌入第1bit数据1后,嵌入第2bit数据1后嵌入后X0,2可按式(22)计算,

第6.1步完成后β已完成嵌入,对应的图像如图9所示。

第6.2步:将附加数据α=(1001)2使用最低3位有效位替换的方式嵌入到X0,0及存在预测差值的像素中,每次嵌入1位数据前,先提取原有数据并同步按式(8)嵌入到第1步完成后的预测差值中,直至α嵌入完成,如步骤6.2.1~6.2.4所示:

第6.2.1步:将嵌入α中的第1bit数据1到X0,0中,由X0,0=190=(10111110)2,因此先备份X0,0的最低有效位0,并将其嵌入到步骤6.1完成后的预测差值中,此时由于e0,2本可嵌入3bit数据,但仅嵌入2bit数据,还剩余1bit嵌入容量,因此该1bit的被备份的数据0将按式(8)嵌入到e0,2中,得到嵌入后X0,2可按式(22)计算,备份数据被嵌入后,将α中的第1bit数据1对X0,0的最低有效位进行替换得到X0,0=(10111111)2=191;

第6.2.2步:将嵌入α中的第2bit数据0到X0,0中,先备份X0,0的次低有效位1,并将其嵌入到步骤6.2.1完成后的预测差值中,此时由于L0,2=0,P0,2=3,e0,3=-15<-T3=-10因此X0,3需按式(8)进行平移操作,对应的预测值与预测差值分别为e0,3=-15得到e0,3=e0,3-(23-1)(T3+1)+Q3=-85,平移后X0,3可按式(22)计算,同样,X1,1,X1,2也将进行平移,分别得到平移后X1,1=185,X1,2=213。由L1,2=0,P1,2=1,-T1=-14≤e1,3=1≤T1=14,即X1,3可按差值扩展嵌入1bit数据,因此被备份的数据1将按式(8)嵌入到X1,3中得到嵌入后X1,3=246,然后将α中的第2bit数据0对X0,0的次低有效位进行替换得到X0,0=(10111101)2=189;

第6.2.3步:将α中的第3bit数据0与第4bit数据1分别对X0,0的第3低有效位及X0,1的最低有效位进行替换,替换前先对被替换的比特位进行备份并分别嵌入到X2,1,X2,2中,得到X0,0=185,X0,1=177,X2,1=215,X2,2=250。

第6.2步完成后完成水印嵌入,嵌入后载体如图10所示。

提取算法是嵌入算法的逆过程,提取算法的具体步骤为:

第1步:按嵌入时相同的扫描序提取用于嵌入附加数据的像素最低3位有效位,得到不含压缩位置图数据的附加数据,即先提取X0,0=185=(11101001)2的最低与次低有效位得到2bit的不含压缩位置图数据的附加数据(10)2,由于该值为例中的假定值,其实际值应为:

可从中得到5bit大小的压缩位置图数据的值len(τ)=2,并继续提取len(τ)=2bit位数据,即提取X0,0的第3低有效位0及X0,1=177=(10110001)2的最低有效位1得到压缩位置图数据(01)2

第2步:从附加数据中得到所有嵌入参数与嵌入结束时位置,即{<Q1,s1,T1>,<Q2,s2,T2>,<Q3,s3,T3>}={<1,700,14>,<3,345,3>,<7,203,10>}与(epi,epj,epk)=(2,2,1),并按位置图压缩方法对应的解压缩方法解压缩位置图得到解压缩位置图数据,由于压缩位置图数据仍为假定值,但提取得到的值和嵌入时相同,因此解压有位置图仍如图8所示;

第3步:从即X2,2开始,按嵌入时逆扫描序进行备份数据提取与像素还原,每提取1位数据即按最低3位有效位替换对嵌入过附加数据的像素最低3位有效位进行恢复,直到嵌入过附加数据的像素最低3位有效位全部恢复,以X2,2=250为例:由L2,1=0可知,X2,2在按式(8)修改后未发生溢出,因此需进行提取与恢复;由式(1)、式(4)与式(3)计算得到其局部复杂度C*(X2,2)=700,预测值与预测差值e2,2=9;结合式(7)对X2,2进行分类,由s2=345<C*(X2,2)=700≤s1=700,因此X2,2对应类别P2,1=1;由式(16),-21T1+Q1-21+1=-28≤e2,2=9≤21T1+Q1=29,因此e2,2是差值扩展后的预测差值,但由于X2,2为嵌入完成时像素,因此其嵌入的数据不能直接通过P2,1=1确定,而应当由epk=1确定,即X2,2中嵌入了1bit数据,即此时可按式(16)、式(17)、式(18)与式(22)进行提取与恢复:即嵌入的数据为得到提取的1bit数据1后,对X0,1的最低有效位进行还原,得到X0,1=177;所有用于嵌入备份数据的像素均完成提取与恢复后,对应的图像如图11所示;

第4步:按嵌入时逆扫描序对剩余像素进行数据提取与还原,即从X0,2,X0,1中分别按式(16)、式(17)、式(18)与式(22)进行提取与还原得到负载数据β=(0111)2并恢复原始载体图像,如图12所示;

至此,负载数据被完整提取且载体图像完整恢复。

去获取专利,查看全文>

相似文献

  • 专利
  • 中文文献
  • 外文文献
获取专利

客服邮箱:kefu@zhangqiaokeyan.com

京公网安备:11010802029741号 ICP备案号:京ICP备15016152号-6 六维联合信息科技 (北京) 有限公司©版权所有
  • 客服微信

  • 服务号