首页> 中国专利> 一种基于稀疏自回归模型建模的多频信号去噪方法

一种基于稀疏自回归模型建模的多频信号去噪方法

摘要

本发明公开了一种基于稀疏自回归模型建模的多频信号去噪方法,其基于稀疏自回归模型,并利用多频信号自身的采样值构建多频信号的自适应过完备稀疏基;然后通过随机抽取自适应过完备稀疏基中不连续的多行构成冗余字典;接着采用正交匹配追踪算法获取多个冗余字典各自对应的向量在对应的冗余字典上的稀疏映射系数向量;之后对这些稀疏映射系数向量求平均向量作为信号复原时所要使用的系数;最后对原多频信号的去噪结果和将原多频信号倒置后的信号的去噪结果合并得到去噪复原信号;优点是计算复杂度低,去噪效果好,而且处理信噪比不同的信号的情况下去噪效果稳定。

著录项

  • 公开/公告号CN104484557A

    专利类型发明专利

  • 公开/公告日2015-04-01

    原文格式PDF

  • 申请/专利权人 宁波大学;

    申请/专利号CN201410719053.9

  • 发明设计人 宋欢欢;叶庆卫;周宇;王晓东;

    申请日2014-12-02

  • 分类号

  • 代理机构宁波奥圣专利代理事务所(普通合伙);

  • 代理人周珏

  • 地址 315211 浙江省宁波市江北区风华路818号

  • 入库时间 2023-12-17 04:31:51

法律信息

  • 法律状态公告日

    法律状态信息

    法律状态

  • 2017-05-03

    授权

    授权

  • 2015-04-29

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

    实质审查的生效

  • 2015-04-01

    公开

    公开

说明书

技术领域

本发明涉及一种信号去噪方法,尤其是涉及一种基于稀疏自回归模型(AR)建模 的多频信号去噪方法。

背景技术

当今,对于大型建筑的健康检查一般都是通过采集建筑上的振动信号,通过分析振 动信号来研究大型建筑的健康状况。然而,由于外界环境的影响和采集设备的局限,会 导致采集到的振动信号含有噪声,因此要先对采集到的振动信号进行降噪处理。

目前,信号降噪处理方法主要有小波去噪法、最小二乘去噪法、基于EMD(Empirical  Mode Decomposition,经验模态分解)阈值降噪法、基于FFT(快速傅里叶变换)降噪 法、中值滤波降噪法、稀疏降噪法等。上述这些降噪方法中小波去噪法是当前最常用的 去噪方法,但是小波去噪过程中,阈值的选择会影响小波去噪结果的好坏,而且在信号 不连续区域还会出现Gibbs现象。此外,上述这些降噪方法存在共同的缺陷,即去噪效 果一般,而且处理信号不同的情况下会导致去噪效果的不稳定。

发明内容

本发明所要解决的技术问题是提供一种基于稀疏自回归模型建模的多频信号去噪 方法,其计算复杂度低,去噪效果好,而且处理信噪比不同的信号的情况下去噪效果稳 定。

本发明解决上述技术问题所采用的技术方案为:一种基于稀疏自回归模型建模的多 频信号去噪方法,其特征在于包括以下步骤:

①将待处理的多频信号以向量形式表示为x=x1x2...xnT,其中, (x1 x2 … xn)T为(x1 x2 … xn)的转置向量,n表示多频信号的采样点数,n≥500, x1表示多频信号的第1个采样值,x2表示多频信号的第2个采样值,xn表示多频信号的 第n个采样值;

②基于稀疏自回归模型,构造的自适应过完备稀疏基,记为Z, Z=x1x2...xpx2x3...xp+1.........xn-pxn-p+1...xn-1,然后根据Z构造Z对应的向量,记为bZ=xp+1xp+2...xn,再令j表示计算次数,并令j的初始值为1,其中,Z的维数为(n-p)×p,的维数为 (n-p)×1,p表示稀疏自回归模型的阶数,x3、xp、xp+1、xp+2、xn-p、 xn-p+1和xn-1对应表示多频信号的第3个采样值、第p个采样值、第p+1个采样值、第 p+2个采样值、第n-p个采样值、第n-p+1个采样值和第n-1个采样值,1≤j≤N, N表示设定的计算总次数;

