高性能PVT解算方法及导航方法与流程

文档序号:30561929发布日期:2022-06-29 03:57阅读:854来源:国知局
导航: X技术> 最新专利> 测量装置的制造及其应用技术
高性能PVT解算方法及导航方法与流程
高性能pvt解算方法及导航方法
技术领域
1.本发明属于导航技术领域,具体涉及一种高性能pvt解算方法及导航方法。


背景技术:

2.随着经济技术的发展和人们生活水平的提高,导航技术已经广泛应用于人们的生产和生活当中,给人们的生产和生活带来了无尽的便利。因此,保证导航过程的精度和实时性,就成为了研究人员的研究重点。
3.gnss导航接收机解算的主要目的是为了求解出接收机的位置(position)、速度(velocity)和时间(time),这个过程称作pvt解算。为了满足多系统、多频点的应用场景,gnss导航接收机需要具备多系统、多频点同时解算的特性。参与pvt解算的有效卫星数目和信号分量越多,就越能提高接收机定位的可靠性和定位精度;但是,相应也会增大pvt解算的运算量和耗时,增加系统负担;这将使得gnss导航接收机的应用很难得到本质的提升和突破。
4.目前,在传统的接收机的解算方法中,pvt解算过程通常包括以下几个步骤:(1)获取观测量数据和导航电文数据;(2)计算卫星状态信息(包括卫星位置、速度和加速度等);(3)最小二乘解算;(4)定位结果处理及输出。该过程的典型方法流程图如图1所示。
5.在上述传统解算的过程中,四个步骤顺序执行,每一个步骤的计算结果均要作为下一个步骤的输入。整个pvt求解运算过程的耗时主要集中在计算卫星状态信息和最小二乘解算两个部分;而其中计算卫星状态信息耗时更多,占整个pvt解算耗时的60~70%。
6.在多系统、多卫星、多频点的应用场景下,为了提高定位精度,同一颗卫星的不同频点往往会使用当前频点的导航电文数据参与解算,这就需要基于当前频点的导航电文数据再次计算该卫星的更高精度的卫星状态信息;这使得运算量成倍增加。另外,高频度解算时,高频的每一拍pvt解算都需要完成所有有效卫星的状态信息计算,这使得高频度解算的接收机压力更大;此外,现有技术在进行pvt解算的过程中,其计算流程繁复,效率也较低。因此,传统的pvt解算过程,解算耗时巨大,限制了高频度、低成本等应用需求。


技术实现要素:

7.本发明的目的之一在于提供一种精度较高、可靠性较好且解算耗时更低的高性能pvt解算方法。
8.本发明的目的之二在于提供一种包括了所述高性能pvt解算方法的导航方法。
9.本发明提供的这种高性能pvt解算方法,包括如下步骤:s1. 获取pvt解算所需要的输入数据;s2. 根据步骤s1获取的输入数据,采用第一线程,以设定的调度频度、卫星号周期性方式和时间间隔抽样方式,将所有卫星状态信息的计算工作分摊到不同的时隙,并逐个卫星进行卫星状态信息的计算;s3. 在步骤s2进行的同时,采用第二线程,结合二次运动模型理论,以二次曲线拟
合的方式进行卫星状态信息外推,并进行最小二乘解算;s4. 根据步骤s3得到的最小二乘解算结果,得到最终的定位结果,完成pvt解算。
10.采用第一线程进行计算的优先级低于采用第二线程进行计算的优先级。
11.步骤s2计算得到卫星状态信息的计算结果后,将计算结果更新到全局卫星状态信息缓存中,供步骤s3进行卫星状态信息外推使用;在对全局卫星状态信息缓存进行任何操作时,均采用互斥量对全局卫星状态信息缓存进行防护,从而保证卫星状态信息的完整性、操作的可靠性和安全性。
12.步骤s1所述的获取pvt解算所需要的输入数据,具体包括obs观测量数据和nav导航电文数据;其中,nav导航电文数据包括eph星历数据,satclk卫星钟差数据,iono电离层数据和tgd群延迟数据;然后将获取的输入数据按照设定的数据结构和格式进行组装。
13.步骤s2所述的设定的调度频度,具体为100hz。
14.步骤s2所述的卫星号周期性方式,具体包括如下步骤:总卫星数n为,其中n1为gps卫星的数量,n2为glo卫星的数量,n3为gal卫星的数量,n4为北斗卫星的数量;在设定的调度频度下,每一拍仅计算一颗卫星的状态信息;当计算完最后一颗卫星的状态信息后,再从第一颗卫星开始周期性的进行计算。
15.步骤s2所述的时间间隔抽样方式,具体包括如下步骤:已经计算过第i颗卫星的状态信息后,当再次轮到对第i颗卫星进行计算时,只有当前时刻与上次进行状态信息计算的时刻的时间间隔大于设定值时,才对第i颗卫星进行状态信息的计算;否则,直接跳过第i颗卫星的状态信息计算过程。
16.步骤s2所述的进行卫星状态信息的计算,具体包括如下步骤:a. 获取卫星的星历数据信息;b. 对全局卫星状态信息缓存中的当前卫星的状态信息有效标志进行判断:若当前卫星的状态信息有效标志无效,则进行后续步骤;若当前卫星的状态信息有效标志有效,则按照时间间隔抽样方式进行判断;c. 对当前卫星的星历信息进行检测:若当前卫星的星历信息有效,则进行后续步骤;若当前卫星的星历信息无效,则设置当前卫星的状态信息有效标志为无效,将结果保存到全局卫星状态信息缓存中,并等待下一拍的调度计算;d. 采用当前卫星的星历信息,计算当前时刻的卫星位置;e. 采用当前卫星的星历信息,计算当前时刻的卫星速度;f. 采用当前卫星的星历信息,计算当前时刻的卫星加速度;g. 计算完成后,设置当前卫星的状态信息有效标志为有效,并将计算结果进行保存,保存到全局卫星状态信息缓存中。
17.步骤a所述的卫星星历的数据信息,具体包括星历参考时间t
oe
、卫星轨道长半轴as的平方根值、轨道偏心率es、星历参考时间t
oe
时刻的轨道倾角i0、周内时为0时的轨道升交点赤经、轨道近地角距、星历参考时间t
oe
时刻的平近点角m0、平均运动角速度校正值、轨道倾角对时间的变化率、轨道升交点赤经对时间的变化率、升交点角距余
弦调和校正振幅c
uc
、升交点角距正弦调和校正振幅c
us
、轨道半径余弦调和校正振幅c
rc
、轨道半径正弦调和校正振幅c
rs
、轨道倾角余弦调和校正振幅c
ic
和轨道倾角正弦调和校正振幅c
is

