CN113971769B - 基于多源大数据的海岸带地域功能长时序识别方法 - Google Patents

Abstract

本发明涉及基于多源大数据的海岸带地域功能长时序识别方法，方法如下：首先，构建以人类生产、生活、生态为基础的海岸带地域功能用地分类，基于多源大数据融合平台采用随机森林算法对长时序影像进行初始分类。其次，利用扫描线种子填充算法和几何特征分析实现不透水面、水体的空间分离，并基于时空变化逻辑规则对耕地转变、近海养殖用地及盐田、围填海涉及的功能类型进行修正。最后，根据分类结果，提取粮食生产、近海养殖、围填海的变化范围和时间阶段。该方法提高了海岸带地域功能长时序识别精度，尤其适用于长时序、大范围、高密度的近海岸人类活动引发的地域功能用地转变与围填海建设的变化检测，能够直接应用于海岸带空间规划与区域政策的辅助决策。

Description

基于多源大数据的海岸带地域功能长时序识别方法

技术领域

本发明涉及基于多源大数据的长时序海岸带地域功能长时序识别方法，特别是涉及一种基于遥感大数据分类的近海岸人类活动引发的地域功能用地转变与围填海建设的长时序识别方法。

背景技术

海岸带区域是人口经济高度集聚且快速增长的区域，也是人地关系作用最为剧烈的区域，突出表现为海岸带陆地部分显著的土地利用变化以及近海海域的围填海活动。受人类日益增长的海产品需求以及临海经济发展的经济收益驱动，无论是海岸带陆域部分还是近海海域，除了传统的海水/混合水养殖、盐田等人类生产活动，围填海用于港口工业园区发展、滨海房地产开发、滨海旅游以及基础设施建设等，同样处于较为显著的扩张状态。随着遥感技术的飞速发展，遥感影像在时间、空间和光谱分辨率上都得到了极大的提高，凭借其广泛的观测范围、高效的数据获取方式和较低的经济成本，在海岸带资源调查、土地利用分类、景观功能监测等方面具有很大的潜力。如何能够快速、准确地对大区域的海岸带区域进行长时序、大范围、高密度的地域功能用地识别，一直是海岸带区域规划与政策研究亟待解决的基础性难题。

以往传统海岸带区域土地利用或地域功能变化检测的研究，普遍采用的是基于分类的变化检测方法。基于分类的变化检测方法，目前主要有两种，一种是较为传统和常用的分类后比较法或者双时态变化检测方法，一种是基于时间轨迹分析的变化检测。前者是先根据单时相影像进行分类，然后比较不同时相的分类结果获取变化检测信息，一般采用周年日期或周年窗口（年度周期或其倍数），多针对多地类要素整体分析。如2020年Liu等在《Ecological Indicators》撰文“Assessing and predicting changes in ecosystemservice values based on land use/cover change in the Bohai Rim coastal zone”，基于多时相土地利用数据（2000、2005、2010、2015年），对环渤海沿岸地区生态系统服务价值进行评估。后者是通过构建时间序列指数重构地物变化过程，利用其季节和周期性特征，通过检测时间突变点获取目标地物转化情况。如2016年Zhang等在《ISPRS Journal ofPhotogrammetry and Remote Sensing》撰文“Annual dynamics of impervious surfacein the Pearl River Delta, China,from 1988 to 2013, using time series Landsatimagery”，开发这一种逐年提取不投水面的有效方法，并应用于1988-2013年中国南部的珠江三角洲。然而，当进行大尺度长时序、高频率的土地覆盖分类与变化检测时，受制于处理平台的计算能力和储存设备，且需要付出巨大的时间成本，传统遥感手段显然不是理想平台。近年来，新开发的地理空间数据分析云平台Google Earth Engine（GEE，https://earthengine.google.com/），改变了传统的遥感处理方法，其庞大的遥感影像数据集和高性能的计算能力，为长时期、大规模的海岸带遥感变化监测分析提供了一种新途径。

借助于GEE平台进行长时序海岸带地域功能或土地利用功能识别与变化检测的研究已取得了一些突破性的进展，主要侧重于基于时间序列特征指数的单一地物或同类地物的研究，通过重构地物生产过程或利用季节性、稳定性等时间序列特征对海岸带特定地物或特定现象进行识别与变化检测，包括沿海滩涂、湿地、红树林、海岸线等特有的海岸带生态景观，以及养殖网箱、水产养殖池塘、围填海等多种沿海人类开发活动。利用GEE平台，对海岸带区域进行全覆盖、长时序土地利用功能谱系分类与变化检测的研究，已有一些初步的探索。如2018年Wang等在《Remote Sensing of Environment》撰文“Tracking annualchanges of coastal tidal flats in China during 1986–2016 through analyses ofLandsat images with Google Earth Engine”，通过计算自然地表水与植被覆盖海岸带的年度频率，得到1986-2020年中国东部沿海滩涂的逐年30m空间分辨率地图；2019年Ma等在《remote sensing》撰文“Change Detection of Mangrove Forests in CoastalGuangdong during the Past Three Decades Based on Remote Sensing Data”，利用基于决策树方法研究广东省1985-2015年30年红树林分布的时空特征。