③在第j次计算过程中,从Z中随机抽取不连续的m行按序构成第j个冗余字典, 记为ΦjΦj=xi1xi1+1...xi1+p-1xi2xi2+1...xi2+p-1.........ximxim+1...xim+p-1,然后根据Φj构造Φj对应的向量,记为bΦj=xi1+pxi2+p...xim+p,其中,15<m<p,Φj的维数为m×p,的维数为m×1,i1,i2,im对应 表示从Z中随机抽取出的m行中的第1行、第2行、第m行在Z中的行序号,且 1≤i1<i2<…<im≤n-p,xi1、xi1+1、xi1+p-1、xi2、xi2+1、xi2+p-1、xim、xim+1、xim+p-1、 xi1+p、xi2+p、xim+p对应表示多频信号的第i1个采样值、第i1+1个采样值、第i1+p-1个 采样值、第i2个采样值、第i2+1个采样值、第i2+p-1个采样值、第im个采样值、第im+1 个采样值、第im+p-1个采样值、第i1+p个采样值、第i2+p个采样值、第im+p个采 样值;

④采用正交匹配追踪算法,计算在Φj上的稀疏映射系数向量,记为其中, 的维数为p×1;

⑤判断j=N是否成立,如果成立,则直接执行步骤⑥,如果不成立,则令j=j+1, 然后返回步骤③继续执行,其中,j=j+1中的“=”为赋值符号;

⑥对得到的N个稀疏映射系数向量做算术平均,得到平均向量,记为将作为信号复原所需的稀疏映射系数向量,其中,表示第1个冗余字典Φ1对应的向量在第1个冗余字典Φ1上的稀疏映射系数向量,表示第2个冗余字典Φ2对应的向量在第2个冗余字典Φ2上的稀疏映射系数向量,表示第N个冗余字典ΦN对应的向量在第N个冗余字典ΦN上的稀疏映射系数向量;

⑦根据Z和计算去噪复原后的后n-p个采样值构成的列向量,记为y1=Za^=yp+1yp+2...yn,其中,yp+1表示去噪复原后的第p+1个的采样值,yp+2表示去噪复原 后的第p+2个的采样值,yn表示去噪复原后的第n个的采样值;

⑧对进行倒置,得到的倒置向量,记为x=xnxn-1...x1T,然后按照 步骤②至步骤⑥的操作过程,以相同的方式获得的自适应过完备稀疏基和平均向量, 对应记为Z'和再根据Z'和计算去噪复原后的前n-p个采样值构成的列向量的 倒置向量,记为y1=Za^=yn-pyn-p-1...y11,其中,(xn xn-1 … x1)T为(xn xn-1 … x1) 的转置向量,yn-p表示去噪复原后的第n-p个采样值,yn-p-1表示去噪复原后的第 n-p-1个采样值,y1表示去噪复原后的第1个采样值;

⑨对进行倒置,得到的倒置向量,记为y1=y1y2...yn-p,然后将中的前 p个采样值和中的所有采样值按序构成的去噪复原信号,记为y=y1y2...ypyp+1yp+2...yn,至 此完成多频信号的去噪过程。

所述的步骤②中N∈[20,40]。

所述的步骤④的具体过程为:

④-1、令t表示迭代的次数,令rt表示第t次迭代的残差,令Λt表示第t次迭代的索 引集,令K表示在Φj上的稀疏映射系数向量的稀疏度,K的值用于代表在Φj上 的稀疏映射系数向量中的非零元素的总个数,其中,t的初始值为1,K≥1;

④-2、计算第t-1次迭代的残差rt-1与Φj中的每列的内积,然后从计算得到的p个 内积值中选出最大值,再将Φj中与该最大值对应的一列的脚注记为λt,其中,rt-1表示 第t-1次迭代的残差,当t=1时rt-1的值为λt∈[1,p];

④-3、令Λt=Λt-1∪{λt},其中,Λt-1表示第t-1次迭代的索引集,当t=1时Λt-1的 值为空集,符号“∪”为并集运算符号,在此符号“{}”表示集合符号;

④-4、根据和Φj,计算第t次迭代时在Φj上的稀疏映射系数向量,记为atj=atgmin||bΦj-Φja||2,其中,符号“|| ||2”为求L2范数符号,表示 取使得的值最小的a,a表示稀疏映射系数向量;

