基于元胞自动机的干旱区自然绿洲空间动态模拟方法技术

本发明专利技术公开了一种上述基于元胞自动机的干旱区自然绿洲空间动态模拟方法，基于元胞自动机的生态输水驱动下干旱区自然绿洲空间动态模拟方法，基于通过遥感影像获取的归一化植被指数、地下水埋深观测数据等多源信息，采用地下水承载力函数和绿洲面积动态函数模拟绿洲总面积，由NDVI序列构建元胞NDVI、邻域植被率、邻域淹没率与元胞状态变化的概率函数，采用t时刻的元胞NDVI、邻域植被率、邻域淹没率计算t+1时刻的覆被转换概率，采用加权平均法计算各元胞的覆被转换得分并从大到小排序，综合绿洲面积模拟结果设定覆被转换动态阈值，并结合覆被转换静态阈值最终判定元胞的状态如何变化，得到绿洲空间动态模拟结果并通过精度评价反映模型模拟效果。

Spatial dynamic simulation method of natural oasis in arid area based on Cellular Automata

【技术实现步骤摘要】
基于元胞自动机的干旱区自然绿洲空间动态模拟方法
本专利技术涉及地球物理下生态水文分支
，尤其涉及一种基于元胞自动机的干旱区自然绿洲空间动态模拟方法。
技术介绍
干旱区河岸带及尾闾绿洲生态系统是典型的地下水依赖型生态系统。流域水资源不合理开发利用过度挤占生态用水，地表水量减少、地下水位降低诱发的水分胁迫导致自然绿洲生态退化。生态输水是干旱区河岸带及尾闾自然绿洲生态系统修复及保护的重要措施之一，掌握绿洲面积变化及空间分布动态是评估生态输水环境效应、优化生态输水策略的必要条件。目前主要通过实地样方调查、遥感影像解译等手段识别生态输水影响范围，评估自然绿洲恢复区的面积；建立绿洲恢复区植被指数与地下水埋深的函数关系，进而模拟生态输水、地下水埋深变化影响下绿洲植被指数的动态变化。这类模拟方法一般先划分绿洲恢复区的空间范围，考虑的绿洲动态变化是恢复区植被指数空间均值的时间动态变化，即模拟静态空间范围内植被指数空间均值的时间动态对生态输水的响应过程，无法反映生态输水累积影响下、地下水位抬升过程中绿洲恢复的空间动态变化。因而，传统的绿洲动态模拟方案效果差。
技术实现思路
针对以上问题，本专利技术提出一种基于元胞自动机的干旱区自然绿洲空间动态模拟方法。为实现本专利技术的目的，提供一种基于元胞自动机的干旱区自然绿洲空间动态模拟方法，包括如下步骤：S10，获取生态输水后研究区的年平均地下水埋深序列、年内最大NDVI序列和年内最小NDVI序列；S20，根据年内最大NDVI数据生成植被分布二值化图，根据年内最小NDVI数据生成淹没区二值化图；所述植被分布二值化图包括植被恢复的元胞和植被尚未恢复的元胞；所述年内最小NDVI数据包括季节性淹没区和季节性非淹没区；S30，基于历年的植被分布二值化图统计生态输水后的绿洲恢复面积，采用Sigmoid方程构建地下水承载力函数，以计算各个地下水埋深条件下研究区所能够承载的最大绿洲面积，采用Verhulst方程构建绿洲面积动态函数，以模拟生态输水驱动下地下水恢复过程中绿洲面积的动态变化；S40，以元胞为中心设置邻域，基于植被分布二值化图统计邻域植被率，基于淹没区二值化图统计邻域淹没率；所述邻域植被率为相应邻域内的植被栅格数目占邻域栅格总数的比例；所述邻域淹没率为相应邻域内的淹没区栅格数目占邻域栅格总数的比例；S50，统计元胞NDVI对应的第一覆被转换概率，邻域植被率对应的第二覆被转换概率，和邻域淹没率对应的第三覆被转换概率；所述第一覆被转换概率包括第一0至1转换概率和第一1至0转换概率；所述第二覆被转换概率包括第二0至1转换概率和第二1至0转换概率；所述第三覆被转换概率包括第三0至1转换概率和第三1至0转换概率；S60，采用t时刻的元胞NDVI、邻域植被率、邻域淹没率作为元胞自动机模型的输入，对于t时刻覆被类型为0的元胞统计第一0至1转换概率、第二0至1转换概率和第三0至1转换概率，对于t时刻覆被类型为1的元胞统计第一1至0转换概率、第二1至0转换概率和第三1至0转换概率，根据归一化后的第一0至1转换概率、第二0至1转换概率和第三0至1转换概率计算0至1覆被转换得分S01，根据归一化后的第一1至0转换概率、第二1至0转换概率和第三1至0转换概率计算1至0覆被转换得分S10；S80，在t时刻0至1覆被转换得分S01、1至0覆被转换得分S10、第一状态转换阈值S01阈值、第二状态转换阈值S10阈值的基础上，生成t+1时刻的元胞状态，得到t+1时刻植被分布二值化图，作为模拟分析的结果。在一个实施例中，步骤S60之后，还包括：S70，设置元胞由0至1覆被转换的静态阈值S01静态阈值，以及元胞由1至0覆被转换的静态阈值S10静态阈值，根据绿洲面积动态函数计算t+1时刻的植被栅格数Vt+1，若Vt+1＝Vt，则绿洲维持动态平衡，S01阈值＝S01静态阈值，S10阈值＝S10静态阈值；若Vt+1>Vt，则S10阈值＝S10静态阈值，统计S10＞S10阈值的栅格数S10＞S10阈值对应的覆被变化由0至1的元胞数目为从S01序列中统计排位第位的特定S01值作为动态阈值S01动态阈值，令S01阈值＝S01动态阈值；若Vt+1<Vt，则S01阈值＝S01静态阈值，统计S01＞S01阈值的栅格数S01＞S01阈值对应的覆被变化由1至0的元胞数目为从S10序列中统计排位第位的S10值作为动态阈值S10动态阈值，令S10阈值＝S10动态阈值；Vt表示研究区在t时刻对应的植被栅格数，S01序列包括按降序排列的各个S01，S10序列包括按降序排列的各个S10。在一个实施例中，基于元胞自动机的干旱区自然绿洲空间动态模拟方法，还包括：获取t+1时刻的实测植被分布二值化图和模拟植被分布二值化图，根据实测植被分布二值化图和模拟植被分布二值化图确定绿洲空间动态模拟的精度。作为一个实施例，根据实测植被分布二值化图和模拟植被分布二值化图确定绿洲空间动态模拟的精度包括：根据实测植被分布二值化图和模拟植被分布二值化图确定正确模拟的状态为0的第一元胞数目n0模拟正确、正确模拟的状态为1的第二元胞数目n1模拟正确、状态为0的实测元胞数目n0实测、状态为0的模拟元胞数目n0模拟，状态为1的实测元胞数目n1实测、以及状态为1的模拟元胞数目n1模拟；计算总体分类精度po和中间变量pe；其中，N表示研究区总栅格数；计算Kappa系数k，根据Kappa系数k确定绿洲空间动态模拟的精度；其中：在一个实施例中，根据年内最大NDVI数据生成植被分布二值化图，根据年内最小NDVI数据生成淹没区二值化图包括：在年内最大NDVI数据中，将NDVI值大于或等于0.2的栅格判定为生态输水影响下植被恢复的元胞，用1表征，将NDVI值小于0.2的栅格判定为植被尚未恢复的元胞，用0表征，根据用1表征的元胞和用0表征的元胞生成植被分布二值化图；在年内最小NDVI数据中，将NDVI值小于0的栅格判定为生态输水影响下的季节性淹没区，用1表征，将NDVI值大于或等于0的栅格判定为非淹没区，用0表征，根据用1表征的季节性淹没区和用0表征的非淹没区生成淹没区二值化图。在一个实施例中，Sigmoid方程包括：Verhulst方程包括：其中，ht为t时刻地下水埋深，VGCC(ht)表示当地下水埋深为ht时，地下水资源对于植被的承载能力，以NDVI≥0.2的植被栅格数目表征；Vmin为地下水对植被的最低承载力，Vmax为地下水对植被的最高承载力；为地下水承载力为(Vmin+Vmax)/2时所对应的地下水埋深；s表征Vmin至Vmax之间Sigmoid方程曲线变化的陡峭程度；V为植被栅格数目；r表征内禀增长率。在一个实施例中，统计元胞NDVI对应的第一覆被转换概率，邻域植被率对应的第二覆被转换概率，和邻域淹没率对应的第三覆被转换概率包括：对于元胞NDVI，基于研究区历年的年内最大NDVI数据，统计元胞本文档来自技高网...