在单一地物或同类地物在利用季节性、稳定性等时间序列特征进行识别与变化检测取得突出进展的同时，全覆盖的海岸带地域功能（或土地利用功能）分类仍侧重于陆域系统，在海陆统筹以及服务于海岸带地区可持续调控等方面需求仍有较大的差距。突出表现为：一是海岸带功能分类体系依然侧重于陆域土地覆盖/土地利用系统，缺乏对海陆统筹的整体考虑，尤其缺少对影响海岸带可持续发展的主要功能变化的关注；二是在解决海岸带这种复杂区域的分类制图与变化检测时，还需要进一步发挥大数据平台的多源数据融合的优势；三是长时序、大范围、高密度的海岸带地域功能识别与变化监测的研究较为薄弱，在已有的这些研究中有极少能够同时考虑到海岸带特有的地域功能类型时序间与空间邻域间的合理变化逻辑，严重制约了其识别与变化检测的精度。

本文提出了基于多源大数据的海岸带地域功能长时序识别方法，采用基于多源大数据融合的长时序、大范围、高密度地域功能用地识别方法，突出了沿海人类生产生活活动的功能特点以及海陆统筹的功能细分特征，实现了天然海水、天然水域与近海养殖用地及盐田，城镇建设用地、农村居民区与工矿及基础设施建设用地的空间分离，能够快速、准确地识别长时间尺度过程中近海岸地域功能转变与围填海建设变化特征。该方法提高了海岸带地域功能长时序识别精度，尤其适用于长时序、大范围、高密度的近海岸人类活动引发的地域功能用地转变与围填海建设的识别与变化检测，能够直接应用于海岸带空间规划与区域政策的辅助决策。

发明内容

本发明要解决技术问题是：如何利用多源大数据更为快速、准确地对大区域的海岸带区域进行长时序、大范围、高密度的地域功能用地识别，进而能够有效识别粮食生产功能、近海养殖用地及盐田、围填海的海岸带地区人类生产生活主要活动引发的地域功能用地变化时间和范围。

为了解决上述技术问题，本发明提出基于多源大数据的海岸带地域功能长时序识别方法，包括以下步骤：

步骤1：海岸带地域功能用地分类——从海岸带人类生产、生活、生态为基础构建海岸带地域功能用地分类，包括：城镇生产生活区、农村生活区、粮食生产功能区、近海生产区域、陆地生态功能区、海洋生态功能区6个一级分类，其中，城镇生产生活区包括城镇建设用地、工矿及基础设施建设用地，农村生活区为农村居民区用地，粮食生产功能区为耕地，近海生产区域为近海养殖用地及盐田，陆地生态功能区包括陆地天然水域、林地、草地、未利用地，海洋生态功能区包括海水水域、沿海滩涂，一共11个二级分类；；

步骤2：多源大数据的收集——借助Google Earth Engine云平台收集指定时间范围和研究范围的云量少于15%的Landsat卫星影像数据，并完成卫星影像拼接和裁剪处理；收集其他的多源辅助大数据包括：VIIRS夜间灯光数据、数字高程数据、最新年度的土地利用调查产品、城市兴趣点；

步骤3：样本点选取——通过实地采样、借助高分辨率影像或城市兴趣点选取最后一年各类地物的样本点，并以此年为基准年，切换上一年参考影像，逐点比对，对地物类别发生改变的样本点进行修改，从而得到上一年的样本点，以此类推，直到完成全部年份的样本点选取；

步骤4：多源大数据融合的初始分类——针对海岸带地区复杂多样的区域特征，采用随机森林算法对融合数据进行长时序影像分类，在此过程，通过不断调整样点分布、特征向量、纹理及最优窗口大小，持续优化初始分类结果，按照以下步骤进行处理：

4.1、以耕地、林地、草地、水体、不透水面、未利用地作为初步分类类型，使用Google Earth Engine内置的随机函数对训练样本点进行10-20次随机位置分布，取精度最高的结果，作为最终训练样本，所述不透水面包括城镇建设用地、工矿及基础设施建设用地和农村居民区；

4.2、计算所有影像中各像元的特征变量及纹理特征，所述特征变量包括：归一化植被指数（NDVI）、归一化建筑指数（NDBI）、归一化差水指数(NDWI)、修正归一化差水指数(MNDWI)，与VIIRS夜间灯光数据、数字高程数据（SRTM）；纹理特征包括差异值（diss）、惯性矩（inertia）、相关度（idm）；

4.3、将上述的特征变量及纹理特征进行组合，分别对纹理特征进行取值为1-9的不同窗口大小设置，利用随机森林算法进行优化，取精度最高的结果，作为最终纹理特征和窗口尺寸；

4.4、使用步骤4.1的最终训练样本、步骤4.2和步骤4.3的最终遥感指数特征、纹理特征和窗口大小，对分类器进行训练，使用训练完成的分类器对对应年份的影像进行分类，获得初始分类结果，并进行滤波处理；

