河口湿地外来种监测方法与流程

文档序号:14773148发布日期:2018-06-23 02:09阅读:394来源:国知局

本发明涉及生态环境评估领域,尤其涉及一种河口湿地外来种监测方法。



背景技术:

研究和监测滨海/河口湿地植被的空间特征及其变化,是掌握湿地生态环境特征的一个主要工具,能为湿地生态系统生物的保护、底栖环境质量的评价、受损生境的生态恢复提供的科学依据和管理建议。随着全球化进程的不断深入,生物入侵造成了滨海/河口湿地生态系统中生物多样性的急剧丧失,监测和管理非本土植物分布状况急需重视。然而,仅依靠传统调查手段已经远远不能在景观尺度上研究入侵植物空间分布。幸而遥感技术的众多优点为解决这一问题提供了新的办法。

遥感技术具有覆盖面广、空间和时间尺度多样、光谱信息丰富、观测灵活及数据获取方便等优势,已成为湿地生态环境监测中最具成本效益和高效率的数据采集手段,在湿地规划、管理和保护中扮演着举足轻重的角色。由于在潮间带受到生物入侵、潮汐作用、放牧等人为活动的影响,容易产生面积较小的植物或地物斑块,要精确监测其空间分布及动态变化,关键需要高空间分辨率的遥感影像数据源。高空间分辨率数据(VHR)越来越多地用于研究和监测入侵植物物种的空间分布,它们包括Pleiades,QuickBird、IKONOS、IRS-P5等航天平台,以及几乎所有航拍影像,它们的空间分辨率有些甚至达到l米以下。另一方面,由于潮间带不同植物种光谱特征相似程度高,有必要在某些物候期以及适当的空间和光谱分辨率下获取数据,通过物候学信息增加可分性。已经表明,数据采集的时机起着重要的作用,因为在植被季节的某些时期植物通常与周围植被更加不同。使用变化监测方法,与背景植被相比,某些物种的生命周期差异可以被识别。由于植物物种在某些物候阶段可能被更好地识别,高空间分辨率遥感数据的采集时机十分重要。然而,VHR卫星数据的获取成本很高,一般情况下,VHR卫星数据的获取受到云和卫星的正常轨迹的约束。卫星空中移动的灵活性不高,并需要提前计划,费用十分昂贵,且无法通过不确定的天气预报确保不受云层遮挡。因此,现阶段进行VHR数据物候阶段的选择是有限的。

近年来我国无人机技术日臻成熟,无人机具有便携性、低成本、低损耗、可重复利用、风险小、易于部署等优势,为地物监测提供了高灵活性。另外,无人机能在中等尺度上获取地物本身的空间特征,或与邻近地物之间纹理特征,或者位置信息,这些信息突破了野外作业的时空限制,为地物的精细分类提供更为全面的信息,进而优化分类算法。然而,在植物监测方面,单单依靠无人机获取的图像进行分类,又面临着无人机图像光谱分辨率较低,由于几何失真和辐射不规则导致需要复杂处理的缺点。



技术实现要素:

针对上述现有技术中的不足,本发明提供一种河口湿地外来种监测方法,利用现有成熟的遥感软件和地理信息系统软件对入侵物种植被进行分类,具有高时效、高分辨率、高信息量、实地调查效率高、准确度高、客观性强和成本低的优点。

为了实现上述目的,本发明提供一种河口湿地外来种监测方法,包括步骤:

S1:根据监测需求选择一无人机平台和一卫星遥感片源;

S2:在一目标监测区域设置至少一无人机的航线并对所述无人机设置对应的飞行参数;

S3:根据目标监测区域的气候数据设置一无人机数据采集任务;

S4:执行所述无人机数据采集任务,获得无人机数据;

S5:将所述无人机数据的地理信息与自所述卫星遥感片源获得的卫星数据进行匹配,获得预处理卫星图片;

S6:提取所述预处理卫星图片中的斑块特征信息;

S7:根据所述斑块特征信息筛选获得最优可分性特征指数;

S8:根据所述最优可分性特征指数对所述卫星数据进行分类并对结果进行精度评价。

优选地,所述S1步骤中,所述卫星遥感片源包括Pleiades遥感高分辨率影像源、QuickBird遥感高分辨率影像源和IKONOS遥感高分辨率影像源。