④-5、令rt=rt-1-Φjatj;

④-6、判断t=K是否成立,如果成立,则结束迭代过程,将作为在Φj上的 稀疏映射系数向量,重新记为如果不成立,则令t=t+1,然后返回步骤④-2继续迭 代,其中,t=t+1中的“=”为赋值符号。

与现有技术相比,本发明的优点在于:

1)本发明方法基于稀疏自回归模型,并利用多频信号自身的采样值构建多频信号 的自适应过完备稀疏基,使得对于不同的多频信号都有自己的自适应过完备稀疏基,有 效地提高了本发明方法的实用性。

2)本发明方法中的冗余字典是通过随机抽取自适应过完备稀疏基中不连续的m行 构成的,因此冗余字典也是自适应的,这样随着多频信号的改变冗余字典也会发生改变, 与使用固定冗余字典相比,本发明方法构造的冗余字典更具灵活性、实用性,更能很好 的捕捉多频信号的结构特征,从而有效地提高了去噪效果,同时使得本发明方法在处理 信噪比不同的信号的情况下去噪效果稳定;此外,因冗余字典是通过随机抽取自适应过 完备稀疏基中不连续的m行构成的,因此计算复杂度低。

3)本发明方法通过重复多次随机抽取自适应过完备稀疏基中不连续的m行构建多 个冗余字典,求得多个对应的稀疏映射系数向量,然后对这些稀疏映射系数向量求算术 平均值,将得到的平均向量作为信号复原时所要使用的系数,这有效地降低了随机抽取 带来的去噪的不稳定性。

4)本发明方法通过对采集到的原多频信号和将原多频信号倒置后的信号进行去噪 处理,并最后合并两次去噪结果,很好地克服了一次去噪中前p个信号未被去噪复原的 缺点。

附图说明

图1为本发明方法的流程框图;

图2为实验仿真中原始信号的图形;

图3为原始的无噪声信号含有信噪比为SNR0=10dB的噪声的信号利用本 发明方法对进行去噪后得到的去噪复原信号的波形对比图;

图4为原始的无噪声信号含有信噪比为SNR0=10dB的噪声的信号利用现 有的小波去噪方法对进行去噪后得到的去噪复原信号的波形对比图;

图5为对原始的无噪声信号加不同大小噪声得到的信号的信噪比SNR0、利用 本发明方法对无噪声信号加不同大小噪声信号后得到的信号(的信噪比为原始信 噪比SNR0)去噪后得到的去噪复原信号的信噪比SNR1和利用现有的小波去噪方法对无 噪声信号加不同大小噪声信号后得到的信号(的信噪比为原始信噪比SNR0)去 噪后得到去噪复原信号的的信噪比SNRw的值对比图。

具体实施方式

以下结合附图实施例对本发明作进一步详细描述。

本发明提出的一种基于稀疏自回归(AR)模型建模的多频信号去噪方法,其流程 框图如图1所示,其包括以下步骤:

①将待处理的多频信号以向量形式表示为x=x1x2...xnT,其中, (x1 x2 … xn)T为(x1 x2 … xn)的转置向量,n表示多频信号的采样点数,n≥500, x1表示多频信号的第1个采样值,x2表示多频信号的第2个采样值,xn表示多频信号的 第n个采样值。

在此,n的取值最好为大于或等于500且小于或等于2000,如取n=1000,这是因 为:如果n的取值过小,则后续构造的自适应过完备稀疏基可能不能很完整的包含多频 信号的结构特征,最终影响多频信号的去噪结果;如果n的取值过大,则将构造出很大 的自适应过完备稀疏基,这样会影响后续去噪处理的计算速度。

②基于稀疏自回归模型,构造的自适应过完备稀疏基,记为Z, Z=x1x2...xpx2x3...xp+1.........xn-pxn-p+1...xn-1,然后根据Z构造Z对应的向量,记为bZ=xp+1xp+2...xn,再令j表示计算次数,并令j的初始值为1,其中,Z的维数为(n-p)×p,的维数为 (n-p)×1,p表示稀疏自回归模型的阶数,p的具体值为在范 围内通过实验选择最佳的值,如当n=1000时取p=300,x3、xp、xp+1、xp+2、xn-p、 xn-p+1和xn-1对应表示多频信号的第3个采样值、第p个采样值、第p+1个采样值、第 p+2个采样值、第n-p个采样值、第n-p+1个采样值和第n-1个采样值,1≤j≤N, N表示设定的计算总次数,N∈[20,40],在本实施例中,对该多频信号去噪方法经过多 次实验仿真,发现N在[20,40]这个范围内取值,去噪效果能够达到最好,在具体操作时 如可取N的值为30。

