基于多源影像的内河滨岸带植被生长适宜区确定方法和装置

本发明属于生态河道治理，具体涉及一种基于多源影像的内河滨岸带植被生长适宜区确定方法和装置。

背景技术：

1、滨岸带绿化是生态河道治理的重要内容，如城市江段两岸的高滩岸线常被开发为覆盖大量植被的休闲景观区，河道和航道整治工程中也多采用生态护坡、护岸技术替代传统护岸模式。这些工程的设计和建设中，滨岸带植被配置需要因地制宜、充分利用现状河道的地形、水文等条件再辅以适当的人工修复，构建出适宜植物生长的稳定生境条件。当前许多生态护岸工程实施后出现植被大面积死亡的现象，一些工程甚至出现“建设施工-效果反馈-修改补救”的模式，造成资源浪费，导致这些现象的原因在于规划设计阶段未能明确滨岸带哪些区域适宜植被生长，以及植被生长的区域需要满足什么水文条件。

2、不同于人工渠道水位较为固定、海滨潮间带滩涂位置较为固定的简单情形，天然河道的岸滩地貌形态复杂且受年内、年际周期性大幅涨落的水位变化制约，其滨岸带植被分布既呈现出沿河道横剖面梯级分布的总体特征，又具有强烈的空间和时间随机性，很难通过量化方法确定出植被适宜区域。以往一些工程建设中，采用实地调查和水动力数值模拟相结合的方法来判定植被适宜区，但该过程需要长期调查资料，建模难度大、成本高，一般情况下不具实施条件。

技术实现思路

1、本发明是为了解决上述问题而进行的，目的在于提供基于多源影像的内河滨岸带植被生长适宜区确定方法和装置，对实测地形数据要求量少，数据互补性强，操作过程简便且定量化，充分考虑了天然河流水位变幅大且随机性强的水文特点，以及岸滩高低起伏且水陆交错的复杂形态特点，淹没频率的计算更加准确，能够为生态河道治理提供更加科学、高效、可靠的实施途径。

2、本发明为了实现上述目的，采用了以下方案：

3、<方法>

4、本发明提供基于多源影像的内河滨岸带植被生长适宜区确定方法，其特征在于，包括以下步骤：

5、步骤1.基于研究河段的影像数据，并逐幅实施水陆区域划分、构建二值化淹没影像数据集；

6、步骤2.基于像元尺度计算滨岸带岸滩长历时平均淹没频率；具体包括如下子步骤：

7、步骤2.1：根据目标河段是否具有长系列水位观测资料来选取滨岸带淹没频率计算方法，若具有水位观测资料则进入步骤2.2，否则进入步骤2.10；

8、步骤2.2：以日期为索引，匹配每幅二值化淹没影像{b11,.......,bij}所对应日期的日均水位，记为{z11,.......,zij}，其中zij表示第i年第j幅影像对应日期的水位；

9、步骤2.3：根据区域淹没频率统计需求对t年内多年日均水位从小到大排序后按算数间隔等距划分m个水位等级区间，25≤m≤35；

10、步骤2.4：统计水位{z11,.......,zij}落在m个水位区间内的频次f，保留频次f≥1的水位区间，对频次f＝0的水位区间与相邻区间进行合并，合并后水位区间级数记为m'，20≤m'；

11、步骤2.5：根据逐年日均水位数据，统计第i年(i＝1,2,……,t)内全年每日水位落在m'个水位区间内的天数nik(k＝1,2,……,m')，天数nik为第i年第k个水位区间的时段长度；

12、步骤2.6：根据时序影像对应水位数据{z11,.......,zij}，统计第i年内第j幅影像对应水位zij落在m'个水位区间内的影像个数c；

13、步骤2.7：按步骤2.5将第i年的各水位区间一一匹配出其代表影像，编码记为bik，bik表示第i年第k个水位区间的二值化淹没影像；

14、步骤2.8：对所有二值化影像bik，按照公式逐像元进行运算，得到第i年各像元的淹没频率，所有像元的淹没频率数值构成第i年目标河段区域的淹没频率空间分布的栅格数据ffi；

15、步骤2.9：在逐年淹没频率基础上，逐像元按照公式进行运算，得到t年内的多年平均淹没频率空间分布栅格数据ff；

