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

一种基于递推最小二乘法的铣削力在线识别方法及系统与流程

2021-10-09 01:49:00 来源:中国专利 TAG:铣削 在线 乘法 小二 机床


1.本发明属于机床铣削加工技术领域,具体涉及一种基于递推最小二乘法的铣削力在线识别方法及系统。


背景技术:

2.机床加工状态的感知与监测是实现加工过程智能化和保障加工质量的重要途径。铣削力作为机床加工过程中的一项重要参数,对于加工过程监测,刀具磨损状态、振动控制以及机床性能具有重要参考价值。因此,对加工过程中的铣削力进行识别研究具有重要意义。
3.铣削力测量与识别方法主要分为以下几种:基于力传感器的直接测量方法;基于铣削力模型的预测方法;基于位移信号的识别方法;基于电流信号的识别方法。有不少国内外学者就此展开了大量的研究,如西北工业大学的万敏利用位移传感器和刀具刚度计算方法进行铣削力识别(wan m,yuan h,feng j,zhang w h&yin w.industry

oriented method for measuring the cutting forces based on the deflections of tool shank[j].international journal of mechanical sciences,2017,130:315

323.),然而该方法需要有足够精确的刀具刚度模型。西安交通大学的王晨希通过加速度传感器结合共轭梯度最小二乘法进行铣削力识别(wang cx,zhang xw,qiao bj,cao hr&chen xf.dynamic force identification in peripheral milling based on cgls using filtered acceleration signals and averaged transfer functions[j].journal of manufacturing science and engineering,2019,141(6),064501.),然而该方法需要截取一段时间响应信号进行识别,无法随时间进行递推。
[0004]
通过文献调研发现,现有方法无法同时保障识别精度,降低设备成本,提高频域带宽,实现在线识别等要求。


技术实现要素:

[0005]
本发明所要解决的技术问题在于针对上述现有技术中的不足,提供一种基于递推最小二乘法的铣削力在线识别方法及系统,采用加速度信号进行铣削力识别,为解决现有铣削力监测的不足提供参考。
[0006]
本发明采用以下技术方案:
[0007]
一种基于递推最小二乘法的铣削力在线识别方法,包括以下步骤:
[0008]
s1、在铣床坐标系下,分别沿x、y进给方向进行力锤敲击实验,获得多组加速度频响函数,将每个方向上的频响进行平均得到跨点加速度频率响应函数;
[0009]
s2、进行铣削加工,选择采样频率,测量加工过程中x、y进给方向的加速度信号,并对加速度信号进行滤波预处理;
[0010]
s3、利用步骤s1获取的跨点加速度频率响应函数和步骤s2滤波预处理后的铣削过程加速度信号,采用递推最小二乘法对铣削力进行在线识别。
[0011]
具体的,步骤s1具体为:
[0012]
s101、分别在主轴箱体x、y方向上,靠近前端轴承处安装加速度传感器;
[0013]
s102、用力锤分别沿x、y方向敲击铣刀刀尖处,获取刀尖到主轴箱体的跨点加速度频率响应函数。
[0014]
进一步的,步骤s102中,刀尖到主轴箱体的跨点加速度频率响应函数表示刀尖处的铣削力与箱体处加速度响应之间的关系,具体为:
[0015][0016]
其中,h(w)为加速度频响函数,x(w)为箱体处的加速度响应,f(w)为刀尖处的铣削力。
[0017]
具体的,步骤s2具体为:
[0018]
保持加速度传感器的位置与步骤s1最终确定的安装位置一致;根据实际工艺需求和铣削稳定性叶瓣图确定加工参数;获取铣削加工过程中的加速度信号;由主轴转速计算刀齿通过频率,并根据刀齿通过频率对加速度信号进行梳状滤波,保留前3~5倍的刀齿通过频率,获取滤波后的加速度信号。
[0019]
进一步的,刀齿通过频率tpf计算如下:
[0020][0021]
其中,ω为主轴转速rpm,n为铣刀齿数。
[0022]
具体的,步骤s3具体为:
[0023]
s301,将步骤s1获取的跨点加速度频率响应函数h(w)
xx
、h(w)
xy
、h(w)
yx
、h(w)
yy
转换为时域上的单位脉冲响应函数h(t)
xx
、h(t)
xy
、h(t)
yx
、h(t)
yy

