太阳城集团

  • / 14
  • 下载费用:30 金币  

一种GNSS和MIMU组合导航中微陀螺G敏感性误差补偿方法.pdf

摘要
申请专利号:

CN201410158231.5

申请日:

2014.04.18

公开号:

CN103954298A

公开日:

2014.07.30

当前法律状态:

授权

有效性:

有权

法律详情: 授权|||实质审查的生效IPC(主分类):G01C 25/00申请日:20140418|||公开
IPC分类号: G01C25/00 主分类号: G01C25/00
申请人: 中国人民解放军国防科学技术大学
发明人: 胡小平; 何晓峰; 范晨; 罗兵; 唐康华; 江明明; 王安成; 李涛; 练军想; 张礼廉
地址: 410073 湖南省长沙市砚瓦池正街47号中国人民解放军国防科学技术大学三院
优先权:
专利代理机构: 湖南兆弘专利事务所 43008 代理人: 周长清
PDF完整版下载: PDF下载
法律状态
申请(专利)号:

CN201410158231.5

授权太阳城集团号:

||||||

法律状态太阳城集团日:

2017.03.01|||2014.08.27|||2014.07.30

法律状态类型:

太阳城集团授权|||实质审查的生效|||公开

摘要

本发明公开了一种GNSS和MIMU组合导航中微陀螺g敏感性误差补偿方法,其步骤为:(1)建立微陀螺的g敏感性误差方程;(2)建立含g敏感性误差的GNSS/MIMU组合导航系统方程;(3)采用非线性全局可观性分析方法,对含g敏感性误差的GNSS/MIMU组合导航系统进行全局可观性分析,设计该组合系统全局可观的运动路径;(4)依据步骤(1)中得到的微陀螺的g敏感性误差方程,构建相应的卡尔曼滤波器,并按照步骤(3)得到的运动路径运动,实现微陀螺g敏感性误差的补偿。本发明具有原理简单、操作简便、适用范围广等优点。

权利要求书

权利要求书
1.  一种GNSS和MIMU组合导航中微陀螺g敏感性误差补偿方法,其特征在于,步骤为:
(1)建立微陀螺的g敏感性误差方程:
Go=Ggfb
式中,Go表示微陀螺的g敏感性误差,Gg表示g敏感性误差的系数矩阵,fb表示加速度计测得的比力太阳城集团;
(2)建立含g敏感性误差的GNSS/MIMU组合导航系统方程;
(3)采用非线性全局可观性分析方法,对含g敏感性误差的GNSS/MIMU组合导航系统进行全局可观性分析,设计该组合系统全局可观的运动路径;
(4)依据步骤(1)中得到的微陀螺的g敏感性误差方程,构建相应的卡尔曼滤波器,并按照步骤(3)得到的运动路径运动,实现微陀螺g敏感性误差的补偿。

2.  根据权利要求1所述的GNSS和MIMU组合导航中微陀螺g敏感性误差补偿方法,其特征在于,在所述步骤(2)中,选取GNSS/MIMU组合系统的状态为位置、速度、姿态、陀螺零偏、加速度计零偏和g敏感性误差系数,则系统状态方程为:
r·e=ve]]>
v·e=Cbe(fb-δfb)-2ωiee×vee+gle]]>
C·be=Cbe(ωebb×)]]>
ωebb=ωibb-δωb-Cebωiee]]>
式中:×表示矢量积运算;re和ve分别是载体在e系的位置和对地速度;为载体坐标系b系到e系的姿态变换矩阵;为在e系中表示的地球自转角速度;为在e系中表示的重力加速度;为在b系中表示的载体相对地球的转动角速度,为其反对称阵;为陀螺测量的载体相对惯性系的转动角速度;δωb和δfb分别为陀螺和加速度计的测量误差;
考虑g敏感性误差后,微陀螺的误差模型可分为两部分,一部分为与载体运动无关的误差项,由陀螺的固定零偏和随机噪声组成,记为另一部分则为g敏感性误差,记为δω2b=Go=Ggfb;]]>即:
δωb=δω1b+δω2b]]>
满足关系:
δω·1b=0]]>
δω·2b=Ggf·b]]>
G·g=03×3]]>
若以卫星接收机输出的位置和速度作为观测量,则系统的观测方程为:
Z=z1z2=reve]]>
满足关系:
z2=z·1.]]>