16、步骤2.10：若不具备水位观测资料，则使用二值化影像集{b11,.......,bij}的时间坐标作为计算淹没历时的依据，将第i年第j-1幅与第j幅影像对应日期的时间中点记为ts，第j幅与第j+1幅影像对应日期的时间中点记为te，取nij＝te-ts作为第j幅影像的代表时段，在此基础上按公式进行运算得到各像元的淹没频率，它们共同形成第i年的淹没频率空间分布栅格数据ffi'，j＝1,2,……,ni，ni为第i年收集到的影像数量；多年平均淹没频率空间分布栅格数据的计算方法同步骤2.8；

17、步骤3.对影像逐幅提取滨岸带植被覆盖区并构建二值化覆被影像集合；

18、步骤4.基于像元尺度计算长历时内滨岸带岸滩植被出现频率；

19、步骤5.将目标河段区域内滨岸带植进行适宜区划分；根据目标河段区域多年平均淹没频率ff和植被出现频率vf的栅格数据得到ff～vf的散点图和核密度云图，进而划分出滨岸带植被生长的适宜区，并推求出滨岸带植被适宜生长的临界淹没频率阈值。

20、优选地，本发明提供的基于多源影像的内河滨岸带植被生长适宜区确定方法，在步骤1中，对二值化淹没影像数据集{mb11,.......,mbij}和{sb11,.......,sbij}中的10m精度的二值化影像，采用重采样技术将像素尺度转换为30m，随后将{mb11,.......,mbij}和{sb11,.......,sbij}合并成30m精度二值化淹没影像数据集{b11,.......,bij}，其中bij表示第i年的第j幅影像。考虑了单一影像数据集在对区域水体植被监测时的时空局限性，引入合成孔径雷达等影像数据集，弥补了多光谱影像易受云雾干扰导致数据量不足等问题，使用多源遥感影像数据集通过差异化水体指数选取和像素尺度的转化等技术手段实现了对数据的扩展，为后续淹没频率和植被出现频率的精细化计算提供了充足的数据量。

21、优选地，本发明提供的基于多源影像的内河滨岸带植被生长适宜区确定方法，步骤1包括如下子步骤：

22、步骤1.1：从数据中心(例如，google earth engine)分别获取历时t年的目标河段影像数据，t≥10，其中以经过几何校正和辐射定标以及大气校正等预处理手段的多光谱影像数据为主，如landsat sr二级数据产品(30m精度)、sentinel-2 2a级数据(10m精度)等，以经过辐射定位、校正、去噪后的合成孔径雷达数据为辅，如sentinel-1sar grd地距影像数据(10m精度)。并筛选出质量、数量满足条件的多光谱影像数据：(1)遥感影像中，目标河段区域应满足不超过5％的云量遮挡，无缺损及条带；(2)每年内应满足每季均有可用影像；

23、步骤1.2：对landsat sr二级数据使用qa波段对云层进行掩膜处理，对sentinel-22a级数据使用slc波段进行掩膜处理；对sar grd数据进行平滑处理和降噪；去云后对目标河段进行裁剪，组成目标河段时序多光谱影像数据集{ms11,.......,msij}，时序合成孔径雷达数据集{sar11,.......，sarij}，其中i表示第i年，j表示第i年的第j幅的影像，i＝1,2,……,t；j＝1,2,……,ni，ni为第i年收集到的影像数量；

24、步骤1.3：选取{ms11,.......,msij}中每幅影像msij的ρgreen绿光波段和ρswir短波红外波段的反射率数据按照公式mndwi＝(ρgreen-ρswir)/(ρgreen+ρswir)进行波段运算，得到改进归一化差异水体指数mndwi影像，从而将{ms11,.......,msij}数据集转化为改进归一化差异水体指数mndwi影像数据集{m11,.......,mij}，其中mij表示第i年的第j幅mndwi影像；针对合成孔径雷达数据集{sar11,.......,sarij}，选取每幅影像sarij中经过投影校正的level 1地面距离探测产品中的双极化数据(vv波段和vh波段)，按照水体提取指数公式sdwi＝ln(10×vv×vh)－8对波段数据进行计算，得到水体提取指数sdwi，形成sdwi影像数据集{s11,.......,sij},其中sij表示第i年的第j幅sdwi影像；