[0024]
s302,利用步骤s301中的单位脉冲响应函数,根据离散化杜哈梅积分构建最小二乘量测模型,将最小二乘量测模型扩展到x、y方向并进行变形;
[0025]
s303,选择数据长度,根据步骤s302中的最小二乘量测模型进行递推最小二乘计算,分别计算增益矩阵,进行状态递推和方差阵递推;
[0026]
s304,将x、y方向的识别铣削力从识别状态向量中分离,得到x、y方向的识别铣削力;
[0027]
s305,引入矩形窗,覆盖前一段无效数据,接收实时加速度数据,对输入信号进行实时更新。
[0028]
进一步的,步骤s302中,将最小二乘量测模型扩展到x、y方向并进行变形具体为:
[0029][0030]
其中,x
n
、y
n
分别为x、y方向第n时刻加速度响应;h
xxn
、h
xyn
、h
yxn
和h
yyn
分别为单位脉冲响应函数h(t)
xx
、h(t)
xy
、h(t)
yx
和h(t)
yy
第n时刻的值,f
xn
、f
yn
分别为x、y方向第n时刻铣削力。
[0031]
进一步的,步骤s303中,增益矩阵为:
[0032][0033]
状态递推为:
[0034][0035]
方差阵递推为:
[0036]
p
k
=(i

k
k
h
k
)p
k
‑1[0037]
其中,表示识别的铣削力,k表示当前时刻;z
k
表示加速度信号,r
k
为测量噪声矩阵,k
k
为增益矩阵,p
k
为方差矩阵,h
k
为k时刻的单位脉冲响应矩阵。
[0038]
进一步的,步骤s305中,根据脉冲响应函数h(t)
xx
、h(t)
xy
、h(t)
yx
、h(t)
yy
衰减接近于零的时间,确定中k值的大小;在k 1时刻,状态递推公式为:
[0039][0040]
k 1时刻用于更新的具体为:
[0041][0042]
其中,z
k 1
=[x
k 1
,y
k 1
]
t
为k 1时刻的加速度信号,k
k 1
为k 1时刻的增益矩阵。
[0043]
本发明的另一技术方案是,一种基于递推最小二乘法的铣削力在线识别系统,包括:
[0044]
函数模块,在铣床坐标系下,分别沿x、y进给方向进行力锤敲击实验,获得多组加速度频响函数,将每个方向上的频响进行平均得到跨点加速度频率响应函数;
[0045]
滤波模块,进行铣削加工,选择采样频率,测量加工过程中x、y进给方向的加速度信号,并对加速度信号进行滤波预处理;
[0046]
识别模块,利用函数模块获取的跨点加速度频率响应函数和滤波模块滤波预处理后的铣削过程加速度信号,采用递推最小二乘法对铣削力进行识别。
[0047]
与现有技术相比,本发明至少具有以下有益效果:
[0048]
本发明一种基于递推最小二乘法的铣削力在线识别方法,在铣床坐标系下,分别沿x、y进给方向进行力锤敲击实验,获得多组加速度频响函数,将每个方向上的频响进行平均得到加速度频响函数;进行铣削加工,选择采样频率,测量加工过程中x、y进给方向的加速度信号,并对加速度信号进行滤波预处理,通过滤波处理可以消除加速度信号中噪声成分的干扰,能够使识别出来的铣削力更加准确;利用加速度频响函数和滤波预处理后的铣削过程加速度信号,采用递推最小二乘法对铣削力进行在线识别,通过对多次模态实验获取的加速度频响函数取平均值能够使获取的频响更加准确,减少频响函数获取过程中随机误差的影响,保障铣削力识别结果的准确性。
[0049]
进一步的,通过将加速度传感器安装在主轴前轴承位置附近,可以保证振动能量衰减最少,减小随机误差的影响。切削力经过主轴支撑轴承将振动传递到箱体,这样可以使切削力整个传递路径最小,能量衰减最小。
[0050]
进一步的,在切削前进行叶瓣图分析确定加工参数,可以避免颤振的产生,提高铣削稳定性,减小非线性因素的影响,保证识别的准确性。
[0051]
进一步的,对采集到的加速度信号进行梳状滤波预处理,可以减小其他频率噪声对于识别过程的影响,进而提高铣削力识别的准确性。
[0052]
进一步的,通过计算刀齿通过频率tpf设置的目的或好处,给出与原理分析说明,可以确定铣削力主要分布频率,同时也给梳状滤波提供了参考频率。
[0053]
进一步的,通过将频域频响函数转换为时域脉冲响应函数,可以将整个识别过程映射到时域,为在线识别提供了识别模型。
[0054]
进一步的,通过将最小二乘模型扩展为x、y两个方向,保证识别的准确性。因为铣削过程中x、y方向的铣削力是相互耦合的,不是相互独立的。
[0055]
进一步的,通过对最小二乘量测模型的变形,将整个模型按照时间顺序重新排序,改变了原来按照方向排序的形式,可以保证对数据的实时更新辨识,为在线识别的基础。
[0056]
综上所述,本发明采用递推最小二乘法并结合矩形窗,对输入加速度信号进行不断的更新,以保证识别切削力的实时性,为铣削力的在线识别提供了新思路;采用加速度信号采集成本低,易于获得,且能够反映切削力各频谱成分。
[0057]
下面通过附图和实施例,对本发明的技术方案做进一步的详细描述。
附图说明
[0058]
图1为机床x、y方向示意图;
[0059]
图2为力锤敲击与加速度传感器安装位置示意图;
[0060]
图3为基于递推最小二乘法的铣削力在线识别方法流程图;
[0061]
图4为本发明用于验证所提出铣削力识别算法的各方向上的模拟仿真频响图,其中,(a)为x

