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

用于航天器损伤红外检测的热扩散效应缺陷定量识别方法与流程

2021-10-24 09:25:00 来源:中国专利 TAG:航天器 缺陷 检测 定量 扩散

用于航天器损伤红外检测的热扩散效应缺陷定量识别方法
1.技术邻域
2.本发明属于缺陷无损检测技术领域,更为具体地讲,涉及一种用于航天器损伤红外检测的热扩散效应缺陷定量识别方法。


背景技术:

3.高性能复合材料在航天领域有着广泛的使用,其安全性对于航天器的正常运行有着重要意义,因此,对其开展缺陷无损检测是一项必要的工作。红外热成像利用被检测材料内部热学性质差异以及热传导的不连续性在被检测材料表面产生的温度差异对缺陷进行成像,红外图像序列在缺陷位置的温度梯度特征信息反映了复合材料的缺陷特征,因此,可以基于红外图像序列构建突出缺陷特征的红外重构图像,实现复合材料的无损检测。
4.在获得突出缺陷特征的红外重构图像后,为了进一步对检测复合材料中的损伤进情况,需对缺陷的定量信息如面积进行一步的获取,才能有效地判断损伤的严重程度。对于红外重构图像中缺陷定量信息的提取,通常的做法是直接对缺陷进行分割并提取缺陷对应的像素点数量。但是,通过红外重构图像获得缺陷定量信息时需要考虑以下问题:对于缺陷面,热能集中在很小的体积内,使得缺陷对应区域的温度迅速升高,而完整区域的温度升温缓慢,因此在缺陷边缘部位就会形成较大的温度梯度,使得缺陷边缘的热量同样升高。并且缺陷边缘上累积的热量受到缺陷内部和完整区域的双重横向热传导的作用会形成边缘热扩散效应反映在热成像中,从而形成热扩散区域。热扩散区域的出现会扩大缺陷在红外重构图像中的显现面积,使得分割后的缺陷特征会包含热扩散区域部分,导致其像素点的统计结果不能准确地反映缺陷的定量信息。
5.有学者在考虑到边缘热扩散效应后,首先,针对样品复合材料,设计一个缺陷进行检测实验,利用温度特征中的温度峰值信息、材料的热传导系数与实际的缺陷大小进行计算拟合,得到拟合公式,然后,利用拟合公式对被检测复合材料的红外重构图像中的缺陷特征做相应的定量研究。由于不同的复合材料具有不同的热传导系数,这种方法需要对同一类别复合材料进行大量实验才能获得单一类型复合材料的拟合公式。


技术实现要素:

6.本发明的目的在于克服现有技术的不足,提供一种用于航天器损伤红外检测的热扩散效应缺陷定量识别方法,不需要进行实验拟合,也能实现缺陷定量的准确识别。
7.为实现上述发明目的,本发明用于航天器损伤红外检测的热扩散效应缺陷定量识别方法,其特征在于,包括:
8.(1)、红外重构图像缺陷特征区域提取
9.1.1)、对于基于红外图像序列进行重构得到的红外重构图像,将其像素值(温度特征值)从rgb颜色空间转换到l*a*b*颜色空间;
10.1.2)、根据红外重构图像的“a*”和“b*”颜色信息转换值,将其像素点聚类为k类,保留符合缺陷对应的高亮颜色信息的一类的像素点,其余类丢弃,保留的则为缺陷特征区
域图像;对缺陷特征区域图像进行二值化,将缺陷特征区域标记为高亮的白色,然后进行形态学开闭运算操作,来连接相邻的像素点,获得包含n个缺陷特征区域de1,de2,...,de
n
的红外重构缺陷分割图像;
11.统计n个缺陷特征区域de1,de2,...,de
n
的像素点数量,分别记为p1,p2,...,p
n

12.(2)、缺陷边缘像素点(边缘点)提取
13.利用边缘提取算法对红外重构缺陷分割图像进行边缘提取,得到边缘轮廓图像,提取每个缺陷特征区域的边缘点像素坐标,得到对应的n个边缘点像素坐标集pi1,pi2,...,pi
n