优选地,所述S2步骤中,所述飞行参数包括航线、拍摄高度、拍摄角度、飞行速度和拍摄间隔;根据所述卫星遥感片源的所述卫星数据设置所述航线、拍摄高度和拍摄角度;根据所述航线设定所述飞行速度;根据所述飞行速度设定所述拍摄间隔,所述无人飞机按照所述拍摄间隔拍摄的相邻照片具有影像重叠区域;获得所述目标监测区域的潮间带横截面距离d。

优选地,所述S3步骤进一步包括步骤:

S31:获取所述目标监测区域的潮汐表;

S32:根据所述潮汐表获取所述目标监测区域大潮期的高潮潮位、低潮潮位和潮时;

S33:根据所述无人机的最远飞行距离D和所述潮间带横截面距离d,确定所述无人机的拍摄执行时间段,形成所述无人机数据采集任务;

当D>d时,所述无人机的拍摄执行时间段位于所述高潮潮位的所述潮时之前,所述无人机的拍摄执行时间段的时长大于等于所述无人机沿所述航线拍摄的往返时间;

当D≤d时,工作人员需携带无人机移动至距离水边的离水距离d’小于所述最远飞行距离D时再进入所述无人机的拍摄执行时间段,所述无人机的拍摄执行时间段的时长大于等于所述无人机沿所述航线拍摄的往返时间与所述工作人员的往返时间之和。

优选地,所述S5步骤进一步包括步骤:

S51:将所述卫星数据进行波段融合和几何校正的预处理,获得预处理遥感数据;

S52:通过地理信息系统软件导入所述无人机数据;

S53:将所述无人机数据矢量化,矢量化后的所述无人机数据包括照片数据、第一GPS信息和植物种类信息,所述照片数据的各第一矢量点与所述第一GPS信息和所述植物种类信息关联;

S54:所述预处理遥感数据包括卫星影像数据、第二GPS信息、未知植物种类光谱信息和空间信息;利用所述第一GPS信息和所述第二GPS信息的一致性,将所述卫星影像数据的各第二矢量点与所述照片数据的所述第一矢量点所述关联的各信息相匹配,使得所述未知植物种类光谱信息和所述植物种类信息相匹配。

优选地,所述S6步骤中,利用eCognitionDeveloper软件的分割功能获得同一种植物在所述预处理卫星图片中所占的范围,按照植物种类在所述预处理卫星图片中形成多个斑块;利用eCognitionDeveloper软件的分类功能为所述斑块赋值关联对应的所述植物种类信息;形成斑块特征信息,所述斑块特征信息包括:斑块光谱信息、斑块形状信息和斑块位置信息。

优选地,所述S7步骤进一步包括步骤:

利用所述eCognitionDeveloper软件中的输出功能,输出所述斑块特征信息的统计数据;

通过多元方差分析所述统计数据,筛选出所述最优可分性特征指数并确定每一类植被的阈值范围。

优选地,所述S8步骤进一步包括步骤:

基于所述最优可分性特征指数通过所述eCognitionDeveloper软件建立分类算法规则集;

利用所述eCognitionDeveloper软件中的算法集的监督分类法对所述卫星数据进行分类,获得一分类像元集合;

将所述分类像元集合与一参考像元集合作比较,评价所述卫星数据的精度。

本发明由于采用了以上技术方案,使其具有以下有益效果:

首先,本发明利用无人机数据的时空便捷性和可准确定位植被物种信息的功能,结合高空间分辨率遥感数据丰富的光谱和空间纹理信息,以及可用成熟分类软件的易处理性,提供了一种适合进行植被物种的精确监测方法,在监测选定的入侵物种具有优势。这一方法克服了在湿地区域仅仅通过高分辨率遥感数据监测,其受天气、潮汐、成本所限制,可用数据少,难以通过足够的物候或光谱信息精确监测入侵植被的难点问题。实现了在空间/光谱/时间分辨率、成本、精度之间的优化组合和平衡。

其次,无人机采样的实施还克服在湿地实地样点采样的困难,无人机易于操作,覆盖范围大,减少人工野外工作的时间和工作量,提高实地调查的效率,不仅可获取人为采样所不能达到的大量数量样本量,并且还能获取其不能获得的中尺度上的斑块信息,提高分类精度。

附图说明

图1为本发明实施例的河口湿地外来种监测方法的流程图。