25、步骤1.4：在步骤1.3的基础上，针对影像数据集{m11,.......,mij}和{s11,.......,sij}，输出目标河段区域每幅影像的灰度直方图，横坐标为灰度值(mndwi或sdwi值)，纵坐标为该灰度值下的像元频数h，该灰度直方图应呈双峰状；

26、步骤1.5：基于每幅影像的灰度直方图双峰状分布特点，采用otsu最大类间方差法逐幅确定水体和非水体的分割阈值；，其实现过程为：(1)对每幅灰度直方图进行分类试算，先假定一个分类的灰度阈值，再按照公式σ2＝ω0×ω1×(μ0-μ1)2计算该灰度阈值下的类间方差σ2，将图幅中小于该灰度值的像素占比记为ω0，其平均灰度为μ0；大于该灰度值占比记为ω1，其平均灰度为μ1；(2)不断重新假设分类灰度阈值，重复步骤(1)，直至找出类间方差σ2取最大值时的灰度阈值threshold，即为水体和非水体的分割阈值。

27、步骤1.6：逐幅计算出水体和非水体分割阈值threshold后，按每幅影像的分割阈值对数据集{m11,.......,mij}及{s11,.......,sij}逐幅二值化，水体区域像素赋值为1，非水体像素区域赋值为0，得到二值化淹没数据集{mb11,.......,mbij}和{sb11,.......,sbij}，其中mbij表示基于landsat sr和sentinel-2 2a数据的第i年第j幅二值化淹没影像；sbij则表示基于合成孔径雷达数据集sentinel-1sar数据的第i年第j幅二值化淹没影像；

28、步骤1.7：对二值化淹没影像数据集{mb11,.......,mbij}和{sb11,.......,sbij}中基于sentinel-2 2a/sentinel-1sar卫星生成的10m精度的二值化影像，采用重采样技术将像素尺度转换为30m，随后将{mb11,.......,mbij}和{sb11,.......,sbij}合并成30m精度二值化淹没影像数据集{b11,.......,bij}，其中bij表示第i年的第j幅影像。

29、优选地，本发明提供的基于多源影像的内河滨岸带植被生长适宜区确定方法，在步骤2.6中，若第k区间内影像个数c≥2，则取对应日均水位最接近该区间中值的影像作为该区间代表影像，例如水位区间划分为(32m～33m)，则选取对应水位最接近32.5m影像作为代表影像；若某年第k区间内影像个数c＝0，则遍历t年内影像对应水位zij，取影像对应水位zij最接近该水位区间中值且时间最接近该年的影像作为该年第k区间代表影像的替代。

30、优选地，本发明提供的基于多源影像的内河滨岸带植被生长适宜区确定方法，步骤3包括如下子步骤：

31、步骤3.1：基于步骤1.1和1.2获取的目标河段时序多光谱影像数据集{ms11,.......,msij}，取其中的每幅影像msij的ρred红光波段和ρnir红外波段的反射率数据按照公式ndvi＝(ρnir-ρred)/(ρnir+ρred)进行波段运算，得到归一化差异植被指数ndvi影像，并形成目标河段区域长时序植被指数ndvi影像数据集{v11,.......,vij}，其中vij表示第i年的第j幅ndvi影像；

32、步骤3.2：针对ndvi影像数据集{v11,.......,vij}，使用波段计算器逐幅掩膜水体区域像元mask(ndvi≤0)，掩膜后生成包含光滩+植被的rndvi影像，形成数据集{rn11,.......,rnij}；

33、步骤3.3：{rn11,.......,rnij}中的rnij影像灰度直方图呈现双峰状，使用步骤1.5中otsu最大类间方差法逐幅计算出分割光滩和植被的阈值threshold'；

34、步骤3.4：按每幅影像的分割阈值threshold'对数据集{rn11,.......,rnij}逐幅二值化，植被区域像素赋值为1，光滩区域像素赋值为0，得到二值化覆被数据集{mvb11,.......,mvbij}，其中mvbij表示第i年第j幅二值化覆被影像；通过简化天然河道地物分类，使用一种较为简单的植被覆盖区提取方法，实现了对区域植被覆盖范围的快速识别和准确提取；

35、步骤3.5：对来源于sentinel-2 2a10m精度的二值化覆被影像执行像素尺度转换，采用重采样技术将其像素由10m转换为30m，形成30m精度二值化覆被数据集{p11,.......,pij}。

36、优选地，本发明提供的基于多源影像的内河滨岸带植被生长适宜区确定方法，步骤4包括如下子步骤：

37、步骤4.1：使用二值化覆被影像集{p11,.......,pij}的时间坐标作为计算植被出现历时的依据，将第i年第j-1幅与第j幅影像对应日期的时间中点记为ts，第j幅与第j+1幅影像对应日期的时间中点记为te，取nij＝te-ts作为第j幅影像的代表时段；在此基础上按公式进行运算计算各像元的植被出现频率，得到第i年的植被出现频率空间分布栅格数据vfi；

38、步骤4.2：在逐年植被出现频率的基础上，逐像元按照公式进行运算，得到t年内的多年平均植被出现频率空间分布的栅格数据。

39、以像元的植被出现频率为植被生长状态和生长范围的判别指标充分考虑了滨岸带植被的多年变化过程，植被出现频率反映了不仅反映出受周期性水淹影响的滨岸带植被生长的区域信息，其频率高低反映植被物种种类，植被出现频率较低的区域为滨岸带新生先锋植物，而植被出现频率较高的区域则代表了已经稳定生长的多年生植物。

40、优选地，本发明提供的基于多源影像的内河滨岸带植被生长适宜区确定方法，步骤5包括如下子步骤：

41、步骤5.1：根据目标河段区域多年平均淹没频率ff和植被出现频率vf的栅格数据，逐像元输出淹没频率和植被出现频率，将其存储到两个数组中；基于两个数组中所有像元的ff～vf关系绘制散点图，其中ff为横坐标，vf为纵坐标；

42、步骤5.2：对ff～vf散点图进行聚类分析，将生成的散点图转换为核密度云图；具体为：利用matlab的ks density函数对散点图数据进行核密度估计，使用高斯核函数对数据进行平滑，并计算估计密度，对其估计概率密度归一化后进行可视化，并绘制ff～vf核密度估计云图。云图中的数值代表了某一淹没频率和某一植被出现频率组合条件下的像元占总像元的比例；

43、步骤5.3：为ff～vf核密度云图添加等值线，根据2％～5％等值线的封闭情况形成三个分区，三个分区两两之间的界限在纵轴上的投影对应了两个植被出现频率阈值cvf1、cvf2，据此将vf栅格数据划分为三类，对应植被出现频率为高、中、低三个等级，其中，vf≥cvf2的像元对应区域为稳定植被区，cvf2＞vf＞cvf1的区域为过渡植被区，cvf1≥vf为植被不适宜生长区；结合地理空间信息对这些植被分区进行可视化，制作出滨岸带植被适宜度分区图；

44、步骤5.4：针对ff～vf散点图，将ff按照步长值1％划分为100个等距区间，各区间的ff均值为ffa，各区间内vf的均值为vfa，通过origin软件根据ffa和vfa拟合logistic s型曲线，其一般形式为y＝a+(b/(1+(x/c)d))，其中a、b、c、d均为常数，a是纵向偏移量，b是曲线的上限，c是中心位置，d是曲线的斜率；

45、步骤5.5：将步骤5.2中界定的两个植被出现频率cvf1、cvf2对应到logistic s型曲线上，得出植被适宜区、过渡区、不适宜区之间的临界淹没频率cff1、cff2，工程实践中可根据这些临界值对滨岸带地形实施适当改造，通过地形高程升降来使某位置的淹没频率大于或小于以上临界值，从而满足生态治理需求。

46、<装置>

47、进一步，本发明还提供能够自动实现以上<方法>的基于多源影像的内河滨岸带植被生长适宜区确定装置，其特征在于，包括：

48、数据集构建部，基于研究河段的影像，逐幅实施水陆区域划分、构建二值化淹没影像数据集；

49、平均淹没频率计算部，采用如下步骤2.1～2.10基于像元尺度计算滨岸带岸滩长历时平均淹没频率；

50、步骤2.1：根据目标河段是否具有长系列水位观测资料来选取滨岸带淹没频率计算方法，若具有水位观测资料则进入步骤2.2，否则进入步骤2.10；

51、步骤2.2：以日期为索引，匹配每幅二值化淹没影像{b11,.......,bij}所对应日期的日均水位，记为{z11,.......,zij}，其中zij表示第i年第j幅影像对应日期的水位；

52、步骤2.3：根据区域淹没频率统计需求对t年内多年日均水位从小到大排序后按算数间隔等距划分m个水位等级区间，25≤m≤35；