18.所述的步骤d,具体包括如下步骤:d1. 计算归一化时间tk为,式中t为当前时刻;计算平均角速度n为,式中为地心引力常数;计算平近点角mk为,式中n为计算得到的平均角速度;d2. 以为初始值,采用算式进行迭代,直至精度满足设定的要求,得到最终的偏近点角e;ej为第j次计算得到的偏近点角;d3. 计算真近点角vk为;计算升交点角距为;d4. 采用如下算式计算摄动校正项:采用如下算式计算摄动校正项:采用如下算式计算摄动校正项:式中为升交点角距修正量;为卫星地心向径修正量;为卫星轨道倾角修正量;d5. 计算升交点角距uk为;计算卫星矢量长度rk为;计算轨道倾角ik为;d6. 采用如下算式计算卫星轨道平面的位置和升交点赤经:::式中为地球自转角速度;d7. 采用如下算式计算地心地固直角坐标系中的卫星坐标为,和。
19.所述的步骤e,具体为将地心地固直角坐标系中的卫星坐标对时间求导,从而得到当前时刻的卫星速度为:为:
式中为升交点赤经对时间的变化率;为轨道倾角ik对时间的变化率;为卫星轨道平面的位置对时间的变化率。
20.所述的步骤f,具体为将步骤e得到的当前时刻的卫星速度对时间求导,从而得到当前时刻的卫星加速度为:为:为:式中为第一系数,且;为第二系数,且;为轨道倾角ik对时间的二次求导;为卫星轨道平面的位置对时间的二次求导。
21.所述的步骤s3,具体包括如下步骤:a. 从全局卫星状态信息缓存中获取当前卫星的状态信息有效标志:若当前卫星的状态信息有效标志为无效,则本步骤结束;并获取下一颗卫星的状态信息有效标志;若当前卫星的状态信息有效标志为有效,则计算并得到当前卫星的外推时间间隔为,式中t
pos
为当前pvt解算时刻,ts为当前卫星的最近一次状态信息计算的时刻;b. 设置时间阈值t
th1
:若当前卫星的外推时间间隔大于t
th1
时,判定当前卫星异常,不再外推当前卫星的状态信息,且不再使用当前卫星参与后续的解算过程;若当前卫星的外推时间间隔小于或等于t
th1
时,采用如下二次曲线模型外推当前卫星在t
pos
时刻的位置和速度:时刻的位置和速度:时刻的位置和速度:时刻的位置和速度:时刻的位置和速度:
式中上标s表示为当前卫星s;为当前卫星s在ts时刻的位置;为当前卫星s在ts时刻的速度;为当前卫星s在ts时刻的加速度;为当前卫星s外推到t
pos
时刻的位置;为当前卫星s外推到t
pos
时刻的速度;c. 重复步骤a~b,直至遍历完成所有卫星;d. 将所有正常卫星的状态信息外推后,采用外推得到的状态信息进行后续的最小二乘解算。
22.本发明还提供了一种包括了所述高性能pvt解算方法的导航方法,包括如下步骤:(1)实时获取导航所需的数据信息;(2)采用上述的高性能pvt解算方法进行实时pvt解算,得到实时的定位结果;(3)根据步骤(2)得到的实时定位结果,完成最终的导航。
23.本发明提供的这种高性能pvt解算方法及导航方法,将运算量较大的卫星状态信息计算放在优先级较低的第一线程中,采用设定的调度频度、按卫星号周期性的方式、按时间间隔抽样方式执行,将原本需要在一拍内完成的所有卫星状态信息的计算工作,分摊到不同的时隙,并逐个卫星开始计算;把运算量较少的卫星状态信息外推和最小二乘解算放在优先级较高的第二线程中执行,结合二次运动模型理论,采用二次曲线拟合方式进行卫星状态信息外推,并采用系统配置的解算频度进行计算;对传统的pvt解算流程进行优化,将卫星状态信息计算过程进行剥离和优化;两个线程独立运算,两者之间进行卫星状态信息的交互,降低了系统负担,减少了解算耗时,提高了解算效率,而且定位精度较高、可靠性较好。
附图说明
24.图1为现有的pvt解算方法的方法流程示意图。
25.图2为本发明的解算方法的方法流程示意图。
26.图3为本发明的解算方法的时序示意图。
27.图4为本发明的导航方法的方法流程示意图。
具体实施方式
28.如图2所示为本发明的解算方法的方法流程示意图,其对应的时序图如图3所示:本发明提供的这种高性能pvt解算方法,包括如下步骤:s1. 获取pvt解算所需要的输入数据;具体包括obs观测量数据和nav导航电文数据;其中,nav导航电文数据包括eph星历数据,satclk卫星钟差数据,iono电离层数据和tgd群延迟数据等;然后将获取的输入数据按照设定的数据结构和格式进行组装;s2. 根据步骤s1获取的输入数据,采用第一线程,以设定的调度频度、卫星号周期性方式和时间间隔抽样方式,将所有卫星状态信息的计算工作分摊到不同的时隙,并逐个卫星进行卫星状态信息的计算;
具体实施时,设定的调度频度优选为100hz,即10ms为一拍,每一拍仅计算一颗卫星的状态信息;卫星号周期性方式具体为:总卫星数n为,其中n1为gps卫星的数量,n2为glo卫星的数量,n3为gal卫星的数量,n4为北斗卫星的数量;在设定的调度频度下,每一拍仅计算一颗卫星的状态信息;当计算完最后一颗卫星的状态信息后,再从第一颗卫星开始周期性的进行计算;具体实施时,n1为32,n2为28,n3为30,n4为37,n为127;时间间隔抽样方式具体为:已经计算过第i颗卫星的状态信息后,当再次轮到对第i颗卫星进行计算时,只有当前时刻与上次进行状态信息计算的时刻的时间间隔大于设定值(优选为5s)时,才对第i颗卫星进行状态信息的计算;否则,直接跳过第i颗卫星的状态信息计算过程;根据以上的设定的调度频度、卫星号周期性方式和时间间隔抽样方式,将原本需要在一拍内完成的所有卫星状态信息的计算工作,分摊到不同的时隙,并逐个卫星开始计算;进行卫星状态信息的计算,则具体包括如下步骤:a. 获取卫星的星历数据信息;具体包括星历参考时间t
oe
、卫星轨道长半轴as的平方根值、轨道偏心率es、星历参考时间t
oe
时刻的轨道倾角i0、周内时为0时的轨道升交点赤经、轨道近地角距、星历参考时间t
oe
时刻的平近点角m0、平均运动角速度校正值、轨道倾角对时间的变化率、轨道升交点赤经对时间的变化率、升交点角距余弦调和校正振幅c
uc
、升交点角距正弦调和校正振幅c
us
、轨道半径余弦调和校正振幅c
rc
、轨道半径正弦调和校正振幅c
rs
、轨道倾角余弦调和校正振幅c
ic
和轨道倾角正弦调和校正振幅c
is
;b. 对全局卫星状态信息缓存中的当前卫星的状态信息有效标志进行判断:若当前卫星的状态信息有效标志无效,则进行后续步骤;若当前卫星的状态信息有效标志有效,则按照时间间隔抽样方式进行判断;c. 对当前卫星的星历信息进行检测:若当前卫星的星历信息有效(包括星历已经收齐、星历健康、星历龄期有效等),则进行后续步骤;若当前卫星的星历信息无效,则设置当前卫星的状态信息有效标志为无效,将结果保存到全局卫星状态信息缓存中,并等待下一拍的调度计算;d. 采用当前卫星的星历信息,计算当前时刻的卫星位置;具体包括如下步骤:d1. 计算归一化时间tk为,式中t为当前时刻;计算平均角速度n为,式中为地心引力常数;计算平近点角mk为,式中n为计算得到的平均角速度;d2. 以为初始值,采用算式进行迭代,直至精度满足设定的要求,得到最终的偏近点角e;ej为第j次计算得到的偏近点角;
d3. 计算真近点角vk为;计算升交点角距为;d4. 采用如下算式计算摄动校正项:采用如下算式计算摄动校正项:采用如下算式计算摄动校正项:式中为升交点角距修正量;为卫星地心向径修正量;为卫星轨道倾角修正量;d5. 计算升交点角距uk为;计算卫星矢量长度rk为;计算轨道倾角ik为;d6. 采用如下算式计算卫星轨道平面的位置和升交点赤经:::式中为地球自转角速度;d7. 采用如下算式计算地心地固直角坐标系中的卫星坐标为,和;e. 采用当前卫星的星历信息,计算当前时刻的卫星速度;具体为将地心地固直角坐标系中的卫星坐标对时间求导,从而得到当前时刻的卫星速度为:为:为:式中为升交点赤经对时间的变化率;为轨道倾角ik对时间的变化率;为卫星轨道平面的位置对时间的变化率;f. 采用当前卫星的星历信息,计算当前时刻的卫星加速度;具体为将步骤e得到的当前时刻的卫星速度对时间求导,从而得到当前时刻的卫星加速度为:为:为:
式中为第一系数,且;为第二系数,且;为轨道倾角ik对时间的二次求导;为卫星轨道平面的位置对时间的二次求导;g. 计算完成后,设置当前卫星的状态信息有效标志为有效,并将计算结果进行保存,保存到全局卫星状态信息缓存中;s3. 在步骤s2进行的同时,采用第二线程,结合二次运动模型理论,以二次曲线拟合的方式进行卫星状态信息外推,并进行最小二乘解算;本步骤所采用的“以二次曲线拟合的方式进行卫星状态信息外推”的技术方案,是一种全新的计算方式;该方案基于卫星轨道理论模型和大量测试数据的统计结果,结合二次运动模型理论,实现了pvt解算的快速进行和优化;具体实施时,步骤s3具体包括如下步骤:a. 从全局卫星状态信息缓存中获取当前卫星的状态信息有效标志:若当前卫星的状态信息有效标志为无效,则本步骤结束;并获取下一颗卫星的状态信息有效标志;若当前卫星的状态信息有效标志为有效,则计算并得到当前卫星的外推时间间隔为,式中t
pos
为当前pvt解算时刻,ts为当前卫星的最近一次状态信息计算的时刻;b. 设置时间阈值t
th1
(优选为10s):若当前卫星的外推时间间隔大于t
th1
时,判定当前卫星异常(包括失锁、星历无效等情况),不再外推当前卫星的状态信息,且不再使用当前卫星参与后续的解算过程;若当前卫星的外推时间间隔小于或等于t
th1
时,采用如下二次曲线模型外推当前卫星在t
pos
时刻的位置和速度:时刻的位置和速度:时刻的位置和速度:时刻的位置和速度:时刻的位置和速度:时刻的位置和速度:式中上标s表示为当前卫星s;为当前卫星s在ts时刻的位置;为当前卫星s在ts时刻的速度;为当前卫星s在ts时刻的加速
度;为当前卫星s外推到t
pos
时刻的位置;为当前卫星s外推到t
pos
时刻的速度;c. 重复步骤a~b,直至遍历完成所有卫星;d. 将所有正常卫星的状态信息外推后,采用外推得到的状态信息进行后续的最小二乘解算;本步骤在进行状态信息外推时,当外推时间间隔在几分钟之内时,采用二次曲线拟合计算所得的卫星位置和速度的误差分别小于10cm和1mm/s,这对定位结果产生的误差极其微小,可以忽略不计;而它的计算比直接用卫星星历参数来计算卫星位置和速度的方法要快约20倍,极大的提升了pvt解算的效率;s4. 根据步骤s3得到的最小二乘解算结果,得到最终的定位结果,完成pvt解算。
29.在具体实施时,采用第一线程进行计算的优先级低于采用第二线程进行计算的优先级;同时,步骤s2计算得到卫星状态信息的计算结果后,将计算结果更新到全局卫星状态信息缓存中,供步骤s3进行卫星状态信息外推使用;在对全局卫星状态信息缓存进行任何操作时,均采用互斥量对全局卫星状态信息缓存进行防护,从而保证卫星状态信息的完整性、操作的可靠性和安全性。
30.以下结合一个实施例,对本发明的解算方法进行进一步说明:导航接收机使用arm a7处理器600mhz主频进行四系统pvt定位解算,四系统分别为北斗b1,gps的l1,glo的l1和gal的e1,第一线程调度频度100hz(要求每10ms计算一颗卫星的状态信息),四系统卫星最大数n为127,pvt定位解算频度为10hz(要求1秒内解算10次,即每100ms解算一次),运行步骤如下:(1) t1+100*k1 ms时刻,接收机跟踪了15颗北斗卫星、10颗gps卫星、7颗glo卫星、8颗gal卫星,这些卫星的观测量数据组装到obs数据结构中;同时,所有卫星的导航电文数据全部接收完成,已组装到nav数据结构中。
31.(2) t1+100*k1 ms时刻,第一线程和第二线程都满足执行时机,各自启动执行。系统优先执行高优先级的第二线程,秒内第一拍定位解算使用之前已经计算并存储好的卫星状态信息集合进行外推,第一拍解算完成耗时3.53ms。
32.(3) t1+100*k1+3.53 ms时刻,开始执行第一线程,运行了0.16ms完成卫星s的状态信息计算。
33.(4) t1+100*k1+10 ms时刻,再次执行第一线程,运行了0.01ms跳过卫星s+1的计算,因为该卫星号没有被跟踪。
34.(5) t1+100*k1+20 ms时刻,再次执行第一线程,运行了0.17ms完成卫星s+2的状态信息计算。
35.……
(6) t1+100*k1+90 ms时刻,再次执行第一线程,运行了0.16ms完成卫星s+9的状态信息计算。
36.(7) t1+100*k2 ms时刻,秒内第二拍定位解算启动执行,再次组装obs数据,并使用已经更新了卫星s、卫星s+2、

、卫星s+9的状态信息的集合进行外推,第二拍解算完成耗时3.54ms。
37.(8)上述步骤(7)中,k2为k1的下一拍解算,把k2当作k1,继续重复步骤(1)。
38.(9) 步骤(2)中,在接收机开机的时候卫星状态信息集合是空的,第二线程定位解算得不到定位结果,导航接收机的第一个定位结果输出会在步骤(2)中k1迭代到某一个kn时给出。
39.(10) 步骤(6)和步骤(7)会存在同时操作全局卫星状态信息缓存的可能性,接收机引入互斥量防护技术消除了多线程同时操作的灾难性后果。
40.以下对传统的pvt解算方法和本发明所述的高性能pvt解算方法进行性能对比:导航接收机使用arm a7处理器600mhz主频进行四系统pvt定位解算,定位解算频度分别配置为10hz、20hz、50hz和100hz,对比传统方法和本发明所述方法的实验结果数据如下表1所示:表1实验数据对比示意表从上表数据可以看出:使用传统解算方法最大只能支持50hz解算频度,对于100hz解算频度已经无法支持了。而使用本发明所述的高性能pvt解算方法,接收机系统的cpu使用率和解算总耗时都大大降低。在20hz解算的情况下,本发明方法系统负担下降为传统方法的60%左右,pvt解算效率提升了62%左右,且能够支持到更高的100hz频度解算,改善效果非常明显。
41.如图4所示为本发明的导航方法的方法流程示意图:本发明提供的这种包括了所述高性能pvt解算方法的导航方法,包括如下步骤:(1)实时获取导航所需的数据信息;(2)采用上述的高性能pvt解算方法进行实时pvt解算,得到实时的定位结果;(3)根据步骤(2)得到的实时定位结果,完成最终的导航。
完整全部详细技术资料下载
当前第1页 1  2 
相关技术
  • 一种基于相关分析的变频电机断...
  • 一种土木工程用差动倾斜测量仪
  • 飞机疲劳强度测试用局部振动载...
  • 一种玻璃耐久性检测设备的制作...
  • 上转换拉曼传感方法及应用与流...
  • 一种适用于低速和高速的潜油电...
  • 载波半周修复方法及其RTK整...
  • 一种潜油电机机组井下工况测试...
  • 一种野外动力取土的装置
  • 一种重新收敛的精密单点定位方...
