首页> 中国专利> 基于HSV颜色协方差特征的目标跟踪方法

基于HSV颜色协方差特征的目标跟踪方法

摘要

本发明公开了一种基于HSV颜色协方差特征的目标跟踪方法,主要解决现有跟踪技术针对彩色图像目标提取的特征存在信息冗余大和特征间独立性差的问题。其实现步骤是:在粒子滤波框架下首先通过目标状态预测得到候选目标,提取候选目标图像的色相、饱和度、明度、拉普拉斯算子响应作为表观特征,融合构建协方差算子;其次通过计算候选目标特征与模板的相似度权值,调整跟踪窗,并更新特征模版;最后根据权值融合候选目标,更新目标状态,实现对目标的有效跟踪。本发明能够有效减少目标特征的信息冗余,增强特征之间的独立性,提高目标特征描述的精确性和鲁棒性,可实现对彩色图像目标的实时精确跟踪。

著录项

  • 公开/公告号CN103985141A

    专利类型发明专利

  • 公开/公告日2014-08-13

    原文格式PDF

  • 申请/专利权人 西安电子科技大学;

    申请/专利号CN201410231201.2

  • 发明设计人 姬红兵;樊振华;刘月;王磊;王磊;

    申请日2014-05-28

  • 分类号G06T7/20;G06T7/40;

  • 代理机构陕西电子工业专利中心;

  • 代理人王品华

  • 地址 710071 陕西省西安市太白南路2号

  • 入库时间 2023-12-17 00:35:36

法律信息

  • 法律状态公告日

    法律状态信息

    法律状态

  • 2022-05-17

    未缴年费专利权终止 IPC(主分类):G06T 7/20 专利号:ZL2014102312012 申请日:20140528 授权公告日:20161109

    专利权的终止

  • 2016-11-09

    授权

    授权

  • 2015-07-08

    著录事项变更 IPC(主分类):G06T7/20 变更前: 变更后: 申请日:20140528

    著录事项变更

  • 2014-09-10

    实质审查的生效 IPC(主分类):G06T7/20 申请日:20140528

    实质审查的生效

  • 2014-08-13

    公开

    公开

说明书

技术领域

本发明属于跟踪监测技术领域,特别涉及一种彩色图像的目标跟踪方法,可用于 视频跟踪、目标监控等系统。

背景技术

彩色图像,由于其图像细节丰富,感光器件成本低,在实际中应用范围广,因此 针对彩色图像的目标跟踪有着广阔的应用前景和实际意义。丰富的图像细节,既给目 标特征描述提供了大量素材,同时也带来了特征甄别选取的难度。如何充分挖掘目标 的图像信息,去除其中的冗余,有效地组织和建立多特征描述模型,成为彩色图像目 标跟踪研究中的重点问题。

目前,针对彩色目标跟踪的特征描述模型主要有:联合特征直方图,协方差算子 等。其中:

联合特征直方图,存在大量的冗余信息,计算复杂度高,对光照变化敏感;

协方差算子,是一个高效精炼的多特征融合描述模型,具有较强的鲁棒性,但其 传统定义中所选取的特征主要针对目标在各个方向上纹理的描述,同方向上不同阶次 的梯度在某些情况下会存在一定程度的信息冗余,而且缺乏对彩色图像中丰富的颜色 信息的挖掘。

传统的彩色图像目标跟踪方法,包括基于联合特征直方图的目标跟踪方法和基于 协方差矩阵的目标跟踪方法。其中,基于联合直方图的跟踪方法选取红色r、绿色g 和蓝色b构成联合颜色直方图;基于协方差矩阵的跟踪方法中对目标所提取的特征包 括灰度值、x方向的一阶梯度、x方向的二阶梯度、y方向的一阶梯度、y方向的二阶 梯度。这些传统方法在实际情况下,由于所提取的特征间具有较高的相似性,特征间 的相互独立性较差,存在信息冗余,导致目标特征描述的准确性和鲁棒性较低,使得 跟踪系统的稳定性较差。

发明内容

本发明的目的在于针对上述已有技术的不足,在粒子滤波框架下,利用协方差算 子特征描述模型,融合了色相H、饱和度S、明度V和拉普拉斯算子响应,提出了一 种基于HSV颜色协方差特征的目标跟踪方法,以提高目标特征描述的准确性和鲁棒 性,实现对彩色图像中运动目标的实时精确跟踪。

实现本发明的技术关键是:在跟踪过程中,通过提取目标图像的色相、饱和度、 明度、拉普拉斯算子响应,获取目标图像的表观特征,通过表观特征对候选目标建立 协方差矩阵,以减少信息冗余,增强特征独立性,从而提升目标特征描述的精确性和 鲁棒性,提高跟踪精度。其实现步骤包括如下:

(1)初始化步骤:

(1a)读入k-1时刻的彩色图Ik-1,根据目标的初始状态,产生k-1时刻的初始 粒子集其中,表示k-1时刻第i个粒子的状态估计值,N表示粒子总 数,i表示粒子的序号,k表示时刻,初始时刻为k=1;

(1b)初始化目标跟踪窗:Bk-1=(rk-1,ck-1)T,其中rk-1和ck-1分别表示k-1时刻目 标跟踪窗的长度和宽度值,T表示向量转置;

(1c)根据目标初始状态与目标跟踪窗Bk-1,初始化目标的特征协方差矩阵M作 为特征模板;

(2)目标状态预测步骤:

(2a)读入k时刻的彩色图Ik,通过对k-1时刻彩色图中粒子集的传递, 得到k时刻彩色图中的预测粒子集其中为k时刻第i个粒子的状态预测值;

(2b)根据k时刻预测粒子集和目标跟踪窗Bk-1,确定k时刻候选目标集 其中为k时刻第i个候选目标,表示以为中心、Bk-1为长宽所界定出的 矩形区域;

(3)候选目标特征提取步骤:

(3a)对于k时刻彩色图Ik,提取其对应的W×H×d维特征图F,其中W表示 彩色图Ik的宽,H表示彩色图Ik的高,d表示特征维数;

(3b)在特征图F的基础上得到特征向量积分图IP和特征向量乘积积分图IQ;

(3c)根据特征向量积分图IP、特征向量乘积积分图IQ和候选目标集提取候选目标的特征集其中Ci表示第i个候选目标的特征协方差矩阵;

(4)计算特征权值步骤:

(4a)求取各候选目标特征集与特征模板M之间的距离集其中ρi表示第i个候选目标的特征协方差矩阵与特征模板之间的距离;

(4b)根据距离集计算候选目标的特征权值集其中ωi表示第i个候 选目标的特征权值;

(5)目标跟踪窗调整步骤:

(5a)对特征权值最大的候选目标的目标跟踪窗进行缩放,分别产生缩小目标 跟踪窗的候选目标和放大目标跟踪窗的候选目标其中β为特征权值最大 粒子的序号,下标small和big无具体的物理含义,仅表示所属变量分别为缩小目标跟 踪窗的候选目标和放大目标跟踪窗的候选目标;

(5b)分别提取缩放目标跟踪窗的候选目标和的特征协方差矩阵和

(5c)根据和计算和对应的特征权值和

(5d)根据与缩放目标跟踪窗前候选目标的特征权值ωβ的对比结果, 更新k时刻的目标跟踪窗Bk,并更新特征权值最大的候选目标的特征协方差矩阵Cβ及最大特征权值ωβ

(6)特征模板更新步骤:

通过计算特征模板M和特征权值最大的粒子对应的特征协方差矩阵Cβ的加权对 数-欧几里得均值,来对特征模板进行更新;

(7)目标状态更新步骤:

(7a)利用重采样算法,根据特征权值集对k时刻预测粒子集进行重 采样,得到k时刻的更新粒子集其中表示k时刻第i个粒子的状态估计值;

(7b)根据k时刻的更新粒子集估计k时刻的目标状态Xk;

(7c)根据k时刻的目标状态Xk和目标跟踪窗Bk,确定出k时刻目标的估计范围 Target;

(8)输出步骤:

输出k时刻的目标的估计范围Target,若下一时刻观测信息到达,令k=k+1, 转到步骤(2)进行迭代,否则,目标跟踪过程结束。

本发明具有以下优点:

(1)本发明由于利用协方差算子对目标进行描述,融合多特征,从而对干扰、 光照变化、以及目标形变具有较强的鲁棒性;