【技术保护点】
1.一种基于元胞自动机的干旱区自然绿洲空间动态模拟方法，其特征在于，包括如下步骤：/nS10，获取生态输水后研究区的年平均地下水埋深序列、年内最大NDVI序列和年内最小NDVI序列；/nS20，根据年内最大NDVI数据生成植被分布二值化图，根据年内最小NDVI数据生成淹没区二值化图；所述植被分布二值化图包括植被恢复的元胞和植被尚未恢复的元胞；所述年内最小NDVI数据包括季节性淹没区和季节性非淹没区；/nS30，基于历年的植被分布二值化图统计生态输水后的绿洲恢复面积，采用Sigmoid方程构建地下水承载力函数，以计算各个地下水埋深条件下研究区所能够承载的最大绿洲面积，采用Verhulst方程构建绿洲面积动态函数，以模拟生态输水驱动下地下水恢复过程中绿洲面积的动态变化；/nS40，以元胞为中心设置邻域，基于植被分布二值化图统计邻域植被率，基于淹没区二值化图统计邻域淹没率；所述邻域植被率为相应邻域内的植被栅格数目占邻域栅格总数的比例；所述邻域淹没率为相应邻域内的淹没区栅格数目占邻域栅格总数的比例；/nS50，统计元胞NDVI对应的第一覆被转换概率，邻域植被率对应的第二覆被转换概率，和邻域淹没率对应的第三覆被转换概率；所述第一覆被转换概率包括第一0至1转换概率和第一1至0转换概率；所述第二覆被转换概率包括第二0至1转换概率和第二1至0转换概率；所述第三覆被转换概率包括第三0至1转换概率和第三1至0转换概率；/nS60，采用t时刻的元胞NDVI、邻域植被率、邻域淹没率作为元胞自动机模型的输入，对于t时刻覆被类型为0的元胞统计第一0至1转换概率、第二0至1转换概率和第三0至1转换概率，对于t时刻覆被类型为1的元胞统计第一1至0转换概率、第二1至0转换概率和第三1至0转换概率，根据归一化后的第一0至1转换概率、第二0至1转换概率和第三0至1转换概率计算0至1覆被转换得分S...