3.  根据权利要求2所述的GNSS和MIMU组合导航中微陀螺g敏感性误差补偿方法,其特征在于,经步骤(3)的全局可观性分析后,对于含g敏感性误差的组合导航系统,以位置和速度为观测量,若载体的运动满足以下两个条件,则系统是全局可观的:
(a)存在一个直线段在这个直线段上存在两个时刻t1≠t2,使比力导数和线性不相关;
(b)存在两个直线运动段,且在这两个路段中存在三个不同时刻t1、t2与t3,使得R(A)=R(A,B)=3或存在直线运动段,使得在此路段中各分量都存在不为零的时刻。

4.  根据权利要求1或2或3所述的GNSS和MIMU组合导航中微陀螺g敏感性误差补偿方法,其特征在于,所述卡尔曼滤波器的状态方程为:
X·(t)=F(t)X(t)+W(t)]]>
式中,系统状态为:
X=(δre,δve,ϵe,δω1b,δfb,G1,G2,G3)T]]>
状态转移矩阵为:
F(t)=03×3I3×303×303×303×303×303×303×3Ne-2Ωiee-Fe03×3Cbe03×303×303×303×303×3-ΩieeCbe03×3F1eF2eF3e03×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×3]]>
系统噪声阵为:
W(t)=[wr wv wε 01×3 01×3 01×3 01×3 01×3]T
观测方程为:
Z(t)=H(t)X(t)+N(t)
式中,系统的观测矩阵H(t)=I3×303×303×303×303×303×303×303×303×3I3×303×303×303×303×303×303×3,]]>观测噪声阵为N(t)=[Nr Nv 01×3 01×3 01×3 01×3 01×3 01×3]T。

说明书

