一种基于高光谱图像比对的字画艺术品鉴定方法与流程

未命名 10-19 阅读:91 评论:0


1.本发明涉及字画鉴定技术,更具体地说,它涉及一种基于高光谱图像比对的字画艺术品鉴定方法。


背景技术:

2.文物艺术品检验与鉴定是经典物证鉴定的延伸扩展,它属于司法物证鉴定和文物保护两个学科的交叉领域。
3.当作赝造假蔚然成风并技艺日显高超之时,对以形为重的文物艺术品传统经验鉴定所构成的威胁当不言自明。在法学领域,物证司法鉴定因其服务于诉讼活动,有严格的科学性、客观性和中立性的要求,用于鉴定的技术方法追求系统性、科学性和客观性基础。同样文物艺术品物证鉴定技术体系必须尊重司法鉴定的本质要求,在借鉴传统经验鉴定方法的基础上,吸收基于现代科学技术手段的客观性判据,从而形成科学合理的文物艺术品物证鉴定技术体系。
4.高光谱成像技术对目标信息的采集是非接触的,满足了文物对无损分析的要求。同时利用光谱图像丰富的空间、光谱信息能够对文物进行更全面的记录与分析。这使得高光谱成像技术能够作为物证鉴定的客观性依据融入到文物鉴定方法之中,并在文物材料检验、痕迹检验、工艺特点和制作方式等细致而微的局部问题的鉴定起着重要的作用。
5.如四川博物院巩梦婷等用高光谱成像系统采集了现代中国画的高光谱图像,利用光谱角填图法对绘画所用颜料进行分类和识别。首都博物馆武望婷利用高光谱仪对《论道图》进行了隐含信息提取和图像增强,隐含信息和图像增强效果对实际工作中的书画鉴定和书画修复提供重要的依据。故宫博物院史宁昌等对故宫博物院的部分馆藏书画文物进行分析,发现高光谱成像技术在文字信息增强、隐藏信息提取、底稿线提取、颜料分析等方面有独特的优势,特别是短波红外波段不仅可以有效增强印章的文字信息,发现书画作品的涂改与涂沫痕迹,而且对于底稿线等线描特征的提取有着非常好的效果。西北大学彭进业等利用vis-nir成像光谱仪提取彩绘文物的图案。中国科学院张立福等对故宫博物院提供的常用国画颜料进行了光谱建库,将颜料分为红色、绿色、青色等6种颜料,光谱波长范围是400-2500nm。北京建筑大学丁新峰等对采集的文物颜料光谱与已知成分的标准颜料光谱进行比较,确定当代颜料与古代颜料最为接近的就是国画颜料,他们将张立福等建立的颜料光谱库与国外矿物颜料光谱库的数据结合起来分析古画中出现的颜料。国外方面,英国思克莱德大学adampolak利用高光谱成像获得绘画作品的高光谱影像,并利用颜料光谱库学习实现艺术品的真伪鉴定。文艺复兴时期画家埃尔
·
格列柯的作品以色彩运用见长,希腊克里特理工大学balas costas等建立了一个文艺复兴颜料光谱库,通过光谱匹配技术识别作品中的颜色。意大利罗马大学pronti等论证了利用vis-nir多光谱成像能够复原十四至十五世纪的羊皮纸上退化的文字。波兰科学院bartosz grabowski等采用端元光谱自动提取和混合光谱分解技术自动的分离画作中的颜料。
6.以高光谱成像为基础的物证鉴定体系包括成像技术、信息提取和鉴定标准与规范
构成。然而,对比于文件检验、违禁毒品鉴定和生物特征识别等发展比较成熟的司法鉴定,文物艺术品的高光谱物证鉴定尚缺乏发展。


技术实现要素:

7.本发明要解决的技术问题是针对现有技术的不足,提供一种基于高光谱图像比对的字画艺术品鉴定方法,该方法检验结果客观,准确度高、实用性好且易操作,能够有效的服务于物证鉴定和艺术品鉴定工作。
8.本发明所述的一种基于高光谱图像比对的字画艺术品鉴定方法,包括以下步骤,
9.步骤一、采集样本字画作品的高光谱图像;
10.步骤二、对所述高光谱图像进行辐射校正和高光谱图像去噪的预处理;
11.步骤三、对经预处理后的所述高光谱图像进行目标区域裁剪,构建具有高光谱子图像参考集和与高光谱子图像参考集相对应的近红外-红绿蓝多光谱子图像参考集的参考图像集b;
12.步骤四、通过自动端元提取和光谱差异比对法提取待检验图像a和所述参考图像集b中高光谱子图像参考集的端元,并进行端元光谱相似性比对;基于词袋模型对待检验图像a和所述参考图像集b中近红外-红绿蓝多光谱子图像参考集进行图像内容相似度对比;
13.步骤五、将端元光谱相似性比对结果和图像内容相似度对比结果进行融合,以对所述待检验字画图像块进行真伪判定。
14.通过以下公式对所述高光谱图像进行辐射校正预处理,
[0015][0016]
式中,r
l
(i,j)表示第l个波段像元位置(i,j)校正后的反射率数据,数值范围介于0-1;dn
l
(i,j)表示第l个波段像元位置(i,j)的字画高光谱图像;a(i,j)为采集像元位置(i,j)的白板图;d
l
(i,j)为自动快门采集像元位置(i,j)的暗信号;l为波段数,l为1-l中的任一数值;i=1,2,...n;j=1,2,...,m;n和m分别为高光谱图像的行数与列数。
[0017]
对所述高光谱图像进行高光谱图像去噪的预处理,具体为,
[0018]
首先利用自适应中值滤波法去除所述高光谱图像中的椒盐噪声,再利用分组主成分变换法对所述高光谱图像进行光谱维噪声去除,最后利用bm4d对所述高光谱图像进行图像维噪声去除。
[0019]
所述自适应中值滤波法具体包括,
[0020]
定义滤波窗口最大尺寸为s
max
,计算所述滤波窗口内的中值r
med
、最大值r
max
和最小值r
min
,并进入第一进程,
[0021]
第一进程:a1=r
med-r
min

[0022]
a2=r
med-r
max

[0023]
如果a1》0且a2《0,则转至第二进程,
[0024]
否则以设定的增加阈值增大滤波窗口尺寸,
[0025]
如果窗口尺寸小于或等于s
max
,则重复第一进程,
[0026]
否则输出中值r
med
替代r
l
(i,j);
[0027]
第一进程:b1=r
l
(i,j)-r
min