【技术特征摘要】
1.一种基于元胞自动机的干旱区自然绿洲空间动态模拟方法，其特征在于，包括如下步骤：
S10，获取生态输水后研究区的年平均地下水埋深序列、年内最大NDVI序列和年内最小NDVI序列；
S20，根据年内最大NDVI数据生成植被分布二值化图，根据年内最小NDVI数据生成淹没区二值化图；所述植被分布二值化图包括植被恢复的元胞和植被尚未恢复的元胞；所述年内最小NDVI数据包括季节性淹没区和季节性非淹没区；
S30，基于历年的植被分布二值化图统计生态输水后的绿洲恢复面积，采用Sigmoid方程构建地下水承载力函数，以计算各个地下水埋深条件下研究区所能够承载的最大绿洲面积，采用Verhulst方程构建绿洲面积动态函数，以模拟生态输水驱动下地下水恢复过程中绿洲面积的动态变化；
S40，以元胞为中心设置邻域，基于植被分布二值化图统计邻域植被率，基于淹没区二值化图统计邻域淹没率；所述邻域植被率为相应邻域内的植被栅格数目占邻域栅格总数的比例；所述邻域淹没率为相应邻域内的淹没区栅格数目占邻域栅格总数的比例；
S50，统计元胞NDVI对应的第一覆被转换概率，邻域植被率对应的第二覆被转换概率，和邻域淹没率对应的第三覆被转换概率；所述第一覆被转换概率包括第一0至1转换概率和第一1至0转换概率；所述第二覆被转换概率包括第二0至1转换概率和第二1至0转换概率；所述第三覆被转换概率包括第三0至1转换概率和第三1至0转换概率；
S60，采用t时刻的元胞NDVI、邻域植被率、邻域淹没率作为元胞自动机模型的输入，对于t时刻覆被类型为0的元胞统计第一0至1转换概率、第二0至1转换概率和第三0至1转换概率，对于t时刻覆被类型为1的元胞统计第一1至0转换概率、第二1至0转换概率和第三1至0转换概率，根据归一化后的第一0至1转换概率、第二0至1转换概率和第三0至1转换概率计算0至1覆被转换得分S01，根据归一化后的第一1至0转换概率、第二1至0转换概率和第三1至0转换概率计算1至0覆被转换得分S10；
S80，在t时刻0至1覆被转换得分S01、1至0覆被转换得分S10、第一状态转换阈值S01阈值、第二状态转换阈值S10阈值的基础上，生成t+1时刻的元胞状态，得到t+1时刻植被分布二值化图，作为模拟分析的结果。