x方向上的频响,(b)为x

y方向上的频响,(c)为y

x方向上的频响,(d)为y

y方向上的频响;
[0062]
图5为本发明所用的x方向仿真铣削力与识别出的铣削力对比图;
[0063]
图6为本发明所用的y方向仿真铣削力与识别出的铣削力对比图。
[0064]
其中,1.加速度传感器;2.力锤。
具体实施方式
[0065]
下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
[0066]
应当理解,当在本说明书和所附权利要求书中使用时,术语“包括”和“包含”指示所描述特征、整体、步骤、操作、元素和/或组件的存在,但并不排除一个或多个其它特征、整体、步骤、操作、元素、组件和/或其集合的存在或添加。
[0067]
还应当理解,在本发明说明书中所使用的术语仅仅是出于描述特定实施例的目的而并不意在限制本发明。如在本发明说明书和所附权利要求书中所使用的那样,除非上下文清楚地指明其它情况,否则单数形式的“一”、“一个”及“该”意在包括复数形式,
[0068]
还应当进一步理解,在本发明说明书和所附权利要求书中使用的术语“和/或”是指相关联列出的项中的一个或多个的任何组合以及所有可能组合,并且包括这些组合。
[0069]
在附图中示出了根据本发明公开实施例的各种结构示意图。这些图并非是按比例绘制的,其中为了清楚表达的目的,放大了某些细节,并且可能省略了某些细节。图中所示出的各种区域、层的形状及它们之间的相对大小、位置关系仅是示例性的,实际中可能由于制造公差或技术限制而有所偏差,并且本领域技术人员根据实际所需可以另外设计具有不同形状、大小、相对位置的区域/层。
[0070]
本发明提供了一种基于递推最小二乘法的铣削力在线识别方法,在机床坐标系下通过多次模态测试取平均值的方法获得x、y方向上的频响函数,并将频响函数转换成时域上的单位脉冲响应函数;通过在x、y合适的位置布置加速度传感器,获取铣削加工过程中加速度响应;根据铣削加工过程中的转速,选择对应的刀齿通过频率及其倍频对加速度信号进行梳状滤波处理;利用获得的单位脉冲响应函数,采用递推最小二乘法配合滑动矩形窗对滤波后的加速度信号进行处理,识别出铣削力。本发明采用递推最小二乘法并结合矩形窗,对输入加速度信号进行不断的更新,以保证识别切削力的实时性,为铣削力的在线识别提供了新思路。
[0071]
本发明一种基于递推最小二乘法的铣削力在线识别方法,包括以下步骤:
[0072]
s1、在如图1和图2所示的铣床坐标系下,分别沿着x、y进给方向进行力锤敲击实验,每个方向上进行多次敲击获得多组跨点加速度频率响应函数,将每个方向上的频响进行平均以获取更加准确的跨点加速度频率响应函数;
[0073]
s101、分别在主轴箱体x、y方向上,靠近前端轴承处安装加速度传感器1;
[0074]
s102、用力锤2分别沿着x、y方向敲击铣刀刀尖处,获取刀尖到主轴箱体的跨点加速度频率响应函数。
[0075]
刀尖到主轴箱体的跨点加速度频率响应函数表示刀尖处的铣削力与箱体处加速度响应之间的关系,具体为:
[0076][0077]
其中,h(w)为加速度频响函数,x(w)为箱体处的加速度响应,f(w)为刀尖处的铣削力。
[0078]
s2、设置合理的参数进行铣削加工,选择合适的采样频率,测量加工过程中x、y进给方向的加速度信号,并对加速度信号进行滤波预处理;
[0079]
s201,保持加速度传感器的位置与s103最终确定的安装位置一致,以确保传递路径的一致性;
[0080]
s202,根据实际工艺需求和铣削稳定性叶瓣图确定合适的加工参数,以避免颤振等非线性效应的影响;
[0081]
s203,利用安装在主轴箱体上的加速度传感器和数据采集系统,获取铣削加工过程中的加速度信号;
[0082]
s204,由主轴转速计算刀齿通过频率(tooth passing frequency,tpf),并根据tpf对加速度信号进行梳状滤波,保留前3至5倍tpf,获取滤波后的加速度信号;
[0083]
刀齿通过频率的计算公式为:
[0084][0085]
其中,ω为主轴转速rpm,n为铣刀齿数。
[0086]
s3、利用步骤s1获取的频响函数和步骤s2获取的铣削过程加速度信号,采用递推最小二乘法对铣削力进行识别。
[0087]
s301,将s102获取的跨点加速度频率响应函数h(w)
xx
、h(w)
xy
、h(w)
yx
、h(w)
yy
转换为时域上的单位脉冲响应函数h(t)
xx
、h(t)
xy
、h(t)
yx
、h(t)
yy