[0028]
b2=r
l
(i,j)-r
max

[0029]
如果b1》0且b2《0,则保持当前像元r
l
(i,j),
[0030]
否则输出输出中值r
med
替代r
l
(i,j)。
[0031]
所述自动端元提取和光谱差异比对法包括,
[0032]
图像端元提取步骤,用于确定所述待检验图像a和参考图像集b的高光谱子图像参考集的端元数p;并将所述待检验图像a和高光谱子图像参考集的每一个波段按逐行排列的方式转换为一个n维的行向量,记为影像x;提取所述影像x中的p个端元作为影像x的端元光谱集e;
[0033]
图像间的端元光谱差异比对步骤,用于根据所述待检验图像a的端元光谱集ea和高光谱子图像参考集的端元光谱集eb内的各端元光谱之间距离值构建矩阵d
a,b
,并根据所述矩阵d
a,b
中的元素获取所述端元距离的欧几里得距离de和余弦角ds;同时计算端元光谱集ea和端元光谱集eb的距离s(ea,eb);
[0034]
其中,p为1以上的自然数。
[0035]
所述提取影像x中的p个端元作为影像x的端元光谱集e,具体包括,
[0036]
s1、所述影像x经pca变换得到p-1维特征影像
[0037]
s2、在所述p-1维特征影像的第一成分中随机选择像元r,并遍历所述p-1维特征影像中所有的像元xi(i=1,2,...,n),计算行列式绝对值将所述行列式值最大的像元作为第一个端元,并记录所述第一个端元的像元坐标为t1;
[0038]
s3、假设(t1,...,tk)(k≥1)是已经提取的端元的像元坐标,从所述p-1维特征影像的前k个主成分特征影像中取出对应的端元主成分特征,对任一像元定义单形体体积,并通过如下单形体体积公式进行计算,
[0039][0040]
采用所述单形体体积公式对每个像元进行计算,以获取相应的单形体体积,并得到单形体体积最大的像元t
k+1
,同时记录端元的像元坐标为(t1,...,tk,t
k+1
);
[0041]
s4、如果k+1《p,则k=k+1,并转到s3;否则(t1,...,tk,t
k+1
)记录的即为最后端元的图像位置;从所述影像x对应位置提取的像元光谱即为端元光谱集e={e1,e2,...,e
p
}。
[0042]
在所述图像间的端元光谱差异比对步骤中,
[0043]
所述矩阵d
a,b
为,
[0044][0045]
其中,d
a,b
的任一元素为d
i,j
,i=1,...,pa;j=1,...,pb;d
i,j
表示端元光谱集ea内端元和eb内端元之间的端元距离;
[0046]
通过以下两个公式分别获取所述欧几里得距离de和余弦角ds,
[0047]de
(e1,e2)=||e
i-e2||,
[0048][0049]
上式中,e1和e2分别表示端元和端元的端元向量,e1·
e2表示e1和e2这两个端元向量之间的点积,|| ||表示2范数;
[0050]
通过以下公式计算端元光谱集ea和端元光谱集eb的距离,
[0051][0052]
式中,式中,
[0053]
基于词袋模型对待检验图像a和所述参考图像集b中近红外-红绿蓝多光谱参考子图像集进行图像内容相似度对比,具体包括,
[0054]
特征提取和描述步骤,用于对所述待检验图像a和近红外-红绿蓝多光谱参考子图像集进行近红外、红色、绿色、蓝色通道的重采样,以获得近重采样多光谱图像;对所述重采样多光谱图像进行超像素分割,以超像素的特征值表达所述重采样多光谱图像的颜色特征;对所述待检验图像a和近红外-红绿蓝多光谱参考子图像集进行尺度空间极值点检测、特征点定位、确定主方向和生成特征点描述处理,以获得surf特征;
[0055]
视觉词典生成步骤,用于对所述重采样多光谱图像的颜色特征和surf特征进行聚类,将聚类的中心定义为视觉单词,并将所有的视觉单词进行组合以构成视觉词典;
[0056]
图像量化步骤,用于将所述近红外-红绿蓝多光谱参考子图像集和待检验图像a中的每一个颜色特征、surf特征与所述视觉词典中的每个特征向量进行比较,并将所述颜色特征、surf特征近似成与其最相似的词典项,然后计算所有所述颜色特征、surf特征在每个词典项出现的次数,得到词典直方图;将所述近红外-红绿蓝多光谱参考子图像集和待检验图像a分别量化成颜色视觉单词直方图和surf视觉单词直方图;
[0057]
词袋特征相似性对比步骤,用于对所述颜色视觉单词直方图和surf视觉单词直方图中的词袋特征进行相似性对比,以获取所述surf视觉单词直方图中surf视觉特征相似值和颜色视觉单词直方图中的颜色视觉特征相似值。
[0058]
在所述词袋特征相似性对比步骤中,视觉特征相似值通过以下公式获取,
[0059][0060]
式中,ha和hb为待比较的两个词袋特征向量。
[0061]
通过下式获取所述surf视觉特征相似值和颜色视觉特征相似值的加权相似性特征,
[0062]
xf=w
shape
*x
shape
+w
color
*x
color

