首页> 中国专利> 基于新息速率优化和抗差估计的GNSS/INS紧组合欺骗检测方法

基于新息速率优化和抗差估计的GNSS/INS紧组合欺骗检测方法

摘要

本发明涉及一种基于新息速率优化和抗差估计的GNSS/INS紧组合欺骗检测方法,属于导航系统欺骗检测技术领域。对于微小突变的阶跃式欺骗或缓慢增长的斜率式欺骗,本发明与现有技术中通过新息向量和协方差直接计算新息速率然后判断的方案相比,先获取归一化新息向量,然后通过卡尔曼滤波器对归一化新息向量实时更新。在时间更新阶段,通过先验估计得到先验估计状态向量;在量测更新阶段,根据归一化新息向量和先验估计状态向量对新息向量进行更新,根据更新后的新息向量和增益矩阵,计算后验估计状态向量,根据后验估计状态向量和新息速率矩阵来计算新息速率,进而结合预设的上、下限检测阈值,判断是否存在欺骗,检测欺骗所用的时间更短,检测效率提高。

著录项

法律信息

  • 法律状态公告日

    法律状态信息

    法律状态

  • 2022-08-26

    实质审查的生效 IPC(主分类):G01S19/21 专利申请号:2022104200514 申请日:20220420

    实质审查的生效

说明书

技术领域

本发明涉及一种基于新息速率优化和抗差估计的GNSS/INS紧组合欺骗检测方法,属于导航系统欺骗检测技术领域。

背景技术

现有技术中,无人机普遍采用全球卫星导航系统和惯性导航系统的组合导航系统进行导航控制。全球卫星导航系统(Global Navigation Satellite System,GNSS)和惯性导航系统(Inertial Navigation System,INS)的组合导航系统具有两者互补误差特性,GNSS能够在全球范围内提供全天候连续性的位置、速度和时间(Position,Velocity&Time,PVT)服务的优点,而INS具有独立自主、连续工作、提供短期抗干扰能力等优点,两者组合增加了系统的冗余度和可靠性。GNSS/INS组合导航可以分为松组合、紧组合和深组合三种组合模式。紧组合是利用GNSS接收机输出的伪距、伪距率等观测量,与INS结合星历反算得到的伪距、伪距率进行组合。

然而,由于GNSS信号功率低、结构公开,使得GNSS服务易受欺骗干扰的影响,第三方可能通过施加欺骗干扰来干涉无人机的运动,此时,无人机的飞行轨迹偏离正常飞行轨迹。然而,当无人机的飞行轨迹存在偏离时,却难以判断是由于无人机本身存在的问题导致的,还是由于第三方施加欺骗干扰导致的。为此,有必要进行实时、快速、准确地欺骗检测,以确保组合导航系统的可靠性和完好性。

欺骗干扰是指干扰源产生与真实信号高度相似的欺骗信号或转发真实信号欺骗目标接收机,迫使其生成错误和可能危险的信息。在GNSS/INS组合导航系统中,GNSS模块锁定欺骗信号并输出错误的信息,进而影响卡尔曼滤波测量更新阶段状态误差的估计值,并输出错误的导航结果,同时错误的状态误差估计值通过信息融合反馈至INS,最终影响GNSS/INS组合导航系统,使得无人控制设备的轨迹发生较大偏移,可能影响设备安全。

GNSS/INS组合导航系统欺骗检测可分两类,第一类是基于测量值的方法,常见的有基于新息或残差的卡方检测法、自主完好性检测外推法(Autonomous IntegrityMonitored Extrapolation,AIME)、扩展接收机自主完好性检测法(Extended ReceiverAutonomous Integrity Monitoring,ERAIM)和速率检测法等;第二类是基于结论的方法,例如多解分离法(Multiple Solution Separation,MSS)。基于新息的卡方检测法使用卡尔曼滤波新息作为检测统计量,具有成本低、效率高和计算量小等优点,是一种应用广泛的检测方法。该方法又可分成“快照法”和“连续法”,“快照法”是以当前时刻的新息向量及其协方差阵构成检验统计量,适合检测GNSS测量值中阶跃式的误差;“连续法”是以过去某时刻到当前时刻所有的新息向量及其协方差阵构成检验统计量,适合检测缓慢增长式的误差。有文献专门介绍了阶跃式和缓慢增长式的误差分类和模型,同时将GNSS欺骗干扰的误差特征分为同步式欺骗和异步式欺骗,同步式欺骗信号在初期与真实信号保持同步,当接收机跟踪环路锁定欺骗信号后,为防止接收机自主完好性监测,采取逐步诱导的方式,使伪距测量值逐渐偏离真实值,相当于对真实伪距的测量值施加了斜率式欺骗;异步欺骗信号的码在相位和载波频率上与真实信号无法保持同步,所以其伪距测量值是任意的,相当于对真实伪距施加了阶跃式欺骗,表明该误差特征类似于完好性分析中的阶跃式和斜坡式欺骗。

