本发明属于遥感影像智能识别，具体涉及一种基于遥感的河湖水域侵占时空回溯分析方法的设计。
背景技术：
1、河湖不仅具有供水、养殖、交通、航运、采砂、景观、娱乐等经济社会功能，还具备提供生物栖息地、维护生物多样性、消纳污染物、调节气候等生态系统服务功能，在经济社会发展和生态文明建设中具有十分重要的地位和作用。河湖水域范围是河湖监管的核心区域，一旦被侵占，将对河道行洪和生态安全造成严重影响。河湖水域侵占主要是人类社会不合理利用河道空间造成的，大致可以分为乱占、乱采、乱堆和乱建四类，如在水域范围内开展滩地利用，非法采砂或取土，倾倒、填埋或堆放固体废物，修建阻碍行洪的建筑物构筑物等。
2、河湖水域侵占范围大、类型多且较为隐蔽，人工巡检需要大量人力物力投入，且巡检的时效性低。卫星遥感技术具有宏观、准确和真实等优势，是开展大范围河湖水域调查的重要手段。随着国产高分遥感技术快速发展，相继发射了高分一号到高分八号等多颗高分卫星，形成了高分辨率多源卫星协同观测系统，为开展大尺度和精细的河湖水域侵占回溯分析提供了丰富且高质量的影像数据。
3、目前已有相关工作开展河湖水域侵占监管的工作，多以高分卫星遥感为数据，采用时序对比和人工目视解译方式，发现河湖水域侵占区域，这种方式存在不足：(1)人工目视解译方法效率较低，大范围河湖水域侵占监管耗时费力，河湖水域侵占遥感监测的时效性低。(2)缺乏河湖水域侵占时空回溯分析，难以展现河湖侵占全过程信息，取证信息匮乏，难以满足精准监管需求。
技术实现思路
1、本发明的目的是为了解决现有河湖水域侵占遥感监测方法时效性较低，且缺乏河湖水域侵占时空回溯分析，难以满足精准监管需求的问题，提出了一种基于遥感的河湖水域侵占时空回溯分析方法。
2、本发明的技术方案为：一种基于遥感的河湖水域侵占时空回溯分析方法，包括以下步骤：
3、s1、采集河湖侵占区域的中分辨率密集时序遥感影像数据，并对其进行几何和辐射校正，得到时序一致的中分辨率卫星遥感数据集。
4、s2、采集河湖侵占区域的高分遥感数据并进行预处理，得到时序一致的高分卫星遥感数据集。
5、s3、根据高分卫星遥感数据集对河湖水域进行智能提取，得到河湖水域变化范围。
6、s4、根据河湖水域变化范围和中分辨率卫星遥感数据集对河湖水域侵占进行时空回溯取证分析。
7、进一步地，步骤s2中河湖侵占区域的高分遥感数据选择标准为：
8、(1)选择晴空无云影像，河湖水域区域无云覆盖；
9、(2)影像地表层次清晰，无明显气溶胶覆盖；
10、(3)影像没有缺失扫描行、条带、噪声和异常像元；
11、(4)影像无冰雪覆盖。
12、进一步地，步骤s2包括以下分步骤：
13、s21、采集河湖侵占区域的30m数字高程模型数据和哨兵2号时序高分影像数据，并以哨兵2号10m空间分辨率的遥感数据作为参考影像采集控制点。
14、s22、将时序高分影像划分为10000×10000个像素的格网，在每个网格中，采用频率相位匹配法采集影像控制点，对误差高于2m的控制点进行剔除，在每个格网中保留20个高精度控制点。
15、s23、利用高精度控制点对时序高分影像中全色影像的有理函数系数进行优化，构建高精度几何校正模型对全色影像进行校正，并通过对时序高分影像中全色和多光谱影像进行配准，实现多光谱影像的几何精确校正。
16、s24、采用pansharp影像融合方法对时序高分影像中全色和多光谱影像进行融合，得到时序一致的高分卫星遥感数据集。
17、进一步地，步骤s3包括以下分步骤：
18、s31、对高分卫星遥感数据集中的遥感影像划分水体和非水体区域，并采用10000×10000个像素网格对遥感影像进行空间划分，针对每一个网格中的水体和非水体区域，分别随机生成500个样本点，并对每个样本点进行水体和非水体划分。
19、s32、构建高分遥感的特征向量。
20、s33、根据高分遥感的特征向量以及每个样本点的水体或非水体类型，采用改进极端梯度提升算法得到河湖水域变化范围。
21、进一步地，步骤s32中高分遥感的特征向量v′表示为：
22、v′＝v(b1,b2,b3,b4,bndvi,bndwi,bndbsi,bslope)
23、其中v(·)表示特征集，b1,b2,b3,b4分别表示高分遥感卫星影像中波段1、波段2、波段3和波段4影像特征波段的像素值，bndvi表示归一化植被指数特征波段的像素值，bndwi表示归一化水体指数特征波段的像素值，bndbsi表示归一化裸土指数特征波段的像素值，bslope表示坡度指数特征波段的像素值。
24、归一化植被指数ndvi的计算公式为：
25、
26、归一化水体指数ndvi的计算公式为：
27、
28、归一化裸土指数ndbsi的计算公式为：
29、
30、坡度指数slope的计算公式为：
31、
32、其中ρnir表示高分遥感多光谱波段4，ρred表示高分遥感多光谱波段3，ρgreen表示高分遥感多光谱波段2，ρblue表示高分遥感多光谱波段1，dz/dx表示高程z与x方向的比值，dx/dy表示高程z与y方向的比值。
33、进一步地，步骤s33包括以下分步骤：
34、s331、根据高分遥感的特征向量v′以及每个样本点的水体或非水体类型yi，构建有k个样本点的集合{t1(v′,yi),t2(v′,yi),…,tk(v′,yi)}，其中yi＝1表示样本点为水体，yi＝0表示样本点为非水体，tk(v′,yi)表示第k个样本点的特征和对应水体或非水体类型。
35、s332、通过在极端梯度提升算法中逐渐加入回归树，得到每个样本点的水体或非水体类型的预测值
36、
37、其中表示第t次生成的第i个水体或非水体类型的预测值，fk(·)表示第k棵回归树的得分，ft(·)表示第t棵回归树的得分。
38、s333、根据水体或非水体类型的预测值和真实值yi构建正则化目标函数l(ft)：
39、
40、其中l(·)表示微分损失函数，ω(ft)表示正则化项，计算公式为：
41、
42、其中γ和λ均为控制正则化程度的常数，t表示回归树叶子的个数，ωj表示第j个叶子节点的权重。
43、s334、计算正则化目标函数l(ft)的一阶统计量gi和二阶统计量hi：
44、
45、其中表示对的偏导。
46、s335、根据一阶统计量gi和二阶统计量hi对正则化目标函数l(ft)进行优化，得到优化目标函数l(t)：
47、
48、其中n表示水体或非水体类型数量。
49、s336、计算优化目标函数的最优解
50、
51、其中q表示样本，ij表示水体或非水体类型集合。
52、s337、将最优解作为回归树的叶子分数，当叶子分数小于预设阈值时，停止回归树增长叶子深度，得到遥感影像的水体和非水体分类结果。
53、s338、根据遥感影像的水体和非水体分类结果构建河湖水域范围集{a(t1,ya,b),a(t2,ya,b),…,a(tm,ya,b)}，其中a(tm,ya,b)表示tm时刻遥感影像第a行第b列的河湖水域类型，a(tm,ya,b)＝1表示河湖水域类型为水体，a(tm,ya,b)＝0表示河湖水域类型为非水体。
54、s339、根据河湖水域范围集获取河湖水域变化范围a(a,b)：
55、
56、其中a(a,b)＝1表示河湖水域侵占区域，a(a,b)＝0表示河湖水域非侵占区域。
57、进一步地，步骤s4包括以下分步骤：
58、s41、根据河湖水域管理范围和河湖水域侵占区域，获取河湖水域回溯取证指标。
59、s42、根据河湖水域变化范围和中分辨率卫星遥感数据集对河湖侵占关键时间节点进行回溯取证，得到河湖水域持续侵占时间。
60、s43、根据河湖水域回溯取证指标和河湖水域持续侵占时间，形成河湖水域侵占时空回溯取证信息，对河湖水域侵占进行时空回溯取证分析。
61、进一步地，步骤s41中获取的河湖水域回溯取证指标包括河湖侵占变化面积自然岸线或堤防破坏长度建筑或养殖变化面积河湖侵占面积年际变化率自然岸线或堤防年际破坏长度和建筑或养殖年际变化面积
62、河湖侵占变化面积的计算公式为：
63、
64、其中st2表示t2时间点河湖水域面积，st1表示t1时间点河湖水域面积。
65、自然岸线或堤防破坏长度的计算公式为：
66、
67、其中lt2表示t2时间点自然岸线或堤防破坏长度，lt1表示t1时间点自然岸线或堤防破坏长度。
68、建筑或养殖变化面积的计算公式为：
69、
70、其中dt2表示t2时间点建筑或养殖面积，dt1表示t1时间点建筑或养殖面积。
71、河湖侵占面积年际变化率的计算公式为：
72、
73、其中sr2表示r2年份河湖水域面积，sr1表示r1年份河湖水域面积，r1表示河湖侵占开始年份，r2表示河湖侵占结束年份。
74、自然岸线或堤防年际破坏长度的计算公式为：
75、
76、其中lr2表示r2年份自然岸线或堤防破坏长度，lr1表示r1年份自然岸线或堤防破坏长度。
77、建筑或养殖年际变化面积的计算公式为：
78、
79、其中dr2表示r2年份建筑或养殖面积，dr1表示r1年份建筑或养殖面积。
80、进一步地，步骤s42包括以下分步骤：
81、s421、根据河湖水域变化范围a(i,j)，确定河湖侵占典型位置，针对河湖侵占典型位置的每一幅中分辨率卫星遥感影像计算归一化水体增强指数endwi：
82、
83、其中green表示绿光波段的反射率，mir表示中红外波段的反射率。
84、s422、以开始监测河湖侵占典型位置的第一幅中分辨率卫星遥感影像的日期为原点，后续每一幅影像的日期距离第一幅影像的天数为x轴，每一幅影像河湖侵占典型位置的endwi值为y轴绘制散点图，并构建关于x和y之间的拟合函数拟合散点图：
85、y＝ax3+bx2+cx+d
86、其中a、b、c、d均为通过多项式拟合得到的系数及常数。
87、s423、根据拟合函数确定当y＝0时x的值，并根据y＝0时x的值选取两个相邻影像日期作为河湖水域侵占开始时间t1和河湖水域侵占结束时间t2。
88、s424、根据河湖水域侵占开始时间t1和河湖水域侵占结束时间t2计算河湖水域持续侵占时间：
89、δt＝t2-t1
90、其中δt表示河湖水域持续侵占时间。
91、本发明的有益效果是：
92、(1)本发明利用多源卫星遥感数据，构建了河湖水域侵占遥感时空回溯分析方法，突破了跨区域、长时序河湖水域侵占回溯取证的技术瓶颈，具备较高的时效性。
93、(2)本发明利用时序高分遥感数据和中分辨率遥感数据，结合改进极端梯度提升算法，进行长时序河湖水域智能化提取，具有较高精度且普适性较强。
94、(3)本发明构建了河湖侵占的遥感回溯取证指标体系，定量分析河湖侵占的时空过程，为河湖侵占监管提供客观依据。
95、(4)本发明结合密集时序的landsat系列和哨兵2号中分辨率密集时序遥感影像数据，实现了河湖侵占关键时间节点判定，为河湖侵占时空取证提供了明确的依据。
1.一种基于遥感的河湖水域侵占时空回溯分析方法，其特征在于，包括以下步骤：
2.根据权利要求1所述的基于遥感的河湖水域侵占时空回溯分析方法，其特征在于，所述步骤s2中河湖侵占区域的高分遥感数据选择标准为：
3.根据权利要求1所述的基于遥感的河湖水域侵占时空回溯分析方法，其特征在于，所述步骤s2包括以下分步骤：
4.根据权利要求1所述的基于遥感的河湖水域侵占时空回溯分析方法，其特征在于，所述步骤s3包括以下分步骤：
5.根据权利要求4所述的基于遥感的河湖水域侵占时空回溯分析方法，其特征在于，所述步骤s32中高分遥感的特征向量v′表示为：
6.根据权利要求4所述的基于遥感的河湖水域侵占时空回溯分析方法，其特征在于，所述步骤s33包括以下分步骤：
7.根据权利要求6所述的基于遥感的河湖水域侵占时空回溯分析方法，其特征在于，所述步骤s4包括以下分步骤：
8.根据权利要求7所述的基于遥感的河湖水域侵占时空回溯分析方法，其特征在于，所述步骤s41中获取的河湖水域回溯取证指标包括河湖侵占变化面积自然岸线或堤防破坏长度vl、建筑或养殖变化面积河湖侵占面积年际变化率自然岸线或堤防年际破坏长度和建筑或养殖年际变化面积
9.根据权利要求7所述的基于遥感的河湖水域侵占时空回溯分析方法，其特征在于，所述步骤s42包括以下分步骤：