53、步骤2.4：统计水位{z11,.......,zij}落在m个水位区间内的频次f，保留频次f≥1的水位区间，对频次f＝0的水位区间与相邻区间进行合并，合并后水位区间级数记为m'，20≤m’；

54、步骤2.5：根据逐年日均水位数据，统计第i年(i＝1,2,……,t)内全年每日水位落在m'个水位区间内的天数nik(k＝1,2,……,m')，天数nik为第i年第k个水位区间的时段长度；

55、步骤2.6：根据时序影像对应水位数据{z11,.......,zij}，统计第i年内第j幅影像对应水位zij落在m'个水位区间内的影像个数c；

56、步骤2.7：按步骤2.5将第i年的各水位区间一一匹配出其代表影像，编码记为bik，bik表示第i年第k个水位区间的二值化淹没影像；

57、步骤2.8：对所有二值化影像bik，按照公式逐像元进行运算，得到第i年各像元的淹没频率，所有像元的淹没频率数值构成第i年目标河段区域的淹没频率空间分布的栅格数据ffi；

58、步骤2.9：在逐年淹没频率基础上，逐像元按照公式进行运算，得到t年内的多年平均淹没频率空间分布栅格数据ff；

59、步骤2.10：若不具备水位观测资料，则使用二值化影像集{b11,.......,bij}的时间坐标作为计算淹没历时的依据，将第i年第j-1幅与第j幅影像对应日期的时间中点记为ts，第j幅与第j+1幅影像对应日期的时间中点记为te，取nij＝te-ts作为第j幅影像的代表时段，在此基础上按公式进行运算得到各像元的淹没频率，它们共同形成第i年的淹没频率空间分布栅格数据ffi'，j＝1,2,......,ni，ni为第i年收集到的影像数量；多年平均淹没频率空间分布栅格数据的计算方法同步骤2.8；

60、影像集合构建部，对影像逐幅提取滨岸带植被覆盖区并构建二值化覆被影像集合；

61、出现频率计算部，基于像元尺度计算长历时内滨岸带岸滩植被出现频率；

62、划分确定部，将目标河段区域内滨岸带植进行适宜区划分并确定临界淹没频率阈值；根据目标河段区域多年平均淹没频率ff和植被出现频率vf的栅格数据得到ff～vf的散点图和核密度云图，进而划分出滨岸带植被生长的适宜区，并推求出滨岸带植被适宜生长的临界淹没频率阈值；

63、控制部，与数据集构建部、平均淹没频率计算部、影像集合构建部、出现频率计算部、划分确定部均通信相连，控制它们的运行。

64、优选地，本发明提供的基于多源影像的内河滨岸带植被生长适宜区确定装置，还可以包括：输入显示部，与控制部通信相连，让操作员输入操作指令，并根据控制指令显示相应的信息。

65、优选地，本发明提供的基于多源影像的内河滨岸带植被生长适宜区确定装置中，输入显示部根据控制指令能够对数据集构建部、平均淹没频率计算部、影像集合构建部、出现频率计算部、划分确定部的输入、输出数据和处理过程以数据表或图像形式进行静态或动态显示。

66、发明的作用与效果

67、(1)基于开放的多源遥感影像，数据之间互补性强，时间和空间覆盖程度高，影像处理基于遥感云平台，操作过程简便且定量化，克服了实地查勘工作量大、数据处理经验性强的不足。

68、(2)计算淹没频率对实测地形数据要求量少，思路清晰，操作简单，充分利用逐年的实测影像和水位资料，考虑了冲淤平衡态下年际间河床演变引起的河道地形微变化，避免使用单一地形对滨岸带淹没频率计算造成的误差，大大提高了研究时长内淹没频率的计算精度，能够较为准确的反映区域的真实淹没频率情况，便于实践中应用。

69、(3)在遥感影像的处理和统计过程中，充分考虑了天然河流水位变幅大且随机性强的水文特点，以及岸滩高低起伏且水陆交错的复杂形态特点，计算更加准确而可靠，并具有较强的普适性。

70、(4)不仅能够得到河滨带内稳定植被区、过渡植被区、不适宜生长区的空间范围，还可得到各区对应的淹没频率阈值，从而可为生态河道治理规划设计中的植被区配置、断面形态改造等提供科学依据，充分契合了工程实践的需求。