|
空-谱融合的条件随机场高光谱影像分类方法 魏立飞1 1. 湖北大学资源环境学院, 湖北 武汉 430062; 收稿日期：2019-01-25；修回日期：2019-10-24 基金项目：国家重点研发计划课题（2017YFB0504202）；国家自然科学基金优秀青年科学基金（41622107）；湖北省技术创新专项重大项目（2018ABA078）；空间数据挖掘与信息共享教育部重点实验室开放基金（2018LSDMIS05）；农业部农业遥感重点实验室开放基金（20170007） 第一作者简介：魏立飞(1983—), 男, 博士, 副教授。研究方向为遥感图像处理、农业遥感、生态遥感。E-mail:weilifei2508@163.com 通信作者：余铭, E-mail:ym17317399819@163.com 摘要：高光谱遥感数据具有光谱信息丰富、图谱合一的特点，目前已经广泛地应用在对地观测中。传统的高光谱分类模型大多过分依赖影像光谱信息，没有充分利用空间特征信息，这使得分类精度还有很大的提升空间。条件随机场是一种概率模型，能够较好地融合空间上下文信息，在高光谱影像分类中已经得到越来越多的关注，但大部分条件随机场模型存在超平滑的现象，会导致影像细节丢失。针对该问题，本文提出了一种优化融合影像空-谱信息的高分辨率/高光谱影像分类方法，该方法将影像的纹理信息与原始光谱信息进行融合，利用SVM分类器对其进行预分类，并将各类概率定义为一元势函数，以融合空间特征信息；然后将空间平滑项和局部类别标签成本项加入二元势函数中，以考虑空间背景信息，并保留各类别中的详细信息。最后，通过两组的高分辨率/高光谱影像数据进行试验。结果表明，与SVM算法、传统的条件随机场方法和面向对象的分类方法相比，本文提出的算法在整体分类精度上分别提高了10%、9%和8%以上，同时在保持地物边缘完整性、避免'同谱异物'与'同物异谱'的现象方面有较明显的优势。 关键词：高光谱遥感影像 条件随机场 空-谱融合 影像分类 Hyperspectral image classification method based on space-spectral fusion conditional random field WEI Lifei1 1. Faculty of Resources and Environmental Science, Hubei University, Wuhan 430062, China; Foundation support: The National Key Research and Development Program of China (No. 2017YFB0504202);The National Natural Science Foundation of China (No. 41622107);The Special Projects for Technological Innovation in Hubei (No. 2018ABA078);The Open Fund of Key Laboratory of Ministry of Education for Spatial Data Mining and Information Sharing (No. 2018LSDMIS05);The Open Fund of Key Laboratory of Agricultural Remote Sensing of the Ministry of Agriculture (No. 20170007) First author: WEI Lifei (1983—), male, PhD, associate professor, majors in remote sensing image processing, agricultural remote sensing and ecological remote sensing.E-mail:weilifei2508@163.com. Corresponding author: YU Ming, E-mail:ym17317399819@163.com. Abstract: Hyperspectral remote sensing image has the characteristics of rich spectral information and combining image with spectrum, which has been widely applied in the earth observation. Most of traditional hyperspectral image classification models don't make fully use of spatial feature information, rely too much on the spectral imformation, making the classification accuracy still have a lot of room to improve. Conditional random field (CRF) is a kind of probability mode that can better integrate spatial context information. It plays a more and more important role in hyperspectral image classification. However, most CRF models have the problem of excess smoothness, which will result in the loss of detail information. Aiming at this problem, this paper proposed a hyperspectral image classification method based on space-spectral fusion conditional random field. The proposed method designs suitable potential functions in a pairwise conditional random field model, fusing the spectral and spatial features to consider the spatial feature information and retain the details in each class. The experiments on two sets of hyperspectral image showed that, compared with the traditional methods, the proposed classification method can effectively improve the classification accuracy, protect the edges and shapes of the features, and relieve excessive smoothing, while retaining detailed information. Key words: hyperspectral remote sensing imagery conditional random field space-spectral fusion image classification 随着空间技术的发展，高光谱遥感技术已经成为遥感领域最重要的研究热点之一。高光谱遥感技术是指具有高光谱分辨率的遥感数据获取、处理、分析及应用的过程。除了具有较高的空间分辨率，高光谱数据还有连续且丰富的光谱波段，波段数多达数十甚至上百个，能够检测到地物细微差别并鉴别出不同的类别[1-3]，因此被广泛地运用于农业生产、植被调查、地质勘探、国防军事、生态保护、海洋研究等领域[4-7]。 利用高光谱遥感数据分类是高光谱遥感应用热点之一，也面临着许多挑战，如维数灾难、非线性数据结构的挑战、空间同质性与异质性问题等。目前针对这些问题，许多从事高光谱影像分类研究工作的国内外学者一方面利用机器学习、模式识别等领域的算法进行高光谱遥感分类，如文献[8]针对高空间分辨率遥感影像的分类问题，提出了基于深度学习的分类方法，利用影像的光谱-纹理特征有效提高了高分辨率遥感影像的分类精度。文献[9]利用主动学习结合空间约束对高光谱影像分类，提升了分类算法的精度和自适应程度。另一方面则挖掘高光谱遥感影像隐藏的特征信息对高光谱分类进行研究，如文献[10]基于波段分组和分类器集成的方法，根据波段之间的相似性信息对光谱波段进行分类。文献[11]提出基于谱聚类和稀疏表示的两级分类算法对高光谱图像进行分类，减少了噪声及区域边界对分类效果的影响。文献[12]提出了一种基于稀疏性的高光谱图像分类算法，将上下文信息结合到稀疏恢复优化问题中，其分类性能优于经典的监督分类器支持向量机。文献[13]使用多尺度自适应稀疏表示有效地利用了空间信息，使得其提出的高光谱分类方法优于众所周知的几种分类方法。文献[14]提出了基于核函数的高光谱图像分类方法框架，采用光谱信息和空间信息相融合的算法对高光谱影像进行分类。这些方法在一定程度上解决了高光谱遥感影像地物分类的问题，但还存在较大的精度提升空间。 随着传感器技术的完善，同时具备高空间分辨率和高光谱分辨率的遥感数据大量产生，与单独的高光谱分辨率遥感数据相比，这类数据除了具有丰富而连续的光谱波段，还具有较高的空间分辨率，但同时也存在影像上地物出现高度的细节化以及“同物异谱”与“同谱异物”的现象更加明显的局限，降低了光谱可分性[15]，因此单纯利用光谱分类无法满足越来越高的空间分辨率。隐藏在高光谱数据中的空间特征逐渐被挖掘，融合光谱-空间特征的方法[16]越来越多地被应用于高光谱影像分类。条件随机场(conditional random field, CRF)由于能够结合空间背景信息而被广泛地运用于图像分割、立体视觉和活动分析等领域[17-19]，条件随机场分类框架不仅可以考虑标签图像的上下文信息，还可以考虑观察数据的空间相互作用。但是传统的条件随机场方法常常会过度平滑，从而导致细节信息的丢失[20]。因此本文提出一种适合高分辨率高光谱影像的条件随机场分类方法：空-谱融合的条件随机场高光谱影像分类方法。该方法将纹理特征与原始光谱融合在一元势函数中以考虑空间特征信息；然后将空间平滑项和局部类别标签成本项加入二元势函数中，考虑空间背景信息以及标签约束，保留各类别中的详细信息。通过对两组高分辨率高光谱数据进行试验，分别计算各类别精度(CA)、总体精度(OA)与Kappa系数，并与传统的分类方法进行对比，结果表明本文提出的分类方法在分类效果与分类精度上表现更加优异。 1 算法模型 本文针对传统条件随机场存在过度平滑的问题，提出了一种空-谱融合的条件随机场高光谱遥感影像分类方法。该方法将由灰度共生矩阵(gray level co-occurrence matrix, GLCM)获得的纹理特征与原始光谱特征组成空-谱融合特征，将利用SVM分类器获得的各类别概率定义为条件随机场的一元势函数，将空间平滑项与局部类别标签项的线性组合定义为二元势函数，从而达到在结合空间上下文信息的同时保留细节信息的分类效果。 1.1 CRF模型 CRF用统一的概率框架模拟随机变量之间的局部邻域相互作用，它直接模拟标签的后验概率，给出观察图像数据为吉布斯(Gibbs)分布
式中，y为输入图像的观察数据；x是整个图像的相应类别标签；Z是分区函数；ψc(xc, y)是势函数，它根据图像中的邻域系统和派系对随机变量的空间交互进行局部建模。根据观测数据中的派系类别及其相应的标签，势函数可以分为一元势函数、二元势函数和高阶势函数。虽然高阶势函数可以模拟更广泛的上下文信息，但是这类势函数的推理过程一般较困难，因此在分类问题中常常选择二阶条件随机场进行研究，其中包括一阶势函数与二阶势函数[20]。 作为判别分类框架，如式(1)所示，CRF直接模拟标签x的后验分布，给定观察值y，相应的吉布斯能为
相应的，分类影像通过贝叶斯最大后验规则(MAP)找到使后验概率P(x|y)最大化的标签图像y。所以随机场的MAP标记xMAP由式(3)给出
因此，当后验概率P(x|y)最大时，能量函数E(x|y)最小。遥感分类问题可以通过设计合适的势函数来描述
式中，ψi(xi, y)和ψij(xi, xj, y)是分别定义在点i局部领域Ni上的一元势函数和二元势函数。 1.2 G-CRF模型 1.2.1 G-CRF方法技术流程 本文提出一种空-谱融合的条件随机场高光谱影像分类方法。该算法的步骤如下： (1) 通过灰度共生矩阵(GLCM)提取待分类影像的纹理特征，联合每个像素的光谱信息，形成纹理-光谱融合特征。 (2) 将融合特征输入SVM分类器，对标签和融合特征的关系进行建模，并根据给出的标签基于特征向量独立的计算每个像素的概率估计，得到概率影像。 (3) 将待分类影像与概率影像输入G-CRF模型中，在空间平滑项与局部类别成本项的作用下，通过标记场和观测场模拟每个像素与其相应邻域的空间上下文信息，考虑空间相关性并保留相邻区域之间的边界和细节信息，得到更优分类结果。 本方法技术流程如图 1所示。
1.2.2 一元势函数 一元势函数ψi(xi, y)对标签和观察图像数据之间的关系进行建模，通过特征向量来计算采用特定类别标签的单个像素的成本。因此，可以用能够给出标签xi的概率估计的判别式分类器单独计算每个像元，给出特征向量。它在分类过程中起主导作用，一般为某监督分类器的后验概率。通常被定义为
式中，f是一个特征映射函数，它将图像中任意一个像元映射到一个特征向量；fi(y)表示位置i处的特征向量，本文在像素原有光谱信息的基础上加入了空间纹理信息，考虑了影像灰度空间中像元点对之间的相关性和空间结构特性，因此，该特征向量在本文中为空间-光谱特征向量。P[xi=lk|fi(y)]是基于特征向量的像素i取的标签lk的概率。由于SVM分类器在遥感图像分类的小训练集中表现良好，本文选取SVM分类器从空间-光谱特征中获取各类别概率作为一元势能。 1.2.2.1 光谱特征 最小噪声分离法(minimum noise fraction, MNF)是一种常用的光谱特征提取方法，其简便，易于实施。MNF变换后各分量按照信噪比大小排列，信息主要集中在第1分量，随着分量增加，影像质量逐渐降低。有研究表明，与原始高维影像数据和PCA变换得到的特征影像相比，MNF变换得到的低维特征影像能够更加有效地提取影像中的光谱信息[21]，因此本文选择该方法提取高分辨率高光谱影像的光谱信息。 1.2.2.2 纹理特征 为了进一步提高分类精度，本文选取稳定性较好、贡献度较大的6种纹理信息进行辅助分类。纹理特征是一种不依赖于亮度和颜色的反映影像中相邻像元相似信息的视觉特征，它是所有物体表面共有的内在特性，包含了物体表面结构组织规律的重要信息以及它们与周围环境的关系[22-23]。灰度共生矩阵是常用的具有较好鉴别能力的提取纹理信息的方法[24-25]，其原理是通过建立图像中存在一定位置关系的两个像元之间的灰度共生矩阵，并从这个矩阵中提取相应的特征量进行纹理分析。 设f(x, y)是一幅大小为M×N，灰度级别为Ng的二维数字图像，那么满足一定空间关系的灰度共生矩阵为
式中，#(x)为集合x中的元素个数；P为Ng×Ng的矩阵，若(x1, y1)与(x2, y2)之间的距离为d，二者与横坐标轴的夹角为θ，那么各种间距及角度的灰度共生矩阵P(i, j, d, θ)为
本文使用到以下几种纹理统计度量指标。 (1) 同质性(homogeneity)：反映影像灰度的均匀性
(2) 角二阶距(angular second moment)：反映影像灰度分布的均匀程度和纹理的粗细度
(3) 对比度(contrast)：反映影像的灰度变化量
(4) 不相似性(dissimilarity)：衡量影像中灰度值不相似的程度
(5) 均值(mean)：表示纹理的规则程度
(6) 熵(entropy)：反映影像纹理的复杂程度或者非均匀度
1.2.3 二元势函数 二元势函数通过考虑标记场和观测场来模拟每个像素与其邻域之间的空间上下文信息。虽然由于光谱变化和噪声的影响，均匀影像上的相邻像素的光谱值可能看起来不相同，但是由于空间相关性的影响，它们很有可能是相同的类别。二元势函数模拟了这种光滑性并且考虑了标签约束，这有利于分布均匀的区域中具有相同地物的像元的分类以及保留两个相邻的区域的边缘。在本文中二元势函数定义为空间平滑项和局部类别标签成本项的线性组合。基于光谱和空间的相关性，空间平滑项和局部类别标签成本项倾向于使相邻像素具有相同的标签，同时局部类别标签成本项在分类过程中考虑了相邻像素的不同类别标签，有利于考虑上下文信息的同时保留大量地物特征的细节信息。 本文二元势函数定义如下
式中，gij(y)表示与数据y有关的平滑项；ΘL(xi, xj|y)是大小为|L|×|L|的局部类标签成本项，它表示邻域内的标签xi和xj之间的成本。参数θ是控制二元势函数中标签成本项的程度的相互作用系数。函数gij(y)模拟相邻像素i和j之间的相互作用，用于衡量相邻像素之间的差异，定义如下
式中，(i, j)是相邻像素的空间位置；函数dist(i, j)是它们的欧几里得距离；yi和yj是表示像素i和j出现的光谱向量，可以使邻域内相互作用的强度与图像数据相关并且在相似区域促进一致性。参数β为图像中所有相邻像素的光谱向量之间的均方差(即，β=(2〈‖yi－yj‖2〉)－1，其中〈‖yi－yj‖2〉是图像上的平均值)。 局部类标签成本项ΘL(xi, xj|y)通过观察到的图像数据来模拟不同邻域类别标签xi和xj之间的空间关系，并且被定义为
式中，P[xi|fi(y)]是特征向量fi(y)的标签概率。局部类别标签成本项考虑了相邻像素的当前类别标签xi，通过xi来衡量图像中相邻元素i和j处标签之间的相互关系，当特征空间中的类存在较强的相似性或者重叠时，它会通过邻域空间标签信息来改变像元的标签。因此，与当前专题标签相关的局部类别标签成本项通过专题类标签的概率分布的估计形式来考虑光谱信息，以在考虑空间上下文信息的同时进行适当的平滑。 2 试验与分析2.1 试验数据 本文使用两组试验以验证本文方法。第1组试验数据为2017年11月利用M600PRO无人机在湖北省洪湖市获取的高分辨率高光谱数据集(WHU-Hi-Farmland-HH UAV dataset)，该数据可用于精准农业制图，由武汉大学智能数据提取与遥感分析组(RSIDEA组)提供。其原始影像如图 2(a)所示，空间分辨率为0.04 m，大小为400×400像素，共有274个波段，包含屋顶、道路、棉花、棉花柴、油菜、白菜、包菜、菜苔、白萝卜、小青菜、青菜、蚕豆、柿子树等13种地物类型，对应的地面真实分布图为图 2(b)。
第2组试验数据是由AVIRIS传感器采集的美国加利福尼亚州萨利纳斯山谷数据，其原始影像如图 3(a)所示，该数据共224个波段，其中去除了20个吸水带，具有3.7 m的空间分辨率，大小为512×217像素，共包含蔬菜、裸地、葡萄园等16类地物类型，对应的地面真实分布图为图 3(b)。
本文使用中国洪湖市与美国加利福尼亚州萨利纳斯山谷高光谱数据集来验证提出的G-CRF方法。同时对比算法包括传统的像素分类法、面向对象分类方法和随机场分类法。像素分类法采用传统的SVM算法。面向对象方法采用均值漂移分割法(mean shift, MS)[26]。随机场方法则采用DPSCRF方法(简称CRF)[27]。 在试验中，每种算法分别随机选择3%、5%、10%的训练样本进行分类，利用剩余97%、95%、90%的样本进行精度验证。使用各类别精度、总体精度(OA)、Kappa系数[28] 3种指标定量评价分类结果。 2.3 结果与分析 本文中各组试验分别随机选取参考样本的3%、5%、10%作为训练数据集，剩余97%、95%、90%样本作为测试样本。 2.3.1 洪湖数据试验结果及分析 在第1组洪湖市的试验中，通过MNF变换将原始影像从270个波段降至17波段，根据影像特征经试验提取同质性、不相似性、均值、熵4类纹理特征。参数θ和β在3%、5%、10%的情况下分别设置为3/1.6、2/0.8、1.5/1。试验结果如下。 从图 4的分类效果对比图可以看出，在训练样本分别为3%、5%、10%时，SVM算法分类效果均较差，整体具有较多的“椒盐”噪声，各部分无法避免分类混淆现象。如图 4(a)、(e)、(i)右上角产生了大片噪声，该区域的部分棉花柴被错分为棉花；中部的柿子树、青菜、油菜、小青菜分类混淆现象严重。图 4(b)、(f)、(j)为CRF的分类图，可以明显看到噪声减少，但总体看来错分现象仍然严重，如中间的棉花柴部分错分为棉花，小青菜错分为菜苔，包菜和蚕豆分类后边界不清晰，同时由于CRF方法过平滑的作用，部分分布零散的作物被当作噪声处理，因此分类效果较SVM算法有一定优化，但仍然有较大的提升空间。MS算法较前两种方法噪声也减少了，但是对于个别作物无法识别，如图 4(c)(g)(k)所示，小青菜完全被错分为菜苔和包菜。图 4(d)、(h)、(l)为G-CRF算法的分类效果图，可以看到“椒盐”噪声现象被优化，各类地物错分现象减少，同时较好地保持了各类地物的界限，尤其是在前两种方法中分类效果较差的青菜、菜苔、油菜等作物有较大的改善。
从表 1对精度的定量分析可以看到，当训练样本分别为3%、5%、10%时，传统的SVM算法的整体精度分别为88.80%、89.88%、90.92%，Kappa系数分别为0.748 3、0.775 0、0.800 2。青菜的精度无论是在训练样本为3%、5%、10%时均较低，分别为1.4%、1.43%、6.02%，分类效果较差。CRF算法的整体精度分别为89.05%、90.98%、92.01%，Kappa系数分别为0.753 3、0.798 6、0.823 5，较SVM算法精度约提高了1%，但由于其过度平滑作用，部分作物分类精度有所下降，尤其是在训练样本为3%时，包菜和青菜的分类精度仅为0.8%与0.84%。面向对象的MS算法的整体精度分别为90.25%、91.74%、93.01%，Kappa系数分别为0.777 9、0.815 5、0.844 7。该方法对于面积较小的包菜和青菜分类效果较差，在训练样本为3%时，精度均为0。G-CRF算法在同时考虑了空间特征信息与光谱信息之后整体分类精度分别达到了98.87%、99.17%、99.34%，Kappa系数分别为0.975 9、0.982 2、0.986 0，与前3种算法相比，精度分别提高约了10%、9%和8%，各类地物的分类精度均有较大的提高，尤其是在前几种算法中分类效果差的包菜、青菜、小青菜，在训练样本为3%的情况下，分类精度分别达到了97.39%、89.64%和92.19%，因此可以看出该方法的分类效果较好。 表 1 洪湖数据不同算法精度对比Tab. 1 The results comparison on Honghu
表选项 2.3.2 Salinas数据试验结果及分析 在第2组Salinas的试验中，通过MNF变换将原始影像从224个波段降至29波段，根据影像特征经试验提取角二阶距、对比度、均值、熵4类纹理特征，参数θ和β在3%、5%、10%的情况下分别设置为2.5/0.4、0.8/0.4、2.8/1.2。试验结果如下： 图 5为Salinas数据在训练样本分别为3%、5%、10%时的分类效果对比图，从图 5(a)、(e)、(i)可以看出SVM算法分类图中整体噪声非常严重，尤其是C8和C15两类地物分类后产生“椒盐”噪声以及分类混淆。图 5(b)、(f)、(j)为CRF的分类图，整体分类效果较SVM有一定改进，但是C8类和C16类仍有错分情况。图 5(c)、(g)、(k)为MS的分类图，C8类和C16类分类效果依然较差，并且各类之间的边缘模糊不清，如C4和C5、C11和C12。从图 5(d)、(h)、(l)G-CRF算法的分类效果图可以看到分类后的噪声在考虑空间信息与光谱信息后得到了优化，各类地物错分现象得到了改善，并且较好地保持了各类地物的界限，尤其是在其他方法中表现较差的C8类和C15类，在该方法中错分现象减少，分类效果较前3种方法有较高的提升。
从表 2的精度评价对比可以看出本文方法具有一定的优越性。在训练样本为3%、5%、10%时，SVM算法和CRF算法的整体精度分别达到87.18%、87.68%、88.86%和88.55%、89.74%、91.00%，Kappa系数分别为0.856 5、0.862 4、0.875 5和0.871 7、0.885 3、0.899 4，但个别地物的分类精度仍然较低，如在训练样本为3%时，地物C15：Vinyard_untrained在SVM方法下精度仅为44.15%。并且在CRF方法下某些地物分类精度低于SVM方法，如C15：Vinyard_untrained在训练样本为3%和5%时，CRF方法的分类精度均低于SVM。这是由于传统的CRF方法将零散分布的地物当作噪声过度平滑从而导致精度降低。MS算法的整体精度分别为86.57%、86.12%、87.58%，Kappa系数分别为0.849 4、0.844 5、0.860 9，对于C15：Vinyard_untrained分类效果依然较差，分别为25.57、27.70、32.92。G-CRF方法在考虑了空间特征信息与空间背景信息后，整体精度分别达到了98.47%、98.57%、99.04%，Kappa系数分别为0.983 0、0.984 1、0.989 3，较前3种方法分别提高了约11%、10%和12%，各类地物的分类精度也有较大提升，尤其是地物C15：Vinyard_untrained的分类精度在融入了空间信息之后分别达到了99.60%、99.00%、99.27%，比SVM算法提高了约45%，分类效果较好。 表 2 Salinas数据不同算法精度对比Tab. 2 The results comparison on Salinas
表选项
由于高分辨率高光谱影像空间分辨率的提高、波段数增加，遥感影像包含的信息更加丰富，地物呈现出高度的细节化，这使得相同地物的光谱信息可能会产生差别，不同地物的光谱信息可能会产生重叠，导致“同谱异物”与“同物异谱”的现象。同时遥感影像不再仅仅反映地物的光谱信息，也包含了大量不同地物的组织结构信息。 因此，本文提出了一种应用于高分辨率高光谱影像分类的G-CRF方法，通过对条件随机场模型中的一元势函数使用SVM分类器对局部邻域和有限的训练样本获得概率估计，在二元势函数中使用空间平滑项和局部类别标签成本项来平滑局部邻域空间以及考虑空间上下文信息。同时利用GLCM提取高分辨率高光谱影像的纹理信息，并融合到原始影像中进行分类。该方法在无人机获取的数据以及公共数据集上的试验表明，融合了纹理信息的条件随机场方法能够有效地减少噪声，保持地物边缘的完整性，减少错分的概率，同时缓解了传统的条件随机场过度平滑的现象，大大地提高了分类精度。 【引文格式】魏立飞, 余铭, 钟燕飞, 等. 空-谱融合的条件随机场高光谱影像分类方法. 测绘学报，2020，49(3)：343-354. DOI: 10.11947/j.AGCS.2020.20190042 |
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
|
, 余铭1
(1)
(2)
(3)
(4)