说明书一种GNSS和MIMU组合导航中微陀螺g敏感性误差补偿方法
技术领域
本发明主要涉及到组合导航技术领域,特指一种适用于高动态条件下GNSS/MIMU组合导航中微陀螺的g敏感性误差误差补偿方法。
背景技术
微惯性测量单元(MIMU)与卫星接收机构成的组合导航系统具有体积小、成本低等优势,应用潜力巨大,成为当前导航领域的一个研究热点。MIMU一般精度较差,尤其是其中的微机械陀螺,多数还处于角速率级的水平。当前多数微机械陀螺都属于哥氏振动陀螺,振动或线加速运动对其性能影响很大。g敏感性(g-sensitivity)误差系数是微陀螺重要的性能参数之一,它表征了载体线加速运动所引起的微陀螺性能恶化程度。
目前,市场上微机械陀螺的g误差系数标定的典型值为0.1°/s/g左右,例如美国Melexis公司的MLX90609和Silicon Sensing公司的CRG20-01,这意味着1g加速度的变化将引起360°/h的零偏变化,在更高动态的应用场合中,微机械陀螺的性能会变得很差甚至不可用。因此高动态下微陀螺g敏感性误差对MIMU/GNSS组合导航性能的影响,是一个值得关注的问题。
总而言之,目前的MEMS精度较低,动态性能较差,难以适应高动态条件下的导航需求,尤其是其表现出显著地g敏感性误差,极大地制约了MEMS的应用,因此需要对高动态条件下微陀螺g敏感性误差分析与补偿技术进行研究。
发明内容
本发明要解决的技术问题就在于:针对现有技术存在的技术问题,本发明提供一种原理简单、操作简便、适用范围广的GNSS和MIMU组合导航中微陀螺g敏感性误差补偿方法。
为解决上述技术问题,本发明采用以下技术方案:
一种GNSS和MIMU组合导航中微陀螺g敏感性误差补偿方法,其步骤为:
(1)建立微陀螺的g敏感性误差方程:
Go=Ggfb
式中,Go表示微陀螺的g敏感性误差,Gg表示g敏感性误差的系数矩阵,fb表示加速度计测得的比力太阳城集团;
(2)建立含g敏感性误差的GNSS/MIMU组合导航系统方程;
(3)采用非线性全局可观性分析方法,对含g敏感性误差的GNSS/MIMU组合导航系统进行全局可观性分析,设计该组合系统全局可观的运动路径;
(4)依据步骤(1)中得到的微陀螺的g敏感性误差方程,构建相应的卡尔曼滤波器,并按照步骤(3)得到的运动路径运动,实现微陀螺g敏感性误差的补偿。
作为本发明的进一步改进:在所述步骤(2)中,选取GNSS/MIMU组合系统的状态为位置、速度、姿态、陀螺零偏、加速度计零偏和g敏感性误差系数,则系统状态方程为:
r·e=ve]]>
v·e=Cbe(fb-δfb)-2ωiee×vee+gle]]>
C·be=Cbe(ωebb×)]]>
ωebb=ωibb-δωb-Cebωiee]]>
式中:×表示矢量积运算;re和ve分别是载体在e系的位置和对地速度;为载体坐标系b系到e系的姿态变换矩阵;为在e系中表示的地球自转角速度;为在e系中表示的重力加速度;为在b系中表示的载体相对地球的转动角速度,为其反对称阵;为陀螺测量的载体相对惯性系的转动角速度;δωb和δfb分别为陀螺和加速度计的测量误差;
考虑g敏感性误差后,微陀螺的误差模型可分为两部分,一部分为与载体运动无关的误差项,由陀螺的固定零偏和随机噪声组成,记为另一部分则为g敏感性误差,记为δω2b=Go=Ggfb;]]>即:
δωb=δω1b+δω2b]]>
满足关系:
δω·1b=0]]>
δω·2b=Ggf·b]]>
G·g=03×3]]>
若以卫星接收机输出的位置和速度作为观测量,则系统的观测方程为:
Z=z1z2=reve]]>
满足关系:
z2=z·1.]]>
作为本发明的进一步改进:经步骤(3)的全局可观性分析后,对于含g敏感性误差的组合导航系统,以位置和速度为观测量,若载体的运动满足以下两个条件,则系统是全局可观的:
(a)存在一个直线段在这个直线段上存在两个时刻t1≠t2,使比力导数和线性不相关;
(b)存在两个直线运动段,且在这两个路段中存在三个不同时刻t1、t2与t3,使得R(A)=R(A,B)=3或存在直线运动段,使得在此路段中各分量都存在不为零的时刻。
作为本发明的进一步改进:所述卡尔曼滤波器的状态方程为:
X·(t)=F(t)X(t)+W(t)]]>
式中,系统状态为:
X=(δre,δve,ϵe,δω1b,δfb,G1,G2,G3)T]]>
状态转移矩阵为:
F(t)=03×3I3×303×303×303×303×303×303×3Ne-2Ωiee-Fe03×3Cbe03×303×303×303×303×3-ΩieeCbe03×3F1eF2eF3e03×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×3]]>
系统噪声阵为:
W(t)=[wr wv wε 01×3 01×3 01×3 01×3 01×3]T
观测方程为:
Z(t)=H(t)X(t)+N(t)
式中,系统的观测矩阵H(t)=I3×303×303×303×303×303×303×303×303×3I3×303×303×303×303×303×303×3,]]>观测噪声阵为N(t)=[Nr Nv 01×3 01×3 01×3 01×3 01×3 01×3]T。
与现有技术相比,本发明的优点在于:本发明的GNSS和MIMU组合导航中微陀螺g敏感性误差补偿方法,原理简单、操作简便、适用范围广,对于解决高动态条件下微陀螺g敏感性误差提供了理论依据,有效地提高了高动态条件下GNSS/MIMU组合系统的精度,对于高动态条件下的GNSS/MIMU组合导航系统设计具有重要的理论指导和工程实践意义。
附图说明
图1是本发明方法的流程示意图。
图2是本发明在具体应用实例中载体全局可观的运动路径示意图。
具体实施方式
以下将结合说明书附图和具体实施例对本发明做进一步详细说明。
如图1所示,本发明的一种GNSS和MIMU组合导航中微陀螺g敏感性误差补偿方法,为:首先,建立了含微陀螺g敏感性误差的GNSS/MIMU组合导航模型;然后,采用非线性全局可观性法对系统进行可观性分析,设计相应的全局可观的路径;最后,依据微陀螺g敏感性误差的模型,设计用于补偿g敏感性误差的卡尔曼滤波器。
在具体应用实例中,本发明方法的具体流程为:
(1)微陀螺g敏感性误差方程为:
Go=Ggfb    (1)
式中,Go表示微陀螺的g敏感性误差,Gg表示g敏感性误差的系数矩阵,fb表示加速度计测得的比力太阳城集团。
等价变换:
Go=G11G12G13G21G22G23G31G32G33fxbfybfzb=G1G2G3Tfb---(2)]]>
式中,G1、G2和G3分别表示x、y和z三个轴向陀螺的g敏感性误差系数矢量。
为了方便设计包含g敏感性误差组合滤波器,则需将g敏感性误差系数作为误差状态。因此利用矩阵变换,将其写成矢量形式:
G1G2G3Tfb=G103×103×1Tfb+03×1G203×1Tfb+03×103×1G3Tfb=fb03×103×1TG1+03×1fb03×1TG2+03×103×1fbTG3---(3)]]>