4.5、对步骤4.4获得的历年分类结果进行合成，构造指定时间范围的时间序列数据集，使用最新年度的土地利用调查产品中的城乡范围栅格数据裁剪不透水面，将工矿及基础设施建设用地从不透水面中分离；

步骤5：基于扫描线种子填充算法的初步分离——利用扫描线种子填充算法将城镇建设用地和农村居民区从不透水面中分离，内陆水域和海水水域从水体中分离，具体按以下步骤进行实施：

5.1、取城市中心兴趣点为种子点，放入堆栈中，作为待填充对象；

5.2、填充当前种子点所处水平扫描线到边界之前的区段，填充完毕后删除堆栈中的该种子点，然后确定与此相邻的上下两条平行扫描线，并在这两者与边界间的区段中，取与种子点上下相邻的点并存入堆栈中，作为下次填充的种子点；

5.3、反复这一过程，直至清空堆栈内所有种子点，则区域填充完成；

5.4、保持种子点不变，遍历所有年度的初步分类结果，对其完成城镇建设用地、农村居民区从不透水面的分离；

5.5、以任一处海水位置为种子点，重复5.1-5.4，完成内陆水域和海水从水体中分离；

步骤6：基于几何特征的生产生态水体分离——将内陆水域分为天然水域和近海养殖用地及盐田，依据两者显著区别的几何特征进行分离，具体步骤如下：

6.1、对所有影像进行二值化处理，内陆水域区域赋值为1，其他为0，并进行对象分割和轮廓提取，获得若干内陆水域对象；

6.2、计算内陆水域对象P _i 的中心线长度L _i 、宽纵比R _i 和凸包Conv _i 三个形态特征，宽纵比R _i 为内陆水域对象P _i 的中心线长度L _i 与内陆水域对象P _i 的像素总数S _i 的比值，凸包Conv _i 为凸周长P _i(c) 与周长P _i(p) 的比值，用来评价对象凸性；

6.3、当内陆水域对象P _i 的中心线长度L _i 、宽纵比R _i 和凸包Conv _i 三个形态特征均达到预设阈值（L _i >30, R _i >0.2, Conv _i <0.45），则判定为天然水域，否则为近海养殖用地及盐田，从而实现天然水域和近海养殖用地及盐田的划分；

步骤7：基于时空变化逻辑规则的海岸带生产生活用地的分类修正——根据海岸带生产生活用地的时间序列变化特征，对不同地域功能用地之间的转变进行参数设定，并基于该参数设定进行双向时空一致性检测方法对耕地、工矿及基础设施建设用地和养殖水域进行修正，所述参数设定如下：一、不可逆型：包括耕地转化、围填海行为，滑动检测窗口W _d 尺寸为3，阈值prob为0.67；二、间断型：水产养殖行为，滑动检测窗口W _d 尺寸为5，阈值prob为0.6。

步骤1中，当多源大数据无法满足大范围海岸带区域某一类功能用地的准确识别时，可考虑对功能类型进行合并，以提高识别的总体精度。

进一步的，步骤7中，所述双向时空一致性检测方法如下：在每个时间序列中分别构造两个种子窗口W _s 和两个滑动检测窗口W _d ，两个种子窗口W _s 的初始位置位于时间序列的首末端，两个滑动检测窗口W _d 的初始位置位于紧邻于对应的种子窗口W _s ，按照以下步骤进行处理：

a）、将种子窗口和滑动检测窗口置于初始位置，并设置滑动检测窗口尺寸；

b）、将滑动窗口置于种子窗口紧邻的内侧；

c）、统计滑动窗口内的主导地类K _d ，窗口内某地类出现频率高于60%则称其为该窗口的主导地类K _d ；若窗口内所有地类出现频率均低于60%，则该窗口没有主导地类；

d）、判断种子窗口W _s 内地物类型K _s 与紧邻的滑动检测窗口W _d 内的主导地类K _d 是否一致，若一致，则将滑动检测窗口W _d 内首尾主导地类栅格之间所有地类均置为K _d ，将种子窗口位置移动到滑动窗口内最后出现主导地类K _d 的栅格位置，转至步骤b)直到两个滑动检测窗口相遇；若不一致，则将种子窗口向内移动一位，并转至步骤b) 直到两个滑动检测窗口相遇；

e）、当两个滑动检测窗口相遇时，两种子窗口之间的所有像素合并为一个大窗口，大窗口内与左侧种子窗口的像素地类一致的像素向左侧种子窗口方向移动，大窗口内与右侧种子窗口的像素地类一致的像素向右侧种子窗口方向移动，若还有其他地类像素，则将其他地类像素集中排列在中间位置，连续地类不一致处，即为地物转化突变年。

更进一步的，本发明基于多源大数据的海岸带地域功能长时序识别方法，还包括步骤8：精度评价——使用验证样本对海岸带地域功能用地分类结果进行精度评价，利用混淆矩阵计算生产者精度和用户精度对各类地物进行评价，并算出总体精度和kappa系数。