GNSS/INS组合导航系统在检测欺骗干扰时存在的难点是:对微小突变的阶跃式欺骗和缓慢增长的斜坡式欺骗进行检测时的时延问题。针对这一难点,WANG ShiZhuang等提出了解决微小突变的阶跃式欺骗的算法是“归一化微小突变幅度”,将实际微小突变的欺骗影响值除以归一化新息协方差的值,在单路通道受5m的微小突变干扰时,采用抗差估计检测时间为28s;解决缓慢增长的斜坡式欺骗的算法是“归一化滑动窗口斜坡幅度”,将在检测窗口内实际缓慢增长的欺骗影响值除以归一化新息协方差的值,在单路通道受0.1m/s的缓慢增长干扰时,采用抗差估计检测时间缩短了约30s。ZHANG Chuang等提出了基于抗差估计和滑动窗口的改进检测算法,对于缓慢增长的欺骗干扰,选择两个合适的阈值计算权重因子,并自适应调整增益矩阵,在单路通道受到0.5m/s斜坡式干扰时漏警率为28%。以上算法存在的不足之处是都需要设置合理的检测窗口,即在一段检测时间内,不包含错误测量的长度成为“浪费窗口”,这对检测时间和性能均有负面影响,因此,检测窗口需要具体情况确定。

针对GNSS/INS组合导航系统欺骗检测需要设置检测窗口的问题,Bhatti等提出了新息速率检测算法,通过归一化新息的变化速率来判断GNSS测量值是否存在异常,而后采用卡尔曼滤波实时估计归一化新息速率,该方法有效提升了检测效率,但受组合导航系统中闭环校正的影响,其余正常通道的新息速率也受到影响,在单路通道受到0.1m/s的缓慢增长式干扰时检测时间为110s。张超等提出了新息速率抗差估计检测算法,该算法能够有效的抑制欺骗干扰对状态向量的影响,提高了数据使用率和算法可靠性,在单路通道0.1m/s的缓慢增长式干扰式检测时,检测时间缩短了60%,漏警率和误警率维持在4%以内。但这两种算法在检测微小突变和缓慢增长的欺骗干扰时,仍然存在检测时间过长,甚至检测不敏感的现象。

综上,现有的欺骗干扰检测算法在检测微小突变和缓慢增长的欺骗干扰时,检测时间过长,甚至存在检测不敏感的现象,导致无人机飞行了一段较长的路程后,才能准确判断是否存在第三方施加的欺骗干扰,严重影响无人机飞行过程的安全。

发明内容

本发明的目的在于提供一种基于新息速率优化和抗差估计的GNSS/INS紧组合欺骗检测方法,用于解决现有欺骗检测干扰算法在检测微小突变和缓慢增长的欺骗干扰时,检测时间过长的问题。

为了实现上述目的,本发明提供了一种基于新息速率优化和抗差估计的GNSS/INS紧组合欺骗检测方法,包括如下步骤:

S1、通过观测模型H

S2、通过卡尔曼滤波器对归一化新息向量ω

在时间更新阶段,通过对GNSS/INS紧组合系统的状态向量进行先验估计,得到k时刻的先验估计状态向量

在量测更新阶段,将归一化新息向量ω

S3、若新息速率向量