(2)本发明由于提取色相H、饱和度S、明度V和拉普拉斯算子响应,充分挖掘 目标彩色图像的丰富信息,并减少信息冗余,增强特征间的相互独立性,从而实现对 目标特征的精确描述,增强了目标跟踪的稳定性。

附图说明

图1是本发明的整体流程图;

图2是本发明方法的特征提取信息比较实验结果示意图;

图3是本发明方法的目标非刚性形变实验结果图;

图4是本发明方法的光照变化实验结果图;

图5是本发明方法的遮挡干扰实验结果图。

具体实施方式

参照图1,本发明的具体实施过程包括以下步骤:

步骤1.初始化粒子、目标跟踪窗和特征模板。

1.1)初始化粒子:令初始时刻k=1,根据目标的初始状态X0,产生N个粒子, 组成k-1时刻的粒子集其中表示k-1时刻第i个粒子的状态估计值, 服从均值为X0方差为Ψ的高斯分布,X0为目标的初始状态,Ψ为过程噪声方 差,i表示粒子的序号,N表示粒子总数;

1.2)初始化目标跟踪窗:Bk-1=(rk-1,ck-1)T,其中rk-1和ck-1分别表示k-1时刻目标 跟踪窗的长度和宽度值,T表示向量转置;

1.3)根据目标初始状态与目标跟踪窗Bk-1,初始化目标的特征协方差矩阵M作 为特征模板。

步骤2.目标状态预测。

2.1)读入k时刻的彩色图Ik,通过对k-1时刻彩色图中粒子集的传递, 得到k时刻彩色图中的预测粒子集

2.1.1)分别得到k时刻每一个粒子的状态预测值:

qki=pk-1i+wi,---1)

其中,为k时刻第i个粒子的状态预测值,为k-1时刻第i个粒子的状态估 计值,wi为状态噪声,服从均值为方差为Ψ的高斯分布;

2.1.2)用步骤2.1.1)中所得的N个预测粒子,组成预测粒子集

{qki}i=1N={qk1,qk2,···,qki,···,qkN};---2)

2.2)根据k时刻的预测粒子集和目标跟踪窗Bk-1,确定k时刻候选目标集

2.2.1)分别计算每一个候选目标区域:

Pki={(x,y)|x-υki|rk-12,|y-vki|ck-12},---3)

其中,表示第i个候选目标区域,和分别表示k时刻第i个预测粒子的 横坐标和纵坐标,rk-1和ck-1分别表示k-1时刻目标跟踪窗的长度和宽度值;

2.2.2)用步骤2.2.1)中所得的N个候选目标,组成候选目标集

{Pki}i=1N={Pk1,Pk2,···,Pki,···,PkN}.---4)

步骤3.提取候选目标特征。

3.1)对于k时刻的W×H维彩色图Ik,提取其W×H×d维特征图F:

F={F(x,y)|F(x,y)=[h(x,y),s(x,y),v(x,y),l(x,y)]T,1≤x≤W,1≤y≤H},5)

其中,d为特征维数,d=4,h(x,y)为彩色图Ik中像素点(x,y)处的色相值, s(x,y)为彩色图Ik中像素点(x,y)处的饱和度值,v(x,y)为彩色图Ik中像素点(x,y) 处的明度值,l(x,y)为彩色图Ik中像素点(x,y)处的拉普拉斯算子响应值,T表示向 量转置;

3.2)根据特征图F计算得到特征向量积分图IP和特征向量乘积积分图IQ:

3.2.1)根据特征图F计算特征向量积分图IP和特征向量乘积积分图IQ中的点:

IP(x',y',a)=Σ1xx',1yy'F(x,y,a),a=1,···,d,---6)

IQ(x',y',a,b)=Σ1xx',1yy'F(x,y,a)F(x,y,b),a,b=1,···,d,---7)

其中,a为第一组特征序号,b为第二组特征序号,d为特征维数,F(x,y,a)为 特征图F中的点,F(x,y,b)为特征图F中的点,IP(x′,y′,a)为特征向量积分图IP中 的点,IQ(x′,y′,a,b)为特征向量乘积积分图IQ中的点;

将积分图IP和IQ中的任意一点IPx,y和IQx,y分别表示如下:

IPx,y=[IP(x,y,1)…IP(x,y,d)]T,       8)

3.2.2)将步骤3.2.1)所得的点组成集合,得到特征向量积分图IP和特征向 量乘积积分图IQ:

IP={IPx,y|1≤x≤W,1≤y≤H},           10)

IQ={IQx,y|1≤x≤W,1≤y≤H},        11)

其中,IPx,y为特征向量积分图IP中的点,IQx,y为特征向量乘积积分图IQ中 的点,W表示IP和IQ的宽,H表示积分图IP和IQ的高;

3.3)根据特征向量积分图IP、特征向量乘积积分图IQ和候选目标集计 算候选目标的特征集

3.3.1)分别计算每一个候选目标的特征协方差矩阵:

Ci=1n-1[IQx'',y''-IQx',y''-IQx'',y'+IQx',y'-1n(IPx'',y''-IPx',y''-IPx'',y'+IPx',y')(IPx'',y''-IPx',y''-IPx'',y'+IPx',y')T],---12)

其中,Ci为第i个候选目标的特征协方差矩阵,(x′,y′)为第i个候选目标对 应矩形区域左上角的顶点坐标,(x′′,y′′)为第i个候选目标对应矩形区域右下角的顶 点坐标,n为区域内像素总数n=(x′′-x′)·(y′′-y′),IQx′′,y′′、IQx′,y′′、IQx′′,y′、IQx′,y′为 特征向量乘积积分图IQ中的不同点,IPx′′,y′′、IPx′,y′′、IPx′′,y′、IPx′,y′为特征向量积分 图IP中的不同点;

3.3.2)用步骤3.3.1)中所得的N个特征协方差矩阵,组成候选目标特征集

{Ci}i=1N={C1,C2,···,Ci,···,CN}.---13)

步骤4.求取粒子特征权值。

4.1)计算候选目标特征集与特征模板M之间的距离集

4.1.1)分别计算每一个候选目标特征与特征模板M的之间距离:

ρi=(Σa=1dln2λa)1/2,---14)

其中,ρi表示第i个候选目标特征与特征模板之间的距离,λa为特征协方差矩阵 Ci和特征模板M的广义特征值:λa=Mxa(Cixa)-1,Ci为第i个候选目标的特征协方 差矩阵,xa为广义特征向量,a表示特征序号,a=1,2,…,d;

4.1.2)用步骤4.1.1)中所得的N个距离,组成距离集

{ρi}i=1N={ρ1,ρ2,···,ρi,···ρN};---15)

4.2)根据距离集计算候选目标特征权值集

4.2.1)分别计算每一个候选目标的特征权值:

ωi=12πRexp{-(ρi)22R},---16)

其中,ωi表示第i个候选目标的特征权值,R为特征观测噪声方差;

4.2.2)用步骤4.2.1)中所求的N个候选目标的特征权值,组成候选目标特征权值 集

{ωi}i=1N={ω1,ω2,···,ωi,···ωN}.---17)

步骤5.目标跟踪窗调整。

5.1)对特征权值最大的候选目标的目标跟踪窗进行缩放,产生缩小目标跟踪 窗的候选目标区域和放大目标跟踪窗的候选目标区域

Pk,smallβ={(x,y)||x-υkβ|0.95rk-12,|y-vkβ|0.95ck-12},---18)

Pk,bigβ={(x,y)||x-υkβ|1.05rk-12,|y-vkβ|1.05ck-12},---19)

其中,β为特征权值最大粒子的序号,下标small和big无具体的物理含义,仅表示 所属变量分别为缩小目标跟踪窗的候选目标和放大目标跟踪窗的候选目标,和分别表示k时刻第β个预测粒子的横坐标和纵坐标,rk-1和ck-1分别表示k-1时刻目 标跟踪窗的长度和宽度值;

5.2)分别计算缩放目标跟踪窗的候选目标和的特征协方差矩阵和

5.2.1)计算缩小目标跟踪窗的候选目标对应的特征协方差矩阵

Csmallβ=1n1-1[IQx2,y2-IQx1,y2-IQx2,y1+IQx1,y1-1n1(IPx2,y2-IPx1,y2-IPx2,y1+IPx1,y1)(IPx2,y2-IPx1,y2-IPx2,y1+IPx1,y1)T],---20)