14.(3)、根据缺陷热扩散区域判断定量识别缺陷
15.3.1)、对于第n个缺陷特征区域de
n
,根据其对应的边缘点像素坐标集pi
n
得到其质心位置的像素点坐标(质心像素点坐标):(x
n_cen
,y
n_cen
);
16.3.2)、根据质心像素点坐标(x
n_cen
,y
n_cen
),在红外图像序列中获取对应的瞬态热响应曲线ttr
n_cen
,根据边缘点像素坐标集pi
n
中的边缘点像素坐标(x
n_m
,y
n_m
),在红外图像序列中获取对应的瞬态热响应曲线ttr
n_m
,m表示边缘点像素坐标集pi
n
的第m个边缘点,m=1,2,

,m
n
,m
n
为第n个缺陷特征区域de
n
对应的边缘点数量;
17.3.3)、对瞬态热响应曲线ttr
n_cen
、第m条瞬态热响应曲线ttr
n_m
分别求各帧(时刻)的单位温度变化率,得到温度变化序列δv
n_cen
、δv
n_m

18.3.4)、对温度变化序列δv
n_cen
、δv
n_m
各个对应时刻的温度变化率进行比较,得到加权因子序列vt
n_m
,其中,对于第t帧的值vt
n_m_t
为:
19.vt
n_m_t
=|δv
n_cen_t
,δv
n_m_t
|
1,2,3
,t=1,2,...,t
‑120.vt
n_m_t
=1,t=t
21.其中,δv
n_cen_t
、δv
n_m_t
分别为温度变化序列δv
n_cen
、δv
n_m
在第t帧的值,t为红外图像序列的帧数,其(第t帧的值vt
n_m_t
)含义为:
22.如果δv
n_cen_t
、δv
n_m_t
的差值小于变化阈值ε
δ
,则第t帧的值vt
n_m_t
为1,如果δv
n_cen_t
、δv
n_m_t
的差值不小于变化阈值ε
δ
,并且同为正或同为负,则第t帧的值vt
n_m_t
为2,如果δv
n_cen_t
、δv
n_m_t
的差值不小于变化阈值ε
δ
,并且正负相反,则第t帧的值vt
n_m_t
为3;
23.3.5)、构建第n个缺陷特征区域de
n
的权重序列ω
n_t
,对于其第t帧的值ω
n_t
,如果在瞬态热响应曲线ttr
n_cen
上,温度值最高的帧,取值为1.5,否则取值为1;
24.3.6)、计算瞬态热响应曲线ttr
n_cen
、第m条瞬态热响应曲线ttr
n_m
的距离:
[0025][0026]
其中,ttr
n_cen_t
、ttr
n_m_t
为瞬态热响应曲线ttr
n_cen
、第m条瞬态热响应曲线ttr
n_m
在第t帧的值;
[0027]
3.7)、在所有的瞬态热响应曲线ttr
n_cen
与m
n
条瞬态热响应曲线ttr
n_m
,m=1,2,...,m
n
的距离d
n_m
,m=1,2,...m
n
中找到最大值,记为d
n_max

[0028]
3.8)、判断瞬态热响应曲线ttr
n_cen
、第m条瞬态热响应曲线ttr
n_m
的距离d
n_m
是否大于ε
ttr
×
d
n_max
,如果大于,则认为第n个缺陷特征区域de的第m边缘点为属于热扩散区域的热扩散点,否则,认为是缺陷点,其中,ε
ttr
为设定的大于0小余1的距离系数;
[0029]
统计热扩散点的数量,得到第n个缺陷特征区域de的热扩散点数量p

n

[0030]
这样得到第n个缺陷特征区域de的实际像素点数量p

n
=p
n

p

n

[0031]
3.9)、计算第n个缺陷特征区域de的面积s
n