③在第j次计算过程中,从Z中随机抽取不连续的m行构成第j个冗余字典,记为 ΦjΦj=xi1xi1+1...xi1+p-1xi2xi2+1...xi2+p-1.........ximxim+1...xim+p-1,然后根据Φj构造Φj对应的向量,记为bΦj=xi1+pxi2+p...xim+p,其中,15<m<p,Φj的维数为m×p,的维数为m×1,i1,i2,im对应 表示从Z中随机抽取出的m行中的第1行、第2行、第m行在Z中的行序号,且 1≤i1<i2<…<im≤n-p,xi1、xi1+1、xi1+p-1、xi2、xi2+1、xi2+p-1、xim、xim+1、xim+p-1、 xi1+p、xi2+p、xim+p对应表示多频信号的第i1个采样值、第i1+1个采样值、第i1+p-1个 采样值、第i2个采样值、第i2+1个采样值、第i2+p-1个采样值、第im个采样值、第im+1 个采样值、第im+p-1个采样值、第i1+p个采样值、第i2+p个采样值、第im+p个采 样值。

④采用现有的正交匹配追踪(OMP)算法,计算在Φj上的稀疏映射系数向量, 记为其中,的维数为p×1。

在此具体实施例中,步骤④的具体过程为:

④-1、令t表示迭代的次数,令rt表示第t次迭代的残差,令Λt表示第t次迭代的索 引集,令K表示在Φj上的稀疏映射系数向量的稀疏度,即K的值用于代表在Φj上的稀疏映射系数向量中的非零元素的总个数,其中,t的初始值为1,K≥1。

④-2、计算第t-1次迭代的残差rt-1与Φj中的每列的内积,然后从计算得到的p个 内积值中选出最大值,再将Φj中与该最大值对应的一列的脚注记为λt,其中,rt-1表示 第t-1次迭代的残差,当t=1时rt-1的值为λt∈[1,p]。

④-3、令Λt=Λt-1∪{λt},其中,Λt-1表示第t-1次迭代的索引集,当t=1时Λt-1的 值为空集,符号“∪”为并集运算符号,在此符号“{}”表示集合符号。

④-4、根据和Φj,计算第t次迭代时在Φj上的稀疏映射系数向量,记为atj=atgmin||bΦj-Φja||2,其中,符号“|| ||2”为求L2范数符号,表示 取使得的值最小的a,a表示稀疏映射系数向量。

④-5、令rt=rt-1-Φjatj.

④-6、判断t=K是否成立,如果成立,则结束迭代过程,将作为在Φj上的 稀疏映射系数向量,重新记为如果不成立,则令t=t+1,然后返回步骤④-2继续 迭代,其中,t=t+1中的“=”为赋值符号。

⑤判断j=N是否成立,如果成立,则直接执行步骤⑥,如果不成立,则令j=j+1, 然后返回步骤③继续执行,其中,j=j+1中的“=”为赋值符号。

⑥对得到的N个稀疏映射系数向量做算术平均,得到平均向量,记为将作为信号复原所需的稀疏映射系数向量,其中,表示第1个冗余字典Φ1对应的向量在第1个冗余字典Φ1上的稀疏映射系数向量,表示第2个冗余字典Φ2对应的向量在第2个冗余字典Φ2上的稀疏映射系数向量,表示第N个冗余字典ΦN对应的向量在第N个冗余字典ΦN上的稀疏映射系数向量。

⑦根据Z和计算去噪复原后的后n-p个采样值构成的列向量,记为y1=Za^=yp+1yp+2...yn,其中,yp+1表示去噪复原后的第p+1个的采样值,yp+2表示去噪复原 后的第p+2个的采样值,yn表示去噪复原后的第n个的采样值。