本发明在对GNSS/INS紧组合导航系统进行观测时,先获取归一化新息向量,然后通过卡尔曼滤波器对归一化新息向量实时更新。在时间更新阶段,通过先验估计得到先验估计状态向量;在量测更新阶段,根据归一化新息向量和先验估计状态向量对新息向量进行更新,根据更新后的新息向量和提前获取的增益矩阵,计算后验估计状态向量,然后根据后验估计状态向量和提前获取的新息速率矩阵来计算新息速率,进而结合预设的上、下限检测阈值,判断是否存在欺骗。对于微小突变的阶跃式欺骗或缓慢增长的斜率式欺骗而言,本发明与现有技术中通过新息向量和协方差直接计算新息速率然后判断的方案相比,先利用归一化新息向量对新息向量进行更新,再利用更新后的新息向量,结合增益矩阵,计算出后验估计状态向量后再以此为基础,计算得到新息速率,使得检测出欺骗所用的时间更短,检测效率提高。

进一步地,在上述方法中,步骤S3中,根据

对于不能够准确判断出是否存在欺骗的情况,引入抗差估计来削弱欺骗干扰的影响,抗差估计时,通过新息速率结合上、下限检测阈值来计算对应第i个测量值的权重,进而通过权重对增益矩阵进行更新,能够提高欺骗检测的敏感能力。

进一步地,在上述方法中,通过如下公式对增益矩阵K

K

式中,K

提供一种较为简单的更新增益矩阵的方式,便于本发明的实施。