[0032][0033]
其中,l为检测区域的长度,b为检测区域的宽度,p
x
为图像长度方向的像素点数量,p
x
为图像宽度方向的像素点数量;
[0034]
3.10)、n个缺陷特征区域de1,de2,...,de
n
均按照步骤3.1)~3.9)进行处理,这样得到n个缺陷特征区域de1,de2,...,de
n
的面积s1,s2,...,s
n
,完成缺陷定量识别。
[0035]
本发明的发明目的是这样实现的:
[0036]
本发明用于航天器损伤红外检测的热扩散效应缺陷定量识别方法,首先,对红外重构图像在l*a*b*颜色空间进行聚类,保留符合缺陷对应的高亮颜色信息的一类的像素点,进行二值化、形态学开闭运算操作,得到红外重构缺陷分割图像以及各个缺陷特征区域的像素点数量;然后提取缺陷边缘像素点;最后对边缘像素点对应的瞬态热响应曲线与质心位置瞬态热响应曲线,基于温度变化序列构建的加权因子序列以及温度值最高的帧构建的权重序列计算其欧式距离,并据此判断边缘点是否是热扩散点,得到各个缺陷特征区域的热扩散点数量,在各个缺陷特征区域的像素点数量中扣除,得到实际像素点数量并转换为面积,完成缺陷定量识别。
[0037]
本发明通过l*a*b*颜色空间进行分类提取缺陷特征区域,能很好地将缺陷与背景区分出来。同时,结合瞬态热响应本身具有的物理属性设计缺陷特征区域的瞬态热响应相似性度量方法,并判断边缘瞬态热响应曲线是否属于热扩散区域,对于属于热扩散区域的瞬态热响应曲线进行剔除得到缺陷特征的实际像素点数量即缺陷定量信息,因而不需要进行实验拟合,并具有更好的定量识别结果。实验结果证明,对热扩散区域进行判定后的缺陷计算面积更接近于真实缺陷大小,验证了本发明的有效性。
附图说明
[0038]
图1是本发明用于航天器损伤红外检测的热扩散效应缺陷定量识别方法一张具体实施方式的流程图;
[0039]
图2是红外重构缺陷分割图像边缘提取示意图;
[0040]
图3是红外重构图像缺陷特征区域提取具体实例图,其中,(a)为作为被检测对象的人工碳纤维板实物图(检测背面),(b)为基于红外图像序列进行重构得到的红外重构图像,(c)提取得到的红外重构缺陷分割图像;
[0041]
图4是缺陷边缘像素点(边缘点)提取示意图,其中,(a)为标记了缺陷的红外重构缺陷分割图像,(b)为边缘轮廓图像;
[0042]
图5是缺陷5边缘点对应的瞬态热响应曲线,
[0043]
图6是缺陷5边缘点对应的瞬态热响应曲线进行比较后得到的加权因子序列图。
具体实施方式
[0044]
下面结合附图对本发明的具体实施方式进行描述,以便本邻域的技术人员更好地理解本发明。需要特别提醒注意的是,在以下的描述中,当已知功能和设计的详细描述也许
会淡化本发明的主要内容时,这些描述在这里将被忽略。
[0045]
在本实施例中,如图1所示,用于航天器损伤红外检测的热扩散效应缺陷定量识别方法包括以下步骤:
[0046]
步骤s1:红外重构图像缺陷特征区域提取
[0047]
基于l*a*b颜色空间进行分类,并将高亮色的缺陷特征区域提取出来作为后续处理对象,可以提高缺陷位置形态信息量化评估的准确度。缺陷特征区域提取的具体步骤为:
[0048]
步骤s1.1:颜色空间转换
[0049]
利用红外重构图像i(x
i
,y
j
),i=1,...,m

,j=1,...,n

((x
1_i
,y
1_j
)表示基准图像i第i列、第j行像素点的坐标,图像宽度为m

个像素点,高度为n

个像素点)中缺陷特征区域与背景颜色不同的特点将红外重构拼接图像进行基于l*a*b颜色空间的缺陷特征提取,将红外重构拼接图像i的像素值(温度特征值)从rgb颜色空间转换到数据处理速度最快的三维l*a*b*颜色空间,将每个红外重构拼接图像中表示温度特征像素点映射得到的“a*”和“b*”颜色信息转换值,构建m

*n

个颜色特征对象放入集合i
lab
(a
i

,b
i

),i

=1,...,m

*n


[0050]
步骤s1.2:聚类、保留丢弃、二值化、形态学开闭运算
[0051]
根据红外重构图像的“a*”和“b*”颜色信息转换值,利用聚类中心
[0052][0053]
将像素点聚类为k类,o
k
是某一聚类簇的聚类中心,其中c
k
是k类中某一聚类簇,n
k
是聚类簇c
k
中红外重构图像颜色特征簇个数。保留符合缺陷对应的高亮颜色信息的一类像素点,其余类丢弃,保留的则为缺陷特征区域图像;对缺陷特征区域图像进行二值化,将缺陷特征区域标记为高亮的白色,然后进行形态学开闭运算操作,来连接相邻的像素点,获得包含n个缺陷特征区域像素点坐标集de1,de2,...,de
n
的红外重构缺陷分割图像。
[0054]
统计n个缺陷特征区域de1,de2,...,de
n
的像素点数量,分别记为p1,p2,...,p
n