[0063]
式中,x
shape
为surf视觉特征相似值;x
color
为颜色视觉特征相似值;
[0064]wshape
+w
color
=1;
[0065]
根据所述加权相似性特征将近红外-红绿蓝多光谱子图像参考集中的子图像进行
排序,并将与所述待检验图像a最相似的子图像选取出来作为最相似图像t;
[0066]
获取所述待检验图像a和最相似图像t的端元光谱集ea和端元光谱集e
t
的距离s(ea,e
t
),将所述距离s(ea,e
t
)与设定的端元光谱相似性阈值进行比较,若所述距离s(ea,e
t
)小于设定的端元光谱相似性阈值,则判定所述待检验图像代表的字画作品为真。
[0067]
有益效果
[0068]
本发明的优点在于:
[0069]
1.本发明涉及了高光谱数据采集、预处理、图像特征提取和字画相似性比对与真伪判定,系统的构建了一种字画艺术品检验的技术。
[0070]
2.本发明采用光谱维、空间维结合的方式有效的去除了采集到的高光谱影像的噪声,计算效率高。
[0071]
3.本发明构建具有辨识意义的字画图像感兴趣区域样本数据集,包括了高光谱子图像和对应的近红外-红绿蓝多光谱子图像参考集,字画图像比对从字画材质的光谱和绘画内容两个方面进行比对检验,检验联合了光谱、视觉颜色和形态特征,不受图像尺寸、旋转、光照和噪声的影响。
[0072]
4.本发明采用自动端元提取的技术获得图像中代表的端元光谱,将这些端元光谱当作反映字画材质的“指纹”,并通过计算待检验图像和参照图像的端元集的光谱距离反映图像间的相似性。
[0073]
5.本发明采用词袋模型描述字画图像内容,即颜色特征词袋和surf特征词袋,通过比较词袋特征的相似性反映图像的相似性,两种词袋特征采用决策级加权的方式进行融合。
[0074]
6.本发明采用决策加权融合的方式融合视觉形态和颜色特征,采用分步判定的方式将视觉特征与光谱特征进行结合得到最终的字画真伪判定。
附图说明
[0075]
图1为本发明的鉴定方法处理与分析流程图;
[0076]
图2为本发明的高光谱图像预处理流程图;
[0077]
图3为本发明的nir-rgb多光谱重采样光谱响应函数图;
[0078]
图4为褚遂良楷书卷.清作品样本图;
[0079]
图5为褚遂良楷书卷.清作品样本基于本发明的鉴定方法采集到的局部图像;
[0080]
图6为褚遂良楷书卷.清作品样本中印章局部基于本发明的鉴定方法去噪处理前后对比示意图;
[0081]
图7为褚遂良楷书卷.清作品样本基于本发明的鉴定方法利用分组pca光谱域变换去噪光谱反射率对比示意图;
[0082]
图8为褚遂良楷书卷.清作品样本基于本发明的鉴定方法得到的三类子图像数据示例图;
[0083]
图9为褚遂良楷书卷.清作品样本基于本发明的鉴定方法获取的部分参考子图像;
[0084]
图10为图9中的参考子图像的端元光谱图;
[0085]
图11为本发明所述的待检验图像及其端元光谱图;
[0086]
图12为褚遂良楷书卷.清作品样本基于本发明的鉴定方法中两种度量方式得到的
前10个比对相似的图像区域及其端元光谱;
[0087]
图13为本发明所述的待检验图像超像素分割与surf特征示意图;
[0088]
图14为褚遂良楷书卷.清作品样本基于本发明的鉴定方法实现的颜色和surf视觉特征相似性比对结果图;
[0089]
图15为褚遂良楷书卷.清作品样本基于本发明的鉴定方法实现的颜色和surf视觉特征相似性比对结果光谱角对比图。
具体实施方式
[0090]
下面结合实施例,对本发明作进一步的描述,但不构成对本发明的任何限制,任何人在本发明权利要求范围所做的有限次的修改,仍在本发明的权利要求范围内。
[0091]
字画检验与鉴定以分析字画的材料组成和书写内容,对字画进行定性或定量的评价为主要内容。
[0092]
本发明提供的一种基于高光谱图像比对的字画艺术品鉴定方法,该方法利用自定义的“可靠”高光谱图像集对待检验的的字画进行鉴定,包括了光谱数据采集、高光谱图像预处理、“可靠”高光谱图像集构建、光谱和书写内容的相似性检验、字画评价。该方法联合了光谱、颜色和形态特征进行图像内容的比对检验,不受图像尺寸、方向、噪声和光照的影响,检验评价结果客观、可靠性好、实用性强。
[0093]
如图1所示,本发明主要包括了五个步骤:
[0094]
步骤一是采集样本字画作品的高光谱图像。步骤二是对高光谱图像进行辐射校正和高光谱图像去噪的预处理。步骤三是构建具有辨识意义的字画图像感兴趣区域数据集,具体是对经预处理后的高光谱图像进行目标区域裁剪,构建具有高光谱子图像参考集和与高光谱子图像参考集相对应的近红外-红绿蓝多光谱子图像参考集的参考图像集b。步骤四是通过自动端元提取和光谱差异比对法提取待检验图像a和参考图像集b中高光谱子图像参考集的端元,并进行端元光谱相似性比对;基于词袋模型对待检验图像a和参考图像集b中近红外-红绿蓝多光谱子图像参考集进行图像内容相似度对比,其实施内容主要是基于参考图像集对待检验字画图像块进行相似性比对,采用联合光谱、颜色和surf词袋特征结合的图像内容相似性比对方法。步骤五是多特征融合的字画真伪判定,具体为将端元光谱相似性比对结果和图像内容相似度对比结果进行融合,以对待检验字画图像块进行真伪判定。
[0095]
各步骤的具体内容如下。
[0096]
步骤一、采集样本字画作品的高光谱图像。其包括两个子步骤,分别是步骤101和步骤102。
[0097]
步骤101、关于数据采集设备与采集参数。本发明采用高光谱凝视相机对字画图像成像,相机采用传动设备内置的方式进行推扫式成像,从而避免磨损字画;相机的光谱波长470-900nm,覆盖可见光和近红外范围;它的平均光谱分辨率为3nm,辐射分辨率为10bits;成像的光源为卤素灯,通过调电压的方式控制灯源功率,为不损伤字画并保证光照均匀,一对卤素灯以50w低功率输出,以左右对称高入射角的方式照射目标物。
[0098]
步骤102、分别采集字画高光谱图像dn
mxnxl
和白板图像a
mxnxl
,图像采集范围约为20cm*20cm,从检验的角度考虑对图幅比较大的物品仅采集局部关键区域。需要说明的是,
在本文中,符号“*”表示数学运算乘以。
[0099]
步骤二、高光谱图像预处理。如图2所示,本步骤包括两方面的处理,即高光谱图像光谱反射率校正和高光谱图像去噪,其中图像去噪采用分步处理的方式完成。
[0100]
步骤201、高光谱图像光谱反射率校正,即辐射校正。反射率校正的目的是对采集得到的数据进行标准化处理,反射率反映字画材料的光学反射特性。本发明利用标准白板做参照对高光谱图像进行光谱反射校正,过程如下:
[0101]
假设成像光谱仪采集的是幅宽mxn,波段数为l的字画高光谱图像,记为dn
mxnxl
,采集的白板图像为a
mxnxl
,自动快门采集的暗信号为d
mxnxl
,则校正后的反射率数据为r
l
(i,j),其表示为式1。
[0102][0103]
式中,r
l
(i,j)表示第l个波段像元位置(i,j)校正后的反射率数据,数值范围介于0-1;l为1-l中任一数值。
[0104]
步骤202、高光谱图像去噪。对图像采集过程中出现的图像噪声进行处理,本发明去噪主要针对图像中包含的椒盐噪声和其它噪声,采用分步处理的方式,即首先利用自适应中值滤波去除椒盐噪声,再利用分组主成分变换进行光谱维噪声去除,最后利用bm4d进行图像维噪声去除。
[0105]
步骤2021、高光谱图像自适应中值滤波法去除椒盐噪声。具体如下。
[0106]
对l个波段的高光谱反射率图像采用逐波段处理的方式,若当前处理波段为第l波段r
l
(l=1,2,...,l),且当前像元r
l
(i,j),定义滤波窗口最大尺寸为s
max
,自适应中值滤波处理如下:
[0107]
计算所述滤波窗口内的中值r
med
、最大值r
max
和最小值r
min
,并进入第一进程,
[0108]
第一进程:a1=r
med-r
min