本发明基于多源大数据的海岸带地域功能长时序识别方法，还包括步骤9：海岸带人类生产生活功能的长时序制图与变化分析；针对粮食生产功能、近海养殖用地及盐田、围填海三类主要的海岸带人类生产生活功能，对海岸带区域地域功能用地长时序制图和变化分析：

9.1、粮食生产功能与近海养殖用地及盐田的变化节点检测与制图，根据海岸带长时序地域功能用地的分类结果和突变年的检测结果，提取粮食生产功能用地、近海养殖用地及盐田的变化时间和范围；

9.2、长时序围填海变化类型与阶段分析；根据海岸带长时序地域功能用地的分类结果，提取天然海水转变为不透水面、近海养殖的空间范围和变化时间，划分长时序围填海变化类型与阶段。

本发明方法采用基于多源大数据融合的长时序、大范围、高密度地域功能用地识别方法，突出了沿海人类生产生活活动的功能特点以及海陆统筹的功能细分特征，能够快速、准确地识别长时间尺度过程中近海岸地域功能用地转变与围填海建设变化特征。该方法适用于长时序、大范围、高密度的海岸带区域地域功能识别，尤其是适用于近海岸人类活动引发的地域功能用地转变与围填海建设的识别与变化检测。

附图说明

下面结合附图对本发明作进一步的说明。

图1是本发明实施例研究区概况图。

图2 本发明方法的总体流程图。

图3是海岸带区域长时序最终分类结果。

图4（a）是1987-2020年粮食生产功能用地变化结果图。

图4（b）是1987-2020年近海养殖用地及盐田功能用地变化结果图。

图5是1987-2020年围填海空间变化结果图。

图6是1987-2020年粮食生产功能、近海养殖用地及盐田和围填海变化阶段划分结果图。

具体实施方式

本实验采用的实验数据包括：1987—2020年间32年的Landsat系列影像，地表反射率产品，30m分辨率；VIIRS夜间灯光数据（Nighttime Day/Night Band CompositesVersion 1），500m分辨率；2020年数字高程数据（SRTM），30m分辨率；2020年度城市兴趣点（POI）数据；2020年度乡镇、行政村的行政区划数据。如图1所示。

研究区选取环渤海海岸带位于我国东北部，北起辽宁省盘山县，南至山东日照市（35°5′-41°27′N，116°42′-125°41′E），岸线总长约为6050千米，占全国岸线总长的三分之一。环渤海海岸带涉及山东、河北、辽宁三省以及天津，共有17个沿海城市。环渤海海岸带地处北温带，受季风环流的影响，大部分地区为暖温带半湿润季风气候，降水存在着明显的季节差异，主要集中在5-9月份，雨热同期。环渤海沿岸大部分地区都是地势低平的平原，如属于黄淮海平原的西南部，因此这里分布着大量耕地，也是中国重要粮仓所在。此外渤海作为我国唯一的内海，沿岸有着黄河、海河、辽河等众多入海河流，这些河流携带大量泥沙堆积在入海口，由此形成面积广阔的滩涂以及绵长的自然岸线，同时渤海海水浅、坡度缓、营养盐含量高，这些都为环渤海地区渔业发展提供了有利条件。2000到2020年环渤海海岸带常住人口数从7717万增长到9107万，增幅达到18%,以全国2%的土地容纳了近6.5%的人口。此外，环渤海海岸带在全国海洋经济中一直占据着重要地位，2020全国海洋生产总值为8.00万亿元，其中环渤海海岸带的海洋生产总值为2.34万亿，占比为29.2%。环渤海海岸带地区是中国北方地区经济发展的重要引擎。

本方法首先借助GEE（Google Earth Engine）平台完成影像数据分类，使用JavaScript语言。其他处理在本地用Python处理，版本3.6.

具体实施步骤如下：

本发明实施例基于多源大数据的海岸带地域功能长时序识别方法，包括以下步骤，如图2所示：

步骤1、海岸带地域功能用地分类。

兼顾海陆统筹特点，着重考虑到海岸带地域功能变化冲突以及制约海岸带可持续发展的主要因素，从人类生产、生活、生态的划分维度出发，划分城镇生产生活功能（城镇建设用地、工矿及基础设施建设用地）、农村生活功能（农村居民区）、食物生产功能（耕地、近海养殖用地及盐田）、生态功能（天然水域、海水水域、林地、草地、沿海滩涂、未利用地）。一共11种地类。

步骤2、多源大数据的收集。

基于Google Earth Engine云平台的多源大数据的收集和预处理。借助GoogleEarth Engine云平台收集从1987-2020年间的云量少于15%的Landsat数据，并完成影像拼接、裁剪等处理。收集辅助数据包括：中科院地理所2020年LULC产品（最新年度的土地利用调查产品）、城市兴趣点（Point of Interest，简称POI）。

步骤3、样本点选取。

