CN111310322A - 基于元胞自动机的干旱区自然绿洲空间动态模拟方法 - Google Patents
基于元胞自动机的干旱区自然绿洲空间动态模拟方法 Download PDFInfo
- Publication number
- CN111310322A CN111310322A CN202010084349.3A CN202010084349A CN111310322A CN 111310322 A CN111310322 A CN 111310322A CN 202010084349 A CN202010084349 A CN 202010084349A CN 111310322 A CN111310322 A CN 111310322A
- Authority
- CN
- China
- Prior art keywords
- neighborhood
- vegetation
- ndvi
- rate
- cellular
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 230000001413 cellular effect Effects 0.000 title claims abstract description 65
- 238000000034 method Methods 0.000 title claims abstract description 53
- 238000005094 computer simulation Methods 0.000 title claims abstract description 29
- 238000006243 chemical reaction Methods 0.000 claims abstract description 160
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims abstract description 90
- 230000003068 static effect Effects 0.000 claims abstract description 40
- 238000004088 simulation Methods 0.000 claims abstract description 37
- 230000008859 change Effects 0.000 claims abstract description 34
- 238000009933 burial Methods 0.000 claims abstract description 28
- 230000007704 transition Effects 0.000 claims description 87
- 238000011084 recovery Methods 0.000 claims description 25
- 239000003673 groundwater Substances 0.000 claims description 24
- 238000011160 research Methods 0.000 claims description 24
- 239000011248 coating agent Substances 0.000 claims description 23
- 238000000576 coating method Methods 0.000 claims description 23
- 238000005259 measurement Methods 0.000 claims description 20
- 238000004458 analytical method Methods 0.000 claims description 15
- 230000008569 process Effects 0.000 claims description 14
- 230000001932 seasonal effect Effects 0.000 claims description 13
- 238000012512 characterization method Methods 0.000 claims description 3
- 239000000126 substance Substances 0.000 claims description 3
- 238000011156 evaluation Methods 0.000 abstract description 7
- 230000000694 effects Effects 0.000 abstract description 6
- 230000006870 function Effects 0.000 description 35
- 230000009466 transformation Effects 0.000 description 13
- 238000010586 diagram Methods 0.000 description 9
- 230000015556 catabolic process Effects 0.000 description 3
- 238000006731 degradation reaction Methods 0.000 description 3
- 230000007613 environmental effect Effects 0.000 description 3
- 238000012886 linear function Methods 0.000 description 3
- 238000012887 quadratic function Methods 0.000 description 3
- 241001494479 Pecora Species 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 230000009467 reduction Effects 0.000 description 2
- 230000004044 response Effects 0.000 description 2
- 239000002689 soil Substances 0.000 description 2
- 208000005156 Dehydration Diseases 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000008020 evaporation Effects 0.000 description 1
- 238000001704 evaporation Methods 0.000 description 1
- 238000009472 formulation Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 238000001556 precipitation Methods 0.000 description 1
- 239000002352 surface water Substances 0.000 description 1
- 230000002194 synthesizing effect Effects 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/002—Biomolecular computers, i.e. using biomolecules, proteins, cells
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
- G06Q50/06—Energy or water supply
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Business, Economics & Management (AREA)
- Theoretical Computer Science (AREA)
- Health & Medical Sciences (AREA)
- General Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Economics (AREA)
- Biophysics (AREA)
- General Physics & Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Organic Chemistry (AREA)
- Evolutionary Computation (AREA)
- Computational Linguistics (AREA)
- Molecular Biology (AREA)
- Computing Systems (AREA)
- Biomedical Technology (AREA)
- Artificial Intelligence (AREA)
- Mathematical Physics (AREA)
- Software Systems (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Data Mining & Analysis (AREA)
- Chemical & Material Sciences (AREA)
- Public Health (AREA)
- Water Supply & Treatment (AREA)
- Human Resources & Organizations (AREA)
- Marketing (AREA)
- Primary Health Care (AREA)
- Strategic Management (AREA)
- Tourism & Hospitality (AREA)
- General Business, Economics & Management (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种上述基于元胞自动机的干旱区自然绿洲空间动态模拟方法，基于元胞自动机的生态输水驱动下干旱区自然绿洲空间动态模拟方法，基于通过遥感影像获取的归一化植被指数、地下水埋深观测数据等多源信息，采用地下水承载力函数和绿洲面积动态函数模拟绿洲总面积，由NDVI序列构建元胞NDVI、邻域植被率、邻域淹没率与元胞状态变化的概率函数，采用t时刻的元胞NDVI、邻域植被率、邻域淹没率计算t+1时刻的覆被转换概率，采用加权平均法计算各元胞的覆被转换得分并从大到小排序，综合绿洲面积模拟结果设定覆被转换动态阈值，并结合覆被转换静态阈值最终判定元胞的状态如何变化，得到绿洲空间动态模拟结果并通过精度评价反映模型模拟效果。
Description
技术领域
本发明涉及地球物理下生态水文分支技术领域，尤其涉及一种基于元胞自动机的干旱区自然绿洲空间动态模拟方法。
背景技术
干旱区河岸带及尾闾绿洲生态系统是典型的地下水依赖型生态系统。流域水资源不合理开发利用过度挤占生态用水，地表水量减少、地下水位降低诱发的水分胁迫导致自然绿洲生态退化。生态输水是干旱区河岸带及尾闾自然绿洲生态系统修复及保护的重要措施之一，掌握绿洲面积变化及空间分布动态是评估生态输水环境效应、优化生态输水策略的必要条件。目前主要通过实地样方调查、遥感影像解译等手段识别生态输水影响范围，评估自然绿洲恢复区的面积；建立绿洲恢复区植被指数与地下水埋深的函数关系，进而模拟生态输水、地下水埋深变化影响下绿洲植被指数的动态变化。这类模拟方法一般先划分绿洲恢复区的空间范围，考虑的绿洲动态变化是恢复区植被指数空间均值的时间动态变化，即模拟静态空间范围内植被指数空间均值的时间动态对生态输水的响应过程，无法反映生态输水累积影响下、地下水位抬升过程中绿洲恢复的空间动态变化。因而，传统的绿洲动态模拟方案效果差。
发明内容
针对以上问题，本发明提出一种基于元胞自动机的干旱区自然绿洲空间动态模拟方法。
为实现本发明的目的，提供一种基于元胞自动机的干旱区自然绿洲空间动态模拟方法，包括如下步骤：
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；其中，
计算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数据，统计元胞NDVI的实际发生值xi-元胞NDVI(i＝1,...,I元胞NDVI)；统计不同元胞NDVIxi-元胞NDVI所对应的0→1事件次数c01-元胞NDVI、0→0事件次数c00-元胞NDVI、1→0事件次数c10-元胞NDVI、1→1事件次数c11-元胞NDVI；第一0至1转换概率P01-元胞NDVI为：P01-元胞NDVI＝c01-元胞NDVI/(c01-元胞NDVI+c00-元胞NDVI)，第一1至0转换概率P10-元胞NDVI为：P10-元胞NDVI＝c10-元胞NDVI/(c10-元胞NDVI+c11-元胞NDVI)；
对于邻域植被率，基于研究区历年植被分布二值化图分析得到的邻域植被率，统计邻域植被率的实际发生值xi-邻域植被率(i＝1,...,I邻域植被率)；统计不同邻域植被率xi-邻域植被率所对应的0→1事件次数c01-邻域植被率、0→0事件次数c00-邻域植被率、1→0事件次数c10-邻域植被率、1→1事件次数c11-邻域植被率；第二0至1转换概率为：P01-邻域植被率＝c01-邻域植被率/(c01-邻域植被率+c00-邻域植被率)，第二1至0转换概率为：P10-邻域植被率＝c10-邻域植被率/(c10-邻域植被率+c11-邻域植被率)；
对于邻域淹没率，基于研究区历年的淹没区分布二值化图分析得到的邻域淹没率，统计邻域淹没率的实际发生值xi-邻域淹没率(i＝1,...,I邻域淹没率)；统计不同邻域淹没率xi-邻域淹没率所对应的0→1事件次数c01-邻域淹没率、0→0事件次数c00-邻域淹没率、1→0事件次数c10-邻域淹没率、1→1事件次数c11-邻域淹没率；第三0至1转换概率为：P01-邻域淹没率＝c01-邻域淹没率/(c01-邻域淹没率+c00-邻域淹没率)，第三1至0转换概率为：P10-邻域淹没率＝c10-邻域淹没率/(c10-邻域淹没率+c11-邻域淹没率)。
在一个实施例中，根据归一化后的第一0至1转换概率、第二0至1转换概率和第三0至1转换概率计算0至1覆被转换得分S01，根据归一化后的第一1至0转换概率、第二1至0转换概率和第三1至0转换概率计算1至0覆被转换得分S10包括：
S01＝w元胞NDVI×P′01-元胞NDVI+w邻域植被率×P′01-邻域植被率+w邻域淹没率×P′01-邻域淹没率，
S10＝w元胞NDVI×P′10-元胞NDVI+w邻域植被率×P′10-邻域植被率+w邻域淹没率×P′10-邻域淹没率，
其中，P′01-元胞NDVI为归一化后的第一0至1转换概率、P′01-邻域植被率为归一化后的第二0至1转换概率、P′01-邻域淹没率为归一化后的第三0至1转换概率，P′10-元胞NDVI为归一化后的第一1至0转换概率、P′10-邻域植被率为归一化后的第二1至0转换概率、P′10-邻域淹没率为归一化后的第三1至0转换概率。
上述基于元胞自动机的干旱区自然绿洲空间动态模拟方法，通过获取生态输水后研究区的年平均地下水埋深序列、年内最大NDVI序列和年内最小NDVI序列，根据年内最大NDVI数据生成植被分布二值化图，根据年内最小NDVI数据生成淹没区二值化图，基于历年的植被分布二值化图统计生态输水后的绿洲恢复面积，采用Sigmoid方程构建地下水承载力函数，以计算各个地下水埋深条件下研究区所能够承载的最大绿洲面积，采用Verhulst方程构建绿洲面积动态函数，以模拟生态输水驱动下地下水恢复过程中绿洲面积的动态变化，以元胞为中心设置邻域，基于植被分布二值化图统计邻域植被率，基于淹没区二值化图统计邻域淹没率，统计元胞NDVI对应的第一覆被转换概率，邻域植被率对应的第二覆被转换概率，和邻域淹没率对应的第三覆被转换概率；所述第一覆被转换概率包括第一0至1转换概率和第一1至0转换概率；所述第二覆被转换概率包括第二0至1转换概率和第二1至0转换概率；所述第三覆被转换概率包括第三0至1转换概率和第三1至0转换概率，采用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，在t时刻0至1覆被转换得分S01、1至0覆被转换得分S10、第一状态转换阈值S01阈值、第二状态转换阈值S10阈值的基础上，生成t+1时刻的元胞状态，得到t+1时刻植被分布二值化图，作为模拟分析的结果，所得到的模拟分析结果具有更高的准确性，可以提升相应的模拟效果，克服现有技术无法模拟自然绿洲恢复面积及空间分布动态变化的不足，提高生态输水、地下水埋深、植被指数观测资料的利用效率，提升干旱区生态输水环境效应的评估能力及评估成果的合理性。
附图说明
图1是一个实施例的基于元胞自动机的干旱区自然绿洲空间动态模拟方法流程图；
图2是一个实施例的基于元胞自动机的生态输水驱动下干旱区自然绿洲空间动态模拟方法的流程图；
图3是一个实施例的甘肃省石羊河流域尾闾青土湖绿洲位置示意图；
图4是一个实施例的元胞NDVI、邻域植被率、邻域淹没率所对应的覆被转换概率散点及概率函数拟合线图；
图5是一个实施例的青土湖绿洲2018年实测植被分布、2019年实测植被分布及2019年模拟植被分布二值化图。
具体实施方式
为了使本申请的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本申请进行进一步详细说明。应当理解，此处描述的具体实施例仅仅用以解释本申请，并不用于限定本申请。
在本文中提及“实施例”意味着，结合实施例描述的特定特征、结构或特性可以包含在本申请的至少一个实施例中。在说明书中的各个位置出现该短语并不一定均是指相同的实施例，也不是与其它实施例互斥的独立的或备选的实施例。本领域技术人员显式地和隐式地理解的是，本文所描述的实施例可以与其它实施例相结合。
参考图1所示，图1为一个实施例的基于元胞自动机的干旱区自然绿洲空间动态模拟方法流程图，包括如下步骤：
S10，获取生态输水后研究区的年平均地下水埋深序列、年内最大NDVI(Normalized Difference Vegetation Index，归一化植被指数)序列和年内最小NDVI序列。
具体地，可以获取生态输水后研究区的地下水埋深观测资料，得到年平均地下水埋深序列；获取同期的遥感影像资料，得到年内最大NDVI序列(考虑植被生长季节，一般按1-12月统计)和年内最小NDVI序列(考虑生态输水时机，例如冬季至次年春季输水，按9月至次年8月统计)。
S20，根据年内最大NDVI数据生成植被分布二值化图，根据年内最小NDVI数据生成淹没区二值化图；所述植被分布二值化图包括植被恢复的元胞和植被尚未恢复的元胞；所述年内最小NDVI数据包括季节性淹没区和季节性非淹没区。
在一个实施例中，根据年内最大NDVI数据生成植被分布二值化图，根据年内最小NDVI数据生成淹没区二值化图包括：
在年内最大NDVI数据中，将NDVI值大于或等于0.2的栅格判定为生态输水影响下植被恢复的元胞，用1表征，将NDVI值小于0.2的栅格判定为植被尚未恢复的元胞，用0表征，根据用1表征的元胞和用0表征的元胞生成植被分布二值化图；
在年内最小NDVI数据中，将NDVI值小于0的栅格判定为生态输水影响下的季节性淹没区，用1表征，将NDVI值大于或等于0的栅格判定为非淹没区，用0表征，根据用1表征的季节性淹没区和用0表征的非淹没区生成淹没区二值化图。
具体地，本实施例可以采用年内最大NDVI数据，基于NDVI的植被阈值(一般取0.2)生成植被分布二值化图，即NDVI≥0.2的栅格判定为生态输水影响下植被恢复的元胞(记为1)，NDVI<0.2的栅格判定为植被尚未恢复的元胞(记为0)；采用年内最小NDVI数据，基于NDVI的水体阈值(一般取0)生成淹没区二值化图，即NDVI<0的栅格判定为生态输水影响下的季节性淹没区(记为1)，NDVI≥0的栅格判定为非淹没区(记为0)。
S30，基于历年的植被分布二值化图统计生态输水后的绿洲恢复面积，采用Sigmoid方程构建地下水承载力函数，以计算各个地下水埋深条件下研究区所能够承载的最大绿洲面积，采用Verhulst方程构建绿洲面积动态函数，以模拟生态输水驱动下地下水恢复过程中绿洲面积的动态变化。
基于历年的植被分布二值化图统计生态输水后恢复的绿洲面积，可用植被栅格数表征；采用Sigmoid方程构建地下水承载力函数，用于计算某一地下水埋深条件下研究区所能够承载的最大绿洲面积，为求解Verhulst方程提供必要的参数；采用Verhulst方程构建绿洲面积动态函数，用于模拟生态输水驱动下地下水恢复过程中绿洲面积的动态变化，可用于预测某一时刻、某一地下水埋深条件下的绿洲面积。
在一个实施例中，Sigmoid方程包括：
Verhulst方程包括：
其中，ht为t时刻地下水埋深，VGCC(ht)表示当地下水埋深为ht时，地下水资源对于植被的承载能力(Groundwater Carrying Capacity,GCC)，以NDVI≥0.2的植被栅格数目表征；Vmin为地下水对植被的最低承载力，Vmax为地下水对植被的最高承载力；h(Vmin+Vmax)/2为地下水承载力为(Vmin+Vmax)/2时所对应的地下水埋深；s表征Vmin至Vmax之间Sigmoid方程曲线变化的陡峭程度；V为植被栅格(元胞)数目；r表征内禀增长率。由年平均地下水埋深和绿洲面积(以植被栅格数目表征)数据率定Sigmoid方程和Verhulst方程的参数。
S40，以元胞为中心设置邻域，基于植被分布二值化图统计邻域植被率，基于淹没区二值化图统计邻域淹没率；所述邻域植被率为相应邻域内的植被栅格数目占邻域栅格总数的比例；所述邻域淹没率为相应邻域内的淹没区栅格数目占邻域栅格总数的比例。
具体地，可以以元胞为中心，设置大小为m×n(一般取m＝n＝3)的邻域。例如当m＝n＝3时，元胞本身及其相邻的左上角、上方、右上角、左侧、右侧、左下角、下方、右下角共计九个栅格即构成邻域。基于植被分布二值化图统计邻域植被率，即邻域内的植被栅格数目占邻域栅格总数的比例；基于淹没区二值化图统计邻域淹没率，即邻域内的淹没区栅格数目占邻域栅格总数的比例。
S50，统计元胞NDVI对应的第一覆被转换概率，邻域植被率对应的第二覆被转换概率，和邻域淹没率对应的第三覆被转换概率；所述第一覆被转换概率包括第一0至1转换概率和第一1至0转换概率；所述第二覆被转换概率包括第二0至1转换概率和第二1至0转换概率；所述第三覆被转换概率包括第三0至1转换概率和第三1至0转换概率。
上述步骤可以采用步骤S40类似的方法分别统计元胞NDVI、邻域植被率、邻域淹没率所对应的覆被转换概率(包括元胞NDVI对应的第一覆被转换概率，邻域植被率对应的第二覆被转换概率，和邻域淹没率对应的第三覆被转换概率)，以经验频次表征概率值，覆被转换类型包括元胞类别由0→1(0至1)的植被恢复变化和元胞类别由1→0(1至0)的植被退化变化。
S60，采用t时刻的元胞NDVI、邻域植被率、邻域淹没率作为元胞自动机模型的输入，对于t时刻覆被类型为0的元胞统计第一0至1转换概率P01-元胞NDVI、第二0至1转换概率P01-邻域植被率和第三0至1转换概率P01-邻域淹没率，对于t时刻覆被类型为1的元胞统计第一1至0转换概率P10-元胞NDVI、第二1至0转换概率P10-邻域植被率和第三1至0转换概率P10-邻域淹没率，根据归一化后的第一0至1转换概率、第二0至1转换概率和第三0至1转换概率计算0至1覆被转换得分S01，根据归一化后的第一1至0转换概率、第二1至0转换概率和第三1至0转换概率计算1至0覆被转换得分S10。
上述步骤可以采用线性变换P′＝(P-P理论最小值)/(P理论最大值-P理论最小值)将P01-元胞NDVI、P01-邻域植被率、P01-邻域淹没率、P10-元胞NDVI、P10-邻域植被率、P10-邻域淹没率归一化到[0,1]范围内，得到归一化后的第一0至1转换概率P′01-元胞NDVI、归一化后的第二0至1转换概率P′01-邻域植被率、归一化后的第三0至1转换概率P′01-邻域淹没率，归一化后的第一1至0转换概率P′10-元胞NDVI、归一化后的第二1至0转换概率P′10-邻域植被率、归一化后的第三1至0转换概率P′10-邻域淹没率。
在一个实施例中，根据归一化后的第一0至1转换概率、第二0至1转换概率和第三0至1转换概率计算0至1覆被转换得分S01，根据归一化后的第一1至0转换概率、第二1至0转换概率和第三1至0转换概率计算1至0覆被转换得分S10包括：
S01＝w元胞NDVI×P′01-元胞NDVI+w邻域植被率×P′01-邻域植被率+w邻域淹没率×P′01-邻域淹没率，
S10＝w元胞NDVI×P′10-元胞NDVI+w邻域植被率×P′10-邻域植被率+w邻域淹没率×P′10-邻域淹没率，
其中，P′01-元胞NDVI为归一化后的第一0至1转换概率、P′01-邻域植被率为归一化后的第二0至1转换概率、P′01-邻域淹没率为归一化后的第三0至1转换概率，P′10-元胞NDVI为归一化后的第一1至0转换概率、P′10-邻域植被率为归一化后的第二1至0转换概率、P′10-邻域淹没率为归一化后的第三1至0转换概率。具体地，P′01-元胞NDVI、P′01-邻域植被率、P′01-邻域淹没率、P′10-元胞NDVI、P′10-邻域植被率、P′10-邻域淹没率分别为归一化后的元胞状态转换概率值，即由归一化后的概率值表征得分值；w元胞NDVI、w邻域植被率、w邻域淹没率可以按元胞NDVI、邻域植被率、邻域淹没率影响的重要程度梯度取值，三者总和为1，例如可分别取0.3、0.2、0.5。将S01、S10从大到小排序，得分高的元胞则倾向于状态变化。
S80，在t时刻0至1覆被转换得分S01、1至0覆被转换得分S10、第一状态转换阈值S01阈值、第二状态转换阈值S10阈值的基础上，生成t+1时刻的元胞状态，得到t+1时刻植被分布二值化图，作为模拟分析的结果。
在t时刻元胞状态转换得分S01、S10计算和状态转换阈值S01阈值、S10阈值分析基础上，依据元胞状态转换规则生成t+1时刻的元胞状态，即可得t+1时刻植被分布二值化图，作为模拟分析的结果。元胞状态转换规则包括：
①对于t时刻的状态为0的元胞，若S01≥S01阈值，则t+1时刻元胞状态发生变化，由0→1，否则t+1时刻元胞维持t时刻的状态0；
②对于t时刻的状态为1的元胞，若S10≥S10阈值，则t+1时刻元胞状态发生变化，由1→0，否则t+1时刻元胞维持t时刻的状态1。
本实施例基于元胞自动机的生态输水驱动下干旱区自然绿洲空间动态模拟方法，基于通过遥感影像获取的归一化植被指数、地下水埋深观测数据等多源信息，采用地下水承载力函数和绿洲面积动态函数模拟绿洲总面积，由NDVI序列构建元胞NDVI、邻域植被率、邻域淹没率与元胞状态变化的概率函数，采用t时刻的元胞NDVI、邻域植被率、邻域淹没率计算t+1时刻的覆被转换概率，即元胞恢复为植被的概率和元胞植被退化的概率，经过线性变换后将取值范围统一为[0,1]，采用加权平均法计算各元胞的覆被转换得分并从大到小排序，综合绿洲面积模拟结果设定覆被转换动态阈值，并结合覆被转换静态阈值最终判定元胞的状态如何变化，得到绿洲空间动态模拟结果并通过精度评价反映模型模拟效果。
上述基于元胞自动机的干旱区自然绿洲空间动态模拟方法，通过获取生态输水后研究区的年平均地下水埋深序列、年内最大NDVI序列和年内最小NDVI序列，根据年内最大NDVI数据生成植被分布二值化图，根据年内最小NDVI数据生成淹没区二值化图，基于历年的植被分布二值化图统计生态输水后的绿洲恢复面积，采用Sigmoid方程构建地下水承载力函数，以计算各个地下水埋深条件下研究区所能够承载的最大绿洲面积，采用Verhulst方程构建绿洲面积动态函数，以模拟生态输水驱动下地下水恢复过程中绿洲面积的动态变化，以元胞为中心设置邻域，基于植被分布二值化图统计邻域植被率，基于淹没区二值化图统计邻域淹没率，统计元胞NDVI对应的第一覆被转换概率，邻域植被率对应的第二覆被转换概率，和邻域淹没率对应的第三覆被转换概率；所述第一覆被转换概率包括第一0至1转换概率和第一1至0转换概率；所述第二覆被转换概率包括第二0至1转换概率和第二1至0转换概率；所述第三覆被转换概率包括第三0至1转换概率和第三1至0转换概率，采用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，在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阈值的栅格数NS10＞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。
上述S01静态阈值、S10静态阈值可以分别取0.7等值。
在一个实施例中，设置元胞由0→1、1→0覆被转换的静态阈值S01静态阈值、S10静态阈值，一般两者均取0.7。由绿洲面积动态函数计算t+1时刻绿洲面积，即植被栅格数Vt+1，依据绿洲面积的动态变化情况设置覆被转换的动态阈值。若Vt+1＝Vt，则绿洲维持动态平衡，S01阈值＝S01静态阈值、S10阈值＝S10静态阈值。若Vt+1>Vt，则绿洲总体倾向于恢复，S10阈值＝S10静态阈值，统计S10＞S10阈值的栅格数
对应的覆被变化由0→1的元胞数目为
从S01序列(降序排列)中统计排位第
位的S01值作为动态阈值，即S01阈值＝S01动态阈值。类似地，若Vt+1<Vt，则绿洲总体倾向于退化，S01阈值＝S01静态阈值，统计S01＞S01阈值的栅格数
对应的覆被变化由1→0的元胞数目为
从S10序列(降序排列)中统计排位第位的S10值作为动态阈值，即S10阈值＝S10动态阈值。
在一个实施例中，上述基于元胞自动机的干旱区自然绿洲空间动态模拟方法，还包括：
获取t+1时刻的实测植被分布二值化图和模拟植被分布二值化图，根据实测植被分布二值化图和模拟植被分布二值化图确定绿洲空间动态模拟的精度。
具体地，确定绿洲空间动态模拟的精度的过程可以包括：定性地，若模拟图和实测图匹配性、一致性较好，则绿洲空间动态模拟的精度较高；定量地，可采用Kappa系数k。
作为一个实施例，根据实测植被分布二值化图和模拟植被分布二值化图确定绿洲空间动态模拟的精度包括：
根据实测植被分布二值化图和模拟植被分布二值化图确定正确模拟的状态为0的第一元胞数目n0模拟正确、正确模拟的状态为1的第二元胞数目n1模拟正确、状态为0的实测元胞数目n0实测、状态为0的模拟元胞数目n0模拟，状态为1的实测元胞数目n1实测、以及状态为1的模拟元胞数目n1模拟；
计算总体分类精度po和中间变量pe；其中，
计算Kappa系数k，根据Kappa系数k确定绿洲空间动态模拟的精度；其中：
具体地，其中，po是正确模拟的栅格数量之和除以研究区总栅格数N，也就是总体分类精度，其本身也可作为精度评价的一种指标。pe不同，仅仅是计算Kappa系数的中间变量；n0模拟正确、n1模拟正确分别为正确模拟的状态为0的元胞数目和状态为1的元胞数目；n0实测、n0模拟分别为状态为0的实测元胞数目和模拟元胞数目，n1实测、n1模拟分别为状态为1的实测元胞数目和模拟元胞数目。Kappa系数k∈[-1,1]，数值越大表明模拟值与实测值的一致性程度越高，即模型效果越好、模拟精度越高；一般认为，k∈(0.2,0.4]表明模拟值与实测值的一致性一般，k∈(0.4,0.6]表明模拟值与实测值的一致性中等，k∈(0.6,0.8]表明模拟值与实测值具有高度的一致性，k∈(0.8,1]表明模拟值与实测值几乎完全一致。
在一个实施例中，统计元胞NDVI对应的第一覆被转换概率，邻域植被率对应的第二覆被转换概率，和邻域淹没率对应的第三覆被转换概率包括：
对于元胞NDVI，基于研究区历年的年内最大NDVI数据，统计元胞NDVI的实际发生值xi-元胞NDVI(i＝1,...,I元胞NDVI)；统计不同元胞NDVIxi-元胞NDVI所对应的0→1(0至1)事件次数c01-元胞NDVI、0→0(0至0)事件次数c00-元胞NDVI、1→0(1至0)事件次数c10-元胞NDVI、1→1(1至1)事件次数c11-元胞NDVI；第一0至1转换概率(0→1覆被转换概率)P01-元胞NDVI为：P01-元胞NDVI＝c01-元胞NDVI/(c01-元胞NDVI+c00-元胞NDVI)，第一1至0转换概率(1→0覆被转换概率)P10-元胞NDVI为：P10-元胞NDVI＝c10-元胞NDVI/(c10-元胞NDVI+c11-元胞NDVI)；以xi-元胞NDVI为横坐标，以不同xi-元胞NDVI所对应的P01-元胞NDVI、P10-元胞NDVI绘制散点图，采用合适的函数(例如二次函数、指数函数、线性函数等)拟合散点图，得P01-元胞NDVI概率函数和P10-元胞NDVI概率函数；计算P01-元胞NDVI、P10-元胞NDVI概率函数理论最大值、理论最小值，用于下一步覆被转换得分的线性变换。
对于邻域植被率，基于研究区历年植被分布二值化图分析得到的邻域植被率，统计邻域植被率的实际发生值xi-邻域植被率(i＝1,...,I邻域植被率)；统计不同邻域植被率xi-邻域植被率所对应的0→1事件次数c01-邻域植被率、0→0事件次数c00-邻域植被率、1→0事件次数c10-邻域植被率、1→1事件次数c11-邻域植被率；第二0至1转换概率为：P01-邻域植被率＝c01-邻域植被率/(c01-邻域植被率+c00-邻域植被率)，第二1至0转换概率为：P10-邻域植被率＝c10-邻域植被率/(c10-邻域植被率+c11-邻域植被率)；以xi-邻域植被率为横坐标，以不同xi-邻域植被率所对应的P01-邻域植被率、P10-邻域植被率绘制散点图，采用合适的函数(例如二次函数、指数函数、线性函数等)拟合散点图，得P01-邻域植被率概率函数和P10-邻域植被率概率函数；计算P01-邻域植被率、P10-邻域植被率概率函数理论最大值、理论最小值，用于下一步覆被转换得分的线性变换。
对于邻域淹没率，基于研究区历年的淹没区分布二值化图分析得到的邻域淹没率，统计邻域淹没率的实际发生值xi-邻域淹没率(i＝1,...,I邻域淹没率)；统计不同邻域淹没率xi-邻域淹没率所对应的0→1事件次数c01-邻域淹没率、0→0事件次数c00-邻域淹没率、1→0事件次数c10-邻域淹没率、1→1事件次数c11-邻域淹没率；第三0至1转换概率为：P01-邻域淹没率＝c01-邻域淹没率/(c01-邻域淹没率+c00-邻域淹没率)，第三1至0转换概率为：P10-邻域淹没率＝c10-邻域淹没率/(c10-邻域淹没率+c11-邻域淹没率)；以xi-邻域淹没率为横坐标，以不同xi-邻域淹没率所对应的P01-邻域淹没率、P10-邻域淹没率绘制散点图，采用合适的函数(例如二次函数、指数函数、线性函数等)拟合散点图，得P01-邻域淹没率概率函数和P10-邻域淹没率概率函数；计算P01-邻域淹没率、P10-邻域淹没率概率函数理论最大值、理论最小值，用于下一步覆被转换得分的线性变换。
在一个实施例中，参考图2所示，图2为本实施例的基于元胞自动机的生态输水驱动下干旱区自然绿洲空间动态模拟方法的流程图，图3为本实施例的甘肃省石羊河流域尾闾青土湖绿洲位置示意图。青土湖绿洲位于我国西北河西走廊石羊河流域尾闾，该区属于温带大陆性干旱荒漠气候，年均气温8.8℃，≥10℃的积温3036℃，无霜期162d，年日照时数3189h，年降水量110mm，蒸发量2640mm。青土湖原为石羊河尾闾湖泊，于1959年干涸。根据《石羊河流域重点治理规划》，从2010年9月以渠道输送的形式开始向青土湖注入生态用水。截至2018年，已连续9年有计划地向青土湖下泄生态用水，累计下泄水量2.5亿立方米。
第一步：从青土湖水文站收集2010-2018年年平均地下水埋深数据资料；从美国国家航空航天局(National Aeronautics and Space Administration)网站下载同期逐月的MODIS TERRA MOD13Q1数据产品，青土湖绿洲所处图幅编号为h26v05；对下载的MODIS数据进行研究区域图幅裁剪、NDVI数值系数校正等常规遥感影像预处理；按1-12月统计年内最大NDVI值，用于后续的植被分布分析；按9月至次年8月统计年内最小NDVI值，用于后续的淹没区分析。
第二步：采用年内最大NDVI数据，以0.2为阈值生成历年植被分布二值化图，用于后续的邻域植被率分析；采用年内最小NDVI数据，以0为阈值生成历年淹没区分布二值化图，用于后续的邻域淹没率分析。
第三步：基于植被分布二值化图统计历年的植被栅格数目，结合地下水埋深数据，以1年为时间步长，率定地下水承载力Sigmoid方程和绿洲面积动态Verhulst方程：
第四步：以元胞为中心，一个栅格即为一个元胞，设置大小为3×3的九宫格邻域；基于植被分布二值化图统计历年的邻域植被率，用于后续的邻域植被率影响下元胞覆被转换概率统计；基于淹没区二值化图统计历年的邻域淹没率，用于后续的邻域淹没率影像下元胞覆被转换概率统计。
第五步：分别统计元胞NDVI、邻域植被率、邻域淹没率所对应的覆被转换概率。图4为元胞NDVI、邻域植被率、邻域淹没率所对应的覆被转换概率散点及概率函数拟合线图，拟合线方程分别为：
P01-元胞NDVI：y＝-0.10+0.01exp(18.62x)
P01-邻域植被率：y＝0.65-1.46(0.67-x)2
P01-邻域淹没率：y＝1.69-1.64exp(-0.26x)
P10-元胞NDVI：y＝-0.03+3.75exp(-10.64x)
P10-邻域植被率：y＝0.13+1.06(x-0.84)2
P10-邻域淹没率：y＝0.03+0.36exp(-3.93x)
第六步：以2018年的元胞NDVI、邻域植被率、邻域淹没率作为元胞自动机模型的输入，利用第五步的概率函数拟合线方程统计各个元胞的P01-元胞NDVI、P01-邻域植被率、P01-邻域淹没率、P10-元胞NDVI、P10-邻域植被率、P10-邻域淹没率；采用线性变换将取值归一化到[0,1]范围内，得到元胞NDVI、邻域植被率、邻域淹没率驱动元胞状态转换的各分项得分值；将元胞NDVI、邻域植被率、邻域淹没率的权重分别设为0.3、0.2、0.5，采用加权平均法计算各元胞状态转换总得分值。
第七步：设置元胞状态转换静态阈值S01静态阈值＝0.7、S10静态阈值＝0.7。根据第三步绿洲面积动态模拟计算2019年植被栅格数V2019＝385，2018年植被栅格数V2018＝359。此种情形下V2019＞V2018，绿洲总体趋向于恢复，S10阈值＝S10静态阈值＝0.7，统计2018年S10＞0.7的栅格数
对应的覆被变化由0→1的元胞数目
在从大到小排序的S01序列中排序第28位的数值为0.4，可得S01阈值＝S01动态阈值＝0.4。
第八步：根据元胞状态转换规则及第六步分析所得的元胞状态转换得分、第七步分析所得的元胞状态转换阈值，生成2019年植被分布二值化模拟图。图5为青土湖绿洲2018年实测植被分布、2019年实测植被分布及2019年模拟植被分布二值化图。
第九步：采用Kappa系数评价2019年实测的植被分布二值化图和模拟的植被分布二值化图的一致性，k＝0.8，说明基于元胞自动机的生态输水驱动下干旱区自然绿洲空间动态模拟方法精度较高。
本实施例进行生态输水驱动下干旱区自然绿洲空间动态模拟的应用是：可以用于预测、评估生态输水的绿洲恢复效应，对于不同生态输水规模下地下水埋深响应的不同工况及淹没区范围变化的不同工况，比对植被空间分布的变化情况，优选生态输水方案，支撑生态输水策略的制定和改进。
以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。
需要说明的是，本申请实施例所涉及的术语“第一\第二\第三”仅仅是区别类似的对象，不代表针对对象的特定排序，可以理解地，“第一\第二\第三”在允许的情况下可以互换特定的顺序或先后次序。应该理解“第一\第二\第三”区分的对象在适当情况下可以互换，以使这里描述的本申请的实施例能够以除了在这里图示或描述的那些以外的顺序实施。
本申请实施例的术语“包括”和“具有”以及它们任何变形，意图在于覆盖不排他的包含。例如包含了一系列步骤或模块的过程、方法、装置、产品或设备没有限定于已列出的步骤或模块，而是可选地还包括没有列出的步骤或模块，或可选地还包括对于这些过程、方法、产品或设备固有的其它步骤或模块。
以上所述实施例仅表达了本申请的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本申请构思的前提下，还可以做出若干变形和改进，这些都属于本申请的保护范围。因此，本申请专利的保护范围应以所附权利要求为准。
Claims (8)
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；其中，
计算Kappa系数k，根据Kappa系数k确定绿洲空间动态模拟的精度；其中：
5.根据权利要求1至4任一项所述的基于元胞自动机的干旱区自然绿洲空间动态模拟方法，其特征在于，根据年内最大NDVI数据生成植被分布二值化图，根据年内最小NDVI数据生成淹没区二值化图包括：
在年内最大NDVI数据中，将NDVI值大于或等于0.2的栅格判定为生态输水影响下植被恢复的元胞，用1表征，将NDVI值小于0.2的栅格判定为植被尚未恢复的元胞，用0表征，根据用1表征的元胞和用0表征的元胞生成植被分布二值化图；
在年内最小NDVI数据中，将NDVI值小于0的栅格判定为生态输水影响下的季节性淹没区，用1表征，将NDVI值大于或等于0的栅格判定为非淹没区，用0表征，根据用1表征的季节性淹没区和用0表征的非淹没区生成淹没区二值化图。
6.根据权利要求1至4任一项所述的基于元胞自动机的干旱区自然绿洲空间动态模拟方法，其特征在于，Sigmoid方程包括：
Verhulst方程包括：
其中，ht为t时刻地下水埋深，VGCC(ht)表示当地下水埋深为ht时，地下水资源对于植被的承载能力，以NDVI≥0.2的植被栅格数目表征；Vmin为地下水对植被的最低承载力，Vmax为地下水对植被的最高承载力；
为地下水承载力为(Vmin+Vmax)/2时所对应的地下水埋深；s表征Vmin至Vmax之间Sigmoid方程曲线变化的陡峭程度；V为植被栅格数目；r表征内禀增长率。
7.根据权利要求1至4任一项所述的基于元胞自动机的干旱区自然绿洲空间动态模拟方法，其特征在于，统计元胞NDVI对应的第一覆被转换概率，邻域植被率对应的第二覆被转换概率，和邻域淹没率对应的第三覆被转换概率包括：
对于元胞NDVI，基于研究区历年的年内最大NDVI数据，统计元胞NDVI的实际发生值xi-元胞NDVI(i＝1,...,I元胞NDVI)；统计不同元胞NDVIxi-元胞NDVI所对应的0→1事件次数c01-元胞NDVI、0→0事件次数c00-元胞NDVI、1→0事件次数c10-元胞NDVI、1→1事件次数c11-元胞NDVI；第一0至1转换概率P01-元胞NDVI为：P01-元胞NDVI＝c01-元胞NDVI/(c01-元胞NDVI+c00-元胞NDVI)，第一1至0转换概率P10-元胞NDVI为：P10-元胞NDVI＝c10-元胞NDVI/(c10-元胞NDVI+c11-元胞NDVI)；
对于邻域植被率，基于研究区历年植被分布二值化图分析得到的邻域植被率，统计邻域植被率的实际发生值xi-邻域植被率(i＝1,...,I邻域植被率)；统计不同邻域植被率xi-邻域植被率所对应的0→1事件次数c01-邻域植被率、0→0事件次数c00-邻域植被率、1→0事件次数c10-邻域植被率、1→1事件次数c11-邻域植被率；第二0至1转换概率为：P01-邻域植被率＝c01-邻域植被率/(c01-邻域植被率+c00-邻域植被率)，第二1至0转换概率为：P10-邻域植被率＝c10-邻域植被率/(c10-邻域植被率+c11-邻域植被率)；
对于邻域淹没率，基于研究区历年的淹没区分布二值化图分析得到的邻域淹没率，统计邻域淹没率的实际发生值xi-邻域淹没率(i＝1,...,I邻域淹没率)；统计不同邻域淹没率xi-邻域淹没率所对应的0→1事件次数c01-邻域淹没率、0→0事件次数c00-邻域淹没率、1→0事件次数c10-邻域淹没率、1→1事件次数c11-邻域淹没率；第三0至1转换概率为：P01-邻域淹没率＝c01-邻域淹没率/(c01-邻域淹没率+c00-邻域淹没率)，第三1至0转换概率为：P10-邻域淹没率＝c10-邻域淹没率/(c10-邻域淹没率+c11-邻域淹没率)。
8.根据权利要求1至4任一项所述的基于元胞自动机的干旱区自然绿洲空间动态模拟方法，其特征在于，根据归一化后的第一0至1转换概率、第二0至1转换概率和第三0至1转换概率计算0至1覆被转换得分S01，根据归一化后的第一1至0转换概率、第二1至0转换概率和第三1至0转换概率计算1至0覆被转换得分S10包括：
S01＝w元胞NDVI×P′01-元胞NDVI+w邻域植被率×P′01-邻域植被率+w邻域淹没率×P′01-邻域淹没率，
S10＝w元胞NDVI×P′10-元胞NDVI+w邻域植被率×P′10-邻域植被率+w邻域淹没率×P′10-邻域淹没率，
其中，P′01-元胞NDVI为归一化后的第一0至1转换概率、P′01-邻域植被率为归一化后的第二0至1转换概率、P′01-邻域淹没率为归一化后的第三0至1转换概率，P′10-元胞NDVI为归一化后的第一1至0转换概率、P′10-邻域植被率为归一化后的第二1至0转换概率、P′10-邻域淹没率为归一化后的第三1至0转换概率。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010084349.3A CN111310322B (zh) | 2020-02-10 | 2020-02-10 | 基于元胞自动机的干旱区自然绿洲空间动态模拟方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010084349.3A CN111310322B (zh) | 2020-02-10 | 2020-02-10 | 基于元胞自动机的干旱区自然绿洲空间动态模拟方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111310322A true CN111310322A (zh) | 2020-06-19 |
| CN111310322B CN111310322B (zh) | 2021-03-16 |
Family
ID=71148273
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010084349.3A Active CN111310322B (zh) | 2020-02-10 | 2020-02-10 | 基于元胞自动机的干旱区自然绿洲空间动态模拟方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111310322B (zh) |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111881879A (zh) * | 2020-08-07 | 2020-11-03 | 宁波市农业科学研究院 | 基于航拍提升机插田杂草取样精度的网格排序取样调查法 |
| CN112270115A (zh) * | 2020-11-25 | 2021-01-26 | 同济大学 | 一种基于元胞自动机的复杂地形洪水淹没进程模拟方法 |
| CN112527442A (zh) * | 2020-12-24 | 2021-03-19 | 中国科学院地理科学与资源研究所 | 一种环境数据多维显示方法、装置、介质及终端设备 |
| CN112668448A (zh) * | 2020-12-24 | 2021-04-16 | 中国科学院地理科学与资源研究所 | 一种生态进程变化分析方法、装置、介质及终端设备 |
| CN114005039A (zh) * | 2021-12-31 | 2022-02-01 | 成都国星宇航科技有限公司 | 基于遥感图像的农作物长势评估方法、装置及电子设备 |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101289251A (zh) * | 2008-05-09 | 2008-10-22 | 北京工业大学 | 基于pca模型的活性污泥吸附和沉降过程的仿真方法 |
| CN102546158A (zh) * | 2011-12-22 | 2012-07-04 | 河海大学 | 一种基于奇偶元胞自动机的分组加密方法 |
| CN105139331A (zh) * | 2015-07-23 | 2015-12-09 | 河海大学 | 一种面向比特的基于仿生元胞自动机的图像置乱方法 |
| US20190220028A1 (en) * | 2019-03-26 | 2019-07-18 | Intel Corporation | Computer-assisted (ca)/autonomous driving (ad) vehicle inference model creation |
| CN110031596A (zh) * | 2019-03-14 | 2019-07-19 | 中国科学院新疆生态与地理研究所 | 一种针对干旱区内陆河流域地下水的实时动态监测方法 |
| KR20190108428A (ko) * | 2018-03-14 | 2019-09-24 | 서울시립대학교 산학협력단 | 셀룰러 오토마타 알고리즘을 이용한 보행 상황 시뮬레이션 방법 및 장치 |
-
2020
- 2020-02-10 CN CN202010084349.3A patent/CN111310322B/zh active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101289251A (zh) * | 2008-05-09 | 2008-10-22 | 北京工业大学 | 基于pca模型的活性污泥吸附和沉降过程的仿真方法 |
| CN102546158A (zh) * | 2011-12-22 | 2012-07-04 | 河海大学 | 一种基于奇偶元胞自动机的分组加密方法 |
| CN105139331A (zh) * | 2015-07-23 | 2015-12-09 | 河海大学 | 一种面向比特的基于仿生元胞自动机的图像置乱方法 |
| KR20190108428A (ko) * | 2018-03-14 | 2019-09-24 | 서울시립대학교 산학협력단 | 셀룰러 오토마타 알고리즘을 이용한 보행 상황 시뮬레이션 방법 및 장치 |
| CN110031596A (zh) * | 2019-03-14 | 2019-07-19 | 中国科学院新疆生态与地理研究所 | 一种针对干旱区内陆河流域地下水的实时动态监测方法 |
| US20190220028A1 (en) * | 2019-03-26 | 2019-07-18 | Intel Corporation | Computer-assisted (ca)/autonomous driving (ad) vehicle inference model creation |
Non-Patent Citations (3)
| Title |
|---|
| YONGJIU FENG 等: "Modeling urban growth with GIS based cellular automata and least squares SVM rules: a case study in Qingpu–Songjiang", 《STOCH ENVIRON RES RISK ASSESS (2016)》 * |
| 井长青 等: "耦合神经网络与元胞自动机的城市土地利用动态演化模型", 《干旱研究》 * |
| 宋冬梅 等: "基于元胞自动机民勤绿洲湖区荒漠化演化预测", 《中国沙漠》 * |
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111881879A (zh) * | 2020-08-07 | 2020-11-03 | 宁波市农业科学研究院 | 基于航拍提升机插田杂草取样精度的网格排序取样调查法 |
| CN111881879B (zh) * | 2020-08-07 | 2022-06-17 | 宁波市农业科学研究院 | 基于航拍提升机插田杂草取样精度的网格排序取样调查法 |
| CN112270115A (zh) * | 2020-11-25 | 2021-01-26 | 同济大学 | 一种基于元胞自动机的复杂地形洪水淹没进程模拟方法 |
| CN112270115B (zh) * | 2020-11-25 | 2022-06-14 | 同济大学 | 一种基于元胞自动机的复杂地形洪水淹没进程模拟方法 |
| CN112527442A (zh) * | 2020-12-24 | 2021-03-19 | 中国科学院地理科学与资源研究所 | 一种环境数据多维显示方法、装置、介质及终端设备 |
| CN112668448A (zh) * | 2020-12-24 | 2021-04-16 | 中国科学院地理科学与资源研究所 | 一种生态进程变化分析方法、装置、介质及终端设备 |
| CN114005039A (zh) * | 2021-12-31 | 2022-02-01 | 成都国星宇航科技有限公司 | 基于遥感图像的农作物长势评估方法、装置及电子设备 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111310322B (zh) | 2021-03-16 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111310322B (zh) | 基于元胞自动机的干旱区自然绿洲空间动态模拟方法 | |
| Lu et al. | Quantifying the impacts of small dam construction on hydrological alterations in the Jiulong River basin of Southeast China | |
| Rocha et al. | Assessing the impacts of sustainable agricultural practices for water quality improvements in the Vouga catchment (Portugal) using the SWAT model | |
| Vanmaercke et al. | A comparison of measured catchment sediment yields with measured and predicted hillslope erosion rates in Europe | |
| CN108875161A (zh) | 基于卷积神经网络深度学习的流量等级预测方法 | |
| Safavi et al. | Conjunctive use of surface and ground water resources using the ant system optimization | |
| Juston et al. | Temporal sampling strategies and uncertainty in calibrating a conceptual hydrological model for a small boreal catchment | |
| Awotwi et al. | Climate change impact on streamflow in a tropical basin of Ghana, West Africa | |
| CN113361742A (zh) | 一种基于水文模拟的区域综合干旱识别方法 | |
| Yasarer et al. | Characterizing ponds in a watershed simulation and evaluating their influence on streamflow in a Mississippi watershed | |
| Hernandez-Suarez et al. | Evaluation of the impacts of hydrologic model calibration methods on predictability of ecologically-relevant hydrologic indices | |
| Shooshtari et al. | Land cover change modelling in Hyrcanian forests, northern Iran: a landscape pattern and transformation analysis perspective | |
| Biggs et al. | The potential for refining nitrogen fertiliser management through accounting for climate impacts: An exploratory study for the Tully region | |
| Ochoa-Zavala et al. | Reduction of genetic variation when far from the niche centroid: Prediction for mangrove species | |
| Xiao et al. | Impact of large-scale tree planting in Yunnan Province, China, on the water supply balance in Southeast Asia | |
| Lo | Quantifying soil erosion for the Shihmen reservoir watershed, Taiwan | |
| Zeng et al. | A land-indicator-based optimization model with trading mechanism in wetland ecosystem under uncertainty | |
| Meena et al. | Information and Communication Technologies for Sustainable Natural Resource Management | |
| CN113850465B (zh) | 一种无资料地区的水文干旱监测系统 | |
| CN112966657A (zh) | 一种大尺度水体覆盖的遥感自动分类方法 | |
| Ceola et al. | Unraveling long-term flood risk dynamics across the murray-darling basin using a large-scale hydraulic model and satellite data | |
| Bromand | Impact assessment of climate change on water resources in the Kabul River Basin, Afghanistan | |
| Jahangir et al. | Application of artificial neural networks to the simulation of climate elements, drought forecast by two indicators of SPI and PNPI, and mapping of drought intensity; case study of Khorasan Razavi | |
| Handayani et al. | Modelling of land use change in indramayu district, west Java Province | |
| Mehrani et al. | Medium-term changes in hydrology of warm Mediterranean wetlands under land use/cover change |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |