一种残膜回收机防缠绕挑膜装置的制 一种秧草收获机用电力驱动行走机构

一种抗洋流扰动的自主水下机器人动力定位方法和系统与流程

2021-10-09 03:12:00 来源:中国专利 TAG:水下 机器人 洋流 自主 扰动


1.本发明涉及自主水下机器人控制技术领域,特别是关于一种抗洋流扰动的自主水下机器人动力定位方法和系统。


背景技术:

2.自主水下机器人(autonomous underwater vehicle,auv)需要在复杂的海洋环境中通过动力定位控制进行定点作业,以完成海洋资源勘测、开发等方面的任务。然而,恶劣海洋环境下动态洋流扰动、强风、海浪等不确定因素作用下自主式水下机器人的受力会发生变化,其动力定位控制系统也面临着控制模型参数摄动问题,对自主式水下机器人定点悬停作业过程中的稳定性造成了较大的影响。
3.现有的自主水下机器人动力定位控制方法存在明显的不足:1)在控制算法的设计上没有有效考虑执行器的饱和约束,会影响控制算法的控制性能;2)没有有效消除洋流扰动和动力学被控模型不匹配等因素所带来不利影响的方法,动力定位方法的鲁棒性较差。


技术实现要素:

4.本发明的目的在于提供一种抗洋流扰动的自主水下机器人动力定位方法和系统来克服或至少减轻现有技术的上述缺陷中的至少一个。
5.为实现上述目的,本发明提供一种抗洋流扰动的自主水下机器人动力定位方法,该方法包括:
6.步骤1,采集自主水下机器人的测量状态量;
7.步骤2,根据设定的自主水下机器人的期望状态信息和所述测量状态量,构建自主水下机器人运动学状态空间模型和依据饱和性预设的约束,在二次规划求解模型预测控制的优化问题后,再通过动力学模型输出自主水下机器人的期望控制量;
8.步骤3,根据所述测量状态量,获取洋流速度估计值;
9.步骤4,根据所述期望控制量、洋流速度估计值和自主水下机器人的速度,通过积分滑模面,获得用于抵消外部扰动的控制量,并通过饱和函数获得滑膜控制量;
10.步骤5,根据所述滑膜控制量和期望控制量,获取实际控制量;
11.步骤6,返回步骤2,直至自主水下机器人的位姿调整至期望位姿。
12.进一步地,所述步骤1中,当前时刻k的所述测量状态量被描述为x(k)=[η(k);v(k

1)],式中,η(k)表示k时刻的惯性坐标系下的位姿测量值,v(k

1)表示k

1时刻的惯性坐标系下的速度测量值;
[0013]
所述步骤3获取所述洋流速度估计值的方法具体包括:
[0014]
通过式(7)描述的洋流观测器,计算所述洋流速度估计值:
[0015]
[0016]
式中,v
f
(k)表示k时刻的洋流观测器计算得到的载体坐标系下的洋流速度估计值,表示k时刻的洋流观测器计算得到的惯性坐标系下的洋流速度估计值,表示的导数,表示k时刻的洋流观测器计算得到的惯性坐标系下的位姿估计值,表示的导数,k
v
和k
η
均表示增益矩阵,j表示惯性坐标系到载体坐标系的坐标转换矩阵。
[0017]
进一步地,所述步骤2中,所述运动学状态空间模型设置为预测时域n
p
和控制时域n
c
下的式(10):
[0018][0019]
所述优化问题被描述为式(12):
[0020][0021]
式中:
[0022][0023][0024]
表示k 1时刻在未来np个预测步长内的预测状态量,x(k|k)表示k时刻预测的第一个状态量,u(k)表示k时刻在未来n
c
个预测步长内的预测控制量,u
*
(k)表示u(k)的最优值,表示k时刻状态量的系数矩阵,表示k时刻控制量的系数矩阵,表示权重矩阵,表示权重矩阵,表示期望状态量,表示k时刻在未来np个预测步长内的预测状态量。
[0025]
进一步地,所述约束设置为s.t.:
[0026]
u(i|k)≤u
max
[0027]

u(i|k)≤

u
min
[0028][0029][0030]
式中,u(i|k)表示k时刻在未来第i时刻的预测控制量,u
max
表示u(k)中的最大值,u
min
表示u(k)中的最小值,v(i|k)表示k时刻在未来第i时刻的预测速度,c表示哥氏力离心矩阵,d表示流体阻尼矩阵,f表示静水恢复力,m表示期望控制量中的力矩控制量在载体坐标系下沿y轴的分量,t为采样时间,τ
max
表示期望控制量中的五自由度力的最大值,τ
min
表示所述五自由度力的最小值。
[0031]
进一步地,所述步骤2中,通过式(9)描述的控制量u(k)、以及式(13),计算期望加速度
[0032]
u(k)=v(k)

v(k

1)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)
[0033][0034]
通过式(14),计算得到期望控制量u
mpc
(k):
[0035]
u
mpc
(k)=mu
*
(k)/t

(c d)v(k)

f
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(14)
[0036]
式中,u
*
(k)是u
*
(k)的第一个元素。
[0037]
进一步地,所述步骤4中,所述积分滑模面被描述为s(v),通过如下式(15)得到所述控制量u
ismc
(t):
[0038][0039]
其中,v(t)表示t时刻的速度,h(v(τ))是状态反馈项。
[0040]
进一步地,所述步骤4中,在载体坐标系下沿x轴方向的所述饱和函数x
ismc
(t)被描述为式(14):
[0041][0042]
其中,ρ
x
为海浪扰动力的上限,s
u
(t)表示载体坐标系相对于惯性坐标系沿x轴方向速度的积分滑模面。
[0043]
本发明实施例提供的抗洋流扰动的自主水下机器人动力定位系统,该系统包括:
[0044]
测量单元,其用于采集自主水下机器人的测量状态量;
[0045]
模型预测控制器,其用于根据设定的自主水下机器人的期望状态信息和所述测量状态量,构建自主水下机器人运动学状态空间模型和依据饱和性预设的约束,在二次规划求解模型预测控制的优化问题后,再通过动力学模型输出自主水下机器人的期望控制量;
[0046]
洋流观测器,其用于根据所述测量状态量,获取洋流速度估计值;
[0047]
滑膜控制器,其用于根据所述期望控制量、洋流速度估计值和自主水下机器人的速度,通过积分滑模面,获得用于抵消外部扰动的控制量,并通过饱和函数获得滑膜控制量;
[0048]
计算单元,其用于根据所述滑膜控制量和期望控制量,获取实际控制量,并监控自主水下机器人的位姿调整至期望位姿。
[0049]
进一步地,所述测量单元获得的当前时刻k的所述测量状态量信息被描述为x(k)=[η(k);v(k

1)],式中,η(k)表示k时刻的惯性坐标系下的位姿测量值,v(k

1)表示k

1时刻的惯性坐标系下的速度测量值;
[0050]
所述洋流观测器获取所述洋流速度估计值的方法具体包括:
[0051]
通过式(7)描述的洋流观测器,计算所述洋流速度估计值:
[0052][0053]
式中,v
f
(k)表示k时刻的洋流观测器计算得到的载体坐标系下的洋流速度估计
值,表示k时刻的洋流观测器计算得到的惯性坐标系下的洋流速度估计值,表示的导数,表示k时刻的洋流观测器计算得到的惯性坐标系下的位姿估计值,表示的导数,k
v
和k
η
均表示增益矩阵,j表示惯性坐标系到载体坐标系的坐标转换矩阵。
[0054]
进一步地,所述运动学状态空间模型设置为预测时域n
p
和控制时域n
c
下的式(9):
[0055][0056]
所述优化问题被描述为式(11):
[0057][0058]
式中:
[0059][0060][0061]
表示k 1时刻在未来np个预测步长内的预测状态量,x(k|k)表示k时刻预测的第一个状态量,u(k)表示k时刻在未来n
c
个预测步长内的控制量,u
*
(k)表示u(k)的最优值,表示k时刻状态量的系数矩阵,表示k时刻控制量的系数矩阵,表示权重矩阵,表示权重矩阵,表示期望状态量,表示k

1时刻在未来np个预测步长内的预测状态量;
[0062]
所述约束设置为s.t.:
[0063]
u(i|k)≤u
max
[0064]

u(i|k)≤

u
min
[0065][0066][0067]
式中,u(i|k)表示k时刻在未来i时刻的预测控制量,u
max
表示u(k)中的最大值,u
min
表示u(k)中的最小值,v(i|k)表示k时刻在未来i时刻的预测速度,c表示哥氏力离心矩阵,d表示流体阻尼矩阵,f表示静水恢复力,m表示期望控制量中的力矩控制量在载体坐标系下沿y轴的分量,t为采样时间,τ
max
表示期望控制量中的五自由度力的最大值,τ
min
表示所述五自由度力的最小值。
[0068]
本发明由于采取以上技术方案,其具有以下优点:
[0069]
在复杂的海洋环境中,自主式水下机器人会受到不规则的洋流干扰力,为此设计了洋流观测器以通过估计实时的洋流速度,以计算洋流对自主式水下机器人的干扰力。而mpc可以通过设计约束考虑执行器的饱和特性。同时,为解决由于控制参考模型与自主式水
下机器人实际动力学模型存在的不匹配问题和洋流估计产生的较小有界误差,本发明设计了ismc以解决模型不匹配的问题,同时消除洋流估计产生的较小有界误差,以提高控制系统的鲁棒性。
附图说明
[0070]
图1为本发明中的描述自主式水下机器人运动的两个参考坐标系。
[0071]
图2是本发明中自主式水下机器人动力定位控制系统结构框图。
具体实施方式
[0072]
下面结合附图和实施例对本发明进行详细的描述。
[0073]
本实施例提供的抗洋流扰动的自主水下机器人动力定位方法包括:
[0074]
步骤1,首先定义两个参考坐标系:惯性坐标系和载体坐标系。其中,惯性坐标系固定于地面。载体坐标系则固定于自主式水下机器人,随自主式水下机器人以任何形式移动。两个参考坐标系如图1所示,图中示意出的e

ξηζ表示固定于地面的惯性坐标系,o

xyz表示固定于自主式水下机器人的载体坐标系,o通常选定在自主式水下机器人的质心,但不限于此。
[0075]
自主式水下机器人的位姿变量被描述为其中,x,y,z为o在惯性坐标系下分别沿x轴、y轴、z轴的坐标,θ,ψ为o分别绕惯性坐标系下的x轴、y轴、z轴的角度,即:纵倾角、航向角、横滚角,这些参数的具体数值可通过自主式水下机器人的导航系统反馈得到。
[0076]
自主式水下机器人的速度变量被描述为v=[u,v,w,p,q,r]
t
。其中,u,v,w为载体坐标系相对于惯性坐标系分别沿x轴、y轴、z轴的速度,p,q,r为载体坐标系相对于惯性坐标系分别围绕x轴、y轴、z轴的角速度。
[0077]
惯性坐标系与载体坐标系之间的转换关系可由自主式水下机器人的运动学方程(1)反映。
[0078]
η
t
=jv
t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0079]
其中:
[0080][0081][0082][0083]
j
11
=cosψcosθ
[0084]
j
12


sinψcosθ cosψsinθsinψ
[0085]
j
13
=sinψsinφ cosψcosφsinθ
[0086]
j
21
=sinψcosφ
[0087]
j
22
=sinψsinθsinφ cosψcosφ
[0088]
j
23


cosψsinψ sinθsinψcosφ
[0089]
j
31


sinθ
[0090]
j
32
=cosψsinφ
[0091]
j
33
=cosψcosφ
[0092]
动力学关系式描述为式(2):
[0093][0094]
其中,m

表示自主式水下机器人惯性矩阵,被描述为式(3);c表示自主式水下机器人的哥氏力离心矩阵,被描述为式(4);d表示自主式水下机器人的流体阻尼矩阵,被描述为式(5);f表示自主式水下机器人的静水恢复力,被描述为式(6);τ=[x;y;z;0;m;n]为自主式水下机器人的控制量,x、y、z表示控制量中的五自由度力控制量分别沿x轴、y轴和z轴的分量,m、n表示控制量中的力矩控制量分别沿y轴和z轴的分量。
[0095][0096]
式(3)中:
[0097][0098][0099][0100][0101]
m表示自主式水下机器人的质量,和表示附加质量,i
xx
表示自主式水下机器人绕x轴的转动惯量,i
yy
表示自主式水下机器人绕y轴的转动惯量,i
zz
表示自主式水下机器人绕z轴的转动惯量,这些动力学特性参数可以通过力学试验获得。
[0102][0103]
式(4)中:
[0104][0105][0106][0107][0108]
y
uv
,x
wq
,x
qq
,x
vr
,x
rr
,y
wp
,y
pq
,y
ur
,z
vp
,z
rp
,z
up
,m
uw
,n
uv
,m
vp
,m
rp
,m
uq
,n
wp
,n
pq
和n
ur
是交叉项附加质量,上述参数的具体数值可以通过力学试验得到。
[0109][0110]
式(5)中:
[0111][0112][0113][0114][0115]
x
uu
,y
vv
,z
ww
,y
rr
,z
qq
,m
ww
,n
vv
,k
pp
,m
qq
和n
rr
表示横流阻力系数,上述参数的具体数值可以通过力学试验得到。
[0116][0117]
w表示自主式水下机器人的重量,b表示自主式水下机器人的浮力,z
b
表示自主式水下机器人的浮心在载体坐标系下对应z轴的坐标。
[0118]
步骤2,洋流观测器用于估计当前时刻的洋流速度,以获得更精准的动力学模型,例如:洋流观测器可以描述为式(7):
[0119][0120]
式中,v
f
(k)表示k时刻的洋流观测器计算得到的载体坐标系下的洋流速度估计值,v
f
(k)=[u
f
(k),v
f
(k),w
f
(k),0,0,0]
t
,u
f
(k),v
f
(k),w
f
(k)分别为v
f
(k)在载体坐标系下分别沿x轴、y轴和z轴的分量,表示k时刻的洋流观测器计算得到的惯性坐标系下的洋流速度估计值,表示的导数,表示k时刻的洋流观测器计算得到的惯性坐标系下的位姿估计值,表示的导数,k
v
和k
η
均表示增益矩阵。
[0121]
本发明实施例通过构建洋流观测器估计洋流的速度以计算洋流扰动力,可实现精准动力学估计。
[0122]
步骤3,基于自主式水下机器人的运动学模型构造,构建由式(8)描述的离散的状态空间模型,用于二次规划的优化求解:
[0123]
x(k 1)=a(k)x(k) b(k)u(k)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0124]
其中:
[0125]
x(k)=[η(k);v(k

1)]
[0126][0127]
b(k)=[j

(k)t;b

]
[0128][0129]
[0130][0131]
u(k)=v(k)

v(k

1)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)
[0132]
x(k 1)表示k 1时刻的测量状态量,a(k)表示状态量的系数矩阵,x(k)表示k时刻的测量状态量,b(k)表示是控制量的系数矩阵,u(k)表示k时刻的控制量,v(k

1)表示k

1时刻的速度信息,t表示采样周期,j(k)表示坐标转换矩阵,j

(k)表示j(k)的逆矩阵,b

表示b的逆矩阵,j1表示线速度坐标转换矩阵,03×2表示3*2的零矩阵,03×3表示是三阶零矩阵,j
′2表示角速度坐标转换矩阵。
[0133]
定义预测时域和控制时域分别为n
p
和n
c
,状态空间模型(8)可转换为式(10):
[0134][0135]
其中:
[0136][0137][0138][0139][0140]
表示k !时刻在未来np个预测步长内的预测状态量,表示k时刻状态量的系数矩阵,表示k时刻控制量的系数矩阵,u(k)表示k时刻在未来n
c
个预测步长内的控制量,x(k|k)表示k时刻预测的第一个状态量,u(k|k)表示k时刻预测的第一个系统控制量,这些参数的具体数值可以通过二次规划算出来的。
[0141]
通过求解一个标准的二次规划优化问题,其中的代价函数描述为式(11)所示:
[0142]
j(k)=0.5*u(k)h(k)u(k) f
t
(k)u(k)
ꢀꢀꢀꢀꢀꢀꢀꢀ
(11)
[0143]
其中:
[0144]
[0145][0146][0147][0148][0149]
表示期望状态量,表示权重矩阵,表示权重矩阵,这些参数的具体数值可以通过试验调试获得。
[0150]
上述优化问题的受到的动力学约束如下:
[0151]
u(k)≤u
max
[0152]

u(k)≤

u
min
[0153][0154][0155]
其中,t为采样时间,u
max
表示控制量的最大值,u
min
表示控制量的最小值,τ
max
表示推力的最大值,τ
min
表示推力的最小值,这些参数可以通过预先设置获得。
[0156]
用于二次规划的优化问题为:
[0157][0158]
u(i|k)≤u
max
[0159]

u(i|k)≤

u
min
[0160][0161][0162]
其中,u
*
(k)=u
*
(k|k)是控制序列u
*
(k)的第一个元素,u(i|k)表示k时刻预测i时刻的控制量,u
max
表示控制序列中的最大值,u
min
表示控制序列中的最小值,这些参数的具体数值可以通过预先设置获得。
[0163]
由于控制量u(k)为当前时刻的速度增量,动力学方程中的加速度可再由下式(11)近似计算:
[0164]
[0165]
t表示采样周期,通过调试选取最佳值。
[0166]
故根据动力学方程(2),可计算得到期望控制量u
mpc
(k),其通过式(14)计算获得,其被描述为u
mpc
=[x
mpc
,y
mpc
,z
mpc
,0,m
mpc
,n
mpc
]
t

[0167]
u
mpc
(k)=mu
*
(k)/t

(c d)v(k)

f
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(14)
[0168]
本发明实施例通过mpc的控制量约束解决执行器的饱和特性问题,通过mpc的滚动优化来提高控制性能。
[0169]
步骤4,通过零阶保持器将离散的mpc控制量u
mpc
(k)转变为连续型信号后,通过如下式(15),得到滑膜控制量u
ismc
(t),其被描述为u
ismc
(t)=[x
ismc
(t),y
ismc
(t),z
ismc
(t),0,m
ismc
(t),n
ismc
(t)]
t