[0055]
本发明基于l*a*b颜色空间的红外重构图像缺陷特征提取,利用红外重构图像突出缺陷的特殊表征形式。将高亮颜色信息的缺陷特征区域分割提取出来作为后续处理对象,可以提高缺陷位置形态信息量化评估的准确度。
[0056]
步骤s2:缺陷边缘像素点提取
[0057]
利用边缘提取算法对红外重构缺陷分割图像进行边缘提取,得到边缘轮廓图像,提取每个缺陷特征区域的边缘点像素坐标,得到对应的n个边缘点像素坐标集pi1,pi2,...,pi
n

[0058]
如图2所示缺陷边界上热流受到缺陷内部和热传导区域的作用会形成边缘热扩散效应反映在热成像(红外图像序列s(i,j,t),i=1,...,m

,j=1,...n

,t=1,...,t,t为红外图像序列帧数)中,使得对使用红外图像序列进行重构后突出缺陷特征的红外重构图像进行直接分割后的缺陷部分,不能准确的反映缺陷的定量信息即分割出来的特征区域与缺陷实际尺寸不能准确对应,于是为了对热扩散区域进行判断,获得正确的缺陷特征区域像素点大小,需要对步骤s1获得的红外重构缺陷分割图像进行边缘提取,得到缺陷特征区域的边缘区域及对应的n个边缘点像素坐标集pi1,pi2,...,pi
n
,其中以特征区域1为例边缘点
像素集为式中(x
1_1
,y
1_1
)表示边缘点的像素坐标,利用边缘点的像素坐标可以提取红外重构图像序列中的瞬态热响应曲线用于相似性度量判断,确定其是否为热扩散点,进而获得准确的缺陷面积,完成缺陷定量的准确识别。
[0059]
步骤s3:根据缺陷热扩散区域判断定量识别缺陷
[0060]
步骤s3.1:求质心位置像素点坐标
[0061]
对于第n个缺陷特征区域de
n
,根据其对应的边缘点像素坐标集pi
n
得到其质心位置的像素点坐标(质心像素点坐标):(x
n_cen
,y
n_cen
)。
[0062]
步骤s3.2:获取质心位置、边缘点对应的瞬态热响应曲线
[0063]
根据质心像素点坐标(x
n_cen
,y
n_cen
),在红外重构图像对应的红外图像序列中获取对应的瞬态热响应曲线ttr
n_cen
,根据边缘点像素坐标集pi
n
中的边缘点像素坐标(x
n_m
,y
n_m
),在红外图像序列中获取对应的瞬态热响应曲线ttr
n_m
,m表示边缘点像素坐标集pi
n
的第m个边缘点,m=1,2,

,m
n
,m
n
为第n个缺陷特征区域de
n
对应的边缘点数量。
[0064]
步骤s3.3:获取温度变化序列
[0065]
对瞬态热响应曲线ttr
n_cen
、第m条瞬态热响应曲线ttr
n_m
分别求各帧(时刻)的单位温度变化率(同一瞬态热响应曲线相邻点温度差值),得到温度变化序列δv
n_cen
、δv
n_m