F1b=fb03×103×1TF2b=03×1fb03×1TF3b=03×103×1fbT---(4)]]>
则可得到:
Go=F1bG1+F2bG2+F3bG3---(5)]]>
(2)选取GNSS/MIMU组合系统的状态为位置、速度、姿态、陀螺零偏、加速度计零偏和g敏感性误差系数,则系统状态方程为:
r·e=ve---(6)]]>
v·e=Cbe(fb-δfb)-2ωiee×vee+gle---(7)]]>
C·be=Cbe(ωebb×)---(8)]]>
ωebb=ωibb-δωb-Cebωiee---(9)]]>
式中:×表示矢量积运算;re和ve分别是载体在e系的位置和对地速度;为载体坐标系b系到e系的姿态变换矩阵;为在e系中表示的地球自转角速度;为在e系中表示的重力加速度;为在b系中表示的载体相对地球的转动角速度,为其反对称阵;为陀螺测量的载体相对惯性系的转动角速度;δωb和δfb分别为陀螺和加速度计的测量误差。
考虑g敏感性误差后,微陀螺的误差模型可分为两部分,一部分为与载体运动无关的误差项,由陀螺的固定零偏和随机噪声组成,记为另一部分则为g敏感性误差,记为δω2b=Go=Ggfb.]]>即:
δωb=δω1b+δω2b---(10)]]>
满足关系:
δω·1b=0---(11)]]>
δω·2b=Ggf·b---(12)]]>
其中,g敏感性误差满足关系
若以卫星接收机输出的位置和速度作为观测量,则系统的观测方程为:
Z=z1z2=reve---(13)]]>
满足关系:
z2=z·1---(14)]]>
(3)采用非线性全局可观性的方法,对含g敏感性误差的GNSS/MIMU组合系统进行了全局可观性分析,得出了系统全局可观的充分条件,对实现微陀螺g敏感性误差的补偿提供了理论依据。该步骤的具体分析过程如下。
系统已知的输入太阳城集团有:fb、和Z,如果在有限太阳城集团内,根据上述输入太阳城集团能够唯一地确定系统的初始状态,则系统是可观的,否则系统不可观。
为了更清楚描述推导过程,先说明以下两个引理:
引理1:对于任意两个线性不相关的向量,如果它们在任意两个坐标系的坐标均已知,则这两个坐标系之间的姿态转换矩阵可以被唯一地确定。
引理2:设A是m×n矩阵,X=(x1,x2,…,xn)T,B=(b1,b2,…,bm)T,其构成的非齐次线性方程组AX=B有唯一解的充要条件是R(A)=R(A,B)=n。
(3.1)加速度计零偏的可观性;
设在载体运行过程中,存在一段直线运动,姿态不发生变化,即:
C·be=0---(15)]]>
由系统的状态方程和观测方程可得:
z·2=Cbe(fb-δfb)-2ωiee×z2+gle---(16)]]>
对上式两边同时求导,并且整理得:
Cbef·b=z··2+2ωiee×z·2---(17)]]>
式中,地球自转角速度已知,故已知;z2为速度观测量,故其各阶导数均已知;fb为加速度计测得的比力太阳城集团,其一阶导数也已知,因此上式中除外,其余状态均为已知。
根据引理1,如果此路段上存在线性不相关的和t1≠t2,则可被唯一地确定,故在此路段上的姿态矩阵已知。
整理得:
δfb=fb-Cbe(z·2+2ωiee×z2-gle)---(18)]]>
根据上式等号右侧各项均已知,故δfb可以被唯一确定。由于δfb为固定常数,因此加速度计零偏是可观的。
(3.2)g敏感性误差系数的可观性;
根据组合导航的系统方程可得:
ωebb=ωibb-Cebωiee-δω1-Ggfb---(19)]]>
对上式两边同时求导数,设则可写为:
Ggf·b=ω·ibb---(20)]]>
微陀螺g敏感性误差是由于载体的大加速度产生,故对于上式,陀螺测量值的一阶导数为因变量,而加速度计测量值的一阶导数为自变量,因此g敏感性系数矩阵Gg是否存在唯一解取决于则上式的求解等价于求解一个非齐次线性方程组。
根据步骤1可得:
G1G2G3Tf·b=ω·ibb---(21)]]>
展开后,写成分量形式:
G11f·xb+G12f·yb+G13f·zb=ω·ibxb---(22)]]>
G21f·xb+G22f·yb+G23f·zb=ω·ibyb---(23)]]>
G31f·xb+G32f·yb+G33f·zb=ω·ibzb---(24)]]>
对于上述,与均可由测量值得到,仅含有三个未知系数G11、G12与G13,因此需得到三个不同时刻的和和其与三个未知系数构成非齐次线性方程组:
G11f·xb(t1)+G12f·yb(t1)+G13f·zb(t1)=ω·ibxb(t1)---(25)]]>
G11f·xb(t2)+G12f·yb(t2)+G13f·zb(t2)=ω·ibxb(t2)---(26)]]>
C11f·xb(t3)+G12f·yb(t3)+G13f·zb(t3)=ω·ibxb(t3)---(27)]]>
记A=(fb(t1) fb(t2) fb(t3))T,则对于上述方程,根据引理2可知,若载体在整个运动过程中,存在两个直线路段,且在这两个路段上存在三个不同时刻t1、t2与t3,使得R(A)=R(A,B)=3,则G11、G12与G13就可被唯一确定。
同理可知,在满足上述条件的同时,其余两个轴的g敏感性误差系数也都可观。
(3.3)的可观性分析;
陀螺的此部分偏差与载体的运动无关,因此单独分开讨论。在直线段满足关系:
δω1b=ωibb-Cebωiee-Ggfb---(28)]]>
易知等号右侧均已知,可以被唯一确定,所以在此路段上可观。
(3.4)姿态的可观性;
由系统的状态方程和观测方程可得:
z·2=Cbe(fb-δfb)-2ωiee×z2+gle---(29)]]>
对上式两边同时求导,整理得:
z··2=Cbef·b+C·be(fb-δfb)-2ωiee×z·2+g·le---(30)]]>
设在载体运行过程中,存在一段直线运动,姿态不发生变化,则:
z··2=Cbef·b-2ωiee×z·2+g·le---(31)]]>
y=z··2+2ωiee×z·2-g·le,]]>则:
y=Cbef·b---(32)]]>
定义惯性坐标系i系,使之与初始时刻t0的e系重合,则当前时刻t的姿态矩阵和初始时刻的姿态矩阵的关系可表示为:
Cbe=Ce(t0)e(t)Cb(t0)e(t0)Cb(t)b(t0)=Ci(t0)e(t)Cb(t0)e(t0)Cb(t)b(t0)=Ci(t)e(t)Cb(t0)e(t0)Cb(t)b(t0)=CieCbe(t0)Cbb(t0)---(33)]]>
其中上标b(t0)表示初始时刻的载体坐标系,表示当前时刻的载体坐标系到初始时刻载体坐标系的坐标转换矩阵,表示惯性系到地固坐标系的坐标转换矩阵。
两边同时左乘得:
Cbe(t0)Cbb(t0)f·b=Ceiy---(34)]]>
是地球自转角速度和太阳城集团的函数,由于地球自转角速度已知,故已知;代表t0时刻的b系与t时刻的b系之间的坐标转换矩阵,此矩阵可通过陀螺输出的角速度(补偿固定零偏与g敏感性误差后)积分得到,故已知;忽略位置误差,故当地的重力加速度矢量已知;可通过加速度计测得的比力求一阶导得到,z2为速度观测量可直接得到,因此上式各项除外均已知的。
令则根据引理1,如果存在和t1≠t2,且与线性不相关,则可以被唯一地确定,因此姿态可观,能够唯一被确定。
(3.5)速度与位置的可观性;
由于系统的观测量为卫星提供的速度和位置,则可知两者一定可观。
综上分析过程可知,对于含g敏感性误差的组合导航系统,以位置和速度为观测量,若载体的运动满足如下两个条件,则系统是全局可观的。
(a)存在一个直线段在这个直线段上存在两个时刻t1≠t2,使比力导数和线性不相关;
(b)存在两个直线运动段,且在这两个路段中存在三个不同时刻t1、t2与t3,使得R(A)=R(A,B)=3或存在直线运动段,使得在此路段中各分量都存在不为零的时刻。
依据上述的全局可观的充分条件,设计了如图2所示的载体运动轨迹示意图。
(4)根据微陀螺g敏感性误差模型,设计了相应的卡尔曼滤波器,以实现g敏感性误差的补偿。
考虑g敏感性误差后的失准角误差方程为:
ϵ·e=-Ωieeϵe+Cbeδω1b+CbeF1bG1+CbeF2bG2+CbeF3bG3---(35)]]>
等价变换:
ϵ·e=-Ωieeϵe+Cbeδω1b+F1eG1+F2eG2+F3eG3---(36)]]>
在组合滤波器中,考虑微陀螺g敏感性误差,增加微陀螺3个轴的g敏感性误差系数作为新的状态变量,则将原有的15维滤波器扩充为24维滤波器,即:
X=(δre,δve,ϵe,δω1b,δfb,G1,G2,G3)T]]>
可得到组合系统的误差状态方程:
X·(t)=δr·eδv·eϵ·eδω·1bδf·bG·1G·2G·3=δve-Feϵe+Cbeδfb-2ωiee×δve+Neδre-Ωieeϵe+Cbeδω1b+F1eG1+F2eG2+F3eG300000---(37)]]>
在此基础上,可得到卡尔曼滤波器的状态方程为:
X·(t)=F(t)X(t)+W(t)---(38)]]>
其中,状态转移矩阵为:
F(t)=03×3I3×303×303×303×303×303×303×3Ne-2Ωiee-Fe03×3Cbe03×303×303×303×303×3-ΩieeCbe03×3F1eF2eF3e03×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×303×3---(39)]]>
系统噪声阵为:
W(t)=[wr wv wε 01×3 01×3 01×3 01×3 01×3]T    (40)
观测方程为:
Z(t)=H(t)X(t)+N(t)    (41)
式中,系统的观测矩阵H(t)=I3×303×303×303×303×303×303×303×303×3I3×303×303×303×303×303×303×3,]]>观测噪声阵N(t)=[Nr Nv 01×3 01×3 01×3 01×3 01×3 01×3]T。
上述卡尔曼滤波器方程,将g敏感性误差系数作为状态进行滤波估计,结合全局可观的充分条件,设计相应的全局可观路径,理论上能够实现g敏感性误差的补偿。
以上仅是本发明的优选实施方式,本发明的保护范围并不仅局限于上述实施例,凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出,对于本技术领域的普通技术人员来说,在不脱离本发明原理前提下的若干改进和润饰,应视为本发明的保护范围。

关 键 词:
一种 GNSS MIMU 组合 导航 陀螺 敏感性 误差 补偿 方法
  专利查询网所有资源均是用户自行上传分享,仅供网友学习交流,未经上传用户书面授权,请勿作他用。
太阳城集团本文
本文标题:一种GNSS和MIMU组合导航中微陀螺G敏感性误差补偿方法.pdf
链接地址:http://zh228.com/p-6140616.html
太阳城集团我们 - 网站声明 - 网站地图 - 资源地图 - 友情链接 - 网站客服客服 - 联系我们

copyright@ 2017-2018 zhuanlichaxun.net网站版权所有
经营许可证编号:粤ICP备17046363号-1 
 


收起
展开
葡京赌场|welcome document.write ('');