⑧对进行倒置,得到的倒置向量,记为x=xnxn-1...x1T,然后按照 步骤②至步骤⑥的操作过程,以相同的方式获得的自适应过完备稀疏基和平均向量, 对应记为Z'和即具体过程为:⑧-1、基于稀疏自回归模型,构造的自适应过完 备稀疏基,记为Z',Z=xnxn-1...xn-p+1xn-1xn-1...xn-p.........xp+1xp...x2,然后根据Z'构造Z'对应的向量, 记为bZ=xn-pxn-p-1...x1,再令j表示计算次数,并令j的初始值为1;其中,xn和xn-2对 应表示多频信号的第n个采样值和第n-2个采样值。⑧-2、在第j次计算过程中,从Z' 中随机抽取不连续的m行按序构成第j个冗余字典,记为Φj', Φj=xi1xi1-1...xi1-p+1xi2xi2-1...xi2-p+1.........ximxim-1...xim-p+1,然后根据Φj'构造Φj'对应的向量,记为bΦj=xi1-pxi2-p...xim-p;其中,n≥i1>i2>…>im≥p+1。⑧-3、采用现有的正交匹配追踪(OMP) 算法,计算在Φj'上的稀疏映射系数向量,记为⑧-4、判断j=N是否成立, 如果成立,则直接执行步骤⑧-5,如果不成立,则令j=j+1,然后返回步骤⑧-2继续 执行;⑧-5、对得到的N个稀疏映射系数向量做算术平均,得到平均向量,记为将作为信号复原所需的稀疏映射系数向量,其中, 表示第1个冗余字典Φ1'对应的向量在第1个冗余字典Φ1'上的稀疏映射系数向 量,表示第2个冗余字典Φ2'对应的向量在第2个冗余字典Φ2'上的稀疏映射系 数向量,表示第N个冗余字典ΦN'对应的向量在第N个冗余字典ΦN'上的稀疏 映射系数向量;再根据Z'和计算去噪复原后的前n-p个采样值构成的列向量的倒 置向量,记为y1=Za^=yn-pyn-p-1...y11,其中,(xn xn-1 … x1)T为(xn xn-1 … x1) 的转置向量,yn-p表示去噪复原后的第n-p个采样值,yn-p-1表示去噪复原后的第 n-p-1个采样值,y1表示去噪复原后的第1个采样值。

⑨对进行倒置,得到的倒置向量,记为y1=y1y2...yn-p,然后将中的前 p个采样值和中的所有采样值按序构成的去噪复原信号,记为y=y1y2...ypyp+1yp+2...yn,至 此完成多频信号的去噪过程。

为说明本发明方法的可行性,针对具体的加有高斯白噪声的多频信号进行去噪处 理。

1)假设待处理的多频信号为x1(t),x1(t)=x0(t)+n(t),其中, x0(t)=sin(2π×100t+0.1)+3cos(2π×130t+0.9)+0.7sin(2π×170t+1.5)+1.2cos(2π×70t),n(t) 是高斯白噪声,在此t为时间变量;然后设置x1(t)的信噪比SNR0=10dB;再对x1(t)进 行奈奎斯特均匀采样,采样频率Fs=1000Hz,采样x1(t)得到的信号写成向量形式为 x1=x1x2...x999x1000T,采样x0(t)得到的信号向量为 x0=x10x20...x9990x10000T,如图2所示。

2)设置需构建的自适应过完备稀疏基有300列,即设置参数p=300,用采样得到 的信号向量构建的自适应过完备稀疏基Z=x1x2...x299x300x2x3...x300x301............x699x700...x997x998x700x701...x998x999,再根 据Z构造Z对应的向量bz=x301x302...x999x1000.

3)设置冗余字典的行数为100,即设置参数m=100,从Z中随机抽取不连续的100 行按序构成第j个冗余字典Φj=xi1xi1+1...xi1+299xi2xi2+1...xi2+299.........xi100xi100+1...xi100+299,然后根据Φj构造Φj对应 的向量ΦΦj=xi1+300xi2+300...xi100+300,其中,1≤i1<i2<…<im≤700,如Φj中的第1行是随机抽取 出的Z中的第5行,即i1=5。

4)计算在Φj上的稀疏映射系数向量采用OMP算法求解稀疏系数向量假设K=15,则得到aj=a1a2...a300,a1,a2,…,a300分别表示中的第1个值、第2个值、…、 第300个值。