[0170][0171]
h(v(i))=m
‑1((c d)v(i) f)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(16)
[0172]
其中,s(v)=(s
u
(t),s
v
(t),s
w
(t),s
p
(t),s
q
(t),s
r
(t))
t
为积分滑模面,s
u
(t)表示u(t)的积分滑模面,s
v
(t)表示v(t)的积分滑模面,s
w
(t)表示w(t)的积分滑模面,s
p
(t)表示p(t)的积分滑模面,s
q
(t)表示q(t)的积分滑模面,s
r
(t)表示r(t)的积分滑模面。文中,t表示连续信号,k表示离散信号。
[0173]
通过引入饱和函数,来缓冲真实系统中抖振现象的不利影响,以x
ismc
(t)为例:
[0174][0175]
其中ρ
x
为海浪扰动力的上限。当然,载体坐标系相对于惯性坐标系沿其它方向的速度和角速度的饱和函数与x
ismc
(t)相同,在此不再展开描述。
[0176]
本发明实施例中的ismc输出的控制量可以在不产生较大系统抖振的前提下减小模型不匹配和扰动误差的影响,以进一步改进控制算法的控制效果。
[0177]
步骤5,根据期望控制量u
mpc
(k)与滑膜控制量u
ismc
(t),通过式(17),获得τ(k)=[x(k),y(k),z(k),0,m(k),n(k)]
t
,其中,x(k)表示k时刻载体坐标系x轴方向期望的力,y(k)表示k时刻载体坐标系y轴方向期望的力,z(k)表示k时刻载体坐标系z轴方向期望的力,m(k)表示k时刻载体坐标系期望的绕y轴的力矩,n(k)表示k时刻载体坐标系期望的绕y轴的力矩。
[0178]
τ(k)=u
mpc
(k) u
ismc
(k)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(17)
[0179]
步骤6,返回步骤2,直至自主水下机器人的位姿调整至期望位姿。
[0180]
本发明通过导航定位系统获得自主式水下机器人的速度和姿态等位姿信息,同时设计了洋流观测器估计洋流速度,为动力定位系统提供信息保障;基于mpc设计了动力定位系统的控制方法,进一步地,结合ismc消除模型不匹配和洋流扰动估计有界误差对动力定位系统控制精度的影响。该方法和系统提高了海洋环境扰动作用下自主式水下机器人动力定位系统的鲁棒性,可以实现自主式水下机器人精准定位作业,具有客观可行的工程应用价值。
[0181]
本发明实施例还提供一种抗洋流扰动的自主水下机器人动力定位系统,该系统包
括测量单元、模型预测控制器、洋流观测器、滑膜控制器和计算单元:
[0182]
测量单元用于采集自主水下机器人的测量状态量。
[0183]
模型预测控制器用于根据设定的自主水下机器人的期望状态信息和所述测量状态量,构建自主水下机器人运动学状态空间模型和依据饱和性预设的约束,在二次规划求解模型预测控制的优化问题后,再通过动力学模型输出自主水下机器人的期望控制量。
[0184]
洋流观测器用于根据所述测量状态量,获取洋流速度估计值;
[0185]
滑膜控制器,其用于根据所述期望控制量、洋流速度估计值和自主水下机器人的速度,通过积分滑模面,获得用于抵消外部扰动的控制量,并通过饱和函数获得滑膜控制量。
[0186]
计算单元用于根据所述滑膜控制量和期望控制量,获取实际控制量,并监控自主水下机器人的位姿调整至期望位姿。
[0187]
所述测量单元获得的当前时刻k的所述测量状态量信息被描述为x(k),
[0188]
x(k)=[η(k);v(k

1)]
[0189]
式中,η(k)表示k时刻的惯性坐标系下的位姿测量值,v(k)表示k时刻的惯性坐标系下的速度测量值;
[0190]
所述洋流观测器获取所述洋流速度估计值的方法具体包括:
[0191]
通过式(7)描述的洋流观测器,计算所述洋流速度估计值:
[0192][0193]
式中,v
f
(k)表示k时刻的洋流观测器计算得到的载体坐标系下的洋流速度估计值,表示k时刻的洋流观测器计算得到的惯性坐标系下的洋流速度估计值,表示的导数,表示k时刻的洋流观测器计算得到的惯性坐标系下的位姿估计值,表示的导数,k
v
和k
η
均表示增益矩阵,j表示惯性坐标系到载体坐标系的坐标转换矩阵。
[0194]
所述运动学状态空间模型设置为预测时域n
p
和控制时域n
c
下的式(9):
[0195][0196]
所述优化问题被描述为式(11):
[0197][0198]
式中:
[0199][0200][0201]
表示k 1时刻在未来np个预测步长内的预测状态量,x(k|k)表示k时刻
预测的第一个状态量,u(k)表示k时刻在未来n
c
个预测步长内的控制量,u
*
(k)表示u(k)的最优值,表示k时刻状态量的系数矩阵,表示k时刻控制量的系数矩阵,表示权重矩阵,表示权重矩阵,表示期望状态量,表示k

1时刻在未来np个预测步长内的预测状态量;
[0202]
所述约束设置为s.t.:
[0203]
u(i|k)≤u
max
[0204]

u(i|k)≤

u
min
[0205][0206][0207]
式中,u(i|k)表示k时刻在未来i时刻的预测控制量,u
max
表示u(k)中的最大值,u
min
表示u(k)中的最小值,v(i|k)表示k时刻在未来i时刻的预测速度,c表示哥氏力离心矩阵,d表示流体阻尼矩阵,f表示静水恢复力,m表示期望控制量中的力矩控制量在载体坐标系下沿y轴的分量,t为采样时间,τ
max
表示期望控制量中的五自由度力的最大值,τ
min
表示所述五自由度力的最小值。
[0208]
上述参数的经验数值可以参考如下表:
[0209]
表1
[0210][0211]
最后需要指出的是:以上实施例仅用以说明本发明的技术方案,而非对其限制。本领域的普通技术人员应当理解:可以对前述各实施例所记载的技术方案进行修改,或者对其中部分技术特征进行等同替换;这些修改或者替换,并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。
再多了解一些

本文用于企业家、创业者技术爱好者查询,结果仅供参考。

发表评论 共有条评论
用户名: 密码:
验证码: 匿名发表

相关文献

  • 日榜
  • 周榜
  • 月榜