具体实施方式

下面根据附图1,给出本发明的较佳实施例,并予以详细描述,使能更好地理解本发明的功能、特点。

请参阅图1,本发明实施例的一种河口湿地外来种监测方法,包括步骤:

S1:根据监测需求选择一无人机平台和一卫星遥感片源。

本实施例中,选择了大疆公司(DJI Corporation)的Phantom 3无人机,这是一种小型的、电池供电的、搭载高清摄像头、相对便宜的无人机,它能自动设置航线规划,拍摄角度、高度、拍摄间隔和飞行速度;获取的照片集成了内置GPS地理信息,与卫星数据可进行无差匹配。

卫星遥感片源包括Pleiades遥感高分辨率影像源、QuickBird遥感高分辨率影像源和IKONOS遥感高分辨率影像源等;也可采用现有的其他遥感高分辨率影像源对其具体种类不做限制。

另外在选择片源时,尽量在可选片源下选择适当的物候期片源,尽量选择监测植被物种与周围植被差异性较大的物候期时间,如出芽期、开花期或枯萎期。

S2:在一目标监测区域设置至少一无人机的航线并对无人机设置对应的飞行参数。

飞行参数包括航线、拍摄高度、拍摄角度、飞行速度和拍摄间隔;根据卫星遥感片源的卫星数据设置航线、拍摄高度和拍摄角度;根据航线设定飞行速度;根据飞行速度设定拍摄间隔,无人飞机按照拍摄间隔拍摄的相邻照片具有影像重叠区域。在设置航线时可获得目标监测区域的潮间带横截面距离d。

本实施例中,可根据卫星遥感片源的卫星数据中的遥感图像景观的异质性特征设置航线;根据物种的可识别程度,图片中斑块大小设定拍摄高度;根据植被冠层特征设定拍摄角度;根据每次航线飞行的距离设定飞行速度;根据飞行速度设定拍摄间隔,使照片重叠和对角线分别在图像高度和宽度的80%到85%之间,保证照片与照片之间有足够的重叠,实现采样区域的均匀覆盖。

S3:根据目标监测区域的气候数据设置一无人机数据采集任务,其进一步包括步骤:

S31:获取目标监测区域的潮汐表;

S32:根据潮汐表获取目标监测区域大潮期的高潮潮位、低潮潮位和潮时,并根据高潮潮位、低潮潮位和潮时获得目标监测区域的潮间带横截面距离d;

S33:根据无人机的最远飞行距离D和潮间带横截面距离d,确定无人机的拍摄执行时间段,形成无人机数据采集任务;

当D>d时,无人机的拍摄执行时间段位于高潮潮位的潮时之前,无人机的拍摄执行时间段的时长大于等于无人机沿航线拍摄的往返时间;

当D≤d时,工作人员需携带无人机移动至距离水边的离水距离d’小于最远飞行距离D时再进入无人机的拍摄执行时间段,无人机的拍摄执行时间段的时长大于等于无人机沿航线拍摄的往返时间与工作人员的往返时间之和。

S4:执行无人机数据采集任务,获得无人机数据。

S5:将无人机数据的地理信息与自卫星遥感片源获得的卫星数据进行匹配,获得预处理卫星图片;其进一步包括步骤:

S51:将卫星数据进行波段融合和几何校正的预处理,获得预处理遥感数据;

S52:通过地理信息系统软件导入无人机数据;本实施例中,地理信息系统软件采用ArcGIS软件;

S53:通过ArcGIS软件中的Data Management的Geotagged photos to point功能模块将无人机数据矢量化,矢量化后的无人机数据包括照片数据、第一GPS信息和植物种类信息,照片数据的各第一矢量点与第一GPS信息和植物种类信息关联。

S54:预处理遥感数据包括卫星影像数据、第二GPS信息、未知植物种类光谱信息和空间信息;利用第一GPS信息和第二GPS信息的一致性,将卫星影像数据的各第二矢量点与照片数据的第一矢量点关联的各信息相匹配,使得未知植物种类光谱信息和植物种类信息相匹配。

通过这一步骤可以获得卫星影像数据上某点的植被名称、光谱信息和空间信息。再通过步骤S7,利用遥感分类软件,将卫星影像数据中已知植被物种样点、样点所在斑块的光谱信息和空间信息特征提取出来,进行分析。