5)设置重复计算次数N=25,重复执行步骤3)和4)共25次,可得到25个冗 余字典、对应的25个向量及25个稀疏映射系数向量。

6)对求得的25个稀疏映射系数向量做算术平均得到平均向量 a^=125(a1+a2+...+a25).

7)根据Z和计算去噪复原后的后700个采样值构成的列向量y1=Za^y301y302...y1000,其中,y301,y302,…,y1000分别是去噪后复原信号的第301个值、第302个值、…、第 1000个值。

8)对进行倒置操作,得到一个新的信号向量x1=x1000x999...x2x1,然后 按照步骤2)至步骤6)的操作过程,以相同的方式获得的自适应过完备稀疏基 Z=x1000x999...x702x701x999x998...x701x700............x302x301...x4x3x301x300...x3x2,对应的向量bZ=x700x699...x2x1,平均向量再根 据Z'和计算去噪复原后的前700个采样值构成的列向量的倒置向量y1y700y699...y2y1.

9)对进行倒置得到y1=y1y2...y699y700,其中,y1,y2,…,y699,y700是最终去噪后复原 信号的第1个值、第2个值、…、第699个值、第700个值;然后取信号的前300 个信号点与的所有信号点合并得到最终完整的去噪复原信号 y=y1y2...y999y1000.

上述得到的去噪复原信号y=y1y2...y999y1000的信噪比为 用MATLAB仿真分别求取的原始信噪比SNR0和的信噪比SNR1,SNR0=9.54dB,SNR1=18.89dB。

设置含有不同信噪比的信号采用本发明方法和现有的小波去噪方法分别对这些 信号进行去噪,其中将原始信噪比SNR0设置在(1dB~20dB)。本发明方法的实验参数 为:Fs=1000Hz,p=300,m=100,N=25。现有的小波去噪方法(小波函数wden()) 的实验参数为:启发式阈值函数形式(heursure),并设定为软门限阈值(s)处理,阈值门 限是由小波去噪函数自动确定的,只根据第一层小波分解(sln)稀疏估计噪声水平,用 sym8小波对信号作5层分解完成对信号的去噪。

上述待处理信号是一多频信号,符合稀疏自回归模型(AR)模型,即具有较强的 时间相关性。图3给出了原始的无噪声信号含有信噪比为SNR0=10dB的噪声的信 号利用本发明方法对进行去噪后得到的去噪复原信号的波形对比图,去噪复原 信号的信噪比SNR1为18.89dB;图4给出了原始的无噪声信号含有信噪比为 SNR0=10dB的噪声的信号利用现有的小波去噪方法对进行去噪后得到的去噪复 原信号的波形对比图,去噪复原信号的信噪比SNRw为12.86dB。从图3和图4中 可以看出,采用本发明方法对信号进行去噪得到的去噪复原信号比采用小波去噪方法 对信号进行去噪得到的去噪复原信号更接近于原始的无噪声信号的波形。通过比 较去噪后得到的去噪复原信号的信噪比,可以更客观的比较出本发明方法比现有的小波 去噪方法的去噪效果好。图5给出了对原始的无噪声信号加噪后得到信号的信噪比 SNR0、采用本发明方法和现有的小波去噪方法分别对无噪声信号加不同大小噪声信 号后得到的信号进行去噪处理得到的去噪复原信号的信噪比对比图。从图5中可以看 出本发明方法的去噪效果比现有的小波去噪方法的去噪效果好,现有的小波去噪方法随 着原始的无噪声信号的信噪比的增大去噪效果相对减弱,即对于原来含噪声信号相对 小的信号去噪效果不是很好。而本发明方法对含不同大小信噪比的信号去噪效果稳定, 即本发明方法对含不同大小噪声的信号都有很好的去噪效果,这是因为本发明方法的自 适应过完备稀疏基Z是用采集到的带噪声信号构建的,冗余字典是随机抽取Z的m行构 建的,即自适应过完备稀疏基和冗余字典具有自适应性,能够很好地利用信号本身特点; 本发明方法中涉及了随机抽取步骤,若只进行一次抽取,仿真结果会不稳定,本发明方 法中采取了重复取平均的操作来稳定去噪效果。

去获取专利,查看全文>

相似文献

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

客服邮箱:kefu@zhangqiaokeyan.com

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

  • 服务号