动态北极湿地测绘:一种基于多任务卫星时间序列的方法
《International Journal of Applied Earth Observation and Geoinformation》:Dynamic Arctic wetland mapping: A multi-mission satellite time series approach
【字体:
大
中
小
】
时间:2026年05月01日
来源:International Journal of Applied Earth Observation and Geoinformation 8.6
编辑推荐:
马苏德·马迪安帕里(Masoud Mahdianpari)| 让·格兰杰(Jean Granger)| 法里巴·穆罕默迪马内什(Fariba Mohammadimanesh)| 奥利弗·松嫩塔格(Oliver Sonnentag)| 穆罕默德·马尔贾尼(Mohammad Marj
马苏德·马迪安帕里(Masoud Mahdianpari)| 让·格兰杰(Jean Granger)| 法里巴·穆罕默迪马内什(Fariba Mohammadimanesh)| 奥利弗·松嫩塔格(Oliver Sonnentag)| 穆罕默德·马尔贾尼(Mohammad Marjani)
加拿大纽芬兰与拉布拉多省圣约翰斯市的C-CORE,A1B 3X5
**摘要**
北极-北方湿地在全球碳循环、生物多样性和水文调节中发挥着核心作用,然而其范围和动态变化仍缺乏准确描述。静态湿地地图往往无法捕捉由融雪、降雨、河流洪水和永久冻土退化引起的短期淹没变化,导致对湿地面积及相关温室气体排放的低估。本研究在两种情景下评估了多传感器框架:(I)基于2021至2025年多传感器中值合成图像的典型静态情景;(II)结合Sentinel-1时间序列反向散射数据和Sentinel-2植被指数的动态情景,以反映季节性和年际性淹没行为。在两个500×500平方公里的研究区域内,动态情景下的效果显著提升:加拿大西北地区的总体准确率达到95.2%,瑞典北部的准确率为96.1%,而静态产品的准确率分别为91%和92%。季节性淹没和植被覆盖区域的改进尤为明显,用户和生产者的准确率相对于静态情景提高了3-7%。动态情景还揭示了14-18%的湿地区域存在季节性或偶发性淹没现象,这是静态方法无法发现的。除了更高的准确率外,动态框架还提供了具有生态意义的输出,如最大淹没范围、淹没频率和季节性水文期地图,对甲烷建模、生物多样性评估和水资源管理具有直接价值。
**1. 引言**
湿地是全球重要的生态系统,因其生态生产力和提供的基本服务(包括水文调节、生物多样性支持及碳储存)而受到重视(Isaac等人,2025年)。在北极-北方地区,湿地的重要性因与永久冻土的紧密联系及其对全球碳循环的贡献而更加突出。快速的气候变暖加速了永久冻土的融化和退化,改变了水文和热力状况,并增加了富含碳的土壤中的温室气体(尤其是甲烷)排放(Zheng等人,2025年)。尽管湿地在地球系统过程中起着关键作用,但在区域和全球数据集中仍被忽视或分类错误,这导致陆地生物圈模型存在持续不确定性(Merchant等人,2024年)。因此,准确且一致地绘制北极-北方湿地地图是一项重要的科学任务(北极动植物保护组织,2021年),但这仍是一项具有挑战性的工作(Mahdianpari等人,2021年)。北极-北方地区的湿地景观本质上是动态的,受复杂且通常短暂的水文周期(几周到几个月)影响,地表水位在季节性和年际尺度上波动(Kreplin等人,2021年;Stuefer等人,2017年)。地貌特征(如冰楔多边形和热喀斯特洼地)(Liljedahl等人,2016年;Merchant等人,2024年)使得基于遥感数据的分类变得复杂(Lara等人,2020年;Olefeldt等人,2016年)。冻融循环进一步影响地表特性(如反照率),从而改变全年光谱和反向散射响应。除了这些自然因素外,气候驱动的过程(如融冻引起的沉降、海岸侵蚀和人类活动的影响)也导致湿地范围和功能的快速变化(Chen等人,2025年)。综上所述,迫切需要先进的地球观测(EO)方法来一致且可扩展地捕捉北极-北方湿地的时空变化。
遥感为解决北极-北方湿地测绘挑战提供了关键框架,能够进行系统化、可重复且空间范围广泛的观测,克服了偏远北部地区实地调查的物流和财务限制(Mahdianpari等人,2017年)。最近卫星任务(如Landsat、Sentinel-1和Sentinel-2)的进步通过提供广泛地理范围内的高质量数据,极大地提升了湿地监测能力。Sentinel-1雷达观测特别宝贵,因为它能够穿透云层并在极夜期间独立运行(Torres等人,2012年)。这些能力与Sentinel-2的多光谱信息相结合,提供了关于植被和地表水动态的详细见解(Drusch等人,2012年)。这些平台支持高频监测,捕捉水文过程、植被物候和冻融转换的季节性和年际变化(Bartsch等人,2017年;Bartsch等人,2020年;Phiri等人,2020年)。在北极-北方环境中,生成密集的时间序列数据尤为重要,因为快速变化需要连续可靠的观测(Mohammadimanesh等人,2018年;Pastick等人,2019年)。此外,这些任务的广泛空间覆盖有助于将局部尺度动态整合到区域和环极评估中,生成用于稳健碳核算和改善气候模型中北极-北方系统表示的统一数据集(Ludwig,2024年)。
最近在EO领域的进展促进了高分辨率湿地清单的开发,例如加拿大湿地清单地图(CWIM),它通过云计算平台整合多传感器卫星图像来提高分类准确性和主题细节(Mahdianpari等人,2021年)。虽然这类产品为生态研究和政策提供了有价值的基准,但它们本质上是静态的,仅提供系统的高度动态性的时间快照。水文周期、植被生长和冻融转换的季节性和年际变化往往未被反映,导致统计分类错误、地图边界简化以及湿地面积估计不准确(Coons等人,2025年)。在气候驱动的水文重组和地貌变化背景下,这些问题日益严重。解决这些问题需要转向动态监测框架,整合时间序列EO数据,以实现更可靠的湿地类型区分及其自然变化和生态过程的准确表示。
最近,湿地测绘研究越来越多地结合来自遥感时间序列数据的动态信息,从而提高了地图的准确性(Marjani等人,2025b)。例如,DeLancey等人(2022年)使用密集的Sentinel-2图像序列捕捉了阿尔伯塔省湿地地区的洪水频率,识别出临时性、季节性或半永久性湿地。这些信息被整合到分类图中,比之前的清单工作识别出更多的湿地面积。类似地,Peng等人(2023年)使用来自Sentinel-2的湿度和绿色度光谱指数时间序列开发了亚洲部分地区的湿地分类框架。通过在时间密集的训练像素上训练随机森林(RF),最终14类分类图的总体准确率为90%。
基于这些研究工作,本研究通过整合多任务卫星数据的时间序列动态分类框架,解决了静态湿地测绘的局限性。主要目标是通过结合时间特征(特别是来自合成孔径雷达(SAR)和光学数据捕获的淹没模式)来提高北极-北方湿地地图的统计准确性。我们假设,结合淹没和物候的时间特征将比静态湿地测绘实现更稳健的湿地类别区分。本研究并不旨在比较不同的分类器或分割范式,而是探讨在保持分类框架其他组件不变的情况下,动态淹没描述符的整合是否能够提高湿地分类性能。为此,我们在相同的实验设置下实施了动态和静态分类情景的对照比较。最终,本研究提供了一种更准确地捕捉北极-北方湿地内在动态特性的湿地监测框架,为区域碳核算和气候影响评估提供了更可靠的数据。
**2. 数据和研究区域**
**2.1. 研究区域**
两个500公里×500公里的研究区域分别位于加拿大西北地区(NWT)(图1a)和瑞典北部的诺尔博滕(Norrbotten)和韦斯特博滕(V?sterbotten)县(图1b)。NWT区域跨越两个生态区:泰加平原和南部北极(生态系统分类组,2009年),以低地、广阔的多边形泥炭高原和连续的永久冻土为特征。起伏的地形中散布着布满巨石的岩石平原。该地区的年平均温度较低,降水量也低于南部相邻地区。马更些三角洲(Mackenzie Delta)位于北部,是加拿大最大的三角洲。数千条河流、小池塘和沼泽湿地遍布三角洲附近区域。三角洲南部则以河流、浅水池和热喀斯特湖泊为主。泥炭地广泛分布,但泥炭层通常较浅。
**2.2. 地理空间数据**
分类框架整合了多传感器卫星观测和辅助产品,以捕捉北极-北方湿地的光谱和物理特征(表1)。Sentinel-2的多光谱信息(MSI)提供了可见光、近红外(NIR)、红边(Red Edge)和短波红外(SWIR)波段的10-20米反射率,从中衍生出光谱指数,以增强对水、饱和植被和高地苔原的区分能力。Sentinel-1的C波段SAR数据通过提供独立于云层和光照条件的观测(包括VV和VH反向散射以及VV/VH比率)补充了光学数据。这些特征能够捕捉土壤湿度、淹没、植被结构和冻融动态。在可用情况下,L波段SAR(ALOS PALSAR-2)提供了对植被覆盖下地下土壤湿度和持续淹没的额外敏感性。来自ArcticDEM(亚5米分辨率)的地形数据支持坡度和朝向的确定,而Multi-Error-Removed Improved-Terrain(MERIT)Hydro数据提供了流量累积、排水密度和流域边界,提供了水文背景。欧洲中期天气预报中心(ECMWF)再分析版本5(ERA5)的气候协变量(温度和降水量)用于表示季节性和年际变化,这对于将湿地状态与甲烷排放联系起来至关重要。训练和验证使用了结合实地观测和系统化桌面划分的湿地和高地多边形的综合数据集(Merchant等人,2022年)。
**表1. 本研究中静态和动态湿地分类工作流程中整合的特征**
| 类型 | 来源数据 | 特征列表 |
| --- | --- | --- |
| 光谱特征 | Sentinel-2(反射率数据) | 10(蓝、绿、红、红边1、红边2、红边3、NIR、红边4、SWIR 1、SWIR 2) |
| Sentinel-2(光谱指数) | 8(归一化差异水指数(NDWI)、归一化差异植被指数(NDVI)、红边NDVI(ReNDVI)、归一化燃烧比(NBR)、增强植被指数(EVI)、土壤调整植被指数(SAVI)、归一化差异雪指数(NDSI)、绿色叶绿素植被指数(GCVI) |
| SAR特征 | Sentinel-1,C波段(VV + VH) | 5(VV极化反向散射系数、VH极化反向散射系数、跨度、VV/VH比率) |
| ALOS PALSAR-2,L波段(HH + HV) | 5(HH极化反向散射系数、HV极化反向散射系数、跨度、HH/HV比率) |
| 地形特征 | ArcticDEM | 3(海拔、坡度、朝向) |
| 环境特征 | ERA5 | 2(温度、降水量) |
| 水文特征 | MERIT Hydro | 2(最近排水口高度(HAND)、流量) |
表1中变量的选择基于先前关于加拿大湿地分类的研究(Mahdianpari等人,2020b;Mahdianpari等人,2020a),其中通过RF进行了变量重要性分析,以确定这些变量对静态分类的相对贡献。
**2.3. 参考数据**
准确的、丰富的测试和训练数据是监督分类方法的必要输入。理想的测试和训练数据通过基于地面的实地数据收集活动获得;然而,在北极-北方地区,由于可访问性、预算和机会等原因,这类数据较为有限。为解决这一数据缺口,在之前的研究中,已在NWT和瑞典研究区域采用了桌面划分方法进行湿地分类(Merchant等人,2021年;Merchant等人,2023年),遵循Merchant等人(2022年)描述的方法。为了增加该项目可用的测试和训练数据的分布和数量,使用相同的桌面划分方法进一步划分了湿地和湖泊多边形。首先,在湿地密集且之前未进行桌面划分的区域分布了样本点。在样本区域上,下载了多季节和无云的Sentinel-2图像,并使用这些图像创建了自然色和假彩色合成图像。还下载了样本区域上的ArcticDEM瓦片,用于计算地形湿度指数(TWI)。Sentinel-2彩色合成图像和TWI被结合使用,以识别、划定和分类每个样本区域内的湿地和湖泊区域。为了防止湿地面积被高估,使用了欧洲航天局(ESA)的WorldCover 2021栅格数据集作为参考,创建了高地(非湿地)多边形。高地区域包括苔原、灌木丛、草地/农田、森林、建筑区/裸地、积雪和烧毁的土地。图2中突出显示了一些划定的湿地和高地多边形的视觉示例。
下载:下载高分辨率图像(979KB)
下载:下载全尺寸图像
图2. 测试和训练网格中划定的湿地和高地多边形。
最终数据集中包含的湿地、水域和高地(非湿地)多边形的总数在表2中列出。空间验证采用了基于网格的分区策略,网格单元为30公里×30公里。整个网格单元被分配到训练子集或测试子集中,以防止空间泄漏并减少相邻样本之间的空间自相关性。70/30的比例划分在网格级别应用,以确保网格单元内的所有样本都被分配到同一个子集中。鉴于研究区域内湿地类型的广泛地理分布,所有主要类别在训练和测试子集中都有代表。
表2. 使用桌面划定方法在瑞典和NWT划定的湿地和水域多边形。
类别/研究区域
瑞典 NWT
沼泽 50 35
湿地 94 53
沼泽地 17 32
湿地 22 52
高地 40 41
总计 63 86
水域 49 64
高地 20 36
3. 方法
本研究评估了两种互补的北极-北方湿地制图方法:静态基线和动态情景。这两种方法都使用了相同的多传感器数据源,包括Sentinel-1 SAR和Sentinel-2光学图像,以进行对照比较。两种情景之间的关键区别不在于输入数据,而在于时间信息的表示和利用方式。静态情景依赖于2021年至2025年的多传感器中值夏季合成图像,将时间变化压缩为单一的时间聚合表示,因此无法捕捉季节内或年际的变化。相比之下,动态情景通过结合Sentinel-1时间序列后向散射和Sentinel-2植被指数来保留和利用时间信息,从而推导出水文特征,如淹没频率、最大淹没范围和每次采集时的淹没状态。这两种情景使得能够评估将淹没动态纳入湿地分类所带来的好处。采用基于对象的图像分析(OBIA)-RF分类方法对两个研究区域进行湿地制图(图3)。工作流中的彩色部分表示静态和动态湿地分类方法之间的不同步骤。
下载:下载高分辨率图像(247KB)
下载:下载全尺寸图像
图3. OBIA-RF分类工作流。蓝色框突出显示了用于生成动态湿地分类结果的额外动态特征信息。(有关此图例中颜色参考的解释,请参阅本文的网络版本。)
3.1. 静态湿地制图
进行静态湿地制图是为了提供一个基准产品,以便评估动态框架。在这种方法中,提取了2021-2025年期间的多传感器卫星特征作为夏季中值合成图像,捕捉生长季节高峰期的植被和地表水状况。本研究中的夏季指的是6月至9月(第6至9个月),这对应于两个研究地点的主要无雪生长季节。使用Sentinel-2光学反射率和衍生的光谱指数来区分湿地植被和高地苔原及水域。Sentinel-1双极化(VV, VH)后向散射和比率指标被用来表征土壤湿度、淹没潜力和植被结构。此外,还结合了从ArcticDEM(高程、坡度、朝向)获得的地形变量,以及来自MERIT Hydro(水流积聚)的水文和上下文属性,以及来自ERA5(温度、降水量)的气候协变量。这些输入共同提供了影响湿地分布的光谱、结构和环境条件的静态表示。
与为每个像素独立分配类别标签的基于像素的方法不同,基于对象的方法首先将图像分割成均匀的单位(对象),然后使用统计和上下文信息对这些对象进行分类。这种策略减少了像素级别的噪声,结合了空间背景,并产生了更具生态意义的分类结果。OBIA与RF的结合为在复杂的北极-北方景观中区分湿地类型提供了一个强大而灵活的框架。
3.1.1. 基于对象的RF分类
OBIA工作流从分割开始,将卫星图像分割成光谱和空间上相似的像素组。在本研究中,使用简单非迭代聚类(SNIC)算法进行分割,这是一种基于超像素的方法,能够高效生成紧凑、均匀的对象,同时保持重要的边界。SNIC从一组均匀分布的种子点开始聚类,并根据组合的光谱相似性和空间邻近性将像素分配到簇中,直到所有像素都被标记。生成的超像素在空间上连贯且计算效率高,非常适合大规模湿地制图(Achanta和Süsstrunk,2017)。为了捕捉生态变化,在具有细粒度特征的区域内生成了较小的超像素,如融池和狭窄的沼泽走廊,而在广阔、更均匀的区域(如沼泽和湿地)则使用较大的超像素。这种自适应缩放确保了分割产生的对象与自然湿地单元紧密对应,从而提高了后续分类的生态可解释性。
分割采用了四连通性(4-connectivity)和紧凑度参数1,这平衡了光谱相似性和空间邻近性,同时保持了湿地边界。而不是使用单一的全局超像素大小,而是应用了自适应的多尺度分割策略来捕捉研究区域内异质北极-北方景观的生态变化。初始种子间距值大约在10到30像素之间。在空间异质性较高的区域(如洪泛区、沼泽-水域交界处、融池、狭窄的沼泽边缘和热喀斯特地形)使用了更细的间距(10到15像素)。在以广阔的沼泽、湿地和高地为主的均匀景观中使用了较粗的间距(25到30像素)。自适应缩放由景观背景指导,包括局部湿地密度、Sentinel-1后向散射变化和从ArcticDEM获得的地形变化,确保分割产生的对象与自然景观单元相对应。
然后通过一套涵盖光谱、雷达、地形、上下文和几何领域的综合属性来描述每个对象(表1)。光谱统计包括Sentinel-2波段及其衍生指数(NDVI、NDWI、EVI、ReNDVI)的平均值、中值和方差。基于雷达的特征捕捉了VV和VH、HH和HV的平均后向散射,以及VV/VH比率和HH/HV比率指标。地形变量包括高程、坡度、朝向和HAND,提供了与湿地形成相关的水文背景,包括与水道的接近程度、到永久水域的距离(PW)、与高地类别的邻近性、局部邻域统计和流域位置指标。几何描述符包括面积、周长、形状指数和紧凑度,有助于区分湖泊、泥炭地多边形和线性沼泽走廊。
对象的分类使用了RF算法,该算法从训练数据的随机子集和预测特征构建决策树集合(Breiman,2001)。每个对象的特征向量作为输入,最终类别分配是通过所有树的多数投票确定的。超参数如树的数量(Ntree)和每次分割考虑的变量数量(Mtry)通过交叉验证进行优化,Ntree设置为500,Mtry设置为预测变量数量的平方根。RF算法还提供了特征重要性的度量,有助于评估哪些对象属性(例如光谱指数、雷达指标、地形导数)对分类准确性贡献最大。
3.2. 动态湿地制图
淹没是区分湿地和湖泊的关键特征,它塑造了水、植被和土壤之间的相互作用,并影响生态系统过程,如栖息地可用性、养分循环和温室气体排放(Marjani等人,2025a)。在北极-北方地区,淹没非常动态,由于融雪、降雨、河流洪水或永久冻土融化,淹没区域会迅速形成和消退。未能捕捉到这些变化可能导致湿地被低估或错误分类,这对甲烷排放估计和生态系统建模有重要影响(Sharma和Naik,2024)。因此,准确的淹没地图对于改进气候预测、生物多样性评估和水资源管理至关重要。
为了解决这一需求,开发了一个多传感器制图框架,将Sentinel-1 SAR时间序列与Sentinel-2光学数据结合起来。虽然静态情景中也使用了这两种数据源,但动态框架在提取和利用时间信息方面有所不同。特别是,Sentinel-1时间序列后向散射被分析以推导出捕捉淹没时间行为的水文特征,包括干涸、暂时淹没和永久淹没状态之间的转换(Zhang等人,2022)。这些特征是通过后向散射动态的时间统计和基于规则的分析获得的,而不是通过聚合合成图像获得的。Sentinel-2数据通过基于NDVI的掩蔽提供了补充的植被信息,支持植被水域和非植被水域的分离。这种组合实现了湿地水文的时间明确表示,其中Sentinel-1驱动淹没动态的检测,而Sentinel-2通过植被背景细化了类别解释。
3.2.1. Otsu阈值分割
在本研究中,首先将Otsu方法应用于最大VH后向散射的直方图,以自动分离暗色、低后向散射的水像素和较高后向散射的土地或植被表面。这是一种自适应的、数据驱动的阈值分割方法,避免了主观性,并确保了场景中PW的一致检测。
Otsu方法通过最大化类别间方差或等效地最小化类别内方差来自动确定直方图中两个类别的最佳阈值。对于给定的阈值t,图像直方图被分为两个类别:C1表示低于t的像素,C2表示高于t的像素。类别概率和均值使用公式(1)-(4)定义:
(1) w1t = ∑i=0tp(i)
(2) w2t = ∑i=t+1L-1p(i)
(3) μ1t = 1/w1(t) ∑i=0ti×p(i)
(4) μ2t = 1/w2(t) ∑i=0L-1i×p(i)
其中p(i)是归一化的直方图,w1和w2表示每个类别的像素计数,μ1t和μ2t表示由阈值候选t分隔的每个类别的像素的平均后向散射值,L是强度级别的数量。最佳阈值t*通过最大化类别间方差来选择,如公式(5)所示:
(5) t* = argmax(t [w1t - w2t(μ1t - μ2t)2]
Otsu阈值一致地识别出暗色(湿润)像素,这些像素用于后续步骤中定义淹没范围。
3.2.2. 最大淹没范围
提出的淹没制图框架的第一步是确定观察期间水存在的最大范围(图4)。此步骤的目标是建立一个掩码,包含至少被淹没过一次的所有像素,从而定义研究区域内的潜在湿地范围。为此,计算了每个像素在时间序列中的Sentinel-1后向散射的时间统计,重点关注VV和VV极化。具体来说,推导出VH后向散射的最大值、最小值和平均值,以及最大VV后向散射值,以表示每个像素随时间经历的所有湿干条件。然后使用三个决策规则将像素分类为淹没类别,这些规则利用VH和VV信号对不同水文和植被条件的敏感性。
下载:下载高分辨率图像(139KB)
下载:下载全尺寸图像
图4. 使用Sentinel-1时间序列数据的提出的淹没制图工作流。
第一条规则将始终具有暗色VH后向散射的像素识别为PW。相应的最大VH(VHmax)阈值来自第3.2.1节中介绍的Otsu方法,满足VHmax3db则k=2,否则k=3。VHmin≤T2的像素被标记为TW。缩放因子k控制了阈值对PW后向散射变化的敏感性。k=2和k=3的值是根据先前的基于SAR的洪水研究通过实证方法确定的,以考虑场景间的噪声和局部异质性(Huang等人,2022年)。当后向散射(PW)的标准差(δ)较大(>3dB)时,较小的乘数(k=2)可以防止将瞬态低后向散射像素错误地分类为水域。相反,当δ较小时,k=3提供了一个更为保守的阈值,能够包含低方差环境中的较弱水信号。这种自适应公式在各种雷达条件下平衡了识别暂时性水域(TW)像素时的遗漏和误分类错误。第三条规则用于隔离浅水上的灌木或木质植被,这些地方会发生强烈的双折射现象。该规则基于方程(7):(7)Δ=VVmax-VH。Δ≥12 dB的像素被标记为暂时性被淹没的植被(TFV)。方程(7)中使用的12 dB偏移量是一个基于实证但具有物理意义的阈值,用于衡量C波段SAR数据中被淹没的木质或灌木植被的强烈双折射散射特征。在洪水条件下,由于树干与水的相互作用,VV极化后的后向散射显著增加,通常会产生10-15 dB的对比度。根据之前的基于SAR的洪水映射研究(Santoro和Wegmüller,2014年;White等人,2015年),采用了12 dB的固定阈值来识别具有这种增强VV-VH差异的像素。虽然12 dB的值并不具有普遍性,但在具有相似植被和雷达采集参数的北方和北极湿地环境中已被证明是可靠的。最后,任何未被归类为上述洪水类别的像素都将自动被归类为干地(DL)。空间清理步骤会移除距离已确认的湿地超过30米的孤立TFV像素,并过滤掉面积小于900平方米的被淹没区域。由此产生的类别,如PW、TW、TFV和DL,代表了在时间序列中至少被淹没过一次的所有位置。这些改进产生了一个稳健的四类洪水掩膜(PW、TW、TFV、DL),为北极-北方湿地的动态场景级分类提供了可靠的基础。
3.2.3 洪水分类
所提出方法的第二步从静态空间范围转向动态监测,通过为每次Sentinel-1采集的数据分类洪水状况(图5)。与第一步定义整个时间序列中的最大洪水范围不同,第二步专注于捕捉湿地和湖泊水文的时间变化。重要的是,这种分析仅限于第一步中确定的潜在湿地像素掩膜,确保分类结果仅限于那些有明显洪水可能性的位置。这种设计大大降低了在高地或永久干燥地形中出现误报的风险。
下载:下载高分辨率图像(301KB)
下载:下载全尺寸图像
图5. 使用组合的Sentinel-1和Sentinel-2时间序列数据进行的动态洪水分类工作流程。
在第二步中,同时使用了雷达后向散射阈值和Sentinel-2的NDVI光学植被信息。从生长季节获取的Sentinel-2场景中创建年度最大NDVI合成图像,并使用0.30的阈值(T4;见图5)来区分植被覆盖表面和非植被表面。这个植被掩膜对于将浮游植物(FPs)和挺水植物(EPs)与PW区分开来至关重要。在第一步的湿地掩膜中,每个新的Sentinel-1场景都是通过将VH后向散射与从PW样本中得出的阈值进行比较来分类的。对于每次采集,计算PW像素中VH的平均值(η)和标准差(δ)。然后在方程(8)中定义一个上限阈值:(8)T5=η+k.δ,其中如果δ>3db则k=2,否则k=3。VH值低于此阈值的像素被视为类似水的区域。低于T5的非植被像素被标记为PW,而低于同一阈值的植被像素被标记为FP。为了进一步区分EPs,定义了第二个阈值,如方程(9)所示:(9)T6=T5+3db。VH值在T5和T6之间的植被像素被分类为EPs,这反映了由于茎秆从浅水中伸出而产生的较强后向散射。TFV是通过基于最大VV和平均VH之间差异的滞后方法检测到的。具体来说,该指标如方程(10)所示:(10)Δ=VV-VHμ。它与一对阈值进行比较,一个高阈值等于固定的12 dB(T3),另一个低阈值设置为T3-δVV,其中δVV是VV的标准差。超过高阈值的像素直接被分类为被淹没的植被(FV),而低于低阈值的连续像素片段被保留下来,以保持FV片段的连续性并最小化由斑点噪声引起的碎片化。任何不符合上述洪水类别的像素在该场景中自动被归类为干地(DL)。
空间清理步骤会移除距离已确认湿地超过30米的孤立TFV像素,并过滤掉面积小于900平方米的被淹没区域。最终得到的类别,如PW、TW、TFV和DL,代表了在时间序列中至少被淹没过一次的所有位置。这些改进产生了一个稳健的四类洪水掩膜(PW、TW、TFV、DL),为北极-北方湿地的动态场景级分类提供了可靠的基础。
3.2.3 洪水分类
所提出方法的第二步通过为每次Sentinel-1采集的数据分类洪水状况,从静态空间范围转向动态监测(图5)。与第一步定义整个时间序列中的最大洪水范围不同,第二步侧重于捕捉湿地和湖泊水文的时间变化。重要的是,这种分析仅限于第一步中确定的潜在湿地像素掩膜,确保分类结果仅限于那些有明显洪水可能性的位置。这种设计大大降低了在高地或永久干燥地形中出现误报的风险。
下载:下载高分辨率图像(301KB)
下载:下载全尺寸图像
图5. 使用组合的Sentinel-1和Sentinel-2时间序列数据进行动态洪水分类的工作流程。
在第二步中,同时使用了雷达后向散射阈值和Sentinel-2的NDVI光学植被信息。从生长季节获取的Sentinel-2场景中创建年度最大NDVI合成图像,并使用0.30的阈值(T4;见图5)来区分植被覆盖表面和非植被表面。这个植被掩膜对于将浮游植物(FPs)和挺水植物(EPs)与PW区分开来至关重要。在第一步的湿地掩膜中,每个新的Sentinel-1场景都是通过将VH后向散射与从PW样本中得出的阈值进行比较来分类的。对于每次采集,计算PW像素中VH的平均值(η)和标准差(δ)。然后在方程(8)中定义一个上限阈值:(8)T5=η+k.δ,其中如果δ>3db则k=2,否则k=3。VH值低于此阈值的像素被视为类似水的区域。低于T5的非植被像素被标记为PW,而低于相同阈值的植被像素被标记为FP。为了进一步区分EPs,定义了第二个阈值,如方程(9)所示:(9)T6=T5+3db。VH值在T5和T6之间的植被像素被分类为EPs,这反映了由于茎秆从浅水中伸出而产生的较强后向散射。TFV是通过基于最大VV和平均VH之间差异的滞后方法检测到的。具体来说,该指标如方程(10)所示:(10)Δ=VV-VHμ。它与一对阈值进行比较,一个高阈值等于固定的12 dB(T3),另一个低阈值设置为T3-δVV,其中δVV是VV的标准差。超过高阈值的像素直接被分类为被淹没的植被(FV),而高于低阈值的连续像素片段被保留下来,以保持FV片段的连续性并最小化由斑点噪声引起的碎片化。任何不符合PW、FP、EP或FV标准的像素在第一步湿地掩膜内的位置都被分配为DL。因此,第二步产生的每次采集的地图包含五个独特的类别:PW、FP、EP、FV和DL。通过随时间聚合场景级分类,该框架生成了季节性和年际产品,包括最大洪水范围和洪水频率地图,这对于理解水文周期、短暂湿地以及北极和北方湿地动态的变化至关重要。所有阈值参数、方程及其确定标准都在表3中进行了总结,以提高分类工作流程的透明度和可重复性。阈值是通过数据驱动算法和实证校准相结合的方法确定的。应用了Otsu方法来自动化识别雷达后向散射直方图中的双峰分离。
表3. 用于基于后向散射的洪水分类的阈值参数总结
阈值 | 描述/类别类型 | 参数值或规则
|-------|-------------|-------------------------|
| T1 | PW | 数据驱动(Otsu) |
| T2 | TW | T2=η+k.δ(如果δ>3 dB,则k=2,否则3) |
| T3 | TFV | Δ≥12 dB |
| T4 | NDVI阈值 | NDVI≥0.30 |
| T5 | FP | T5=η+k.δ |
| T6 | EP | T6=T5+3dB |
| T7 | FV | 滞后规则中的下限阈值 |
需要注意的是,第二步针对每次Sentinel-1场景独立进行洪水分类。时间变化并没有直接在单个分类中建模,而是通过将这些场景级输出聚合到后续阶段的多时相描述符中来引入的。设置的阈值包括固定和自适应两部分。PW阈值(T1)遵循Otsu对最大VH的优化,因此适应每个站点的直方图。TW阈值(T2)和每个场景的类似水阈值(T5)使用相同的自适应形式,其中η和δ来自每个场景中的PW像素,因此这两个值会根据当地的后向散射统计数据进行调整。相比之下,NDVI阈值(T4=0.30)和FV的双折射阈值(T3=12 dB)在各个站点保持固定,这是基于之前的SAR洪水研究和北方及北极湿地中常见的植被-水对比。这种自适应和固定阈值的组合在保持局部灵活性的同时,提供了物理上可解释的决策规则。
3.2.4 在OBIA-RF湿地分类中整合动态特征
洪水输出作为额外的对象级属性输入到分类系统中(图3)。重要的是,动态框架并不修改底层的OBIA-RF分类算法,而是通过从Sentinel-1时间序列分析中派生的时间描述符来扩展特征空间。第二步独立地为每次采集生成洪水地图,然后通过聚合这些场景级分类来引入时间变化。对于每个分割的对象,从2021年5月到2025年9月的完整Sentinel-1时间序列中计算水文描述符。洪水频率定义为在某个对象内,像素被分类为PW、FP、EP或FV的Sentinel-1采集次数,该次数除以该对象可用的有效采集总数(不包括受缺失或无效观测影响的日期)。最大洪水比例表示在整个时间序列中,任何单次采集中该对象区域被分类为被淹没(PW、FP、EP或FV)的最大比例。PW、TW和FV的存在量量化为在第一步最大洪水掩膜内分配给每个类别的对象像素的比例,反映了长期的水文行为,而不是特定场景的条件。与依赖时间聚合的光谱和雷达合成图像的静态场景不同,这些动态属性明确编码了洪水特征的持久性、重复性和变化性。通过将这些描述符添加到统一特征集中的静态光谱、雷达和地形特征中,OBIA-RF模型将时间行为纳入对象级决策空间,而不会改变分类算法本身。
3.3 湿地分类准确性评估
使用来自实地调查的独立验证数据和先进的桌面方法评估了OBIA-RF分类的准确性。为了避免空间自相关性,验证样本是从未用于训练的区域中选取的。计算了Kappa系数、类别特定的用户(UA)和生产者(PA)准确性来评估性能。准确性评估按北极-北方地区分层,以识别性能的空间变化,并确保在不同环境和生态条件下的分类结果的稳健性。
4. 结果与讨论
4.1 北极-北方湿地中的洪水空间和时间模式
洪水分类产生了NWT和瑞典北极-北方湿地中永久性和季节性被淹没区域的空间明确地图(图6),区分了四个类别:PW、TW、TFV和DL。在两个研究地点,大型湖泊和河流通道都被映射为PW,形成了每个景观的长期水文特征。在两个地区,TW和TFV区域环绕着永久性水体的边缘和低洼地区,代表了从永久性水体到更动态湿地带的特征性过渡。特别是TFV像素在湖泊边缘、浅洼地和河流通道周围较为普遍,表明这些区域是周期性洪水的植被带,通常对应于沼泽和湿地湿地的位置。
下载:下载高分辨率图像(1003KB)
下载:下载全尺寸图像
图6. 加拿大NWT(a)和瑞典(b)湿地的第一级洪水分类。
两个地点之间的一个关键区别在于洪水模式的空间范围和异质性。在瑞典,PW是水文景观的主要组成部分,代表了研究区域内的许多稳定湖泊。通常,TW和TFV在空间上受到限制,局限于泥炭地边缘和低洼盆地,反映了研究区域内以aapa和palsa沼泽(沼泽和泥炭湿地)为主的地形。沼泽和泥炭湿地的水文条件相对稳定,季节性波动较小,受降水、逐渐的水位变化、土壤饱和度和地下水排放控制。这些湿地通常不位于水体边缘,洪水较少,水文变化也不如沼泽或湿地湿地那么剧烈。一些TFV像素出现在湖泊边缘和浅洼地,表明这些区域是周期性洪水的植被带。
相比之下,NWT研究区域表现出由马更些河三角洲、马更些河及其在研究区域内众多分支河流驱动的相对异质性和广泛的洪水动态。在这里,辫状河流网络和牛轭湖占据了大部分PW区域,而相邻的洪泛平原支持广泛的TW和TFV区域。TFV类别在河道边缘和河岸洪泛平原尤为突出,那里的沼泽和湿地湿地经常发生洪水。与瑞典相比,NWT有更多的大型、广阔的沼泽和湿地湿地,由于三角洲的存在,TW和TFV的洪水更为普遍。为了评估不同生态单元随时间受洪水影响的行为,分析了5年内PW、水生植被(AV;包括FP和EP)和FV的检测频率。在这种情况下,AV对应于在第二步分类中被标记为FP或EP的植被覆盖水面。此外,频率指的是每个像素在分析期间所有可用的Sentinel-1采集中被分配给给定动态洪水类别(例如PW、FP、EP、FV)的次数。这种基于频率的分析突出了多年生水体、植被-水相互作用以及塑造北极和北方湿地景观结构和功能的短暂洪水过程之间的相互作用。
PW频率(图7)在NWT和瑞典都显示出稳定的模式。高频(持久性)水体包括大型河流和湖泊。持久性较低的区域显示了季节性被淹没的洪泛平原和较小的通道。在NWT,辫状通道和洪泛平原湖泊的异质性更大,而在瑞典,PW区域的异质性较小。
下载:下载高分辨率图像(938KB)
下载:下载全尺寸图像
图7. 2021年至2025年间NWT和瑞典地点的PW检测频率。
AV的变异性高于PW,并且对水位波动更敏感。在NWT,AV频繁出现在近岸和洼地。然而,在NWT的洪泛平原,AV频率处于中间水平,表明植被受到一定程度的但不是连续的洪水影响。在瑞典,AV在空间上局限于沿海地带和低洼盆地。对AV的频率分析显示,这些区域受到水文变化的影响非常大,尤其是在沼泽和湿地类别中。FP在空间上受到限制,在两个研究区域都表现出高度局部化的频率模式。在NWT,FP集中在受保护的回水区和牛轭湖中。在易受季节性侵蚀和临时淹没影响的洪泛区边缘,FP的频率适中。在瑞典,FP在平静的水域中频繁出现。尽管FP的存在范围有限,但它们在生态上非常重要,因为它们创造了缺氧环境,从而支持甲烷的产生和水生生物的生产力。FV是本节描述的湿地类型中分布最广的。在NWT,AV频率揭示了反复受到淹没和干燥影响的洪泛区植被,这表明流动的水道与湿地之间有着密切的关系。在瑞典,FV存在于泥炭地的边缘和低洼盆地中。频率分析揭示了FV在维持水文变化、支持适应交替有氧和厌氧条件的植被群落以及调节北极-北方湿地碳循环方面所起的生态重要作用。这些变化多端的模式展示了北极-北方洪泛区和泥炭地的不同水文行为。NWT的洪泛区以永久性和临时性淹没的高变异性为特征,这种变异性受到河流连通性和季节性洪水的影响。在瑞典,FV占主导地位,反映了泥炭地地形和土壤饱和度对淹没动态的影响。这些差异说明了淹没类别的持续性和变异性如何有助于理解生态系统功能、恢复力以及对温室气体通量的潜在贡献。
4.2. 静态湿地分类与动态湿地分类
在本文中,为了公平比较静态和动态湿地分类方案,除了动态淹没描述符之外的所有实验组件都保持不变。两种模型都使用了相同的输入数据、OBIA分割、特征提取程序、分类器配置和空间独立的验证策略。因此,任何在OA、类别F1分数或类别混淆度上的改进都可以直接归因于动态淹没信息的纳入,而不是分类器设计、分割方法或采样协议的差异。变量重要性分析(图3)还揭示了纳入静态和动态分类工作流的特征的相对贡献(表1)。变量重要性是通过计算所有随机森林模型中每个变量对节点纯度(基尼指数)减少的平均值来确定的,该指数量化了每个变量对节点纯度减少的总体贡献。在节点分割过程中导致基尼纯度较大减少的变量被认为在提高分类性能方面更具影响力。如图8所示,来自Sentinel-2的光学特征被证明是最有影响力的,特别是红边和红外波段,它们对植被特征(如叶绿素浓度)非常敏感(Kaplan和Avdan,2019)。因此,光学波段有效地捕捉了湿地植被的特征,这些特征是湿地类别的定义组成部分(国家湿地工作组,1997)。动态信息是湿地类别的另一个定义特征,在特征重要性方面同样很高,与光学波段一起排名在前半部分。相比之下,来自SAR和水文及气候数据集的变量影响相对中等至较低。地形特征的影响介于光学和SAR及环境特征之间。这些变量通过提供有关表面结构、水文条件和景观梯度的额外信息,补充了光学数据,从而增强了湿地分类的OA和鲁棒性。
鉴于动态信息的重要性,将这些特征整合到湿地分类工作流中应该能够更准确地、生态上有意义地描述湿地类别,捕捉到Sentinel-1和Sentinel-2检测到的湿地行为的全谱。为了比较使用静态湿地分类工作流和整合动态信息的工作流产生的结果,对两种结果进行了准确性评估。固定的分割、分类器配置和空间验证协议确保报告的改进具体归因于添加的动态淹没描述符,而不是模型结构或采样设计的差异。表4、图9和图10分别报告了静态和动态结果的混淆矩阵和OA。
表4. 静态和动态湿地分类结果的OA。
分类结果
OAKappa系数
静态NWT 92.71 10.86
静态瑞典 91.44 0.86
动态NWT 95.42 0.91
动态瑞典 96.38 0.92
下载:下载高分辨率图像(245KB)
下载:下载全尺寸图像
图9. NWT和瑞典静态和动态分类结果的生产者准确性。
下载:下载高分辨率图像(441KB)
下载:下载全尺寸图像
图10. 静态NWT(a)、瑞典(b)分类结果以及动态NWT(c)和瑞典(d)分类结果的混淆矩阵。预测类别列在Y轴上,实际类别(由测试数据定义)列在X轴上。图例描述了最终分类中正确识别的测试像素数量。
NWT的OA从静态结果中的92.71%增加到动态结果中的95.42%。在瑞典,OA从91.44%增加到96.38%(表4)。混淆矩阵的类别级分析(图9和图10)显示,动态框架带来的准确性提升在很大程度上依赖于类别。在NWT,植被湿地(沼泽、湿地、沼泽和湿地)的PA增加了7-18%,其中沼泽和湿地类别的改进最为显著。在瑞典,动态映射将生产者准确性提高了9-10%,而开放水域和高地类别在两种情况下都保持较高水平。
在静态NWT产品中,植被湿地像素的正确分类率为70-75%。在动态产品中,相同的植被湿地类别的正确分类率增加到91%。在瑞典,大多数静态植被湿地像素的正确分类率在84-92%之间,而在动态结果中,所有类别的正确分类率都超过了92%。在瑞典,测试数据与分类结果之间的一致性高于NWT。这可能部分是由于瑞典研究区域位于更南的位置,那里的平均温度较高,夏季的高温天数更多,这意味着植被有更多的时间生长并展现出EO数据可以检测到的独特特征。
在两个研究区域中,类别对之间的混淆现象都很明显,包括沼泽和湿地、湿地和水域以及高地和沼泽。鉴于这些类别对之间的植被和水文模式的相似性(国家湿地工作组,1997),这些错误是预期的。沼泽和湿地之间的混淆尤其常见,因为它们都是泥炭地类型,共享许多功能性的植被群组和物种(特别是贫瘠的沼泽和湿地),并且经常在泥炭复合体中紧密相关,边界模糊,特别是在瑞典的高地边缘(Granlund等人,2022;Kolari等人,2022)。类似的混淆也发生在以树木和灌木为主的沼泽湿地和森林高地之间。这些类别的植被结构相似,特别是在沼泽经历季节性干燥或高地森林出现在地形低洼处时。湿地和水域之间的混淆是由于它们之间的自然邻近性,其中挺水植被逐渐过渡到没有明显边界的开阔水域。尽管如此,动态结果中的误分类率低于静态结果。
重要的是,动态映射大大减少了这些光谱相似类别之间的系统性混淆。例如,在NWT,静态情景下的沼泽到水域的误分类率从19%降低到动态情景下的5%,而沼泽-湿地之间的混淆减少了25%以上。这些改进主要发生在水文动态环境中,静态光学特征在这里趋于一致。通过纳入淹没频率、最大淹没比例和动态淹没状态,分类器获得了明确的水文周期信息,从而能够更稳健地区分以植被为主的湿地,并显著减少分类错误。
虽然OA和类别准确性的提高展示了动态框架的总体益处,但仔细检查混淆矩阵表明,最显著的改进发生在水文相似的类别之间。特别是,在动态情景中,开放水域和季节性淹没湿地(例如,湿地和沼泽)之间的误分类,以及沼泽和湿地之间的误分类显著减少。这些改进是因为动态描述符明确编码了水文周期特征,如持续性和重复性,这些特征在静态光谱或雷达合成图像中并未体现。
从空间上看,这些改进在洪泛区边缘、湖边和过渡湿地区域最为明显,在这些区域,静态反射特征通常不明确,但水文行为存在显著差异。因此,动态淹没描述符通过将基于过程的信息引入对象特征空间,增强了类别的可分离性。虽然变量重要性指标证实了这些特征的贡献,但它们在减少具有强烈时间变异性的区域的类别混淆方面的作用更加明显。动态情景下OA的增加(高达+5%)以及植被主导湿地类别的特定一致性明显改进反映了动态淹没特征的附加值。它们在光谱植被特征单独不明确的情况下提供了对抗误分类的鲁棒性,例如在湿地-水域边界或沼泽-湿地复合体中,静态光学反射率经常趋于一致。这些结果表明,明确表示水文周期变异性增强了分类器的信心,并支持更生态一致的湿地区分。通过将动态淹没作为已建立的湿地分类工作流中的一个特征,分类的准确性显著提高,改善了OA并增强了湿地分类的生态意义。
最近的湿地时间序列研究报道了不同的准确性结果,这些结果取决于地区、类别图例和传感器数据。对于Sentinel-2水文周期工作流,(DeLancey等人,2022)报告在多日期水检测和水文周期指标后,七个湿地类别的OA约为95-96%,并增加了冠层高度线索。对于基于物候的时间序列方法,(Peng等人,2024)报告了10种湿地和6种非湿地类型的OA为91.6% ± 1.2%,并且比基于像素的RF(88.3% OA)和简单的时间合成图像(85.7% OA)有所改进。对于基于SAR规则的动态淹没产品,(Huang等人,2022)报告使用Sentinel-1时间序列和NDVI掩码得到的PW和植被水状态图的OA为75-86%。相比之下,本研究中得到的动态OBIA-RF湿地分类产品在NWT的OA为95.42%,在瑞典为96.38%(表4),这与使用光学物候特征的最近时间序列基准相当或超过了它们,同时提供了明确的水文周期描述符(最大淹没、淹没频率)作为最终湿地类别的直接预测因子,而不是作为单独的后续产品。这一外部背景使得所提出的框架与当前的最佳实践时间序列湿地方法具有竞争力,特别是在植被主导的湿地类别方面,其中静态合成图像在沼泽-湿地和湿地-水域类型之间的混淆较高。
虽然本研究中使用的动态淹没阈值和决策规则在两个研究地点都展示了改进的性能,但一些参数(如Δ ≥ 12 dB的双重反射阈值和自适应k因子)可能仍然反映了当地的水文生态条件。这些参数在不同植被结构、土壤湿度制度或雷达采集几何形状的区域中的可转移性应在未来的工作中进行检验。本研究中使用的所有输入数据集都是公开可用的(Sentinel-1、Sentinel-2和表1中的辅助层),并且该工作流支持在云或高性能计算平台上进行区域到环极地的处理,通过可扩展的并行执行实现。这使得在大面积上进行可转移性测试成为可能,而无需新的野外活动。将验证扩展到其他北极-北方环境以及温带或热带气候的湿地,将有助于评估在更广泛环境变异性下的鲁棒性。此外,结合自动阈值学习,其中SAR-光学可分离性指导参数优化,将减少手动调整并支持完全可扩展的部署。因此,尽管当前框架提供了一种具有强站点级性能的可重复方法,但需要进一步的研究来评估参数的通用性,并加强动态湿地映射在多样化景观中的操作适用性。
5. 结论
湿地是水文、生物多样性和碳循环的关键调节器,然而它们的范围和动态性仍然是北极-北方系统科学中最大的不确定性来源之一。本研究展示了整合时间序列遥感数据在提高北极-北方湿地测绘的准确性和生态相关性方面的价值。通过比较两种方案——一种是基于2021至2025年多传感器夏季合成数据的静态方案,另一种是结合Sentinel-1后向散射时间序列和Sentinel-2植被指数来表征淹没行为的动态方案——研究发现,动态方法在多种北极-北方景观中表现出更优的性能。结果显示,动态产品的总体准确率在加拿大西北地区达到了95.42%,在瑞典北部达到了96.38%,而静态地图的准确率仅为91-92%。对于季节性被淹没和植被覆盖的地类,分类特异性准确率和精度分别提高了3-7%,这凸显了捕捉淹没过程及植被-水相互作用时间变化的优势。除了提高准确性外,动态框架还产生了具有生态意义的输出结果,包括最大淹没范围和频率地图,这些成果有助于更好地理解湿地水文和植被动态。这类产品对于碳和甲烷建模、生物多样性评估以及水资源管理尤为重要,因为准确反映湿地范围及其变化至关重要。同时,静态地图在缺乏密集时间序列数据的情况下,仍可用于建立一致的基准并支持长期监测。这两种方法共同为了解北极-北方湿地的分布和功能提供了互补的视角。将该框架扩展到环极地区,并引入额外的传感器(如L波段SAR时间序列),将进一步提升大规模湿地测绘的能力。将这些产品直接整合到生态系统过程模型和自下而上的甲烷建模中,有助于减少北极-北方碳预算的不确定性,并改善地球系统模型中对湿地动态的描述。在气候变化加速的背景下,动态且大规模地监测湿地的能力对于推进科学理解和政策决策至关重要。
**作者贡献声明:**
- Masoud Mahdianpari:撰写——审阅与编辑、撰写——初稿、项目监督、方法论设计、数据收集、数据分析、概念构建。
- Jean Granger:撰写——审阅与编辑、数据分析。
- Fariba Mohammadimanesh:撰写——审阅与编辑、方法论设计、数据分析、概念构建。
- Oliver Sonnentag:撰写——审阅与编辑。
- Mohammad Marjani:撰写——审阅与编辑、方法论设计。
生物通微信公众号
生物通新浪微博
今日动态 |
人才市场 |
新技术专栏 |
中国科学人 |
云展台 |
BioHot |
云讲堂直播 |
会展中心 |
特价专栏 |
技术快讯 |
免费试用
版权所有 生物通
Copyright© eBiotrade.com, All Rights Reserved
联系信箱:
粤ICP备09063491号