确定地物类型和分类体系，通过实地采样、借助高分辨率影像或城市兴趣点选取2020年不同地物的样本点，并以此年为基准年，切换2019年参考影像（Google Earth中的影像），逐点比对，对地物类别发生改变的样本点进行修改，也可以视情况进行删除或新增样本点（通过目视识别进行样本点新增），从而得到2019年的样本点，以此类推，直到完成全部年份的样本点选取。

步骤4、多源大数据的初始分类和预处理。

以耕地、林地、草地、水体、不透水面（包括城镇建设用地、工矿及基础设施建设用地、农村居民区）、未利用地等作为初步分类类型，将80%的样本点作为训练样本点，20%作为验证样本点。采用随机森林算法对时间序列影像进行初始分类，在此过程中，通过不断调整样点分布、纹理及最优窗口大小等参数不断优化初始分类结果。具体步骤如下：

4.1、训练样本点分布优化。以耕地、林地、草地、水体、不透水面、未利用地作为初步分类类型，使用Google Earth Engine内置的随机函数对训练样本点进行10-20次随机位置分布，取精度最高的结果，作为最终训练样本；其中，不透水面包括城镇建设用地、工矿及基础设施建设用地和农村居民区。

4.2、计算所有影像中各像元的特征变量及纹理特征，所述特征变量包括：归一化植被指数（NDVI）、归一化建筑指数（NDBI）、归一化差水指数(NDWI)、修正归一化差水指数(MNDWI)，与VIIRS夜间灯光数据、数字高程数据（SRTM）；纹理特征包括差异值（diss）、惯性矩（inertia）、相关度（idm）。

4.3、将上述的特征变量及纹理特征进行组合，分别对其进行取值为1-9的不同窗口大小设置，利用随机森林算法进行优化，取精度最高的结果，作为最终纹理特征和窗口尺寸。最佳的窗口尺寸范围大致在3-5，每幅影像的最佳窗口尺寸不尽相同。

4.4、进行使用步骤4.1的最终训练样本、步骤4.2和步骤4.3的最终遥感指数特征、纹理特征和窗口大小，对分类器进行训练。使用训练完成的分类器对对应年份的影像进行分类，获得初始分类结果，从GEE云平台下载到本地，进行后续处理。由于初步分类结果会存在“椒盐”现象，使用ArcGIS软件中众数滤波工具对其进行迭代处理，去除碎小图斑，直至结果不再变化，保持相对稳定。

4.5、使用ArcGIS软件中波段合成工具对预处理后历年分类结果进行合成，1987为第一波段，2020为最后一个波段，构造1987-2020时间序列数据集。可以理解为将所有年份的影像分类结果叠加，构建分类结果的时间序列栅格数据集，时间序列即为取该栅格数据集中的某一行某一列像素的多年集合。使用中科院地理所2020年LULC产品的城乡范围栅格数据裁剪不透水面，得到2020年的城乡和建筑用地区域范围。

步骤5、基于扫描线种子填充算法的初步分离方法。

本实施例中，利用扫描线种子填充算法将城市和农村从不透水面中分离，内陆水域和海水从水域中分离，该算法的基本思想是非同类对象间不存在空间连通，实施的基础是城乡之间有空间间隔、内陆水域与海水间有人工坝体相隔。具体步骤如下：

5.1、取城市中心兴趣点（本例中以市政府作为中心兴趣点）为种子点，放入堆栈中，作为待填充对象。

5.2、对m行n列图像内的城镇用地范围、农村居民区范围取i个起始种子点seed _i ，分别赋值图像坐标(x _i ,y _i )，填入空栈stack，对各目标对象进行区域填充，遍历x _i 行所有列像素A(x _i ,y _i ±n)，若不是边界，则为真（同类对象），反之，该行填充完毕，填充完毕后删除堆栈中的该种子点，新起始点(x _i ±1,y _i )填入栈stack。新起始点(x _i ±1,y _i )为种子点(x _i ,y _i )上下相邻的点。

5.3、反复这一过程，直至清空堆栈内所有种子点，则区域填充完成。

5.4、保持种子点不变，遍历1987-2020年的初步分类结果，对其完成城镇建设用地、农村居民区从不透水面的分离。

5.5、以任一处海水位置为种子点，重复5.1-5.4，完成内陆水域和海水水域从水体中分离。

步骤6：基于几何特征的生产生态水体分离方法。将内陆水域分为天然水体和近海养殖用地及盐田，通过提取两者显著区别的集合特征进行分离，具体步骤：

6.1、对所有影像进行二值化、对象分割和轮廓提取等预处理。包括将步骤5结果中内陆水域二值化为1，再通过对象分割得到各目标内陆水域对象B _(1,2,...,i)，最后通过八邻域消除得到轮廓P _(1,2,...,i)。

6.2、计算内陆水域对象P _i （i为内陆水域对象P _i 的序号）的中心线长度L _i 、宽纵比R _i 和凸包Conv _i 三个形态特征，宽纵比R _i 为内陆水域对象P _i 的中心线长度L _i 与内陆水域对象P _i 的像素总数S _i 的比值，凸包Conv _i 为凸周长P _i(c) 与周长P _i(p) 的比值，用来评价对象凸性，可调用opencv3.0中convexityDefects函数计算。