[0109]
a2=r
med-r
max

[0110]
如果a1》0且a2《0,则转至第二进程,
[0111]
否则以设定的增加阈值增大滤波窗口尺寸,
[0112]
如果窗口尺寸小于或等于s
max
,则重复第一进程,
[0113]
否则输出中值r
med
替代r
l
(i,j);
[0114]
第一进程:b1=r
l
(i,j)-r
min

[0115]
b2=r
l
(i,j)-r
max

[0116]
如果b1》0且b2《0,则保持当前像元r
l
(i,j),
[0117]
否则输出输出中值r
med
替代r
l
(i,j)。
[0118]
以上过程中局部的最大值r
max
和最小值r
min
是可能的椒盐噪声,处理方式是逐个波段图像逐个像元执行该过程,其优点是在去除椒盐噪声的同时保留原光谱数据。
[0119]
步骤2022、分组主成分变换重构光谱维去噪。具体如下。
[0120]
主成分分析(pca-principle component analysis)变换后各主成分分量彼此不相关,且随着主成分编号的增加该分量包含的信息量减小。对高光谱图像进行pca变换后,大部分信息集中在前面的主成分分量中,其它的主成分分量以噪声为主,称之为噪声成分分量。取前若干信息量大的主成分分量而丢掉信息量小的噪声成分分量,进行pca逆变换重
构高光谱图像。以下为采用分组pca变换分析的方法对高光谱图像进行光谱维去噪的具体步骤。
[0121]
第一步、计算波段分组。假设含噪声的高光谱图像空间维大小为mxn,波段数为l,将其表达为lxn(n=mxn)的矩阵x。其中,x=[x1,x2,...,xn]。计算相关矩阵r=(x*x
t
)/n,并根据相关矩阵的相关系数确定波段分组。假设波段分组断点为[b1,b2,...,bt],表示对高光谱图像在光谱维按波段顺序分成t-1组。
[0122]
第二步、分组光谱为pca变换。记xc=[x
1c
,x
2c
,...,x
nc
]为第c个波段分组,波段范围是[bc,bc+1],波段数为lc,xc为均值向量,∑=e((x
ic-xc)*(x
ic-xc)
t
)为协方差矩阵,对协方差矩阵进行特征值分解∑=vλv
t
,a=diag(λ1,λ2,...,λ
lc
)是协方差矩阵∑的非0特征值按从大到小排列组成的对角矩阵,特征向量矩阵v=[v1,v2,...,v
lc
]顺序对应特征值,对xc的pca变换如式2所示。
[0123][0124]
y的各特征之间是不相关的,大部分信息集中在前面的主成分分量中,其它的分量以噪声为主,设保留的前l个特征值的权重占比等于取前l信息量大的主成分分量而丢掉信息量小的噪声成分分量,v
*
=[v1,v2,...,v
l
],通过下式重构,重构的均方误差为
[0125][0126]
第三步、重组分组重构后的各波段数据。重构后的高光谱图像矩阵为
[0127]
步骤2023、基于bm4d的图像域去噪。
[0128]
四维块匹配滤波(block-matching 4d filtering,bm4d)是三维图像去噪算法,由于高光谱图像是由二维空间信息和一维光谱信息组成的三维图像,因此bm4d算法适用于高光谱图像去噪处理。本发明利用开源bm4d进行去噪,但并不是直接利用,而是对经光谱域去噪并经光谱域重采样得到近红外-红绿蓝四通道多光谱图像进行去噪,并将去噪后的图像应用于突出空间维特征的图像比对中。
[0129]
步骤三、构建具有辨识意义的字画图像感兴趣区域数据集。
[0130]
高光谱图像经预处理后裁剪感兴趣区域子图像构建一个“可靠”高光谱图像参考图像集b。该图像集由具有辨识字画作者或绘画时代信息作用的子图像组成。本发明将子图像分为三类,即反映字画本底(纸或绢等)的留白区域,反映字画色彩和墨质的绘制区域(特别是反映作者习惯的用色和用笔区域),和作者印章子图像。根据检验内容的差异,数据集由两部分构成,分别是高光谱子图像参考集和对应的近红外-红绿蓝多光谱子图像参考集。前者用于图像块端元光谱比对,后者则是用于基于词袋的比对。
[0131]
步骤四、通过自动端元提取和光谱差异比对法提取待检验图像a和所述参考图像集b中高光谱子图像参考集的端元,并进行端元光谱相似性比对;基于词袋模型对待检验图像a和所述参考图像集b中近红外-红绿蓝多光谱子图像参考集进行图像内容相似度对比。
[0132]
为了对字画进行有效的检验与鉴定,并同时提高检验的效率,图像比对的方式采用对比关键的图像区域而非成像仪采集的整个书画图像。因而,在具体实现时由用户选择具有辨识意义的感兴趣区域,从字画材质的光谱和绘画(书写)内容两个方面与参考图像集进行比对。以下分别说明基于图像端元光谱和基于词袋模型的图像比对方法。
[0133]
步骤401、字画端元光谱相似性比对。端元光谱是高光谱图像中代表纯物质的光谱形式,它被看成是反映字画材质的“指纹”,因而比对待检验图像与参照图像间的端元光谱的相似性成为材质比对的实现的基本途径。本发明采用自动端元提取技术提取待检验图像和参照图像集的端元,并定义了反映图像块差异的端元光谱距离。其中自动端元提取包含端元数确定和端元光谱自动提取两步。
[0134]
步骤4011、图像端元提取。
[0135]
第一步、端元数的确定。端元数p是根据字与画的图像的复杂性来确定。对于字帖作品而言,其选取的感兴趣图像中的内容主要是字或印章,以及字帖的纸张或裱材,因而绘制区域子图像的端元数通常定义为2;而留白区域的端元数为1;对于绘画作品而言,绘画区域的端元数目则根据作者用色的习惯来确定,留白区域的端元数为1。裱材区则根据一般的材料组成来定义。
[0136]
第二步、端元光谱自动提取。留白区域取图像的平均光谱作为留白端元,字画绘制区域的端元光谱自动提取自动提取过程如下:
[0137]
假设高光谱图像(此处记为r)幅宽为m行n列,波段数为l。将高光谱图像的每一个波段按逐行排列的方式转换为一个n(n=mxn)维的行向量,即提取x中p个端元的步骤如下:
[0138]
s1、影像x经pca变换得到p-1维特征影像
[0139]
s2、在所述p-1维特征影像的第一成分中随机选择像元r,并遍历所述p-1维特征影像中所有的像元xi(i=1,2,...,n),计算行列式绝对值将所述行列式值最大的像元作为第一个端元,并记录所述第一个端元的像元坐标为t1;
[0140]
s3、不失一般性地假设(t1,...,tk)(k≥1)是已经提取的端元的像元坐标,从所述p-1维特征影像的前k个主成分特征影像中取出对应的端元主成分特征对任一像元定义单形体体积,并通过式4进行计算单形体体积。
[0141][0142]
采用式4对每个像元进行计算,以获取相应的单形体体积,并得到单形体体积最大的像元t
k+1
,即记录端元的像元坐标为(t1,...,tk,t
k+1
)。
[0143]
s4、如果k+1《p,则k=k+1,并转到s3;否则(t1,...,tk,t
k+1
)记录的即为最后端元的图像位置;从所述影像x对应位置提取的像元光谱即为端元光谱集e={e1,e2,...,e
p
}。
[0144]
步骤4012、图像间的端元光谱差异比对。
[0145]
假设待检验图像a和参照图像集b提取的端元光谱集分别为ea和eb,它们的端元数分别为pa和pb,定义端元光谱集ea和eb的内各端元光谱之间距离值构成矩阵d
a,b
。矩阵d
a,b
表示为式5。
[0146][0147]
其中,d
a,b
的任一元素为d
i,j
,i=1,...,pa;j=1,...,pb;d
i,j
表示端元光谱集ea内端元和eb内端元之间的端元距离。
[0148]
用两种形式定义该端元距离,分别是欧几里得距离de和余弦角ds。
[0149]de
(e1,e2)=||e
1-e2||
ꢀꢀ
(式6)
[0150][0151]
上式中,e1和e2分别表示端元和端元的端元向量,e1·
e2表示e1和e2这两个端元向量之间的点积,|| ||表示2范数。
[0152]
欧几里得距离de和余弦角ds都具有相同的意义,即数值越小越表明两个端元相似。因而,定义端元光谱集ea和eb的距离如式8所示。
[0153][0154]
式8中,
[0155][0156][0157]
式9和式10分别代表端元光谱集ea内端元与端元光谱集eb内所有端元之间距离的最小值,和端元光谱集eb内端元与端元光谱集ea内所有端元之间距离的最小值。
[0158]
以s(ea,eb)表达待检验图像a和参照图像集b中子图像的相似性,该值越小越表明两者相似。在参照图像集b所有图像中,以比对的端元光谱相似性结果解释待检验图像a。
[0159]
步骤402、基于词袋描述图像内容的图像比对。词袋模型最初应用于文本处理领域,用来对文档进行分类和识别。其基本原理为将文档看作是无序的关键词的集合,通过统计每个关键词在单个文档中出现的频率来对文档进行向量表示,从而进行文档识别。对应到本发明中,应用词袋模型进行近红外-红绿蓝多光谱子图像集比对共包括以下4个步骤,分别为特征提取和描述、生成视觉词典、图像量化、词袋特征相似性比对,具体如下。
[0160]
步骤4021、特征提取和描述。特征提取和描述环节的主要任务是从图像中抽取具有描述图像内容的特征,本发明中采用包含空间信息的颜色和surf(speeded up robust features)特征描述图像内容。其中颜色特征具有全局性质,surf特征形态描述能力强,有局部性质。
[0161]
第一步、颜色特征提取。
[0162]
颜色特征组成表达为(nirn,rn,gn,bn,xn,yn),它们表示归一化的四波段多光谱值与归一化的坐标值结合的特征。根据采集的字画高光谱图像特点,这里的颜色图像特征是由高光谱反射率图像的重采样获得。引用landsat 8 oli多光谱响应函数中的nir、red、green和blue四个光谱响应函数进行光谱重采样,它们的中心波长分别是860nm、650nm、560nm、480nm,如附图3所示。假设red通道的光谱响应函数为f
red
(x),对470nm-900nm的任一像元光谱r={r
470
,...,r
900
}重采样的红色通道值为:
[0163][0164]
类似的绿色、蓝色和近红外通道的图像可按此方法重采样得到,从而得到近红外-红绿蓝(重采样)多光谱图像。
[0165]
接着对光谱重采样得到的重采样多光谱图像进行超像素分割,每一幅图像的颜色特征用超像素的特征值表达,即颜色特征的数目等于超像素数,每一个归一化颜色特征(nirn,rn,gn,bn,xn,yn)的中的颜色值是超像素内像元的颜色平均值,由于颜色值是光谱反射率重采样得到,所以其值介于0-1之间;坐标也是超像素内像元的横纵坐标平均值除以图像的行数和列数。
[0166]
第二步、surf特征提取。
[0167]
surf特征是一种快速特征点检测算法,该算子采用卷积方式提取特征,在保持sift算子缩放、旋转不变和亮度保持不变性等优势的同时,surf算子具有仿射不变和运算性能高等优点。其特征提取过程主要分为尺度空间极值点检测、特征点定位、确定主方向和生成特征点描述符四个步骤。
[0168]
step 1、极值点检测。surf选择图像尺度空间中的极值点作为候选特征点。为了具有尺度不变性,surf算法用不同尺寸的方框滤波器对原始图像进行滤波处理,组成图像金字塔。对图像金字塔中的每一层使用hessian矩阵进行极值点检测。图像i(x,y)中的点x在尺度σ处的hessian矩阵定义如式(12)所示。
[0169][0170]
其中,l
xx
(x,σ)是高斯函数二阶偏导数和图像的二维卷积,l
xy
(x,σ)和l
yy
(x,σ)的含义与之类似。
[0171]
step 2、特征点定位。根据hessian矩阵求出尺度图像在(x,y,σ)处的极值后,先在极值点的3*3*3立体邻域内进行非极大值抑制,只有上下尺度及本尺度周围的26个邻域值都大或者都小的极值点,才能作为候选特征点,然后在尺度空间和图像空间中进行插值运算,得到稳定特征点位置及所在尺度值。
[0172]
step 3、为了保持旋转不变性。首先以特征点为中心,计算半径为6s(s为特征点的尺度值)邻域内的点在水平和垂直方向上的harr小波响应,然后给这些响应值赋予高斯权重系数,接着将60
°
范围内的响应累加形成新的矢量,最后遍历整个圆形区域,选择最长矢量方向作为特征点的主方向。
[0173]
step 4、生成特征描述符。以特征点为中心,将坐标轴旋转到主方向,按主方向选
取边长为20s(采样步长)的正方形区域,将该窗口区域划分成4*4的子区域,计算5s*5s范围内的小波响应。相对于主方向的水平、垂直方向的haar小波响应为d
x
、dy,同样赋予响应值系数,然后将每个子区域的响应系数及其绝对值相加形成矢量v=(∑d
x
,∑dy,∑|d
x
|,∑|dy|)。因此,对每一特征点生成4*(4*4)=64维的特征描述向量,再进行向量的归一化,从而对光照具有一定的鲁棒性。
[0174]
步骤4022、生成视觉词典。
[0175]
采用k-means聚类方法对由“特征提取和描述”步骤得到的大量单词
‑‑
颜色特征和surf特征分别聚类。将聚类的中心定义为视觉单词,所有的视觉单词进行组合则构成视觉词典,词典分为颜色特征视觉词典和surf视觉词典。两个视觉词典的大小分别为k1和k2,即视觉单词的个数,需要根据经验提前定义。
[0176]
步骤4023、图像量化。
[0177]
将待检验图像a和参考图像集b都表示成词典的直方图分布,即一个长度为k的向量。具体的,将所述近红外-红绿蓝多光谱参考子图像集和待检验图像a中的每一个颜色特征、surf特征与所述视觉词典中的每个特征向量进行比较,并将所述颜色特征、surf特征近似成与其最相似的词典项,然后计算所有所述颜色特征、surf特征在每个词典项出现的次数,得到词典直方图;将所述近红外-红绿蓝多光谱参考子图像集和待检验图像a分别量化成颜色视觉单词直方图和surf视觉单词直方图。
[0178]
步骤4024、词袋特征相似性比对。
[0179]
对所述颜色视觉单词直方图和surf视觉单词直方图中的词袋特征进行相似性对比,以获取所述surf视觉单词直方图中surf视觉特征相似值和颜色视觉单词直方图中的颜色视觉特征相似值。其中,视觉特征相似值通过式13获取。
[0180][0181]
式中,ha和hb是待比较的两个词袋特征向量。待检验图像a和参考图像集b内图像的距离越小,表明两者越相似。
[0182]
步骤五、多特征融合的字画真伪判定
[0183]
图像端元光谱反映字画材质,surf视觉特征反映图像的关键结构特征,颜色视觉特征反映空间布局特征,这三种特征互补。本发明结合这三种特征的比对结果评价图像被比对的图像真伪性质。具体方法如下。
[0184]
通过式14获取所述surf视觉特征相似值和颜色视觉特征相似值的加权相似性特征。
[0185]
xf=w
shape
*x
shape
+w
color
*x
color
ꢀꢀ
(式14)
[0186]
式中,x
shape
为surf视觉特征相似值;x
color
为颜色视觉特征相似值;
[0187]wshape
+w
color
=1。
[0188]
根据所述加权相似性特征将近红外-红绿蓝多光谱子图像参考集中的子图像进行排序,并将与所述待检验图像a最相似的子图像选取出来作为最相似图像t。
[0189]
获取所述待检验图像a和最相似图像t的端元光谱集ea和端元光谱集e
t
的距离s(ea,e
t
)。s(ea,e
t
)表示待检验图像a和图像t的端元光谱相似性,设定σ为端元光谱相似性阈
值。若式15成立,则可以判定待比对图像代表的字画作品为真。
[0190]
s(ea,e
t
)《σ
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(式15)
[0191]
为了验证上述字画艺术品鉴定方法的可行性,本实施例基于上述五大步骤,以某地的艺术博物馆中的褚遂良楷书卷.清作品,该作品纵:23cm,横:159cm,样品质地为纸,作为样本字画作品,进行鉴定。
[0192]
步骤一、利用snapscan vnir高光谱成像系统(470-900nm)对样品采用分块采集的方式,采集像素分辨率为3600*2048,在470-900nm范围内含有连续的150个光谱通道,样品如图4,分块采集部分图如图5。
[0193]
步骤二、对采集得到的图像进行预处理,即图像反射率校正和去噪。图6为经自适应中值滤波处理和bm4d去噪处理前后的图像,图7为利用分组pca光谱域变换去噪光谱反射率对比。从图6中可以明显的看到采集的图像存在噪声,经去噪处理后图像的同秩区域变得更光滑;而光谱随着保留的特征值权重占比减小而变得更光滑,但细节保留的则会减少。
[0194]
步骤三、对样本字画作品构建具有辨识意义的字画图像数据集,具体为在样本字画作品的成像光谱图像中利用roi工具选取子图像,子图像分为三类,即字画留白类、墨质文字类和印章类,如图8所示。图9为部分参考子图像,这些子图像有高光谱和近红外-红绿蓝多光谱图像两种形式。
[0195]
步骤四、基于光谱和书写内容相似性检验的图像块比对。具体为:
[0196]
步骤401、字画端元光谱相似性比对。
[0197]
步骤4011、首先获取得到高光谱参考图像集和待检验图像的端元光谱。每幅图像的端元光谱数目1-4不等,具体根据图像内容进行自定义,图10是图9对应的高光谱图像端元光谱。每一个子图像的端元数有区别,同类图像的端元数接近。
[0198]
步骤4012、图像间的端元光谱差异比对。以图8中的印章例作为待检验图像,该图像的端元光谱如图11,分别利用欧式距离和光谱余弦角进行图像间端元光谱集比对。实施例结果如图12,两种度量方式均得到了前10个比对相似的端元光谱及其图像,这些图像以印章图像为主,其中包括目标参考图像。
[0199]
步骤402、基于词袋描述图像内容的图像比对。
[0200]
步骤4021、特征提取和描述。具体为,对图像的进行超像素分割的方法是slic(simple linear iterative clustering),每一个颜色特征(nirn,rn,gn,bn,xn,yn)均由超像素内所有像元的均值计算得到,超像素分割使得各个子图像间具有相对稳定的特征数目;而形态特征以surf特征点表示,这类特征点在保持sift算子缩放、旋转不变和亮度保持不变性等优势的同时,具有仿射不变和运算性能高等优点。图13为超像素分割和surf特征提取的示例。
[0201]
步骤4022-步骤4024是以颜色和surf词袋特征进行图像比对的,图14是基于两种特征将待检验图像与参考图像进行比对后分别得到的前16幅相似的图像,可以明显的看到比对后的结果中正确的找到了合适的参考图像子图。
[0202]
步骤五、多特征融合的字画真伪判定。图14中对比了等权相加比对与单特征的比对结果,加权融合的方式对单特征的结果进行了适当调整。为进一步验证,将待检验图像进行仿射变换,并计算三种方式比对得到的相似子图像与待检验图像之间的光谱相似性,按相似性距离平均值进行对比,如图15。从图中可以看到几何变换前后,三种方式都能正确比
对到参照图像,加权方式得到的前10个相似子图像的光谱角平均值总体要小于其它两种方式,且角度值平均值小于3
°
,比较符合预期。
[0203]
以上所述的仅是本发明的优选实施方式,应当指出对于本领域的技术人员来说,在不脱离本发明结构的前提下,还可以作出若干变形和改进,这些都不会影响本发明实施的效果和专利的实用性。