网友询问留言 已有0条留言
  • 还没有人留言评论。精彩留言会获得点赞!
1

玻璃钢生产厂家自贡玻璃钢花盆花器玻璃钢传统人物雕塑哪里实惠台州玻璃钢海豚雕塑价格山西步行街玻璃钢雕塑销售电话清远欧式玻璃钢人物雕塑杭州商场美陈布置价格怎么算郑州玻璃钢浮雕铜雕塑定制厂家内蒙古玻璃钢卡通雕塑晋江玻璃钢雕塑厂家徐州玻璃钢雕塑厂家供应福建孙思邈人物玻璃钢雕塑厂商场美陈构思常州玻璃钢海豚雕塑定制焦作玻璃钢铜雕塑南阳铜玻璃钢卡通雕塑厂家山西抽象玻璃钢雕塑图片玻璃钢花盆模具制作流程定西广场玻璃钢雕塑定制新乡广场玻璃钢雕塑制作湛江玻璃钢雕塑厂家直销新郑玻璃钢仿铜雕塑定制厂家玻璃钢动物雕塑报价明细表栖霞玻璃钢卡通雕塑玻璃钢雕塑贴金箔台湾玻璃钢花盆昆明火烈鸟玻璃钢雕塑公司达川玻璃钢花盆花器江西玻璃钢雕塑人像定做厂家山东农场迎宾玻璃钢雕塑双层玻璃钢花盆图片价格香港通过《维护国家安全条例》两大学生合买彩票中奖一人不认账让美丽中国“从细节出发”19岁小伙救下5人后溺亡 多方发声单亲妈妈陷入热恋 14岁儿子报警汪小菲曝离婚始末遭遇山火的松茸之乡雅江山火三名扑火人员牺牲系谣言何赛飞追着代拍打萧美琴窜访捷克 外交部回应卫健委通报少年有偿捐血浆16次猝死手机成瘾是影响睡眠质量重要因素高校汽车撞人致3死16伤 司机系学生315晚会后胖东来又人满为患了小米汽车超级工厂正式揭幕中国拥有亿元资产的家庭达13.3万户周杰伦一审败诉网易男孩8年未见母亲被告知被遗忘许家印被限制高消费饲养员用铁锨驱打大熊猫被辞退男子被猫抓伤后确诊“猫抓病”特朗普无法缴纳4.54亿美元罚金倪萍分享减重40斤方法联合利华开始重组张家界的山上“长”满了韩国人?张立群任西安交通大学校长杨倩无缘巴黎奥运“重生之我在北大当嫡校长”黑马情侣提车了专访95后高颜值猪保姆考生莫言也上北大硕士复试名单了网友洛杉矶偶遇贾玲专家建议不必谈骨泥色变沉迷短剧的人就像掉进了杀猪盘奥巴马现身唐宁街 黑色着装引猜测七年后宇文玥被薅头发捞上岸事业单位女子向同事水杯投不明物质凯特王妃现身!外出购物视频曝光河南驻马店通报西平中学跳楼事件王树国卸任西安交大校长 师生送别恒大被罚41.75亿到底怎么缴男子被流浪猫绊倒 投喂者赔24万房客欠租失踪 房东直发愁西双版纳热带植物园回应蜉蝣大爆发钱人豪晒法院裁定实锤抄袭外国人感慨凌晨的中国很安全胖东来员工每周单休无小长假白宫:哈马斯三号人物被杀测试车高速逃费 小米:已补缴老人退休金被冒领16年 金额超20万

玻璃钢生产厂家 XML地图 TXT地图 虚拟主机 SEO 网站制作 网站优化