6.3、当内陆水域对象P _i 的中心线长度L _i 、宽纵比R _i 和凸包Conv _i 三个形态特征均达到预设阈值（L _i >30, R _i >0.2, Conv _i <0.45），则判定为天然水域，否则为近海养殖用地及盐田；从而实现天然水域和近海养殖用地及盐田的划分。

步骤7：基于时空变化逻辑规则的海岸带生产生活用地的分类修正。利用双向时空变化逻辑一致性算法对耕地、工矿及基础设施建设用地和养殖水域进行修正。

7.1、规则设定：根据海岸带区域特征，对不同演变模式地类进行参数设定：

一、不可逆型：包括耕地转化、围填海行为，设置小滑动检测窗口W _d 和高阈值prob=0.67。

二、间断型：水产养殖行为，设置大滑动检测窗口W _d ，和低阈值prob=0.60。

7.2、设置滑动窗口和地物分布概率阈值

在每个时间序列中分别构造两个种子窗口Ws和两个滑动检测窗口W _d ，两个种子窗口Ws的初始位置位于时间序列的首末端，两个滑动检测窗口W _d 的初始位置位于紧邻于对应的种子窗口Ws，按照以下步骤进行处理：

a）、将种子窗口和滑动检测窗口置于初始位置，并设置滑动检测窗口尺寸；

b）、将滑动窗口置于种子窗口紧邻的内侧；

c）、计算滑动窗口内地物分布概率prob并同时确定主导地类K _d ，

d）、判断种子窗口Ws内地物类型Ks与紧邻的滑动检测窗口W _d 内的主导地类K _d 是否一致，若一致，则将滑动检测窗口Wd内首尾主导地类栅格之间所有地类均置为K _d ，将种子窗口位置移动到滑动窗口内最后出现主导地类K _d 的栅格位置，转至步骤b)直到两个滑动检测窗口相遇；若不一致，则将种子窗口向内移动一位，并转至步骤b)直到两个滑动检测窗口相遇；

7.3、检测生产生活用地的突变年。

当两个滑动检测窗口相遇时，两种子窗口之间的所有像素合并为一个大窗口，大窗口内与左侧种子窗口的像素地类一致的像素向左侧种子窗口方向移动，大窗口内与右侧种子窗口的像素地类一致的像素向右侧种子窗口方向移动，若还有其他地类像素，则将其他地类像素集中排列在中间位置，连续地类不一致处，即为地物转化突变年。

示例：

一、不可逆型：

测试序列：[1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 1, 1 ]