[0088]
s302,利用步骤s301中的单位脉冲响应函数,根据离散化杜哈梅积分构建最小二乘量测模型,将其扩展到x、y方向并进行一定的变形,以便进行在线识别;
[0089]
杜哈梅积分公式为:
[0090][0091]
其中,x(t)为加速度响应,h
xx
(t)为加速度单位脉冲响应函数,f
x
(t)为铣削力;
[0092]
离散化杜哈梅积分公式为:
[0093][0094]
等效为最小二乘量测模型:
[0095]
x=h
xx
f
x
v
[0096]
扩展到x、y方向为:
[0097][0098]
其中,δt为离散时间区间;n为采样数;v为测量噪声;x为x方向加速度响应;y为y方向加速度响应;f
x
为x方向铣削力;f
y
为y方向铣削力;h
xx
表示频响函数,下标表示不同方向的频响函数,如h
xx
表示xx方向的频响函数;
[0099]
变形后的扩展最小二乘量测模型为:
[0100][0101]
其中,x
n
、y
n
分别为x、y方向第n时刻加速度响应;h
xxn
、h
xyn
、h
yxn
和h
yyn
分别为单位脉冲响应函数h(t)
xx
、h(t)
xy
、h(t)
yx
和h(t)
yy
第n时刻的值,f
xn
、f
yn
分别为x、y方向第n时刻铣削力。
[0102]
s303,选择一定的数据长度,根据步骤s302中的最小二乘量测模型进行递推最小二乘计算,分别计算增益矩阵,进行状态递推和方差阵递推;
[0103]
增益矩阵公式为:
[0104][0105]
状态递推公式为:
[0106][0107]
方差阵递推公式为:
[0108]
p
k
=(i

k
k
h
k
)p
k
‑1[0109]
其中,表示识别的铣削力,k表示当前时刻,z
k
=[x
k
,y
k
]
t
表示加速度信号,r
k
为测量噪声矩阵,k
k
为增益矩阵,p
k
为方差矩阵,h
k
为k时刻的单位脉冲响应矩阵。
[0110]
s304,将x、y方向的识别铣削力从识别状态向量中分离出来,得到x、y方向的识别铣削力;
[0111]
s305,根据脉冲响应函数h(t)
xx
、h(t)
xy
、h(t)
yx
、h(t)
yy
衰减接近于零的时间,确定中k值的大小。
[0112]
当时间超过k时,如当前处于k 1时刻,则状态递推公式为:
[0113][0114]
该时刻用于更新的需要去除第一时刻识别的铣削力,具体形式为需要去除第一时刻识别的铣削力,具体形式为z
k 1
=[x
k 1
,y
k 1
]
t
为k 1时刻的加速度信号,k
k 1
为k 1时刻的增益矩阵,此处需要注意,当时间超过k时h
k
始终保持不变。
[0115]
通过覆盖前一段无效数据,接收实时加速度数据,对输入信号进行实时更新,以保证识别切削力的实时性。
[0116]
本发明再一个实施例中,提供一种基于递推最小二乘法的铣削力在线识别系统,该系统能够用于实现上述基于递推最小二乘法的铣削力在线识别方法,具体的,该基于递
推最小二乘法的铣削力在线识别系统包括函数模块、滤波模块以及识别模块。
[0117]
其中,函数模块,在铣床坐标系下,分别沿x、y进给方向进行力锤敲击实验,获得多组加速度频响函数,将每个方向上的频响进行平均得到跨点加速度频率响应函数;
[0118]
滤波模块,进行铣削加工,选择采样频率,测量加工过程中x、y进给方向的加速度信号,并对加速度信号进行滤波预处理;
[0119]
识别模块,利用函数模块获取的跨点加速度频率响应函数和滤波模块滤波预处理后的铣削过程加速度信号,采用递推最小二乘法对铣削力进行识别。
[0120]
本发明再一个实施例中,提供了一种终端设备,该终端设备包括处理器以及存储器,所述存储器用于存储计算机程序,所述计算机程序包括程序指令,所述处理器用于执行所述计算机存储介质存储的程序指令。处理器可能是中央处理单元(central processing unit,cpu),还可以是其他通用处理器、数字信号处理器(digital signal processor、dsp)、专用集成电路(application specific integrated circuit,asic)、现成可编程门阵列(field

programmable gatearray,fpga)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等,其是终端的计算核心以及控制核心,其适于实现一条或一条以上指令,具体适于加载并执行一条或一条以上指令从而实现相应方法流程或相应功能;本发明实施例所述的处理器可以用于基于递推最小二乘法的铣削力在线识别方法的操作,包括:
[0121]
在铣床坐标系下,分别沿x、y进给方向进行力锤敲击实验,获得多组加速度频响函数,将每个方向上的频响进行平均得到跨点加速度频率响应函数;进行铣削加工,选择采样频率,测量加工过程中x、y进给方向的加速度信号,并对加速度信号进行滤波预处理;利用跨点加速度频率响应函数和滤波预处理后的铣削过程加速度信号,采用递推最小二乘法对铣削力进行在线识别。
[0122]
本发明再一个实施例中,本发明还提供了一种存储介质,具体为计算机可读存储介质(memory),所述计算机可读存储介质是终端设备中的记忆设备,用于存放程序和数据。可以理解的是,此处的计算机可读存储介质既可以包括终端设备中的内置存储介质,当然也可以包括终端设备所支持的扩展存储介质。计算机可读存储介质提供存储空间,该存储空间存储了终端的操作系统。并且,在该存储空间中还存放了适于被处理器加载并执行的一条或一条以上的指令,这些指令可以是一个或一个以上的计算机程序(包括程序代码)。需要说明的是,此处的计算机可读存储介质可以是高速ram存储器,也可以是非不稳定的存储器(non

volatile memory),例如至少一个磁盘存储器。
[0123]
可由处理器加载并执行计算机可读存储介质中存放的一条或一条以上指令,以实现上述实施例中有关基于递推最小二乘法的铣削力在线识别方法的相应步骤;计算机可读存储介质中的一条或一条以上指令由处理器加载并执行如下步骤:
[0124]
在铣床坐标系下,分别沿x、y进给方向进行力锤敲击实验,获得多组加速度频响函数,将每个方向上的频响进行平均得到跨点加速度频率响应函数;进行铣削加工,选择采样频率,测量加工过程中x、y进给方向的加速度信号,并对加速度信号进行滤波预处理;利用跨点加速度频率响应函数和滤波预处理后的铣削过程加速度信号,采用递推最小二乘法对铣削力进行在线识别。
[0125]
为使本发明实施例的目的、技术方案和优点更加清楚,下面将结合本发明实施例
中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例是本发明一部分实施例,而不是全部的实施例。通常在此处附图中的描述和所示的本发明实施例的组件可以通过各种不同的配置来布置和设计。因此,以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围,而是仅仅表示本发明的选定实施例。基于本发明中的实施例,本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
[0126]
请参阅图1,将本发明所采用的数控铣床设置一个直角坐标系。将两个加速度传感器分别布置在x、y轴靠近前段轴承的位置。
[0127]
请参阅图3,用力锤分别沿x、y方向敲击刀尖位置,获取刀尖到箱体的频响函数,且需要经过多次模态实验求频响函数的平均值。
[0128]
按照上述方法,本实施例获取的频响函数如图4所示,分别为h(w)
xx
、h(w)
xy
、h(w)
yx
、h(w)
yy

[0129]
本实施例所采用的铣削参数为:
[0130]
转速8000rpm,进给480mm,切宽1mm,切深1mm。通过安装在箱体上的两个加速度传感器获取铣削过程中的加速度,并利用本发明所提出的方法识别出铣削力。
[0131]
识别结果如图5和图6所示,分别为x、y方向识别铣削力与测量铣削力对比图。
[0132]
综上所述,本发明一种基于递推最小二乘法的铣削力在线识别方法,操作简单,易于实现,能够准确的识别铣削过程中的铣削力,且通过数据不断更新,为铣削力在线识别提供一种新思路。
[0133]
本领域内的技术人员应明白,本技术的实施例可提供为方法、系统、或计算机程序产品。因此,本技术可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且,本技术可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd

rom、光学存储器等)上实施的计算机程序产品的形式。
[0134]
本技术是参照根据本技术实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器,使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
[0135]
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中,使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品,该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
[0136]
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上,使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理,从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
[0137]
以上内容仅为说明本发明的技术思想,不能以此限定本发明的保护范围,凡是按
照本发明提出的技术思想,在技术方案基础上所做的任何改动,均落入本发明权利要求书的保护范围之内。
再多了解一些

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

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

相关文献

  • 日榜
  • 周榜
  • 月榜