2.根据权利要求1所述的基于元胞自动机的干旱区自然绿洲空间动态模拟方法，其特征在于，步骤S60之后，还包括：
S70，设置元胞由0至1覆被转换的静态阈值S01静态阈值，以及元胞由1至0覆被转换的静态阈值S10静态阈值，根据绿洲面积动态函数计算t+1时刻的植被栅格数Vt+1，若Vt+1＝Vt，则绿洲维持动态平衡，S01阈值＝S01静态阈值，S10阈值＝S10静态阈值；若Vt+1>Vt，则S10阈值＝S10静态阈值，统计S10＞S10阈值的栅格数S10＞S10阈值对应的覆被变化由0至1的元胞数目为从S01序列中统计排位第位的特定S01值作为动态阈值S01动态阈值，令S01阈值＝S01动态阈值；若Vt+1<Vt，则S01阈值＝S01静态阈值，统计S01＞S01阈值的栅格数S01＞S01阈值对应的覆被变化由1至0的元胞数目为从S10序列中统计排位第位的S10值作为动态阈值S10动态阈值，令S10阈值＝S10动态阈值；Vt表示研究区在t时刻对应的植被栅格数，S01序列包括按降序排列的各个S01，S10序列包括按降序排列的各个S10。


3.根据权利要求1所述的基于元胞自动机的干旱区自然绿洲空间动态模拟方法，其特征在于，还包括：
获取t+1时刻的实测植被分布二值化图和模拟植被分布二值化图，根据实测植被分布二值化图和模拟植被分布二值化图确定绿洲空间动态模拟的精度。


4.根据权利要求3所述的基于元胞自动机的干旱区自然绿洲空间动态模拟方法，其特征在于，根据实测植被分布二值化图和模拟植被分布二值化图确定绿洲空间动态模拟的精度包括：
根据实测植被分布二值化图和模拟植被分布二值化图确定正确模拟的状态为0的第一元胞数目n0模拟正确、正确模拟的状态为1的第二元胞数目n1模拟正确、状态为0的实测元胞数目n0实测、状态为0的模拟元胞数目n0模拟，状态为1的实测元胞数目n1实测、以及状态为1的模拟元胞数目n1模拟；
计算总体分类精度po和中间变量pe；其中，

N表示研究区总栅格数；
计算Kappa系数k，根据Kappa系数k确定绿洲空间动态模拟的精度；其中：





5.根据权利要求1至4任一项所述的基于元胞自动机的干旱区自然绿洲空间动态模拟方法，其特征在于，根据年内最大NDVI数据...

【专利技术属性】
技术研发人员：黄峰，潘若云，杨丽娟，张清清，
申请(专利权)人：河海大学，
类型：发明
国别省市：江苏;32