输出序列：[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

其中，1代表耕地/围填海地物，2代表其他地物，可见在测试序列中,出现逆向转换的情况，这与实际情况是不符的，而经过双向时空检测滤波处理后，消除了某些年份孤立存在的地物点，重新梳理了地物类型转换趋势。

二、间断型：

测试序列：[1, 1, 1, 1, 1, 2, 1, 2, 1, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2 ]

输出序列：[1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]

其中，2代表近海养殖用地及盐田，1代表其他地物，可见在测试序列中,在近海养殖用地及盐田为主导地物的阶段，中间会间断型出现其他地物类型，对地物转换趋势造成干扰。

经过双向时空检测滤波处理后，重新梳理了地物类型转换趋势，确定了地物转换的年份。

至此可得到该地区长时序最终分类结果，如图3所示。

步骤8、精度评价。使用验证样本对结果进行精度评价，利用混淆矩阵计算生产者精度和用户精度对各类地物进行评价，并算出总体精度和kappa系数。

选取完整率和正确率两个指标，其中完整率为某地类分类所得像元数与地类实际像元总数之比，对应的是漏分；正确率为正确分类像元数与地类实际像元总数之比，对应的是错分。

研究发现，初步分类结果的总体准确率为70.97%-78.74%，总体平均准确率为74.27%；双向检测结果的总准确率为76.70%-85.60%，总平均准确率为82.30%。很明显，双向检测方法的总体变化检测精度明显优于初始分类检测,总体平均精度提高了8.03%，主要关注功能用地中，内陆天然水体修正后的漏分率降低了7.23%，城镇建设用地、农村居民区、近海养殖用地及盐田两个结果均保持了较高的精度和稳定性，错分率差分别为：2.53%、0.72%、1.65%，漏分率差分别为：2.16%、0.07%、3.65%。另外，城镇建设用地、耕地、工矿及基础设施建设用地以及天然水域、近海养殖用地及盐田、海水水域的精度识别精度都较好，漏分率分别为：13.04%、8.17%、14.91%和0.07%、23.18%、5.46%，错分率分别为：11.38%、9.06%、17.76%和19.93%、19.27%、0.81%。

步骤9：耕地与近海养殖用地及盐田的变化节点检测与制图。通过对海岸带长时序地域功能用地的分类结果和突变年的检测结果分析，计算出该区域耕地从1987年的105463.5 km2锐减到2020年的73674.16 km2。而该时期近海养殖用地及盐田累计增加了7457.76 km2的，且存在两个明显的时间节点：1995和2007。如图4（a）、图4（b）、图6所示。

步骤10：长时序围填海用地变化类型与阶段分析。利用1987-2020年海岸带地域功能用地的分类结果，计算1987-2020减少的海水水域总面积，即围填海总面积，达到4344.73km2。根据对所有地物突变年的汇总，计算环渤海围填海变化存在着两个重要的时间节点：2003和2013年。据此将其划分成了三个阶段：第一阶段是1987-2003，在这一阶段围填海每年的增量基本维持在一个较低的水平，增幅都在150 km2以下。第二阶段是2004-2013，2004年环渤海的围填海数量出现激增，仅一年增幅就达到454.14 km2，随后几年的平均增幅也达到271.24 km2.第三阶段是2014-2020，经过上一阶段跳跃式的递增，这一时期的围填海活动逐渐得到控制。如图5、图6所示。

除上述实施例外，本发明还可以有其他实施方式。凡采用等同替换或等效变换形成的技术方案，均落在本发明要求的保护范围。

Claims (5)

1.一种基于多源大数据的海岸带地域功能长时序识别方法，包括以下步骤：

步骤1：海岸带地域功能用地分类——从海岸带人类生产、生活、生态为基础构建海岸带地域功能用地分类，包括：城镇生产生活区、农村生活区、粮食生产功能区、近海生产区域、陆地生态功能区、海洋生态功能区6个一级分类，其中，城镇生产生活区包括城镇建设用地、工矿及基础设施建设用地，农村生活区为农村居民区用地，粮食生产功能区为耕地，近海生产区域为近海养殖用地及盐田，陆地生态功能区包括陆地天然水域、林地、草地、未利用地，海洋生态功能区包括海水水域、沿海滩涂，一共11个二级分类；

步骤2：多源大数据的收集——借助Google Earth Engine云平台收集指定时间范围和研究范围的云量少于15%的Landsat卫星影像数据，并完成卫星影像拼接和裁剪处理；收集其他的多源辅助大数据包括：VIIRS夜间灯光数据、数字高程数据、最新年度的土地利用调查产品、城市兴趣点；

步骤3：样本点选取——通过实地采样、借助高分辨率影像或城市兴趣点选取最后一年各类地物的样本点，并以此年为基准年，切换上一年参考影像，逐点比对，对地物类别发生改变的样本点进行修改，从而得到上一年的样本点，以此类推，直到完成全部年份的样本点选取；

步骤4：多源大数据融合的初始分类——针对海岸带地区复杂多样的区域特征，采用随机森林算法对融合数据进行长时序影像分类，在此过程，通过不断调整样点分布、特征向量、纹理及最优窗口大小，持续优化初始分类结果，按照以下步骤进行处理：

4.1、以耕地、林地、草地、水体、不透水面、未利用地作为初步分类类型，使用GoogleEarth Engine内置的随机函数对训练样本点进行10-20次随机位置分布，取精度最高的结果，作为最终训练样本，所述不透水面包括城镇建设用地、工矿及基础设施建设用地和农村居民区；

4.2、计算所有影像中各像元的特征变量及纹理特征，所述特征变量包括：归一化植被指数（NDVI）、归一化建筑指数（NDBI）、归一化差水指数(NDWI)、修正归一化差水指数(MNDWI)，与VIIRS夜间灯光数据、数字高程数据（SRTM）；纹理特征包括差异值（diss）、惯性矩（inertia）、相关度（idm）；

4.3、将上述的特征变量及纹理特征进行组合，分别对纹理特征进行取值为1-9的不同窗口大小设置，利用随机森林算法进行优化，取精度最高的结果，作为最终纹理特征和窗口尺寸；

4.4、使用步骤4.1的最终训练样本、步骤4.2和步骤4.3的最终遥感指数特征、纹理特征和窗口大小，对分类器进行训练，使用训练完成的分类器对对应年份的影像进行分类，获得初始分类结果，并进行滤波处理；

4.5、对步骤4.4获得的历年分类结果进行合成，构造指定时间范围的时间序列数据集，使用最新年度的土地利用调查产品中的城乡范围栅格数据裁剪不透水面，将工矿及基础设施建设用地从不透水面中分离；

步骤5：基于扫描线种子填充算法的初步分离——利用扫描线种子填充算法将城镇建设用地和农村居民区从不透水面中分离，内陆水域和海水水域从水体中分离，具体按以下步骤进行实施：

5.1、取城市中心兴趣点为种子点，放入堆栈中，作为待填充对象；

5.2、填充当前种子点所处水平扫描线到边界之前的区段，填充完毕后删除堆栈中的该种子点，然后确定与此相邻的上下两条平行扫描线，并在这两者与边界间的区段中，取与种子点上下相邻的点并存入堆栈中，作为下次填充的种子点；

5.3、反复这一过程，直至清空堆栈内所有种子点，则区域填充完成；

5.4、保持种子点不变，遍历所有年度的初步分类结果，对其完成城镇建设用地、农村居民区从不透水面的分离；

5.5、以任一处海水位置为种子点，重复5.1-5.4，完成内陆水域和海水水域从水体中分离；

步骤6：基于几何特征的生产生态水体分离——将内陆水域分为天然水域和近海养殖用地及盐田，依据两者显著区别的几何特征进行分离，具体步骤如下：

6.1、对所有影像进行二值化处理，内陆水域区域赋值为1，其他为0，并进行对象分割和轮廓提取，获得若干内陆水域对象；

6.2、计算内陆水域对象P _i 的中心线长度L _i 、宽纵比R _i 和凸包Conv _i 三个形态特征，宽纵比R _i 为内陆水域对象P _i 的中心线长度L _i 与内陆水域对象P _i 的像素总数S _i 的比值，凸包Conv _i 为凸周长P _i(c) 与周长P _i(p) 的比值，用来评价对象凸性；

6.3、当内陆水域对象P _i 的中心线长度L _i 、宽纵比R _i 和凸包Conv _i 三个形态特征均达到预设阈值，则判定为天然水域，否则为近海养殖用地及盐田，从而实现天然水域和近海养殖用地及盐田的划分；

步骤7：基于时空变化逻辑规则的海岸带生产生活用地的分类修正——根据海岸带生产生活用地的时间序列变化特征，对不同地域功能用地之间的转变进行参数设定，并基于该参数设定进行双向时空一致性检测方法对耕地、工矿及基础设施建设用地和养殖水域进行修正，所述参数设定如下：一、不可逆型：包括耕地转化、围填海行为，滑动检测窗口W _d 尺寸为3，阈值prob为0.67；二、间断型：水产养殖行为，滑动检测窗口W _d 尺寸为5，阈值prob为0.6。

2.根据权利要求1所述的基于多源大数据的海岸带地域功能长时序识别方法，其特征在于：步骤7中，所述双向时空一致性检测方法如下：在每个时间序列中分别构造两个种子窗口W _s 和两个滑动检测窗口W _d ，两个种子窗口W _s 的初始位置位于时间序列的首末端，两个滑动检测窗口W _d 的初始位置位于紧邻于对应的种子窗口W _s ，按照以下步骤进行处理：

a）、将种子窗口和滑动检测窗口置于初始位置，并设置滑动检测窗口尺寸；

b）、将滑动窗口置于种子窗口紧邻的内侧；

c）、统计滑动窗口内的主导地类K _d ，窗口内某地类出现频率高于60%则称其为该窗口的主导地类K _d ；若窗口内所有地类出现频率均低于60%，则该窗口没有主导地类；

d）、判断种子窗口W _s 内地物类型K _s 与紧邻的滑动检测窗口W _d 内的主导地类K _d 是否一致，若一致，则将滑动检测窗口W _d 内首尾主导地类栅格之间所有地类均置为K _d ，将种子窗口位置移动到滑动窗口内最后出现主导地类K _d 的栅格位置，转至步骤b)直到两个滑动检测窗口相遇；若不一致，则将种子窗口向内移动一位，并转至步骤b) 直到两个滑动检测窗口相遇；