进一步地,在上述方法中,等价权矩阵表示为W=diag(w

式中,

提供一组较为简单的公式来计算等价权矩阵中的元素,便于本发明的实施。

进一步地,在上述方法中,步骤S2中,还在时间更新阶段得到先验协方差估计向量

式中,H

通过在时间更新阶段对协方差进行先验估计,得到先验协方差估计向量,在此基础上,提供一种较为简单的计算增益矩阵的方式,便于本发明的实施。

进一步地,在上述方法中,步骤S2中,通过如下公式计算后验估计状态向量

式中,

通过增益矩阵对更新的新息向量进行增益,然后对先验估计状态向量进行优化,从而得到后验估计状态向量,计算简单,便于实施。

进一步地,在上述方法中,步骤S2中,通过如下公式计算得到先验协方差估计向量

式中,

进一步地,在上述方法中,步骤S2中,通过如下公式计算得到先验估计状态向量

式中,

提供两个具体的公式,用于在时间更新阶段,计算先验估计状态向量和先验协方差估计向量,计算简单,便于实施。

附图说明

图1为本发明方法实施例中基于新息速率优化和抗差估计的GNSS/INS紧组合欺骗检测方法的流程框图;

图2为本发明方法实施例中目标接收机受欺骗干扰示意图;

图3为本发明方法实施例中欺骗干扰在卡尔曼滤波时的影响分析示意图;

图4为本发明方法实施例中模拟的飞行轨迹示意图;

图5(a)为本发明方法实施例中阶跃式欺骗下新息速率欺骗检测算法的仿真结果;

图5(b)为本发明方法实施例中阶跃式欺骗下新息速率抗差估计欺骗检测算法的仿真结果;

图6(a)为本发明方法实施例中斜坡式欺骗下新息速率欺骗检测算法的仿真结果;

图6(b)为本发明方法实施例中斜坡式欺骗下新息速率抗差估计欺骗检测算法的仿真结果;

图7(a)为本发明方法实施例中5m微小突变的阶跃式欺骗下新息速率抗差估计欺骗检测算法的仿真结果;

图7(b)为本发明方法实施例中0.05m/s缓慢增长的斜坡式欺骗下新息速率抗差估计欺骗检测算法的仿真结果;

图8(a)为本发明方法实施例中微小突变的阶跃式欺骗下新息速率抗差估计欺骗检测算法和微小突变的阶跃式欺骗检测算法的仿真结果对比示意图;

图8(b)为本发明方法实施例中微小突变的阶跃式欺骗对位置误差的影响效果示意图;

图9(a)为本发明方法实施例中缓慢增长的斜坡式欺骗下新息速率抗差估计欺骗检测算法和缓慢增长的斜坡式欺骗检测算法的仿真结果对比示意图;

图9(b)为本发明方法实施例中缓慢增长的斜坡式欺骗对位置误差的影响效果示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明了,以下结合附图及实施例,对本发明进行进一步详细说明。

方法实施例:

无人机在飞行过程中,当判断到存在第三方施加的欺骗干扰时,需要快速、准确地进行验证。本发明首先建立GNSS/INS组合导航系统的欺骗干扰模型,并对欺骗干扰带来的影响进行分析。

为了构建GNSS/INS组合导航系统欺骗干扰仿真环境,建立GNSS测量层面的欺骗干扰模型,先对复杂信号层面的欺骗模型进行分析。如图1所示,表示目标接收机受到的欺骗干扰。

第i颗卫星R

R

式中,c是光速,τ

对于真实信号,传播时延

考虑到接收机噪声

对于欺骗信号,传播时延

式中,

式中,▽τ

由此,根据第i路通道的欺骗伪距测量值和真实伪距测量值分别为

式中,t

基于GNSS/INS的紧组合导航系统使用GNSS伪距和伪距率作为输入,在闭环校正中,每次滤波迭代将估计的位置、速度和姿态误差反馈给INS处理器,对INS的解进行校正。误差状态扩展卡尔曼滤波(Extended Kalman Filter,EKF)的17维状态向量X表示为:

其中,

式中,

式中,

假设在k时刻发生第i次测量值伪距附加值为μ的欺骗干扰,其向量Λ表示为:

Λ=[0 0 … μ … 0]

除第i个元素外,其余元素均为0。根据卡尔曼滤波理论,可得:

其中,

式中,Φ

r

从公式中可以看出,欺骗新息r

Δr=-H

由于欺骗新息减少了Δr,那么其它新息值增加了|Δr|,因此欺骗干扰会造成真实新息与欺骗新息之间的差异,从而降低了欺骗干扰检测性能。对于微小突变和缓慢增长的欺骗干扰,其影响效果是不明显且较为隐蔽的,但因为检测有时间延迟,在这段检测时间内,新息也受到了一定的影响,最终影响GNSS/INS组合导航系统。

为此,本发明提供的基于新息速率优化和抗差估计的GNSS/INS紧组合欺骗检测方法如图3所示,包括如下步骤:

1、在k时刻,获取由GNSS观测值和INS预测值作差的观测向量Z

然后,根据新息向量r

式中,

2、采用卡尔曼滤波实时更新,更新过程包括时间更新和量测更新。

定义状态向量为

卡尔曼滤波实时更新过程中,需要考虑噪声,因此定义系统模型为:

其中,定义改进的新息速率是与时间相关的随机过程,β的取值不同,噪声的值也会有所不同。β为相关系数,一般取值为0.5~0.9;ε为噪声,一般取10

定义观测模型为:

式中,υ

将步骤1得到的归一化新息向量ω

在时间更新阶段,通过对GNSS/INS紧组合系统的状态向量和协方差进行先验估计,得到k时刻的先验估计状态向量

式中,符号“∧”表示估计值,上标“-”表示先验估计,“+”表示后验估计,上标“T”表示矩阵的转置。

在量测更新阶段,将归一化新息向量ω

根据观测矩阵H

根据观测矩阵H

还计算后验估计协方差向量

式中,R

根据后验估计状态向量

式中,C=[0 1 0 0]。

3、根据改进的新息速率,判断GNSS测量值中是否存在欺骗。

假设欺骗不存在的情况下,原假设为

则对应

式中,N

因此,欺骗干扰的检测标准为:

式中,

4、根据新息速率,利用IGG-3等价权函数,计算等价权矩阵W。等价权矩阵可表示为:W=diag(w

式中,

5、根据等价权矩阵W计算新的增益矩阵K

K

采用新的增益矩阵K

根据前面对欺骗干扰造成影响进行分析的过程判断出:欺骗新息与真实新息之间存在Δr,是导致欺骗新息下降的直接原因。因此降低Δr是提高欺骗检测能力的有效途径。通过自适应调整等价权来降低欺骗干扰的权重,便是通过将计算Δr的增益矩阵K

Δr=-H

接下来以缓慢增长的斜坡式欺骗和微小突变的阶跃式欺骗为例,对本发明的效果进行验证。通过构建GNSS/INS组合导航系统的欺骗干扰模型,模拟缓慢增长的斜坡式欺骗和微小突变的阶跃式欺骗,然后采用本发明进行欺骗检测,以验证本发明的检测效果。

1、根据欺骗类型选择欺骗检测算法,计算改进的归一化新息。若为缓慢增长的斜坡式欺骗,则选择缓慢增长的斜坡式欺骗检测算法;若为微小突变的阶跃式欺骗,则选择微小突变的阶跃式欺骗检测算法。

缓慢增长的斜坡式欺骗检测算法中,假设施加缓慢增长的伪距附加值为A的斜坡式欺骗干扰,进行第i次测量时,改进归一化新息

式中,A为缓慢增长的伪距附加值,

微小突变的阶跃式欺骗检测算法中,假设施加微小突变的伪距附加值为B的阶跃式欺骗干扰,进行第i次测量时,改进归一化新息

式中,B为微小突变的伪距附加值,

2、在缓慢增长的斜坡式欺骗检测算法中,通过卡尔曼滤波器实时更新改进归一化新息

以微小突变的阶跃式欺骗检测算法为例,对更新过程进行说明:

定义状态向量为

定义系统模型为:

定义观测模型为:

在时间更新阶段,得到k时刻的先验估计状态向量

在量测更新阶段,将归一化新息向量

根据观测矩阵H

根据观测矩阵H

还计算后验估计协方差向量

然后,根据后验估计状态向量

3、根据欺骗干扰的检测标准进行欺骗检测。

在新息速率满足

在缓慢增长的斜坡式欺骗检测算法中,采用同样的方法,计算缓慢增长新息速率

为了验证伪距附加值优化的欺骗检测算法的有益效果,对比新息速率欺骗检测算法(M1)、新息速率抗差估计欺骗检测算法(M2)、微小突变的阶跃式欺骗检测算法(M3)和缓慢增长的斜坡式欺骗检测算法(M4)四种算法。

设计4个GNSS/INS紧组合欺骗干扰检测实验场景,以比较不同算法的检测能力。(1)在3路通道受到相同伪距附加值的阶跃式欺骗和斜坡式欺骗时,比较M1与M2的检测能力;(2)在1路通道受到不同伪距附加值的阶跃式欺骗和斜坡式欺骗时,验证M2的检测能力;(3)在1路通道受到微小突变的阶跃式欺骗时,比较M2与M3的检测能力;(4)在1路通道受到缓慢增长的斜坡式欺骗时,比较M2与M4的检测能力。

1、参数设置。GNSS/INS紧组合的实验参数设置包括:模拟的飞行轨迹数据、欺骗场景和GNSS、IMU传感器的参数,分别如图4和表1、2所示。

表1欺骗场景设置

表2仿真参数

2、结果分析

(1)场景1。M2能够抑制新息速率偏离正常值的效果,即正常通道不受欺骗通道的影响而异常偏离,比较M2与M1的检测能力。

阶跃式欺骗仿真结果如图5所示,在350s至550s时间内,对通道1、2和3施加伪距偏差为30m的阶跃式欺骗。数据表明:(1)图5(a)是M1仿真结果图,显示了通道1、2和3的检测时间分别为10s、9s和10s,但其余5个通道均受到了不同程度的欺骗干扰影响,导致其对应新息速率偏离正常值,出现误警情况。(2)图5(b)是M2仿真结果图,显示了通道1、2和3的检测时间分别为9s、8s和9s,相比图5(a)受欺骗3路通道的检测时间各缩短了1s、2s和1s,且其余通道新息速率正常。由此,可以看出在施加相对较大伪距偏差的阶跃式欺骗时,M2相对M1的检测效率提升不明显,但M2能够抑制新息速率偏离正常值,增加了系统的容错率和鲁棒性。

斜坡式欺骗仿真结果如图6所示,在350s至550s时间内,对通道1、2和3施加斜率为0.3m/s的斜坡式欺骗。数据表明:(1)图6(a)是M1仿真结果,显示了通道1、2和3的检测时间分别为59s、50s和64s,与阶跃式欺骗类似,其余通道均受到了不同程度的欺骗干扰影响,出现误警情况。(2)图6(b)是M2仿真结果,显示了通道1、2和3的检测时间分别为47s、43s和49s,相比图6(a)受欺骗3路通道的检测时间各缩短了12s、7s和15s,且其余通道新息速率正常。

结合图5和图6数据得出:受到阶跃式欺骗时,M2比M1的检测时间分别缩短了10%、22.2%和10%,平均缩短时间了35.5%。受到斜坡式欺骗时,M2比M1的检测时间分别缩短了20.3%、14%和64.2%,平均缩短时间了32.8%。可以得出,M2比M1的平均检测效率提升了三分之一左右。说明M2的抗差估计很好的抑制了正常通道的新息速率偏离正常值,降低了误警率,且检测效率和检测性能相对M1更优。

(2)场景2。场景2设置了不同欺骗类型对应不同伪距附加值的场景,验证M2对微小突变(5m)的阶跃式欺骗或缓慢增长型(0.05m/s)的斜坡式欺骗有局限性,即检测时间过长甚至检测不敏感。

其仿真结果如图7所示,数据表明:(1)图7(a)是M2阶跃式欺骗仿真结果,在350s至550s时间内,对通道3施加伪距偏差分别为40m、30m、20m、10m和5m的阶跃式欺骗;显示了施加对应伪距偏差检测时间分别为4s、6s、10s、25s,其中5m的阶跃式欺骗检测无效。(2)图7(b)是M2斜坡式欺骗仿真结果,在350s至550s时间内,对通道3施加伪距偏差分别为0.4m/s、0.3m/s、0.2m/s、0.1m/s和0.05m/s的斜坡式欺骗;显示了施加对应附加值的斜坡式欺骗检测时间分别为36s、44s、57s、92s和157s。由此可知,验证了M2随着阶跃式施加伪距偏差(由40m至5m)和斜坡式施加斜率(由0.4m/s至0.05m/s)的减小,检测时间增加甚至检测不敏感。其中,对斜坡式欺骗(0.05m/s)最大检测时间为157s;对微小突变欺骗(5m)检测无效。

(3)场景3。针对场景2中M2对微小突变的阶跃式欺骗检测的局限性,设置了对微小突变的阶跃式欺骗场景3,比较M2与M3的检测能力。

其仿真结果如图8所示,在350s至550s时间内,对通道3施加伪距偏差为10m和5m的阶跃式欺骗。数据表明:(1)图8(a)是M2与M3对微小突变欺骗仿真结果,显示了对10m的伪距偏差,M2和M3检测时间分别为25s和9s;对5m的伪距偏差,M2检测无效,而M3检测有效其检测时间为22s,由此可知,在微小突变的欺骗检测中,M3比M2更加敏感。(2)图8(b)是微小突变欺骗对位置误差的影响,在施加10m和5m欺骗时,北向最大误差分别为1.041m和1.042m,东向最大误差分别为1.305m和1.398m,高度误差有微小变化,说明微小突变的欺骗很隐蔽且误差满足精度要求,一般很难检测出来,但对于高精度定位的应用危害很大,如导弹精确制导、智能驾驶和无人机等。

结合场景2和场景3可以得出:通道3受到10m的欺骗干扰时,M3比M2的检测时间缩短了64%;受到5m的欺骗干扰时,M3检测有效,M2检测不敏感。由此得出,M3可以检测隐蔽性强的微小突变的阶跃式欺骗。

(4)场景4。针对场景2中M2对缓慢增长的斜坡式欺骗检测的局限性,设置了对缓慢增长的斜坡式欺骗场景4,比较M2与M4的检测能力。

其仿真结果如图9所示,在350s至550s时间内,对通道3施加斜率为0.1m/s和0.05m/s的斜坡式欺骗。数据表明:(1)图9(a)是M2与M4对缓慢增长欺骗仿真结果,显示了对0.1m/s的斜坡式欺骗,M2和M4检测时间分别为92s和57s;对0.05m/s,M2h和M3检测时间分别为147s和86s,由此可知,在缓慢增长的斜坡式欺骗检测中,M3比M2检测效率更高。(2)图9(b)是缓慢增长欺骗对位置误差的影响,在施加0.1m/s和0.05m/s欺骗时,北向最大误差分别为1.213m和1.173m,东向最大误差分别为1.316m和1.344m,高度误差有微小变化,说明缓慢增长的斜率越小,达到最大误差的时间越长。

结合场景2和场景4可以得出:通道3受到缓慢增长0.1m/s和0.05m/s的欺骗时,M4比M2的检测时间分别缩短了38%和41.5%,平均缩短39.8%。对此,验证了M4可以检测隐蔽性强、危害性大的缓慢增长的斜坡式欺骗。

去获取专利,查看全文>

相似文献

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

客服邮箱:kefu@zhangqiaokeyan.com

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

  • 服务号