其中,(x1,y1)为缩小目标跟踪窗的候选目标对应矩形区域左上角的顶点坐 标,(x2,y2)为缩小目标跟踪窗的候选目标对应矩形区域右下角的顶点坐标,n1为区域内像素总数n1=(x2-x1)·(y2-y1),为特征向量乘 积积分图IQ中的不同点,为特征向量积分图IP中的不同 点;

5.2.2)计算放大目标跟踪窗的候选目标对应的特征协方差矩

Cbigβ=1n2-1[IQx4,y4-IQx4,y4-IQx2,y1+IQx3,y3-1n2(IPx4,y4-IPx3,y4+IPx4,y3+IPx3,y3)(IPx4,y4-IPx3,y4-IPx4,y4+IPx3,y3)T],---21)

其中,(x3,y3)为缩小目标跟踪窗的候选目标对应矩形区域左上角的顶点坐 标,(x4,y4)为缩小目标跟踪窗的候选目标对应矩形区域右下角的顶点坐标,n2为 区域内像素总数n2=(x4-x3)·(y4-y3),为特征向量乘积 积分图IQ中的不同点,为特征向量积分图IP中的不同 点;

5.3)根据特征协方差矩阵和计算对应的缩放目标跟踪窗的候选目标 和的特征权值和

5.3.1)计算特征协方差矩阵和与特征模板M之间的距离ρsmall和ρbig

ρsmall=(Σa=1dln2λa,small)1/2,---22)

ρbig=(Σa=1dln2λa,big)1/2,---23)

其中,λa,small为特征协方差矩阵和特征模板M的广义特征值: λa,big为特征协方差矩阵和特征模板M的广义特征 值:xa,small为λa,small对应的广义特征向量,xa,big为λa,big对应 的广义特征向量,a表示特征序号,a=1,2,…,d,下标small和big无具体的物理含 义;

5.3.2)根据距离ρsmall和ρbig计算对应的缩放目标跟踪窗的候选目标和的特征权值和

ωsmallβ=12πRexp{-(ρsmall)22R},---24)

ωbigβ=12πRexp{-(ρsmall)22R}---(25)

其中,R为特征观测噪声方差,ρsmall为特征协方差矩阵与特征模板M之间 的距离,ρbig为特征协方差矩阵与特征模板M之间的距离;

5.4)根据和与缩放目标跟踪窗前候选目标的特征权值ωβ比较的结果, 更新k时刻的目标跟踪窗Bk、特征协方差矩阵Cβ及最大特征权值ωβ

Bk=0.95Bk-1,ωsmallβ=ωmaxBk-1,ωβ=ωmax1.05Bk-1,ωbigβ=ωmax,---26)

Cβ=Csmallβ,ωsmallβ=ωmaxCβ,ωβ=ωmaxCbigβ,ωbigβ=ωmax,---27)

ωβ=ωsmallβ,ωsmallβ=ωmaxωβ,ωβ=ωmaxωbigβ,ωbigβ=ωmax,---28)

其中,ωmax=max{ωsmallβ,ωβ,ωbigβ}.

步骤6.特征模板更新。

6.1)计算特征模板M和特征权值最大粒子的协方差矩阵Cβ的加权对数-欧几里得 均值M′:

M′=exp[γlog(M)+(1-γ)log(Cβ)],       29)

其中,β为特征权值最大的粒子序号,γ为更新速率因子,0<γ<1。

6.2)用步骤6.1)所得均值M′对特征模板M进行更新:

M=M′。       30)

步骤7.目标状态更新。

7.1)利用重采样算法,根据特征权值集对k时刻预测粒子集重采 样,得到k时刻的更新粒子集

7.2)根据k时刻的更新粒子集估计k时刻的目标状态Xk:

Xk=1NΣi=11pki---31)

其中,N表示粒子总数,表示k时刻第i个更新粒子的状态估计值;

7.3)根据k时刻的目标状态Xk和目标跟踪窗Bk,确定k时刻的目标估计范围 Target:

Target={(x,y)||x-xk|rk2,|y-yk|ck2},---32)

其中,xk和yk分别表示k时刻目标状态Xk的横坐标和纵坐标,rk和ck分别表示k 时刻目标跟踪窗的长度和宽度值。

步骤8输出和迭代判断。

8.1)输出步骤7.3)所得的目标估计范围Target;

8.2)检查下一时刻的观测信息是否到达,若是,令k=k+1,转到步骤2进行迭 代,否则,目标跟踪过程结束。