[0066]
瞬态热响应曲线反映了检测对象不同区域温度随着时间的变化特征。于是考虑温度根据时间的变化趋势,本发明对于瞬态热响应曲线通过温度变化率构建温度变化序列。
[0067]
步骤s3.4:获取加权因子序列
[0068]
对瞬态热响应曲线对应的温度变化序列δv
n_cen
、δv
n_m
各个对应时刻的温度变化率进行比较离散化比较结果,得到加权因子序列vt
n_m
,其中,对于第t帧的值vt
n_m_t
(即t帧的加权因子)为:
[0069]
vt
n_m_t
=|δv
n_cen_t
,δv
n_m_t
|
1,2,3
,t=1,2,...,t
‑1[0070]
vt
n_m_t
=1,t=t
[0071]
其中,δv
n_cen_t
、δv
n_m_t
分别为温度变化序列在第t帧的温度变化率值,t为红外图像序列的帧数,其(第t帧的值vt
n_m_t
)含义为:
[0072]
如果δv
n_cen_t
、δv
n_m_t
的差值小于变化阈值ε
δ
,即缺陷边缘与特征实际区域参考点瞬态热响应曲线相比当前温度点变换趋势相似,则第t帧的值vt
n_m_t
为1,如果δv
n_cen_t
、δv
n_m_t
的差值不小于变化阈值ε
δ
,并且同为正或同为负,即缺陷边缘与特征实际区域参考点瞬态热响应曲线相比在当前温度点变化趋势相同数值不同且差距较大,则第t帧的值vt
n_m_t
为2,如果δv
n_cen_t
、δv
n_m_t
的差值不小于变化阈值ε
δ
,并且正负相反,即缺陷边缘与特征实际区域参考点瞬态热响应曲线相比在当前温度点变化趋势相反,则第t帧的值vt
n_m_t
为3。
[0073]
本发明将质心像素点对应的瞬态热响应曲线与边缘点对应的瞬态热响应曲线进行比较,得到相应加权因子序列。然而,对于t维瞬态热响应曲线在求取单位温度变化率的条件下有t

1个对应计算温度变化率,因此把加权因子序列vt
n_m
最后一维vt
n_m_t
的设定为1。
[0074]
加权因子序列vt
n_m
的加权因子在一定程度上放大了瞬态热响应曲线变化的形态特征,当两条瞬态热响应曲线变化趋势相差越大对应的加权因子也会因此增大进而放大两者之间的距离,更好地对热扩散点进行了区分,得到相应的热扩散区域大小。
[0075]
步骤s3.5:构建权重序列
[0076]
构建第n个缺陷特征区域de
n
的权重序列ω
n_t
,对于其第t帧的值ω
n_t
,如果在瞬态热响应曲线ttr
n_cen
上,温度值最高的帧,取值为1.5,否则取值为1。
[0077]
对于瞬态热响应曲线,如果每一维(每一帧的值)都是平等看待,则不能准确区分质心像素点对应的瞬态热响应曲线与边缘点对应的瞬态热响应曲线,进而影响热扩散点的判断。
[0078]
对于不同的缺陷特征区域其峰值温度值得重点关注,于是对于瞬态热响应曲线ttr
n_cen
上,对温度值最高的帧,权重ω
n_t
取值为1.5,否则取值为1。于是对于加权欧式距离的权值有:
[0079][0080]
步骤s3.6:计算质心位置与边缘点对应的瞬态热响应曲线之间的距离
[0081]
计算瞬态热响应曲线ttr
n_cen
、第m条瞬态热响应曲线ttr
n_m
的加权欧式距离:
[0082][0083]
其中,ttr
n_cen_t
、ttr
n_m_t
为瞬态热响应曲线ttr
n_cen
、第m条瞬态热响应曲线ttr
n_m
在第t帧的值。
[0084]
步骤s3.7:获取距离最大值
[0085]
在所有的瞬态热响应曲线ttr
n_cen
与m
n
条瞬态热响应曲线ttr
n_m
,m=1,2,...,m
n
的加权欧式距离d
n_m
,m=1,2,...m
n
中找到最大值,记为d
n_max

[0086]
根据边缘点对应瞬态热响应曲线与质心像素点对应的瞬态热响应曲线之间的权重欧式距离的大小来确定是否为热扩散点。由于考虑每一个缺陷特征区域和热扩散区域的变化差值可能不同,无法确定一个统一的阈值,因此,本发明采用最大距离d
n_max
乘以一个距离系数来确定每一个缺陷特征区域的热扩散区域判断的具体阈值。
[0087]
步骤s3.8:判断边缘点是否是热扩散点
[0088]
判断瞬态热响应曲线ttr
n_cen
、第m条瞬态热响应曲线ttr
n_m
的距离d
n_m
是否大于ε
ttr
×
d
n_max
,如果大于,则认为第n个缺陷特征区域de的第m边缘点为属于热扩散区域的热扩散点,否则,认为是缺陷点,其中,ε
ttr
为设定的大于0小余1的距离系数;
[0089]
统计热扩散点的数量,得到第n个缺陷特征区域de的热扩散点数量p