技术特征:
1.一种基于高光谱图像比对的字画艺术品鉴定方法,其特征在于,包括以下步骤,步骤一、采集样本字画作品的高光谱图像;步骤二、对所述高光谱图像进行辐射校正和高光谱图像去噪的预处理;步骤三、对经预处理后的所述高光谱图像进行目标区域裁剪,构建具有高光谱子图像参考集和与高光谱子图像参考集相对应的近红外-红绿蓝多光谱子图像参考集的参考图像集b;步骤四、通过自动端元提取和光谱差异比对法提取待检验图像a和所述参考图像集b中高光谱子图像参考集的端元,并进行端元光谱相似性比对;基于词袋模型对待检验图像a和所述参考图像集b中近红外-红绿蓝多光谱子图像参考集进行图像内容相似度对比;步骤五、将端元光谱相似性比对结果和图像内容相似度对比结果进行融合,以对所述待检验字画图像块进行真伪判定。2.根据权利要求1所述的一种基于高光谱图像比对的字画艺术品鉴定方法,其特征在于,通过以下公式对所述高光谱图像进行辐射校正预处理,式中,r
l
(i,j)表示第l个波段像元位置(i,j)校正后的反射率数据,数值范围介于0-1;dn
l
(i,j)表示第l个波段像元位置(i,j)的字画高光谱图像;a(i,j)为采集像元位置(i,j)的白板图;d
l
(i,j)为自动快门采集像元位置(i,j)的暗信号;l为波段数,l为1-l中的任一数值;i=1,2,...n;j=1,2,...,m;n和m分别为高光谱图像的行数与列数。3.根据权利要求2所述的一种基于高光谱图像比对的字画艺术品鉴定方法,其特征在于,对所述高光谱图像进行高光谱图像去噪的预处理,具体为,首先利用自适应中值滤波法去除所述高光谱图像中的椒盐噪声,再利用分组主成分变换法对所述高光谱图像进行光谱维噪声去除,最后利用bm4d对所述高光谱图像进行图像维噪声去除。4.根据权利要求3所述的一种基于高光谱图像比对的字画艺术品鉴定方法,其特征在于,所述自适应中值滤波法具体包括,定义滤波窗口最大尺寸为s
max
,计算所述滤波窗口内的中值r
med
、最大值r
max
和最小值r
min
,并进入第一进程,第一进程:a1=r
med-r
min
,a2=r
med-r
max
,如果a1>0且a2<0,则转至第二进程,否则以设定的增加阈值增大滤波窗口尺寸,如果窗口尺寸小于或等于s
max
,则重复第一进程,否则输出中值r
med
替代r
l
(i,j);第一进程:b1=r
l
(i,j)-r
min
,b2=r
l
(i,j)-r
max
,如果b1>0且b2<0,则保持当前像元r
l
(i,j),否则输出输出中值r
med
替代r
l
(i,j)。5.根据权利要求1所述的一种基于高光谱图像比对的字画艺术品鉴定方法,其特征在
于,所述自动端元提取和光谱差异比对法包括,图像端元提取步骤,用于确定所述待检验图像a和参考图像集b的高光谱子图像参考集的端元数p;并将所述待检验图像a和高光谱子图像参考集的每一个波段按逐行排列的方式转换为一个n维的行向量,记为影像x;提取所述影像x中的p个端元作为影像x的端元光谱集e;图像间的端元光谱差异比对步骤,用于根据所述待检验图像a的端元光谱集e
a
和高光谱子图像参考集的端元光谱集e
b
内的各端元光谱之间距离值构建矩阵d
a,b
,并根据所述矩阵d
a,b
中的元素获取所述端元距离的欧几里得距离d
e
和余弦角d
s
;同时计算端元光谱集e
a
和端元光谱集e
b
的距离s(e
a
,e
b
);其中,p为1以上的自然数。6.根据权利要求5所述的一种基于高光谱图像比对的字画艺术品鉴定方法,其特征在于,所述提取影像x中的p个端元作为影像x的端元光谱集e,具体包括,s1、所述影像x经pca变换得到p-1维特征影像s2、在所述p-1维特征影像的第一成分中随机选择像元r,并遍历所述p-1维特征影像中所有的像元x
i
(i=1,2,...,n),计算行列式绝对值将所述行列式值最大的像元作为第一个端元并记录所述第一个端元的像元坐标为t1;s3、假设(t1,...,t
k
)(k≥1)是已经提取的端元的像元坐标,从所述p-1维特征影像的前k个主成分特征影像中取出对应的端元主成分特征对任一像元定义单形体体积,并通过如下单形体体积公式进行计算,采用所述单形体体积公式对每个像元进行计算,以获取相应的单形体体积,并得到单形体体积最大的像元t
k+1
,同时记录端元的像元坐标为(t1,...,t
k
,t
k+1
);s4、如果k+1<p,则k=k+1,并转到s3;否则(t1,...,t
k
,t
k+1
)记录的即为最后端元的图像位置;从所述影像x对应位置提取的像元光谱即为端元光谱集e={e1,e2,...,e
p
}。7.根据权利要求5所述的一种基于高光谱图像比对的字画艺术品鉴定方法,其特征在于,在所述图像间的端元光谱差异比对步骤中,所述矩阵d
a,b
为,其中,d
a,b
的任一元素为d
i,j
,i=1,...,p
a
;j=1,...,p
b
;d
i,j
表示端元光谱集e
a
内端元和e
b
内端元之间的端元距离;通过以下两个公式分别获取所述欧几里得距离d
e
和余弦角d
s
,d
e
(e1,e2)=||e
1-e2||,
上式中,e1和e2分别表示端元和端元的端元向量,e1·
e2表示e1和e2这两个端元向量之间的点积,|| ||表示2范数;通过以下公式计算端元光谱集e
a
和端元光谱集e
b
的距离,式中,式中,8.根据权利要求5所述的一种基于高光谱图像比对的字画艺术品鉴定方法,其特征在于,基于词袋模型对待检验图像a和所述参考图像集b中近红外-红绿蓝多光谱参考子图像集进行图像内容相似度对比,具体包括,特征提取和描述步骤,用于对所述待检验图像a和近红外-红绿蓝多光谱参考子图像集进行近红外、红色、绿色、蓝色通道的重采样,以获得近重采样多光谱图像;对所述重采样多光谱图像进行超像素分割,以超像素的特征值表达所述重采样多光谱图像的颜色特征;对所述待检验图像a和近红外-红绿蓝多光谱参考子图像集进行尺度空间极值点检测、特征点定位、确定主方向和生成特征点描述处理,以获得surf特征;视觉词典生成步骤,用于对所述重采样多光谱图像的颜色特征和surf特征进行聚类,将聚类的中心定义为视觉单词,并将所有的视觉单词进行组合以构成视觉词典;图像量化步骤,用于将所述近红外-红绿蓝多光谱参考子图像集和待检验图像a中的每一个颜色特征、surf特征与所述视觉词典中的每个特征向量进行比较,并将所述颜色特征、surf特征近似成与其最相似的词典项,然后计算所有所述颜色特征、surf特征在每个词典项出现的次数,得到词典直方图;将所述近红外-红绿蓝多光谱参考子图像集和待检验图像a分别量化成颜色视觉单词直方图和surf视觉单词直方图;词袋特征相似性对比步骤,用于对所述颜色视觉单词直方图和surf视觉单词直方图中的词袋特征进行相似性对比,以获取所述surf视觉单词直方图中surf视觉特征相似值和颜色视觉单词直方图中的颜色视觉特征相似值。9.根据权利要求8所述的一种基于高光谱图像比对的字画艺术品鉴定方法,其特征在于,在所述词袋特征相似性对比步骤中,视觉特征相似值通过以下公式获取,式中,h
a
和h
b
为待比较的两个词袋特征向量。10.根据权利要求8所述的一种基于高光谱图像比对的字画艺术品鉴定方法,其特征在于,通过下式获取所述surf视觉特征相似值和颜色视觉特征相似值的加权相似性特征,x
f
=w
shape
*x
shape
+w
color
*x
color
;式中,x
shape
为surf视觉特征相似值;x
color
为颜色视觉特征相似值;w
shape
+w
color
=1;
根据所述加权相似性特征将近红外-红绿蓝多光谱子图像参考集中的子图像进行排序,并将与所述待检验图像a最相似的子图像选取出来作为最相似图像t;获取所述待检验图像a和最相似图像t的端元光谱集e
a
和端元光谱集e
t
的距离s(e
a
,e
t
),将所述距离s(e
a
,e
t
)与设定的端元光谱相似性阈值进行比较,若所述距离s(e
a
,e
t
)小于设定的端元光谱相似性阈值,则判定所述待检验图像代表的字画作品为真。