本步骤中包含对高空间分辨率遥感数据,即卫星数据的预处理,原始数据需要进行波段融合和几何校正,才能变成高清晰和位置准确的可用数据。这些步骤都可以通过现有遥感软件的波段融合和几何校正功能模块实现。

由于高空间分辨率遥感数据包括多光谱波段(如P星为2m分辨率)和全色波段(如P星为0.5m分辨率)。多光谱波段产品提供比全色波段产品更为丰富的地物光谱信息,能形成彩色合成像片,有利用地物的识别,但空间分辨率较全色波段低。全色波段单波段,在图上显示是灰度图片,一般空间分辨率高,但无法显示地物色彩。将全色波段与多光谱波段影像融合处理,得到既有全色影像的高分辨率,又有多波段影像的彩色信息的影像。

几何校正一般是指通过一系列的数学模型来改正和消除遥感影像成像时,由于飞行器的姿态、高度、速度以及地球自转等因素的影响,造成图像相对于地面目标发生几何畸变,这种畸变表现为像元相对于地面目标的实际位置发生挤压、扭曲、拉伸和偏移等,针对几何畸变进行的误差校正就叫几何校正。几何精校正:利用控制点进行的几何校正,它是用一种数学模型来近似描述遥感图像的几何畸变过程,并利用畸变的遥感图像与标准地图之间的一些对应点(即控制点数据对)求得这个几何畸变模型,然后利用此模型进行几何畸变的校正。

无人机内置了GPS模块,只要能同时收到3颗以上卫星信号后,就能越快越准确地获得GPS定位数据,并记录在照片的EXIF信息中。

S6:按照植物种类在预处理卫星图片中形成多个斑块,并形成斑块的斑块特征信息。

利用eCognitionDeveloper软件的分割功能获得同一种植物在预处理卫星图片中所占的范围,按照植物种类在预处理卫星图片中形成多个斑块;利用eCognitionDeveloper软件的分类功能为斑块赋值关联对应的植物种类信息;形成斑块特征信息,斑块特征信息包括:斑块光谱信息、斑块形状信息和斑块位置信息。

S7:根据斑块特征信息筛选获得最优可分性特征指数。

利用eCognitionDeveloper软件中的输出(Export)功能,输出斑块特征信息的统计数据;通过多元方差分析统计数据,筛选出最优可分性特征指数并确定每一类植被的阈值范围。

本实施例中,为了判别不同植被类别斑块的不同特征信息是否有显著性差异,通过多元方差分析统计数据,可以通过各种统计软件的该功能模块实现运算,本实施例中,是通过SPSS软件进行多元方差分析(分析>一般线性模型>多变量)来实现,对不同类型植物斑块属性的样本均数差异性进行显著性检验,筛选出最优可分性特征指数。

在统计分析中,会得到不同植被种类的不同特征参数下的均值、标准误差以及95%置信区间范围,这些范围构成了具有显著性差异的特征参数下的每一类植被的阈值范围。

S8:根据最优可分性特征指数对卫星数据进行分类并对结果进行精度评价。

首先,基于最优可分性特征指数通过eCognitionDeveloper软件建立分类算法规则集;本实施例中,利用eCognition Developer软件建立分类算法规则集,仅需要先建立所有的地物类别名称,然后在每一类地物下使用的监督分类算法;例如最小邻近分类算法;再将确定的优化特征参数选定加入每一类地物类别算法下即可。

其次,利用eCognitionDeveloper软件中的算法集的监督分类法对卫星数据进行分类,获得一分类像元集合;监督分类是以建立统计识别函数为理论基础、依据典型样本训练方法进行分类的技术,即根据已知训练区提供的样本,通过选择特征参数,求出特征参数作为决策规则,建立判别函数以对各待分类影像进行的图像分类,是模式识别的一种方法。在监督分类中,特征参数的选择非常重要,因此,本步骤的目的就是分析确定最优的可区分不同植被类型的特征参数;

最后,将分类像元集合与一参考像元集合作比较,评价卫星数据的精度。

现以本发明所采用的方法在长江河口湿地监测入侵植物互花米草应用为例:

研究地位于长江口最大的岛屿—崇明岛最东端的东滩,该区域地处崇明东滩国家级自然保护区,是一个典型的滨海/河口湿地。为了土地的需要,互花米草1979年开始引入我国,是一种适应海滩潮间带生长的耐盐、耐淹植物,主要用于保滩护堤、促淤造陆和改良土壤等。互花米草隶属禾本科米草属,为多年生高大草本植物,植株健壮而挺拔,平均株高约1.5m。由于互花米草具有一系列有利于其种群生存和扩散的机制,抑制其它植物生长,也使贝类在密集的米草草滩中活动困难,甚至会窒息死亡,威胁到鱼类、鸟类的食物来源,降低滩涂的生物多样性,严重破坏崇明东滩生态敏感区的自然平衡,对当地产生了严重的生态与进化后果改变了整个生态系统的食物网与营养结构,致使原有的很多鱼类及候鸟失去营养支撑,面临危险。

崇明东滩植物群落具有优势植物种类少,群落类型简单等特点,有利于遥感监测研究。虽然崇明东滩共有高等植物约96种,但优势种主要为芦苇(Phragmitesaustralis)、互花米草(Spartinaalterniflora)、糙叶苔草(Carexscabrifolia)、海三棱藨草(Scirpusmariqueter)。在地理分布上,植被具有一定的分带现象,即首先由低潮位的光滩逐渐演变为以海三棱藨草为优势种,间有藨草、糙叶苔草的群落,其次在中潮位和高潮位演替为以芦苇或互花米草为主的群落,并形成相互竞争的局面。芦苇和互花米草通常可以生长到1.5~2.0m,而藨草、糙叶苔草、海三棱藨草的生长高度约为0.3~0.7m。

第一步:选择无人机平台和卫星遥感片源。

本实施例中,选择的是大疆公司(DJI Corporation)的Phantom 3无人机,对于卫星遥感片源,选择了主流的Pleiades星遥感高分辨率影像源,空间分辨率为0.5m,根据实地观察,崇明东滩的芦苇、互花米草、海三棱藨草、糙叶苔草等优势物种表现出一定的物候差异,尤其是春季出苗期和秋季枯萎期有前后顺序,根据可选影像的质量(角度、云量),最终选择了2015年10月中的遥感影像作为目标解译影像。该时间段,芦苇10月中旬为其盛花期;而互花米草7月末为盛花期,10月中旬为结种子期;糙叶苔草、海三棱藨草在10月中旬已经开始枯萎。

第二步:在目标监测区域设置无人机的航线并对无人机设置对应的飞行参数。

在无人机飞行任务期间,根据遥感图像景观的异质性特征设置航线,由南向北设置三个飞行站位,每个站位向三个方向设置三条航线,一共9条航线,使所有航线较为均匀地分布整个研究区范围;根据物种的可识别程度(斑块大小)设定拍摄高度,根据野外实地调查经验,大部分中小植物斑块直径范围约为10~30m左右,根据实际飞行测试,飞行高度和斑块直径范围相当,拍摄照片的可识别程度较高15~20m;根据植被冠层高度和植被种类设定拍摄角度,东滩湿地优势植物为草本植物,叶片冠层细窄,株高大多范围在0.5~2m左右,若正射俯视拍摄难以区分植被种类,拍摄角度设为仰角30°~45°,能较好地分辨照片中植物种类;根据每次航线飞行的距离设定飞行速度,由于每条航线距离在1km~2km左右,根据测试,飞行速度设为10m/s,拍摄间隔2s/张,能使照片重叠和对角线分别在图像高度和宽度的80%到85%之间,保证照片与照片之间有足够的重叠,实现采样区域的均匀覆盖。

第三步:根据目标监测区域的气候数据设置一无人机数据采集任务。

查询所调查的湿地区域的潮汐表,获悉该区域的大潮期的高潮潮位、低潮潮位和潮时;根据大潮期的高潮潮位、低潮潮位和潮时以及登陆和撤离区域的高程,确定作业人员登陆和撤离的时间以及无人机低空飞机器航拍的可执行时间段。

第四步:执行无人机数据采集任务,获得无人机数据。

第五步:将无人机数据的地理信息与来自卫星遥感片源获得的卫星数据进行匹配,获得预处理卫星图片。