本发明的效果可通过以下仿真实验进一步说明:

1.仿真条件。

仿真环境:计算机采用Intel Core i5-2400CPU3.1Ghz,4GB内存,软件采用 Matlab R2011a仿真实验平台。

仿真参数:粒子数N=100,过程噪声方差Ψ=diag([16,16]),特征观测噪声方差 R=0.1,更新速率因子γ=0.8。

2.仿真方法。

方法1:本发明方法;

方法2:现有基于联合特征直方图特征提取的粒子滤波跟踪方法,该方法选取r (红色)、g(绿色)和b(蓝色)构成联合颜色直方图;

方法3:现有基于传统协方差算子特征提取的粒子滤波跟踪方法,该方法选取灰 度值、x方向的一阶和二阶梯度、y方向的一阶和二阶梯度等5维特征融合构成协方 差算子。

3.仿真内容与结果。

仿真1:用所述三种方法,对目标特征进行提取,结果如图2所示,其中:

图2(a)为目标的原始彩色图像;

图2(b)为用方法1提取的一组特征:色相、饱和度、明度和拉普拉斯算子响应;

图2(c)为用方法2提取的一组特征:红色、绿色和蓝色;

图2(d)为用方法3提取的一组特征:灰度、x方向一阶梯度、x方向二阶梯度、y 方向一阶梯度和y方向二阶梯度。

从图2可以看出,方法2提取的特征中,红色、绿色和蓝色之间存在较强的相似 性,信息冗余大,特征间独立性差;方法3提取的特征中,x方向一阶梯度和二阶梯 度,y方向一阶梯度和二阶梯度等特征之间也存在相似性,且特征能量较弱,易受干 扰;而本发明提取的特征中,色相、饱和度、明度和拉普拉斯算子响应等特征之间有 显著差异,独立性强,信息冗余小,且特征能量较强,抗干扰能力强。

仿真2:用所述三种方法,在目标非刚性形变情况下进行跟踪,结果如图3所 示,其中:

图3(a)为用方法1对第1、26、34和67帧的跟踪结果;

图3(b)为用方法2对第1、26、34和67帧的跟踪结果;

图3(c)为用方法3对第1、26、34和67帧的跟踪结果。

从图3可以看出,方法2和方法3对目标的非刚性形变的适应性较差,跟踪结果 更容易发生偏离;而本发明的跟踪精度更高。

仿真3:用所述三种方法,在光照变化情况下进行跟踪,结果如图4所示,其 中:

图4(a)为用方法1对第1、21、61和81帧的跟踪结果;

图4(b)为用方法2对第1、21、61和81帧的跟踪结果;

图4(c)为用方法3对第1、21、61和81帧的跟踪结果。

从图4可以看出,方法2对光照变化造成目标特征整体发生偏移的适应性较差; 现有方法3在光照变暗时的跟踪结果较好,但在光照变亮时跟踪结果发生了较大的偏 离;而本发明对光照变化具有较强的鲁棒性。

仿真4:用所述三种方法,在遮挡干扰情况下进行跟踪,结果如图5所示,其 中:

图5(a)为用方法1对第1、27、48和63帧的跟踪结果;

图5(b)为用方法2对第1、27、48和63帧的跟踪结果;

图5(c)为用方法3对第1、27、48和63帧的跟踪结果。

从图5可以看出,方法2和方法3在遮挡干扰情况下,均会发生较大程度的偏 离;而本发明对遮挡干扰具有较强的鲁棒性。

对图3~图5中的场景分别进行100次蒙特卡罗实验,统计平均跟踪误差Err和 平均每帧运行时间RT,结果如表1所示。

表1

由表1中的统计数据可以看出:在跟踪误差方面,本发明在目标非刚性形变情况 下,比方法2和方法3降低了30%;在光照变化情况下,本发明比方法2和方法3降 低了60%~80%;在遮挡干扰情况下,本发明比方法2和方法3降低了15%~50%; 在运行时间方面,本发明与方法2、方法3相当。综上可以得出,在用时相当的情况 下,本发明在跟踪精度和鲁棒性方面均优于方法2和方法3。

去获取专利,查看全文>

相似文献

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

客服邮箱:kefu@zhangqiaokeyan.com

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

  • 服务号