技术领域
[0001] 本发明属于生产建设管理技术领域，涉及生产建设项目建设全周期变化检测方法。
相关背景技术
[0002] 生产建设项目在施工建设前需要取得立项、用地许可、水行政许可、施工许可等审批手续，在建设过程中需要接受环境、水利等多个部门的监督管理并根据意见整改，施工结束后需完成验收和备案等工作。而在实际工作中生产建设项目常存在未取得完善审批手续的情况下开工建设，或当前建设阶段未及时履行相关手续和整改，或建设状态长期没变化建设单位成为僵尸企业等情况，这些情况直接影响资源禀赋的高效利用、监督管理工作的全面落实以及生态环境的安全保障，制约着社会经济发展和美好生活建设。根据统计，全国每年有超过100万个项目开工或正在建设，管理部门现常用的方式是通过实地调查取证查明建设状态和履行手续情况，难以全覆盖、全周期监管生产建设项目，难以及时发现僵尸企业予以指导，影响社会经济生态的高质量发展。因此，急需研发生产建设项目建设全周期变化检测方法及介质，自动识别项目状态、动态跟踪建设变化，及时发现并预警僵尸项目，实现生产建设项目建设过程全覆盖、全过程、全链路的信息化管理，为管理部门提供技术支撑。
具体实施方式
[0069] 下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。
[0070] 实施例1
[0071] 本发明公开了一种生产建设项目建设全周期变化检测方法，该方法能够根据多期遥感影像数据，自动识别生产建设项目建设状态、变化情况，并能够自动更新项目建设状态和建设日期，实现生产建设项目建设周期全过程自动化智慧化管理，如图1‑3所示，该方法具体步骤包括：
[0072] S1：获取监管区域卫星遥感影像及数字高程模型DEM，以及生产建设项目防治责任范围矢量文件；上述监测数据作为本发明检测方法的基础数据；
[0073] S2：确定项目建设阶段变化监测指标，从遥感影像和DEM两个方面建立生产建设项目阶段变化监测方法判断项目建设阶段是否发生变化；
[0074] S3：确定项目建设阶段判别指标，建立项目建设阶段判定样本库，从而构建项目建设阶段判别模型以判定当前项目建设阶段；
[0075] S4：将样本库输入项目建设阶段判别模型进行训练，直至建模精度不低于95％、验证精度不低于90％认为建模成功，然后将步骤1采集的数据利用生产建设项目阶段变化监测方法和模型进行周期变化检测。
[0076] 实施例2
[0077] 与实施例1的区别在于，
[0078] S1中卫星遥感影像为红绿蓝RGB三波段栅格数据(GeoTIFF等格式)，分辨率优于2m；也可以采用低空无人机遥感影像。DEM生产时间应与遥感影像总体保持一致，分辨率为优于10m。生产建设项目防治责任范围矢量文件属性信息应包含项目名称、责任面积、建设状态、当前建设状态时间等必要信息。
[0079] S2、确定项目建设阶段变化监测指标，从遥感影像和DEM两个方面建立生产建设项目阶段变化监测方法；
[0080] S21、项目建设阶段变化监测方法
[0081] 本步骤中以某研究区A的某生产建设项目P第i期和第i+1期变化检测计算为例说明项目建设阶段变化监测方法。
[0082] 已知项目P从第i期到第i+1期项目建设内容和阶段均发生了显著变化，现在通过这两期的遥感影像和DEM数据计算其变化监测指标CHi_i+1。
[0083] 收集到项目P第i期遥感影像数据Image_i和DEM数据Dem_i；第i+1期遥感影像数据Image_i+1和DEM数据Dem_i+1。那么CHi_i+1计算方法如下：
[0084] S22读取数据
[0085] 读取遥感影像Image_i和Image_i+1红绿蓝三个波段的原始矩阵数组，分别为ri、gi、bi和ri+1、gi+1、bi+1；读取DEM数据Dem_i和Dem_i+1单波段的原始矩阵数组，分别为di和di+1。其中ri、gi、bi和ri+1、gi+1、bi+1六个矩阵维度均完全相同；di和di+1矩阵维度相同。
[0086] S23差异度指标CHi_i+1衡量标准
[0087] CHi_i+1计算实际上是第i期和第i+1期遥感数据(包括遥感影像和数字高程模型DEM)在空间上的距离，如果差异度越小，两期间空间距离越小，反之越大。本申请采用统计学中常用的5种衡量空间距离大小的计算方法，分别为欧式空间距离(euc)、曼哈顿距离(man)、切比雪夫距离(che)、闵可夫斯基距离(mik)和余弦相似度(cos)。当然也可以采用其他空间距离指标用于衡量项目P第i期和第i+1期遥感数据的差异度。
[0088] S24项目P遥感影像差异度CHi_i+1_image
[0089] 计算遥感影像Image_i和Image_i+1红绿蓝三个波段上述5种方法的波段差异值CHi_i+1_image。此外，与红光波段、蓝光波段相比，绿光波段在项目建设阶段变化中具有更为重要的敏感性，因此，本申请中明确遥感影像红绿蓝三个波段差异度CH权重值确定为wr:wg:wb＝0.3:0.4:0.3。
[0090] 采用欧式空间距离方法计算遥感影像Image_i和Image_i+1红光波段差异值CHi_i+1_r_euc＝feuc(ri,ri+1)，绿光波段差异值CHi_i+1_g_euc＝feuc(gi,gi+1)和蓝光波段差异度值CHi_i+1_b_euc＝feuc(bi,bi+1)。最终第i期到第i+1期遥感影像差异度值为：
[0091] CHi_i+1_euc＝0.3*CHi_i+1_r_euc+0.4*CHi_i+1_g_euc+0.3*CHi_i+1_b_euc (3)；
[0092] 采用曼哈顿距离方法计算遥感影像Image_i和Image_i+1红光波段差异值CHi_i+1_r_man＝fman(ri,ri+1)，绿光波段差异值CHi_i+1_g_man＝fman(gi,gi+1)和蓝光波段差异度值CHi_i+1_b_man＝fman(bi,bi+1)。最终第i期到第i+1期遥感影像差异度值为：
[0093] CHi_i+1_man＝0.3*CHi_i+1_r_man+0.4*CHi_i+1_g_man+0.3*CHi_i+1_b_man (4)；
[0094] 采用切比雪夫距离方法计算遥感影像Image_i和Image_i+1红光波段差异值CHi_i+1_r_che＝fche(ri,ri+1)，绿光波段差异值CHi_i+1_g_che＝fche(gi,gi+1)和蓝光波段差异度值CHi_i+1_b_che＝fche(bi,bi+1)。最终第i期到第i+1期遥感影像差异度值为：
[0095] CHi_i+1_che＝0.3*CHi_i+1_r_che+0.4*CHi_i+1_g_che+0.3*CHi_i+1_b_che (5)；
[0096] 采用闵可夫斯基距离方法计算遥感影像Image_i和Image_i+1红光波段差异值CHi_i+1_r_mik＝fmik(ri,ri+1)，绿光波段差异值CHi_i+1_g_mik＝fmik(gi,gi+1)和蓝光波段差异度值CHi_i+1_b_mik＝fmik(bi,bi+1)。最终第i期到第i+1期遥感影像差异度值为：
[0097] CHi_i+1_mik＝0.3*CHi_i+1_r_mik+0.4*CHi_i+1_g_mik+0.3*CHi_i+1_b_mik (6)；
[0098] 采用余弦相似度方法计算遥感影像Image_i和Image_i+1红光波段差异值CHi_i+1r_cos＝fcos(ri,ri+1)，绿光波段差异值CHi_i+1_g_cos＝fcos(gi,gi+1)和蓝光波段差异度值CHi_i+1_b_cos＝fcos(bi,bi+1)。最终第i期到第i+1期遥感影像差异度值为：
[0099] CHi_i+1_cos＝0.3*CHi_i+1_r_cos+0.4*CHi_i+1_g_cos+0.3*CHi_i+1_b_cos (7)。
[0100] 最后将上述5种方法计算的结果求平均值得到第i期到第i+1期两期遥感影像差异度值：
[0101] CHi_i+1_image＝(CHi_i+1_euc+CHi_i+1_man+CHi_i+1_che+CHi_i+1_mik+CHi_i+1_cos)/5。
[0102] 上述5种计算方法的计算函数feuc、fman、fche、fmik和fcos均采用公开算法即可。
[0103] S25项目P数字高程模型DEM差异度CHi_i+1_dem
[0104] 基于项目P第i期和第i+1期数字高程模型DEM数据Dem_i和Dem_i+1计算项目P第i期到第i+1期数字高程模型DEM差异度CHi_i+1_dem指标。
[0105] 采用欧式空间距离方法计算结果为CHi_i+1_dem_euc＝feuc(di,di+1)。
[0106] 采用曼哈顿距离方法计算结果为CHi_i+1_dem_man＝fman(di,di+1)。
[0107] 采用切比雪夫距离方法计算结果为CHi_i+1_dem_che＝fche(di,di+1)。
[0108] 采用闵可夫斯基距离方法计算结果为CHi_i+1_dem_mik＝fmik(di,di+1)。
[0109] 采用余弦相似度方法计算结果为CHi_i+1_dem_cos＝fcos(di,di+1)。
[0110] 最后将上述5种方法计算的结果求平均值得到第i期到第i+1期两期DEM数据差异度
[0111] CHi_i+1_dem＝(CHi_i+1_dem_euc+CHi_i+1_dem_man+CHi_i+1_dem_che+CHi_i+1_dem_mik+CHi_i+1_dem_cos)/5。
[0112] S26项目从第i期到第i+1期差异度CHi_i+1
[0113] 基于S24和S25的计算结果，采用算数平均法的方法计算得到项目从第i期到第i+1期差异度CHi_i+1值。
[0114] CHi_i+1＝0.5*CHi_i+1_imge+0.5*CHi_i+1_dem
[0115] 对于特定区域而言，如果遥感影像在变化监测中发挥的作用更大、提供了更多关于变化监测的信息，可以将CHi_i+1_imge的计算权重适当增加，反之依然。
[0116] S27、项目建设阶段变化检测阈值
[0117] 上述S21‑S26提供了某研究区A的某个项目P从第i期到第i+1期变化监测差异度的计算方法，对于该研究区A的项目建设阶段变化检测阈值按照如下方法计算得到：
[0118] 收集研究区A的N个项目P1、P2、...、PN，进而分别收集各项目不同建设阶段的遥感影像和数字高程模型DEM数据，例如得到项目Pi的第j期和j+1期遥感影像和DEM数据分别为Imagei_j、Imagei_j+1、Demi_j和Demi_j+1，其中0≤i≤N。按照S21‑S26计算方法得到如下表格：
[0119] 表1项目建设阶段变化检测阈值计算表
[0120] CHi_j_j+1_image CHi_j_j+1_dem CHi_j_j+1项目P1第1期到第2期差异度 0.7727 0.6583 0.7155
项目P1第2期到第3期差异度 0.7738 0.5545 0.6642
项目P1第3期到第4期差异度 0.7210 0.6019 0.6615
项目P1第j期到第j+1期差异度 0.7531 0.6052 0.6792
项目P2第1期到第2期差异度 0.7556 0.6166 0.6861
项目P2第2期到第3期差异度 0.9469 0.6866 0.8168
项目P2第3期到第4期差异度 0.9760 0.7099 0.8430
项目P2第j期到第j+1期差异度 1.0550 0.6972 0.8761
项目P3第1期到第2期差异度 1.1482 0.7452 0.9467
项目P3第2期到第3期差异度 0.6926 0.7508 0.7217
项目P3第3期到第4期差异度 0.7298 0.4604 0.5951
项目P3第j期到第j+1期差异度 0.8243 0.6122 0.7183
项目P4第1期到第2期差异度 0.7662 0.6692 0.7177
项目P4第2期到第3期差异度 0.7901 0.5641 0.6771
项目P4第j期到第j+1期差异度 0.8291 0.6118 0.7205
项目Pi第1期到第2期差异度 0.6256 0.6162 0.6209
项目Pi第2期到第3期差异度 0.6663 0.4622 0.5643
项目Pi第3期到第4期差异度 0.8502 0.5383 0.6943
项目Pi第j期到第j+1期差异度 0.8415 0.5657 0.7036
[0121] 对于选定的N个项目而言，尽量保证项目类别的多样性，最大可能覆盖如下项目类别：交通工程、机场工程、输变电及电力工程、水利水电工程、采矿工程、油气开采与存储加工工程、工业园区工程、房地产工程、农业开发工程等，且要求N≥100。此外，对于某个项目Pi而言，其项目变化阶段应尽量多样化，包括如建设早期、建设高峰期、建设后期等。
[0122] 对表1中的最后一列取最小值即认为是某研究区A的项目建设阶段变化检测阈值CHthreshold。对于该研究区A的项目Pk而言，收集得到该项目两期(ti和ti+1)遥感数据(包括遥感影像和数字高程模型DEM数据)，计算得到该项目从ti到ti+1的变化差异度CHk_t_t+1，如果CHk_t_t+1≥CHthreshold，则表明项目Pk从ti到ti+1的建设阶段已经发生了变化；反之亦然。
[0123] 若检测结果表明项目B建设阶段发生了变化，则项目B当前遥感影像和DEM作为“待判定建设阶段的遥感数据”输入至建设阶段识别深度学习模型中，用以识别判定项目B的当前建设状态和建设日期。
[0124] S3、项目建设阶段判别
[0125] S31、项目建设阶段判别指标
[0126] S311项目建设阶段划分
[0127] 根据水土保持监督管理需求和工程建设管理行业规定，本申请将项目建设阶段划分为5个，分别为未开工、施工早期、施工高峰期、施工后期和完工，分别用0、1、2、3、4编码表示，记为Y。如Y＝1表示施工早期。
[0128] S312项目建设阶段判定指标确定与计算方法
[0129] 不同项目建设阶段最直接的表现有以下4个方面，分别是地表扰动面积、项目区高程复杂度、项目区林草植被覆盖情况以及项目区地表硬化面积。因此，本申请采用4个指标判定项目建设阶段，分别为地表扰动面积占比(X1)、项目区绝对高程均方差(X2)、项目区林草植被覆盖恢复率(X3)和项目区硬化面积占比(X4)。
[0130] 以项目i的建设阶段j为例说明各判定指标含义与计算方法如下：
[0131] 1)地表扰动面积占比X1_i_j＝S/ST，其中S表示i项目j建设阶段地表扰动面积，ST表示i项目设计占地总面积。ST数据来源于“建设项目信息库”，S数据基于i项目j建设阶段高分辨率遥感影像计算获取。
[0132] 2)项目区绝对高程均方差(X2_i_j)。基础数据源为项目i的建设阶段j时期的数字高程模型(DEM)，分辨率优于5m。假设项目i范围边界覆盖的DEM数据有N个像元，各像元值为dk(也是绝对高程值)，那么项目i的建设阶段j绝对高程均值(DEMmean)可由下式计算得到：
[0133] DEMmean＝(d1+d2+d3+...+dk+…+dN)/N (8)
[0134] 则项目项目i的建设阶段j绝对高程均方差X2_i_j可表示为：
[0135] X2_i_j＝([(d1‑DEMmean)2+(d2‑DEMmean)2+(d3‑DEMmean)2+…+(dN‑DEMmean)2+]/N)0.5[0136] (9)
[0137] 3)项目区林草植被覆盖恢复率(X3_i_j)。从“建设项目信息库”获取项目i水土保持方案水行政管理部门行政许可批复中的“林草植被覆盖率”，该值记为Vc。基于i项目j建设阶段高分辨率遥感影像，计算提取得到i项目当前阶段林草植被覆盖率，该值记为Vd。那么X3_i_j＝Vd/Vc。
[0138] 4)项目区硬化面积占比(X4_i_j)。基于i项目j建设阶段高分辨率遥感影像获取当前阶段项目区硬化面积，该值记为E。那么X4_i_j＝E/ST。
[0139] S32、项目建设阶段判定样本数据库建立
[0140] 收集研究区A的N个项目P1、P2、...、PN，同时收集每个项目5个建设阶段的高分遥感影像和数字高程模型DEM数据；如图4(a)‑图4(e)所示为未开工、施工早期、施工高峰期、施工后期和完工5个建设阶段。基于S312计算方法分别计算N各项目5个建设阶段判断的4个指标(X1、X2、X3、X4)，汇总得到下表：
[0141] 表2研究区A的N个项目5个建设阶段判定样本数据
[0142] 序号 X1_i_j X2_i_j X3_i_j X4_i_j Y 备注1 0.0514 22.0805 1.0000 0.0722 0 项目P1建设阶段1
2 0.1202 37.0211 1.0577 0.0401 1 项目P1建设阶段2
3 0.6893 78.0088 0.2644 0.1042 2 项目P1建设阶段3
4 0.2485 58.1761 0.7933 0.4809 3 项目P1建设阶段4
5 0.0481 21.1549 1.0000 0.7133 4 项目P1建设阶段5
...... ...... ...... ...... ...... ...... ......
...... ...... ...... ...... ...... ...... ......
...... ...... ...... ...... ...... ...... ......
5N ...... ...... ...... ...... ...... ......
[0143] S33、项目建设阶段判别模型建立
[0144] 基于上述4个判定指标确定生产建设项目所属建设阶段，实质上是基于已有指标进行分类。基于建立好的样本数据库，采用人工神经网络模型中的多层感知机(Multilayer Perceptron，MLP)用于构建项目建设阶段判别模型。
[0145] S331多层感知机模型搭建
[0146] 基于Tensorflow 2.0框架，借助机器学习中的经典模型—多层感知机(Multilayer Perceptron，MLP)搭建“项目建设阶段判别模型”。模型包含3层全连接层(Dense)，其中前2层神经元节点个数分别为16和8，激活函数均采用Relu，第三层神经元节点个数与目标多分类个数相同为4个，激活函数采用Softmax。模型训练优化器采用“自适应梯度优化算法(Adam)”，相对于传统的梯度优化算法(如基于梯度下降的SGD或SGD+Momentum等)，Adam可以自适应地调整每个参数的学习率，在训练过程中逐渐减小学习率，提高模型的训练效率和泛化能力。采用多分类交叉熵损失函数(categorical_
crossentropy)表征和记录模型训练过程中的损失。所搭建完成的“项目建设阶段判别模型”记录为MLP_model。
[0147] 除本申请提到的上述多层感知机模型外，诸如卷积神经网络、支持向量机、随机森林等均可以实现上述目标。
[0148] S332多层感知机模型训练与验证
[0149] 基于S32所构建的项目建设阶段判别样本数据库，对S331构建的多层感知机模型进行训练，得到最终的项目建设阶段判别多层感知机模型(记为MLPtrained)。图5(a)和图5(b)为模型训练过程的精度和损失的变化。
[0150] 使用与上述模型训练不同的样本数据对得到的MLPtrained模型进行验证(图6)，20组数据中仅有1组发生了预测错误，总体预测精度达95％，这表明所构建的项目建设阶段判别多层感知机模型MLPtrained具有较强的应用泛化能力。
[0151] S333建设阶段识别结果应用
[0152] 对于研究区A的某个项目U，得到其某个时期的遥感影像和正射高程模型DEM数据，采用S312方法计算得到项目建设阶段判断的4个指标向量：X＝[X1X2X3X4]，将向量X输入到训练好的模型MLPtrained中，得到项目建设阶段预测编码Y，若Y＝0，则项目当前建设状态是“未开工”；若Y＝1，则项目当前建设状态是“施工早期”；若Y＝2，则项目当前建设状态是“施工高峰期”；若Y＝3，则项目当前建设状态是“施工后期”；若Y＝4，则项目当前建设状态是“完工”，同时更新项目当前建设状态的日期。
[0153] S4、生产建设项目建设全周期变化检测方法总流程
[0154] 获取第i个时期某项目F遥感影像和数字高程模型DEM数据，以及F项目第i‑1个时期的项目建设阶段(记为Status)；
[0155] 基于S21项目建设阶段变化监测技术，判断F项目第i个时期的建设状态是否发生改变，如果没有发生改变，则更新F项目第i个时期的建设阶段仍旧为Status；如果建设状态发生改变，则基于S3项目建设阶段判别技术获取项目F第i个时期项目的建设阶段(记为Status')，并更新F项目第i个时期的建设阶段为Status'。
[0156] 实施例3
[0157] 本申请所述介质包括数据端、存储器等。数据端以软件管理界面体现，包含访问数据库接口、自动更新数据库模块、建设阶段变化检测技术模块、建设阶段识别深度学习模型模块、产出成果编辑模块等。存储器主要具有数据运算、超大空间存储功能，是数据库和软件的载体。
[0158] 本实施例公开了生产建设项目建设全周期变化检测方法，上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。
[0159] 本申请优点：通过现有的工作成果基础数据，能够获取特定区域项目建设阶段变化检测阈值，得到项目建设阶段判别模型，进而通过非接触非现场的遥感数据即可了解特定区域某生产建设项目当前的建设状态，可为水行政管理部门对生产建设项目监督管理提供强有力的技术支撑，避免对建设周期较长的生产建设项目的监管遗漏等问题。