首先,对遥感影像数据(卫星数据)的多波段数据和全色波段数据进行融合,使其分辨为0.5米。利用国家1:2000的地形图进行多项式系数几何校正。具体地:对购买的原始高分辨率遥感数据产品进行几何校正,找到5~10对地形图和遥感影像数据上同一位置上的控制点(一般是识别度高的道路交叉点、桥梁或标志房屋等),通过这几对控制点的各自坐标(XY值),利用多项式系数算法,算出误差最低时的公式系数,通过该公式,对遥感影像上所有像素点进行转换,完成整图的几何校正。

通过几何校正,能够确保卫星数据和无人机照片中研究区域和验证点匹配在相同区域。接着,通过ArcGIS软件导入无人机照片GPS信息与高空间分辨率遥感数据(VHR数据)匹配,利用ArcGIS软件中的Data Management的Geotaggedphotos to point工具,将照片中基于EXIF的GPS信息转为矢量点格式,并通过identify查询工具可溯源原始照片植被信息,据此为矢量点数据添加植被物种信息,物种属性信息名称与最后在eCognitionDeveloper软件分类过程类别名称要保持一致。

第六步:提取预处理卫星图片中的斑块特征信息。

利用eCognitionDeveloper软件中的分割(Segmentation)、分类(Classification)功能,提取无人机样点所在的入侵物种和其他物种斑块的光谱、形状、位置信息。

分割分类算法中,先对VHR遥感影像(高空间分辨率遥感数据)进行多尺度分割(Multiresolution segmentation),确定符合大部分植物斑块特征的最佳尺度(Scale);在此最佳尺度下,通过棋盘分割法(Chessboard segmentation),得到与样本点相同并带有类别信息的最小分割对象,即为后述步骤中分析所使用的样点斑块。利用分类算法(Assign class by thematic layer)通过样本点的属性(Class_name)对样点图斑进行分类;继续利用分类算法(Assign class)将和样本点对象重叠的样本斑块赋予和样本点同样的类别(通过添加属性限制条件:Class-related features中的border to class属性),此时获得所有已知植被物种的样点所在匀质斑块。

通过在eCognitionDeveloper软件中的斑块特征选择功能,选择斑块光谱、空间特征:在Feature View界面添加创建需要的指数特征指令,例如:NDVI,layer mean/standard deviation values,geometry-extent-area/length/width,position-distance to the scene border,geometry-shape-asymmetry/border index/compactness/main direction/shape index,texture-layer value texture based on sub-object,texture-shape value texture based on sub-object,来获得已知的样点斑块的相关空间特征值。

第七步:根据斑块特征信息筛选获得最优可分性特征指数。

利用eCognitionDeveloper软件中的使用输出对象统计算法(Export object statistics),设置level,class filter,Features,输出各植物斑块各个属性信息(光谱、形状、位置)的统计数据,然后进行统计再分析对比出对不同类别分类的有效特征。

通过SPSS软件进行多元方差分析(分析>一般线性模型>多变量),对不同类型植物斑块属性的样本均数差异性进行显著性检验,筛选出最优的可分性特征指数,并确定这些指数用于区分每一类植被的阈值范围。

第八步:根据最优可分性特征指数对卫星数据的精度进行评价。

先将第四步提取的影像对象,即已经分类的样本斑块转化成分类的子样本(Classified image objects to samples);再基于得到的最优可分性特征指数,在eCognition Developer软件中建立分类算法规则集,根据算法集对整幅VHR遥感图像进行监督分类,本实施例中,选择最近邻近分类法(Nearest neighbor classification)对VHR遥感数据进行分类,并通过误差矩阵法(Error matrix based on TTA mask)进行精度评价。

误差矩阵(error matrix)法是一类精度评价的方法,是适用于利用样本特征进行监督分类后的精度评价方法。是一个用于表示分为某一类别的像元个数与地面检验为该类别数的比较阵列。通常,阵列中的列代表参考数据,行代表由遥感数据分类得到的类别数据。有像元数和百分比表示两种。从误差矩阵总体分类精度、Kappa系数、错分误差、漏分误差、每一类的制图精度和用户精度,通过这些精度值来确定分类的精度和可靠性。

以上结合附图实施例对本发明进行了详细说明,本领域中普通技术人员可根据上述说明对本发明做出种种变化例。因而,实施例中的某些细节不应构成对本发明的限定,本发明将以所附权利要求书界定的范围作为本发明的保护范围。

当前第1页1 2 3 
网友询问留言 已有0条留言
  • 还没有人留言评论。精彩留言会获得点赞!
1