n

[0090]
这样得到第n个缺陷特征区域de的实际像素点数量p

n
=p
n

p

n

[0091]
步骤s3.9:计算缺陷特征区域面积
[0092]
计算第n个缺陷特征区域de的面积s
n

[0093][0094]
其中,l为检测区域的长度,b为检测区域的宽度,p
x
为图像长度方向的像素点数量,p
x
为图像宽度方向的像素点数量;
[0095]
步骤s3.10:重复完成所有缺陷特征区域的缺陷定量识别
[0096]
n个缺陷特征区域de1,de2,...,de
n
均按照步骤3.1)~3.9)进行处理,这样得到n个
缺陷特征区域de1,de2,...,de
n
的面积s1,s2,...,s
n
,完成缺陷定量识别。
[0097]
实验验证
[0098]
为了说明本发明的可行性和正确性,下面分别通过实验验证来对本发明的技术方案进行说明。
[0099]
图3是红外重构图像缺陷特征区域提取具体实例图,其中,(a)为作为被检测对象的人工碳纤维板实物图(检测背面),(b)为基于红外图像序列进行重构得到的红外重构图像,(c)提取得到的红外重构缺陷分割图像。
[0100]
如图3(a)所示,在被检测对象的人工碳纤维板的背面做了8个缺陷。如图3(b)所示,红外重构图像中缺陷与背景比较呈现高亮,图3(b)所示,提取得到的红外重构缺陷分割图像显示缺陷特征区域被良好地提取出来,缺陷特征区域完整保留了缺陷特征,缺陷特征区域与检测背景区域实现了分离,准确地提取到了检测对象中缺陷对应的特征区域。通过对红外重构图像缺陷特征区域提取,实现了对红外重构图像缺陷特征区域的判断,去除了背景噪音,分割后得到的缺陷特征区域和边缘热扩散区域能反映缺陷区域的轮廓和形态分布,为接下来更好地进行的缺陷特征区域量化做好铺垫。将图3(c)所示的红外重构缺陷分割图像作为接下来的处理对象,经过二值化、形态学处理,形成连通域。
[0101]
如图4(a)所示。利用双阈值canny算法对红外重构缺陷分割图像进行边缘提取,获得如图4(b)所示的边缘轮廓图像以及对应的8邻域边缘点像素坐标集pi1,pi2,...,pi8。
[0102]
以缺陷5为例,加权因子序列构建过程中,边缘点像素坐标集pi5中边缘点对应的瞬态热响应曲线如图5所示,对应的加权因子如图6所示。
[0103]
瞬态热响应曲线ttr
5_cen
上在第95帧,温度值最高,ω
5_95
=1.5,其余的ω
5_t
为1。
[0104]
在所有的瞬态热响应曲线ttr
6_cen
与m
n
条瞬态热响应曲线ttr
5_m
,m=1,2,...,m5的距离d
5_m
,m=1,2,...m5中找到最大值,记为d
5_max
,进而得到比较阈值ε
ttr
×
d
5_max
=9.8955。
[0105]
对于缺陷5的热扩散点数量为53,于是对于缺陷5实际像素点数量为161。
[0106]
具体各个缺陷的热扩散点数量如表1所示,
[0107][0108][0109]
表1
[0110]
具体各个缺陷特征区域面积(去除热扩散点前)如表2所示
[0111][0112]
表2
[0113]
具体各个缺陷特征区域面积(去除热扩散点后)如表3所示
[0114][0115]
表3
[0116]
从表2、3可以看出,通过本发明用于航天器损伤红外检测的热扩散效应缺陷定量识别方法,定量识别的像素点数量更加逼近真实值,减少了横向热扩散带来的模糊边缘效应干扰,证明本发明具有更好的定量识别结果。
[0117]
尽管上面对本发明说明性的具体实施方式进行了描述,以便于本技术邻域的技术人员理解本发明,但应该清楚,本发明不限于具体实施方式的范围,对本技术邻域的普通技术人员来讲,只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内,这些变化是显而易见的,一切利用本发明构思的发明创造均在保护之列。
再多了解一些

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

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

相关文献

  • 日榜
  • 周榜
  • 月榜