e）、当两个滑动检测窗口相遇时，两种子窗口之间的所有像素合并为一个大窗口，大窗口内与左侧种子窗口的像素地类一致的像素向左侧种子窗口方向移动，大窗口内与右侧种子窗口的像素地类一致的像素向右侧种子窗口方向移动，若还有其他地类像素，则将其他地类像素集中排列在中间位置，连续地类不一致处，即为地物转化突变年。

3.根据权利要求1所述的基于多源大数据的海岸带地域功能长时序识别方法，其特征在于：步骤6.3中，若以下三个不等式均成立，则内陆水域对象P _i 判定为天然水域，否则为近海养殖用地及盐田：L _i >30, R _i >0.2, Conv _i <0.45。

4.根据权利要求1所述的基于多源大数据的海岸带地域功能长时序识别方法，其特征在于还包括步骤8：精度评价——使用验证样本对海岸带地域功能用地分类结果进行精度评价，利用混淆矩阵计算生产者精度和用户精度对各类地物进行评价，并算出总体精度和kappa系数。

5.根据权利要求1所述的基于多源大数据的海岸带地域功能长时序识别方法，其特征在于还包括步骤9：海岸带人类生产生活功能的长时序制图与变化分析；针对粮食生产功能、近海养殖用地及盐田、围填海三类主要的海岸带人类生产生活功能，对海岸带区域地域功能用地长时序制图和变化分析：

9.1、粮食生产功能与近海养殖用地及盐田的变化节点检测与制图，根据海岸带长时序地域功能用地的分类结果和突变年的检测结果，提取粮食生产功能用地、近海养殖用地及盐田的变化时间和范围；

9.2、长时序围填海变化类型与阶段分析；根据海岸带长时序地域功能用地的分类结果，提取海水水域转变为其他功能类型的空间范围和变化时间，划分长时序围填海变化类型与阶段。

Images