技术总结
本发明公开了一种基于高光谱图像比对的字画艺术品鉴定方法,涉及字画鉴定技术。采集样本字画作品的高光谱图像;对所述高光谱图像进行预处理;构建具有辨识意义的字画图像感兴趣区域数据集;通过自动端元提取和光谱差异比对法提取待检验图像和所述参考图像集中高光谱子图像参考集的端元,并进行端元光谱相似性比对;基于词袋模型对待检验图像和所述参考图像集中近红外-红绿蓝多光谱子图像参考集进行图像内容相似度对比;将端元光谱相似性比对结果和图像内容相似度对比结果进行融合,以对所述待检验字画图像块进行真伪判定。本发明检验结果客观,准确度高、实用性好且易操作,能够有效的服务于物证鉴定和艺术品鉴定工作。效的服务于物证鉴定和艺术品鉴定工作。效的服务于物证鉴定和艺术品鉴定工作。


技术研发人员:黄威 张旭毅 黄远程 钟燕飞 张洪艳 李志刚 许小京 王桂强 张良培 罗旭东 刘栋卓
受保护的技术使用者:广州星博科仪有限公司
技术研发日:2023.07.21
技术公布日:2023/10/15
版权声明

本文仅代表作者观点,不代表航家之家立场。
本文系作者授权航家号发表,未经原创作者书面授权,任何单位或个人不得引用、复制、转载、摘编、链接或以其他任何方式复制发表。任何单位或个人在获得书面授权使用航空之家内容时,须注明作者及来源 “航空之家”。如非法使用航空之家的部分或全部内容的,航空之家将依法追究其法律责任。(航空之家官方QQ:2926969996)

航空之家 https://www.aerohome.com.cn/

航空商城 https://mall.aerohome.com.cn/

航空资讯 https://news.aerohome.com.cn/

分享:

扫一扫在手机阅读、分享本文

评论

相关推荐