CN113361774B - 一种生态空间集约利用优化方法及系统 - Google Patents

一种生态空间集约利用优化方法及系统 Download PDF

Info

Publication number
CN113361774B
CN113361774B CN202110633616.2A CN202110633616A CN113361774B CN 113361774 B CN113361774 B CN 113361774B CN 202110633616 A CN202110633616 A CN 202110633616A CN 113361774 B CN113361774 B CN 113361774B
Authority
CN
China
Prior art keywords
ecological
space
service function
ecosystem
year
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.)
Active
Application number
CN202110633616.2A
Other languages
English (en)
Other versions
CN113361774A (zh
Inventor
欧定华
秦景
高雪松
夏建国
周伟
吴能君
李媛茜
邓欧平
刘涛
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Sichuan Agricultural University
Original Assignee
Sichuan Agricultural University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Sichuan Agricultural University filed Critical Sichuan Agricultural University
Priority to CN202110633616.2A priority Critical patent/CN113361774B/zh
Publication of CN113361774A publication Critical patent/CN113361774A/zh
Application granted granted Critical
Publication of CN113361774B publication Critical patent/CN113361774B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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
    • G06Q10/00—Administration; Management
    • G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
    • 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
    • G06Q10/00—Administration; Management
    • G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
    • G06Q10/063—Operations research, analysis or management
    • G06Q10/0639—Performance analysis of employees; Performance analysis of enterprise or organisation operations
    • G06Q10/06393—Score-carding, benchmarking or key performance indicator [KPI] analysis
    • 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—Systems or methods specially adapted for specific business sectors, e.g. utilities or tourism
    • G06Q50/10—Services
    • G06Q50/16—Real estate
    • G06Q50/165—Land development
    • 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—Systems or methods specially adapted for specific business sectors, e.g. utilities or tourism
    • G06Q50/10—Services
    • G06Q50/26—Government or public services

Abstract

本发明公开了一种生态空间集约利用优化方法及系统，该方法包括：S1：根据多年生态系统服务功能综合值及其变化特征，动态获取生态空间；S2：根据各单项生态系统服务功能值，识别生态空间主导生态系统服务功能；模拟基期年至目标年生态空间主导生态系统服务功能格局，进行主导生态系统服务功能演变分析；S3：采用LP‑Markov‑CA复合模型，进行目标年主导生态系统服务功能数量结构和空间布局的耦合优化，实现生态空间主导生态系统服务功能的集约利用优化。本发明实现了生态空间数量结构、空间布局与服务功能的协同优化，成功获取了生态系统服务功能最大的生态空间集约利用布局，为国土空间集约利用优化配置提供了新方法。

Description

一种生态空间集约利用优化方法及系统
技术领域
本发明涉及国土空间规划技术领域，具体涉及一种生态空间集约利用优化方法及系统。
背景技术
生态空间是具有自然属性、以提供生态产品或生态服务为主体功能的国土空间。生态空间集约利用优化是通过对一定时期生态系统服务功能数量结构和空间布局的优化配置，实现有限生态空间功能效益最大化的一种空间治理行为，是优化国土空间开发保护格局，实现空间资源集约利用的有效途径。生态空间与其他国土空间一样具备面积有限性和供给稀缺性，同样面临空间集约利用问题。因此，如何实现生态空间集约利用优化，让有限生态空间发挥更大的生态系统服务效能，对提高生态系统服务功能、保障区域生态安全，优化国土空间开发格局、促进资源节约集约利用，缓和发展与保护的矛盾、实现可持续发展具有重要理论和现实意义。
随着国家生态文明建设实践的深入推进，许多理论研究者也围绕生态空间识别、功能演变、集约利用开展了大量有益探索。在生态空间识别方面，大量研究集中在方法层面，主要包括基于土地利用/覆被类型直接划定生态空间的定性研究和综合生态系统服务功能重要性和生态环境敏感性划定生态空间的定量研究，而且这些研究多以静态为主，很少有考虑到生态系统服务功能的动态演变特征。尽管国土空间集约利用优化是一个复杂而具有挑战性的资源空间优化分配问题，但仍然有许多研究者就该主题相关领域进行了广泛研究，出现了大量可用于生态空间集约利用优化的数学模型，主要包括数量结构优化模型、空间布局优化模型和集成优化模型。
总体而言，集成优化模型为解决复杂空间优化决策问题提供了有用的解决方案，亦为国土空间集约利用优化提供了方法参考与借鉴。但是，当前国土空间集约利用优化研究仍然存在一些需要解决的问题：(1)生态空间识别多以当前生态系统结构和服务功能信息为参照，而忽略了功能结构演变带来的生态空间时空异质性；(2)生态空间演变分析多集中于数量与空间变化，而忽略了生态空间主导生态系统服务功能演变特征研究；(3)国土空间优化多集中在空间布局优化上，常忽略数量结构优化与空间布局优化的耦合，所得方案虽能实现国土资源空间配置上的优化，但不能确保实现空间利用效益的最大化；(4)在整合元胞自动机与其他优化方法进行空间优化时，大多采用MC等模型预测所得空间需求作为数量转化规则，未能真正实现目标年国土空间集约利用优化，亦即同时实现目标年空间布局最优化和空间利用效益的最大化。因此，在生态空间动态识别和主导生态系统服务功能演变分析的基础上，开展生态空间集约利用优化研究显得非常有必要。
发明内容
针对上述问题，本发明目的在于提供一种生态空间集约利用优化方法及系统，本发明利用土地利用、卫星影像等数据，构建LP-Markov-CA模型进行生态空间集约利用优化研究，其目的在于实现生态空间数量结构、空间布局与服务功能的协同优化，获取生态系统服务功能最大的生态空间集约利用布局，为国土空间集约利用优化配置提供新方法。
本发明通过下述技术方案实现：
一方面，本发明提供了一种生态空间集约利用优化方法，该方法包括：
S1：根据多年生态系统服务功能综合值及其变化特征，动态获取生态空间；所述生态系统服务功能类型包括初级产品生产、气体调节、气候调节、环境净化、水文调节、土壤保持、生物多样性、美学景观功能；
S2：根据各单项生态系统服务功能值，识别生态空间主导生态系统服务功能；模拟基期年至目标年生态空间主导生态系统服务功能格局，进行主导生态系统服务功能演变分析；
S3：采用LP-Markov-CA复合模型，进行目标年主导生态系统服务功能数量结构和空间布局耦合优化，实现生态空间主导生态系统服务功能的集约利用优化。
本发明采用了上述技术方案，与现有技术相比，本发明实现了主导生态系统服务功能数量结构、空间布局的耦合优化，成功获取了生态系统服务功能最大的生态空间集约利用优化布局，为国土空间集约利用优化配置提供了新方法、新技术。
进一步地，步骤S1具体包括以下子步骤：
S11：建立生态空间生态系统服务功能评价单元，其中以网格为生态空间生态系统服务功能评价单元，其大小与评价尺度密切相关；
S12：计算评价单元的生态系统服务功能价值：获取生态系统服务价值当量(根据已知当量表查询得到)，将生态系统服务价值当量作为生态系统服务功能大小的度量衡，计算评价单元的生态系统服务功能价值；采用式(1)计算评价单元的单项生态系统服务功能值，根据所述单项生态系统服务功能值，采用式(2)汇总得到评价单元的生态系统服务功能综合值；
Figure GDA0003567359330000021
Figure GDA0003567359330000022
式中：υlj表示第l类生态系统第j种生态系统服务功能当量；sk′l表示第k′个评价单元第l类生态系统面积；Vk′j,Vk′分别表示第k′个评价单元第j种生态系统服务功能值和功能综合值；M,n,K分别表示生态系统类型数、生态系统服务功能类型数和评价单元个数；l,j,k′分别表示生态系统类型、生态系统服务功能类型和评价单元编号，l＝1,2,…,M；j＝1,2,…,n；k′＝1,2,…,K；
S13：识别生态空间：在长时序生态系统服务综合功能变化幅度分析基础上，结合其功能综合值大小合理划定生态空间，即将生态系统服务功能综合值大于等于某临界值且其平均变化率小于某临界值的单元划为生态空间，其数学表达式为：
Figure GDA0003567359330000031
式中：E表示生态空间；ek表示组成生态空间E的矢量网格(评价单元)；
Figure GDA0003567359330000032
表示t年单元k的生态系统服务功能综合值，计算方法见式(2)；k,m分别表示组成生态空间E的网格编号和个数；t,T分别表示任一时间断面和时间断面个数；δ,ε分别表示划分生态空间与非生态空间的生态系统服务功能综合值、生态系统服务综合功能值年均变化率临界值，采用标准差分级法分析确定，将t年生态系统服务功能综合值不小于δ且年均变化率不大于ε的网格识别为生态空间，否则为非生态空间。
进一步地，步骤S2具体包括以下子步骤：
S21：根据各单项生态系统服务功能值，识别生态空间主导生态系统服务功能：以矢量网格为识别单元，将网格单元对应的功能值最大的生态系统服务类型作为该矢量网格的主导功能，各网格单元主导功能集合即为生态空间主导功能格局，其数学表达式为：
Figure GDA0003567359330000033
式中：D表示生态空间主导生态系统服务功能格局；
Figure GDA0003567359330000034
表示第k个网格的主导生态系统服务功能，其编号为jk∈[1,n]，编码值1、2、3、4、5、6、7、8分别代表初级产品生产、气体调节、气候调节、环境净化、水文调节、土壤保持、生物多样性、美学景观功能；
Figure GDA0003567359330000035
表示第k个网格主导功能jk的生态系统服务功能值；Vkj表示第k个网格第j类生态系统服务功能值。
S22：进行生态空间主导生态系统服务功能格局模拟：应用Markov-CA模型模拟基期年至目标年主导生态空间生态系统服务功能格局；
S23：进行生态空间主导生态系统服务功能演变分析：基于T个时间断面主导生态空间生态系统服务功能分布图，统计得到t年j项主导生态系统服务功能面积atj，分析得到主导生态空间生态系统服务功能区最大面积max1≤t≤Tatj和最小面积min1≤t≤Tatj。
进一步地，步骤S3中的所述LP-Markov-CA复合模型的数学表达式为：
U＝(L,S,N,f) (5)
式中：U表示LP-Markov-CA模型框架；L表示元胞空间；S表示元胞状态集合；N表示领域内元胞状态空间向量；f表示元胞转换规则。
利用LP-Markov-CA模型进行生态空间集约利用优化的技术关键是制定主导功能面积转移矩阵(元胞数量转化规则)、主导功能适宜性图集(元胞空间转化规则)和主导功能初始分布(元胞空间初始状态)。
进一步地，步骤S3具体包括以下子步骤：
S31，将元胞C定义为主导生态系统服务功能分布图的栅格，其大小为栅格图分辨率(如：100m)；
S32，将元胞空间L定义为生态空间E；
S33，将元胞状态S定义为初级产品生产、气体调节、气候调节、环境净化、水文调节、水土保持、生物多样性维持、美学景观8种生态系统服务功能，各功能依次编码为1,2,3,4,5,6,7,8，则S＝{1,2,3,4,5,6,7,8}；
S34，确定元胞领域，采用标准5×5邻近滤波器；
S35，计算元胞转换面积矩阵；
S36，制作元胞转换适宜性图集，将n种经标准化处理的基期年生态空间单项生态系统服务功能值分布栅格图，按主导功能类型编码顺序组合成适宜性图集，作为元胞空间转化规则；
S37，将调整后的基期年主导生态系统服务功能分布图作为元胞空间初始状态，将步骤S35确定的元胞转换面积矩阵和步骤S36生成的元胞转换适宜性图集分别作为元胞数量和空间转换规则，利用Markov-CA模型即可模拟得到目标年生态空间集约利用优化布局。
进一步地，步骤S35计算元胞转换面积矩阵的具体步骤为：
S351：进行目标年生态空间主导生态系统服务功能数量结构优化；
以初级产品生产、气体调节、气候调节、环境净化、水文调节、土壤保持、生物多样性、美学景观8类主导生态系统服务功能面积为决策变量，以主导生态系统服务功能面积条件为约束，以在有限生态空间实现最大生态系统服务功能效用为目标，建立线性规划模型，对目标年生态空间范围内不同主导生态系统服务功能面积进行优化；模型表达式为：
Figure GDA0003567359330000041
式中：f(a)为模型优化目标函数，表示生态空间各主导生态系统服务功能总值(当量)；υtj表示t年主导生态系统服务功能j单位面积功能值(当量/hm2)；aj表示j类主导生态系统服务功能面积(hm2)；aE表示生态空间E总面积(hm2)；atj表示t年主导生态系统服务功能j的面积(hm2)；j，n分别表示生态空间主导生态系统服务功能类型编号和类型数，其中j＝1，2，3，4，5，6，7，8，n＝8；t，T分别表示时间断面编号和时间断面数；
S352：进行基期年至目标年生态空间主导生态系统服务功能转移概率矩阵预测；
利用Markov模型计算得到1，2，3，…，t，…，T-1逐时段主导生态系统服务功能转移概率方阵(pij)n×n(i，j＝1，2，…，n)；转移概率不但会随时段发生变化，而且变化率非常量；假设转移概率pij是随时间t连续变化的函数pij(t)，单位时间变化率r也是随时间变化的函数rij(t)，则可得到描述转移概率变化规律的微分方程及其定解条件：
Figure GDA0003567359330000051
求解定解问题，可得到：
Figure GDA0003567359330000052
式中：B为待定常数，初值为(pij)0；
利用公式(8)分别计算基期年(第t年)至目标年逐时段主导生态系统服务功能转移概率
Figure GDA0003567359330000053
在此基础上采用式(9)计算得到基期年至目标年逐时段转移概率矩阵
Figure GDA0003567359330000054
Figure GDA0003567359330000055
式中：
Figure GDA0003567359330000056
表示第t′个时段n类主导生态系统服务功能转移概率之和；Pt′＝(pij)n×n，表示t时段主导生态系统服务功能转移概率方阵；
最后采用式(10)计算确定基期年至目标年生态空间主导生态系统服务功能转移概率矩阵P；
Figure GDA0003567359330000057
S353：计算基期年至目标年生态空间主导生态系统服务功能转移面积矩阵；
设基期年、目标年主导生态系统服务功能优化面积矩阵分别为
Figure GDA0003567359330000058
则基期年至目标年主导生态系统服务功能转移面积矩阵A＝(aij)n×n为：
A＝diag(x1，x2，…，xn)P (11)
式中：x1，x2，…，xn表示基期年主导生态系统服务功能优化面积矩阵x的元素，其计算公式为：
x＝[PT]-1A*(|PT|≠0) (12)。
进一步地，步骤S36制作元胞转换适宜性图集的具体步骤为：
生态系统服务功能值越大说明越适宜布局该类生态系统服务功能区，反之布局适宜度则越低。为此，将基期年生态空间单项生态系统服务功能值分布栅格图，采用公式(13)标准化后作为该项功能空间分布适宜性图，按主导功能类型编码顺序组合成适宜性图集；并以此为依据对基期年主导生态系统服务功能实际分布进行调整，得到数量结构与面积矩阵x一致的空间分布，如下：
Figure GDA0003567359330000061
式中：Vkj，V′kj表示第j种生态系统服务功能在评价单元(栅格)k的生态系统服务功能值和标准化值；k，m分别表示组成生态空间E的网格编号和个数。
另一方面，本发明还提供了一种生态空间集约利用优化系统，该系统支持所述的一种生态空间集约利用优化方法，该系统包括：
生态空间动态获取模块，用于根据多年生态系统服务功能综合值及其变化特征，动态获取生态空间；所述生态系统服务功能类型包括初级产品生产、气体调节、气候调节、环境净化、水文调节、土壤保持、生物多样性、美学景观功能；
主导生态系统服务功能模拟及演变模块：用于根据各单项生态系统服务功能值，识别生态空间主导生态系统服务功能；模拟基期年至目标年生态空间主导生态系统服务功能格局，进行主导生态系统服务功能演变分析；
集约利用优化模块，用于采用LP-Markov-CA复合模型，进行目标年主导生态系统服务功能数量结构和空间布局耦合优化，实现生态空间主导生态系统服务功能的集约利用优化，得到空间利用效益最大化的生态空间集约利用优化布局；
输出模块，用于输出所述空间利用效益最大化的生态空间集约利用优化布局。
进一步地，所述集约利用优化模块中的所述LP-Markov-CA复合模型的数学表达式为：
U＝(L，S，N，f) (5)
式中：U表示LP-Markov-CA模型框架；L表示元胞空间；S表示元胞状态集合；N表示领域内元胞状态空间向量；f表示元胞转换规则。
利用LP-Markov-CA模型进行生态空间集约利用优化的技术关键是制定主导功能面积转移矩阵(元胞数量转化规则)、主导功能适宜性图集(元胞空间转化规则)和主导功能初始分布(元胞空间初始状态)。
进一步地，所述集约利用优化模块的执行过程如下：
将元胞C定义为主导生态系统服务功能分布图的栅格，其大小为栅格图分辨率(如：100m)；
将元胞空间L定义为生态空间E；
将元胞状态S定义为初级产品生产、气体调节、气候调节、环境净化、水文调节、水土保持、生物多样性维持、美学景观8种生态系统服务功能，各功能依次编码为1,2,3,4,5,6,7,8，则S＝{1,2,3,4,5,6,7,8}；
确定元胞领域，采用标准5×5邻近滤波器；
计算元胞转换面积矩阵；
制作元胞转换适宜性图集，将n种经标准化处理的基期年生态空间单项生态系统服务功能值分布栅格图，按主导功能类型编码顺序组合成适宜性图集，作为元胞空间转化规则；
将调整后的基期年主导生态系统服务功能分布图作为元胞空间初始状态，将确定的元胞转换面积矩阵和生成的元胞转换适宜性图集分别作为元胞数量和空间转换规则，利用Markov-CA模型即可模拟得到目标年生态空间集约利用优化布局。
本发明与现有技术相比，具有如下的优点和创新：
1、本发明方法及系统提出立足于生态系统服务功能动态性的生态空间识别方法，增强了生态空间划定的合理性和可靠性。
2、本发明方法及系统进一步证实了主导生态系统服务功能的非平稳性，并在其数量结构优化模型中以不等式约束方式表征了生态空间各主导生态系统服务功能的动态性，增强了目标年生态空间数量结构优化的合理性。
3、本发明方法及系统实现了目标年生态空间主导生态系统服务功能最优数量结构与最佳空间布局的完美耦合，在实现生态空间布局优化的同时实现了数量结构的优化配置，获得了空间利用效益最大化的生态空间集约利用优化布局。
4、本发明方法及系统关键技术创新是，基于马尔科夫过程的非时齐性和线性规划模型获取的最优主导生态系统服务功能数量结构，制定了基期年至目标年主导生态系统服务功能转移概率矩阵和转移面积矩阵，借助Markov-CA模型框架完美获取了目标年空间利用效益最大化的生态空间集约利用布局，具体为：
(1)现行的应用Markov模型进行空间变化模拟的技术普遍基于马尔科夫过程时齐性，将过去某时段的转移概率作为未来时段的转移概率，确定基期年到目标年的转移面积矩阵，这类技术主要存在两点缺陷：一是本发明发现生态系统服务功能转移概率具备明显非时齐性，亦即各时段转移概率会随时间变化，直接应用现行技术难以准确制定未来一定时期主导生态系统服务功能转移概率矩阵；二是现行模拟技术得到的目标年主导生态系统服务功能最优布局对应的数量结构往往与预期数量结构(如：利用线性规划模型制定的主导生态系统服务功能最优数量结构)不一致，实际上所得到的空间布局并非为空间利用效益最大的空间集约利用优化布局，无法应用服务于国土空间规划；
(2)本发明方法及系统根据过去一定时期生态空间逐时段主导生态系统服务功能转移概率，假设转移概率及其变化率为连续变化函数，建立微分方程模型确定基期年至目标年主导生态系统服务功能转移概率矩阵，解决了现行技术未考虑马尔科夫过程非时齐性不能直接应用于主导生态系统服务功能转移概率矩阵制定的技术缺陷；
(3)本发明方法及系统根据马尔科夫转移概率矩阵特性，建立基期年主导生态系统服务功能优化面积、基期年至目标年主导生态系统服务功能转移概率矩阵与面积矩阵之间的非齐次线性方程组，求解出基期年主导生态系统服务功能优化面积，然后根据基期年主导生态系统服务功能空间分布适宜性，对其实际分布进行调整获得生态空间主导生态系统服务功能初始分布，在此基础上借助Markov-CA模型即可获取与目标年主导生态系统服务功能最优数量结构完美匹配的最优空间布局，解决了现行技术难以实现空间利用效益最大化的生态空间集约利用布局的技术缺陷。
附图说明
此处所说明的附图用来提供对本发明实施例的进一步理解，构成本申请的一部分，并不构成对本发明实施例的限定。在附图中：
图1为本发明一种生态空间集约利用优化方法流程图。
图2为本发明一种生态空间集约利用优化方法的框架图。
图3为本发明实施例2003-2019年生态系统服务综合功能值标准差分级图。
图4为本发明实施例生态空间识别过程图。
图5为本发明实施例生态空间(综合功能值≥6.79且其年平均变化率≤5％)识别结果图。
图6为本发明实施例主导生态系统服务功能适宜性空间分布图。
图7为本发明实施例2025年主导生态系统服务功能空间优化初始分布与生态空间集约利用优化布局图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚明白，下面结合实施例和附图，对本发明作进一步的详细说明，本发明的示意性实施方式及其说明仅用于解释本发明，并不作为对本发明的限定。
实施例1
如图1至图7所示，本发明一种生态空间集约利用优化方法，如图1所示，该方法包括：
S1：根据多年生态系统服务功能综合值及其变化特征，动态获取生态空间；所述生态系统服务功能类型包括初级产品生产、气体调节、气候调节、环境净化、水文调节、土壤保持、生物多样性、美学景观功能；
S2：根据各单项生态系统服务功能值，识别生态空间主导生态系统服务功能；模拟基期年至目标年生态空间主导生态系统服务功能格局，进行主导生态系统服务功能演变分析；
S3：采用LP-Markov-CA复合模型，进行目标年主导生态系统服务功能数量结构和空间布局耦合优化，实现生态空间主导生态系统服务功能的集约利用优化。
本发明采用了上述技术方案，与现有技术相比，本发明实现了主导生态系统服务功能数量结构、空间布局的耦合优化，成功获取了生态系统服务功能最大的生态空间集约利用优化布局，为国土空间集约利用优化配置提供了新方法、新技术。
本发明最主要的创新点在于基于马尔科夫过程的非时齐性和线性规划模型获取的最优主导生态系统服务功能数量结构，制定了基期年至目标年主导生态系统服务功能转移概率矩阵和转移面积矩阵，借助Markov-CA模型框架完美获取了目标年空间利用效益最大化的生态空间集约利用布局。
本发明一种生态空间集约利用优化方法及系统的框架图如图2所示，本实施例以四川省邛崃市为例，以2025年为时间节点，根据上述方法进行生态空间生态系统服务功能集约利用优化模拟。
具体实施如下：
S1：根据多年生态系统服务功能综合值及其变化特征，动态获取生态空间；所述生态系统服务功能类型包括初级产品生产、气体调节、气候调节、环境净化、水文调节、土壤保持、生物多样性、美学景观功能；具体如下：
综合考虑研究区范围大小和数据运算量，确定评价单元网格大小为100m×100m，应用ESRI ArcGIS软件Create Fishnet工具构建生态空间生态系统服务功能评价单元网格矢量数据，研究区共涉及139738个单元网格。
利用ESRI ArcGIS软件TabulateArea工具统计评价单元内生态系统类型面积，采用式(1)计算评价单元初级产品生产、气体调节等单项生态系统服务功能值，最后采用式(2)汇总得到评价单元生态系统服务功能综合值。
Figure GDA0003567359330000091
Figure GDA0003567359330000092
式中：υlj表示第l类生态系统第j种生态系统服务功能当量；sk′l表示第k′个评价单元第l类生态系统面积；Vk′j，Vk′分别表示第k′个评价单元第j种生态系统服务功能值和功能综合值；M，n，K分别表示生态系统类型数、生态系统服务功能类型数和评价单元个数；l，j，k′分别表示生态系统类型、生态系统服务功能类型和评价单元编号，l＝1，2，...，M；j＝1，2，...，n；k′＝1，2，...，K；
将生态系统服务功能综合值大于等于某临界值且其平均变化率小于某临界值的单元划为生态空间，其数学表达式为：
Figure GDA0003567359330000101
式中：E表示生态空间；ek表示组成生态空间E的矢量网格(评价单元)；
Figure GDA0003567359330000102
表示t年单元k的生态系统服务功能综合值，计算方法见式(2)；k，m分别表示组成生态空间E的网格编号和个数；t，T分别表示任一时间断面和时间断面个数；δ，ε分别表示划分生态空间与非生态空间的生态系统服务综合功能值、生态系统服务综合功能值年均变化率临界值，采用标准差分级法分析确定，将t年生态系统服务功能综合值不小于δ且年均变化率不大于ε的网格识别为生态空间，否则为非生态空间。
研究测算了邛崃市2003、2005、2007、2009、2011、2013、2015、2017、2019年共9个时间断面各网格生态系统服务功能综合值，采用标准差分级法对生态系统服务综合功能值进行分级，结果显示各年份生态系统服务综合功能值分级临界值表现出较强的一致性(图2)。其中，第一个分级临界值分别为6.88、6.84、6.86、6.96、6.84、6.78、6.67、6.66、6.65，标准差为0.1042，表明该临界值波动不大，比较稳定。因此，将各年生态系统服务综合功能值的第一个分级临界值的均值6.79作为邛崃市划分生态区和非生态区的分界值。
但是，各年生态系统服务综合功能值大于6.79的空间面积波动较大，呈现先增后减的变化趋势(图4-a)。由此可见，在划定生态空间时，还应当考虑生态系统服务综合功能值的变化率特征，以确保识别出的生态空间稳定可靠。为此，应用标准差分级法对生态系统服务综合功能值年均变化率进行分级，结果显示生态系统服务综合功能值年均变化率共分为4级，分级临界值分别为5％、13％、22％、34.4％(图4-b)。理论上，生态系统服务综合功能值未变(即年均变化率为0)的区域为稳定生态空间，但由于生态系统本身的复杂性和不可避免的数据分析误差，实际中很少有数值绝对不变的区域，所以将第一个生态系统服务综合功能值年均变化率分级临界值5％作为邛崃市划分稳定生态空间和非稳定生态空间的分界值。
综合生态系统服务功能及其变化率对生态空间的影响，将逐年生态系统服务综合功能值大于6.79(图4-c)且年均变化率小于5％(图4-d)的区域识别为生态空间(图5)。
S2：根据各单项生态系统服务功能值，识别生态空间主导生态系统服务功能；模拟基期年至目标年生态空间主导生态系统服务功能格局，进行主导生态系统服务功能演变分析；具体如下：
S21：根据各单项生态系统服务功能值，识别生态空间主导生态系统服务功能：以矢量网格为识别单元，将网格单元对应的功能值最大的生态系统服务类型作为该矢量网格的主导功能，各网格单元主导功能集合即为生态空间主导功能格局，其数学表达式为：
Figure GDA0003567359330000111
式中：D表示生态空间主导生态系统服务功能格局；
Figure GDA0003567359330000112
表示第k个网格的主导生态系统服务功能，其编号为jk∈[1,n]，编码值1、2、3、4、5、6、7、8分别代表初级产品生产、气体调节、气候调节、环境净化、水文调节、土壤保持、生物多样性、美学景观功能；
Figure GDA0003567359330000113
表示第k个网格主导功能jk的生态系统服务功能值；Vkj表示第k个网格第j类生态系统服务功能值。
S22：进行生态空间主导生态系统服务功能格局模拟：应用Markov-CA模型模拟2021、2023、2025年主导生态空间生态系统服务功能格局；主要步骤为：
第1步，定义模型参数。将元胞定义为主导生态系统服务功能分布图的栅格(大小为栅格图分辨率100m)，元胞空间定义为主导生态系统服务功能空间格局，元胞状态定义为初级产品生产、气体调节、气候调节、环境净化、水文调节、水土保持、生物多样性维持、美学景观8种生态系统服务功能类型；采用标准5×5邻近滤波器。
第2步，检验模型精度。首先，根据生态空间主导生态系统服务功能演变特征和模拟数据需求，选择2015、2017年，2011、2015年，2007、2013年主导生态系统服务功能分布栅格图，应用IDRISI软件Markov模块分别计算得到2015-2017年、2011-2015年、2007-2013年主导生态系统服务功能转移概率矩阵和条件概率图，分别将其作为CA模型元胞数量和空间转换规则；然后，基于第1步的参数定义，利用IDRISI软件Markov-CA模块分别模拟得到2019年主导生态系统服务功能分布栅格图；最后，分别计算模拟得到的2019年主导生态系统服务功能分布图与实际分布图的Kappa系数，评测模型精度，若模拟分布图与实际分布图一致性程度达到适中及其以上，则可用其进行相应空间变化模拟，否则应调整模型参数或模拟数据，重新进行精度检验至到符合要求。
从Kappa系数计算结果(表1)可知，时间梯度为2、4、6年的Kappa系数分别为0.9737、0.9730、0.9689，表明随梯度增加模型模拟结果与实际分布的一致性程度在缓慢减弱，但都达到了最佳水平，说明模型具有很强的可信度，可用其进行生态空间主导生态系统服务功能空间变化模拟。
表1 Markov-CA模型模拟精度检验结果
Figure GDA0003567359330000121
第3步，应用模型模拟。若模型精度评测合格，则可参照精度评测数据时间间隔，选择2017、2019年，2015、2019年，2013、2019年主导生态系统服务功能分布栅格图，参照(2)中相应步骤方法构建元胞数量和空间转化规则，设置模型参数，利用IDRISI软件Markov-CA模块模拟得到2021、2023、2025年生态空间主导生态系统服务功能格局图。
S23：进行生态空间主导生态系统服务功能演变分析：基于12个时间断面(2003、2005、2007、2009、2011、2013、2015、2017、2019、2021、2023、2025年)生态空间主导生态系统服务功能分布图，利用ESRI ArcGIS统计得到t年j项主导生态系统服务功能面积atj，分析得到2003～2025年主导生态系统服务功能区最大面积max1≤t≤12atj和最小面积min1≤t≤12atj。
2003-2025年间，气候调节功能面积减少量最大，初级产品生产功能面积增加量最大，水文调节功能面积略微减少。其中，气候调节功能区的面积从2009年的最大面积55444hm2降至2023年的52356hm2；初级产品生产功能面积总体呈现上升的趋势，从2009年的29297hm2增加至2025年的32793hm2；水文调节功能面积变化不明显，在2003-2025年间的最大面积为13650hm2、最小面积为13024hm2。2003-2025年生态空间主导生态系统服务功能面积见表2：
表2 2003-2025年生态空间主导生态系统服务功能面积统计
年份 初级产品生产(hm<sup>2</sup>) 气候调节功能(hm<sup>2</sup>) 水文调节功能(hm<sup>2</sup>)
2003 29327 55416 13564
2005 29302 55440 13565
2007 29342 55418 13547
2009 29297 55444 13566
2011 29463 55287 13557
2013 29747 55023 13537
2015 29882 54907 13518
2017 29906 54874 13527
2019 29937 54860 13510
2021 32221 52512 13574
2023 32301 52356 13650
2025 32793 52490 13024
S3：采用LP-Markov-CA复合模型，进行目标年主导生态系统服务功能数量结构和空间布局耦合优化，实现生态空间主导生态系统服务功能的集约利用优化。
所述LP-Markov-CA复合模型的数学表达式为：
U＝(L，S，N，f) (5)
式中：U表示LP-Markov-CA模型框架；L表示元胞空间；S表示元胞状态集合；N表示领域内元胞状态空间向量；f表示元胞转换规则。
利用LP-Markov-CA模型进行生态空间集约利用优化的技术关键是制定主导功能面积转移矩阵(元胞数量转化规则)、主导功能适宜性图集(元胞空间转化规则)和主导功能初始分布(元胞空间初始状态)。
步骤S3具体包括以下子步骤：
S31，将元胞C定义为主导生态系统服务功能分布图的栅格，其大小为栅格图分辨率(如：100m)；
S32，将元胞空间L定义为生态空间E；
S33，将元胞状态S定义为初级产品生产、气体调节、气候调节、环境净化、水文调节、水土保持、生物多样性维持、美学景观8种生态系统服务功能，各功能依次编码为1，2，3，4，5，6，7，8，则S＝{1，2，3，4，5，6，7，8}；
S34，确定元胞领域，采用标准5×5邻近滤波器；
S35，计算元胞转换面积矩阵；首先，利用MATLAB约束最优化问题求解函数FMINCON，求解主导生态系统服务功能数量结构优化线性规划模型，得到2025年初级产品生产、气候调节、水文调节等主导生态系统服务功能区优化面积
Figure GDA0003567359330000131
然后，基于2003、2005、2007、2009、2011、2013、2015、2017、2019年主导生态系统服务功能图，利用IDRISI软件Markov模块计算得到0，1，2，…，t，…，T-1等时段各主导生态系统服务功能间转移概率pij，在此基础上利用MATLAB拟合得到函数pij(t)解析式，利用该函数解析式即可得到2019-2021年、2021-2023年、2023-2025年三个时段的主导生态系统服务功能转移概率
Figure GDA0003567359330000132
最后，利用公式(9-12)计算出2019-2025年主导生态系统服务功能转移面积矩阵Aij，将其作为元胞转换面积矩阵。
S36，制作元胞转换适宜性图集。利用IDRISI软件Collection Editor工具将n种经标准化处理的2019年生态空间单项生态系统服务功能值分布栅格图，按主导功能类型编码顺序组合成适宜性图集，作为元胞空间转化规则。图6为主导生态系统服务功能适宜性空间分布图。
S37，将调整后的2019年主导生态系统服务功能分布图作为元胞空间初始状态，将步骤S35确定的元胞转换面积矩阵和步骤S36生成的元胞转换适宜性图集分别作为元胞数量和空间转换规则，设置循环次数为6，利用IDRISI软件Markov-CA模型即可模拟得到2025年生态空间集约利用优化布局。图7为2025年主导生态系统服务功能空间优化初始分布与生态空间集约利用优化布局图。
具体地，S35中计算元胞面积转换矩阵的具体步骤为：
(1)2025年生态空间主导生态系统服务功能数量结构优化。
以初级产品生产、气体调节、气候调节、环境净化、水文调节、土壤保持、生物多样性、美学景观8类主导生态系统服务功能面积为决策变量，以主导生态系统服务功能面积条件为约束，以在有限生态空间实现最大生态系统服务功能效用为目标，建立线性规划模型，对2025年生态空间范围内不同主导生态系统服务功能面积进行优化。模型数学表达式为：
Figure GDA0003567359330000141
式中：f(a)为模型优化目标函数，表示生态空间各主导生态系统服务功能总值(当量)；υtj表示t年主导生态系统服务功能j单位面积功能值(当量/hm2)；aj表示j类主导生态系统服务功能面积(hm2)；aE表示生态空间E总面积(hm2)；atj表示t年主导生态系统服务功能j的面积(hm2)；j，n分别表示生态空间主导生态系统服务功能类型编号和类型数，其中j＝1，2，3，4，5，6，7，8，n＝8；t，T分别表示时间断面编号和时间断面数，其中t＝2003，2005，2007，2009，2011，2013，2015，2017，2019，2021，2023，2025，T＝12。
(2)2019-2025年生态空间主导生态系统服务功能转移概率矩阵预测。
利用Markov模型计算得到第0(2003-2005)、1(2005-2007)、2(2007-2009)、3(2009-2011)、4(2011-2013)、5(2013-2015)、6(2015-2017)、7(2017-2019)个时段的逐时段主导生态系统服务功能转移概率方阵(pij)n×n(i，j＝1，2，…，n)。转移概率不但会随时段发生变化，而且变化率非常量。假设转移概率pij是随时间t连续变化的函数pij(t)，单位时间变化率r也是随时间变化的函数rij(t)，则可得到描述转移概率变化规律的微分方程及其定解条件：
Figure GDA0003567359330000142
求解定解问题，可得：
Figure GDA0003567359330000143
式中：B为待定常数，初值为(pij)0；
主导生态系统服务功能转移概率预测模型如表3：
表3主导生态系统服务功能转移概率预测模型
Figure GDA0003567359330000151
利用公式(8)分别计算第8个时段(2019-2021年)、第9个时段(2021-2023年)和第10个时段(2023-2025年)主导生态系统服务功能转移概率
Figure GDA0003567359330000152
在此基础上按式(9)计算得到2019-2025年间逐时段转移概率矩阵
Figure GDA0003567359330000153
Figure GDA0003567359330000154
式中：
Figure GDA0003567359330000155
表示第t个时段n类主导生态系统服务功能转移概率之和；Pt′＝(pij)n×n，表示t时段主导生态系统服务功能转移概率方阵。
最后按式(10)计算确定基期年至目标年生态空间主导生态系统服务功能转移概率矩阵P。
Figure GDA0003567359330000156
(3)2019-2025年生态空间主导生态系统服务功能转移面积矩阵计算。
设2019、2025年主导生态系统服务功能优化面积矩阵分别为
Figure GDA0003567359330000157
则2019-2025年主导生态系统服务功能转移面积矩阵A＝(aij)n×n为：
A＝diag(x1，x2，…，xn)P (11)
式中：x1，x2，…，xn表示2019年主导生态系统服务功能优化面积矩阵x的元素，其计算公式为：
x＝[PT]-1A*(|PT|≠0) (12)
计算出的2019-2025主导生态系统服务功能转移概率矩阵和转移面积矩阵如表4：
表4 2019-2025主导生态系统服务功能转移概率矩阵和转移面积矩阵表
Figure GDA0003567359330000161
具体地，S36中制作元胞转换适宜性图集的具体步骤为：
生态系统服务功能值越大说明越适宜布局该类生态系统服务功能区，反之布局适宜度则越低。为此，将2019年生态空间单项生态系统服务功能值分布栅格图，按公式(13)标准化后作为该项功能空间分布适宜性图，并以此为依据对2019年主导生态系统服务功能实际分布进行调整，得到数量结构与面积矩阵x一致的空间分布。
Figure GDA0003567359330000162
式中：Vkj,V′kj表示第j种生态系统服务功能在评价单元(栅格)k的生态系统服务功能值和标准化值；k,m分别表示组成生态空间E的网格编号和个数。
具体地，利用Markov-CA模型即可模拟得到2025年生态空间集约利用优化布局，具体为：
将调整后的2019年主导生态系统服务功能分布图作为元胞空间初始状态，将元胞转换面积矩阵和元胞转换适宜性图集分别作为元胞数量和空间转换规则，设置循环次数为6，利用Markov-CA模型即可模拟得到目标年生态空间集约利用优化。
利用Markov-CA模型模拟得到的生态空间集约利用优化的数量结构与LP模型确定的生态空间主导生态系统服务功能最优数量结构基本一致(表5)。其中，水文调节功能空间集约利用优化布局误差最大，为0.04％，初级产品生产和气候调节功能空间集约利用优化误差最小，为-0.01％，这主要是由于模型运行过程中小数位数保留所致，表明LP-Markov-CA复合模型很好地解决了生态空间数量结构优化与空间布局的耦合优化问题，同时实现了目标年生态空间布局最优化和空间利用效益最大化，真正实现了目标年生态空间集约利用优化。
2025年生态空间集约利用优化数量结构与主导生态系统服务功能最优数量结构如表5：
表5 2025年生态空间集约利用优化数量结构与主导生态系统服务功能最优数量结构对比表
Figure GDA0003567359330000171
由此可知，LP-Markov-CA模型所得优化布局数量结构与LP制定的最优数量结构基本一致，表明LP-Markov-CA复合模型很好地解决了生态空间数量结构优化与空间布局的耦合优化问题，同时实现了目标年生态空间布局最优化和空间利用效益最大化，真正实现了目标年生态空间集约利用优化。
实施例2
如图1至图7所示，本实施例与实施例1的区别在于，本实施例提供了一种生态空间集约利用优化系统，该系统支持实施例1所述的一种生态空间集约利用优化方法，该系统包括：
生态空间动态识别单元，用于根据多年生态系统服务功能综合值及其变化特征，动态获取生态空间；所述生态系统服务功能类型包括初级产品生产、气体调节、气候调节、环境净化、水文调节、土壤保持、生物多样性、美学景观功能；
主导生态系统服务功能模拟及演变单元：用于根据各单项生态系统服务功能值，识别生态空间主导生态系统服务功能；模拟基期年至目标年生态空间主导生态系统服务功能格局，进行主导生态系统服务功能演变分析；
生态空间集约利用优化单元，用于采用LP-Markov-CA复合模型，进行目标年主导生态系统服务功能数量结构和空间布局耦合优化，实现生态空间主导生态系统服务功能的集约利用优化，得到空间利用效益最大化的生态空间集约利用优化布局；
输出单元，用于输出空间利用效益最大化的生态空间集约利用优化布局。
进一步地，所述生态空间集约利用优化单元中所述LP-Markov-CA复合模型的数学表达式为：
U＝(L,S,N,f) (5)
式中：U表示LP-Markov-CA模型框架；L表示元胞空间；S表示元胞状态集合；N表示领域内元胞状态空间向量；f表示元胞转换规则。
利用LP-Markov-CA模型进行生态空间集约利用优化的技术关键是制定主导功能面积转移矩阵(元胞数量转化规则)、主导功能适宜性图集(元胞空间转化规则)和主导功能初始分布(元胞空间初始状态)。
具体实施时，所述集约利用优化单元的执行过程如下：
将元胞C定义为主导生态系统服务功能分布图的栅格，其大小为栅格图分辨率(如：100m)；
将元胞空间L定义为生态空间E；
将元胞状态S定义为初级产品生产、气体调节、气候调节、环境净化、水文调节、水土保持、生物多样性维持、美学景观8种生态系统服务功能，各功能依次编码为1,2,3,4,5,6,7,8，则S＝{1,2,3,4,5,6,7,8}；
确定元胞领域，采用标准5×5邻近滤波器；
计算元胞转换面积矩阵；
制作元胞转换适宜性图集，将n种经标准化处理的基期年生态空间单项生态系统服务功能值分布栅格图，按主导功能类型编码顺序组合成适宜性图集，作为元胞空间转化规则；
将调整后的基期年主导生态系统服务功能分布图作为元胞空间初始状态，将确定的元胞转换面积矩阵和生成的元胞转换适宜性图集分别作为元胞数量和空间转换规则，利用Markov-CA模型即可模拟得到目标年生态空间集约利用优化布局。
其他各个模块的执行过程按照实施例1中的一种生态空间集约利用优化方法的相应流程步骤执行即可，此实施例中不再一一赘述。
本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种生态空间集约利用优化方法，其特征在于，该方法包括：
S1：根据多年生态系统服务功能综合值及其变化特征，动态获取生态空间；所述生态系统服务功能类型包括初级产品生产、气体调节、气候调节、环境净化、水文调节、土壤保持、生物多样性、美学景观功能；
S2：根据各单项生态系统服务功能值，识别生态空间主导生态系统服务功能；模拟基期年至目标年生态空间主导生态系统服务功能格局，进行主导生态系统服务功能演变分析；
S3：采用LP-Markov-CA复合模型，进行目标年主导生态系统服务功能数量结构和空间布局耦合优化，实现生态空间主导生态系统服务功能的集约利用优化；进而指导国土空间集约利用优化配置；
步骤S3具体包括以下子步骤：
S31，将元胞C定义为主导生态系统服务功能分布图的栅格，其大小为栅格图分辨率；
S32，将元胞空间L定义为生态空间E；
S33，将元胞状态S定义为初级产品生产、气体调节、气候调节、环境净化、水文调节、水土保持、生物多样性维持、美学景观8种生态系统服务功能，各功能依次编码为1,2,3,4,5,6,7,8，则S＝{1,2,3,4,5,6,7,8}；
S34，确定元胞领域，采用标准5×5邻近滤波器；
S35，计算元胞转换面积矩阵；
S36，制作元胞转换适宜性图集，将n种经标准化处理的基期年生态空间单项生态系统服务功能值分布栅格图，按主导功能类型编码顺序组合成适宜性图集，作为元胞空间转化规则；
S37，将调整后的基期年主导生态系统服务功能分布图作为元胞空间初始状态，将步骤S35确定的元胞转换面积矩阵和步骤S36生成的元胞转换适宜性图集分别作为元胞数量和空间转换规则，利用Markov-CA模型即可模拟得到目标年生态空间集约利用优化布局；
步骤S35计算元胞转换面积矩阵的具体步骤为：
S351：进行目标年生态空间主导生态系统服务功能数量结构优化；
以初级产品生产、气体调节、气候调节、环境净化、水文调节、土壤保持、生物多样性、美学景观8类主导生态系统服务功能面积为决策变量，以主导生态系统服务功能面积条件为约束，以在有限生态空间实现最大生态系统服务功能效用为目标，建立线性规划模型，对目标年生态空间范围内不同主导生态系统服务功能面积进行优化；模型表达式为：
Figure FDA0003567359320000021
式中：f(a)为模型优化目标函数，表示生态空间各主导生态系统服务功能总值；υtj表示t年主导生态系统服务功能j单位面积功能值；aj表示j类主导生态系统服务功能面积；aE表示生态空间E总面积；atj表示t年主导生态系统服务功能j的面积；j,n分别表示生态空间主导生态系统服务功能类型编号和类型数，其中j＝1,2,3,4,5,6,7,8，n＝8；t,T分别表示时间断面编号和时间断面数；
S352：进行基期年至目标年生态空间主导生态系统服务功能转移概率矩阵预测；
利用Markov模型计算得到逐时段(0,1,2,…,t,…,T-1)主导生态系统服务功能转移概率方阵(pij)n×n(i,j＝1,2,…,n)；假设转移概率pij是随时间t连续变化的函数pij(t)，单位时间变化率r也是随时间t变化的函数rij(t)，则得到描述转移概率变化规律的微分方程及其定解条件：
Figure FDA0003567359320000022
求解定解问题，得到：
Figure FDA0003567359320000023
式中：B为待定常数，初值为(pij)0；
利用公式(8)分别计算基期年(第t年)至目标年逐时段主导生态系统服务功能转移概率
Figure FDA0003567359320000024
在此基础上采用式(9)计算得到基期年至目标年逐时段转移概率矩阵
Figure FDA0003567359320000025
Figure FDA0003567359320000026
式中：
Figure FDA0003567359320000027
表示t′时段第n类主导生态系统服务功能转移概率之和；Pt′＝(pij)n×n，表示t′时段主导生态系统服务功能转移概率方阵；
最后采用式(10)计算确定基期年至目标年生态空间主导生态系统服务功能转移概率矩阵P；
Figure FDA0003567359320000028
S353：计算基期年至目标年生态空间主导生态系统服务功能转移面积矩阵；
设基期年、目标年主导生态系统服务功能优化面积矩阵分别为x＝(xij)n×1、
Figure FDA0003567359320000029
则基期年至目标年主导生态系统服务功能转移面积矩阵A＝(aij)n×n为：
A＝diag(x1,x2,…,xn)P (11)
式中：x1,x2,…,xn表示基期年主导生态系统服务功能优化面积矩阵x的元素，其计算公式为:
x＝[PT]-lA*(|PT|≠0) (12)；
步骤S37中将调整后的基期年主导生态系统服务功能分布图作为元胞空间初始状态，根据建立基期年主导生态系统服务功能优化面积、基期年至目标年主导生态系统服务功能转移概率矩阵与面积矩阵之间的非齐次线性方程组，求解基期年主导生态系统服务功能分布图作为元胞空间初始状态。
2.根据权利要求1所述的一种生态空间集约利用优化方法，其特征在于，步骤S1具体包括以下子步骤：
S11：建立生态空间生态系统服务功能评价单元，其中以网格为生态空间生态系统服务功能评价单元；
S12：计算评价单元的生态系统服务功能价值：获取生态系统服务价值当量，将生态系统服务价值当量作为生态系统服务功能大小的度量衡，计算评价单元的生态系统服务功能价值；采用式(1)计算评价单元的单项生态系统服务功能值，根据所述单项生态系统服务功能值，采用式(2)汇总得到评价单元的生态系统服务功能综合值；
Figure FDA0003567359320000031
Figure FDA0003567359320000032
式中：υlj表示第l类生态系统第j种生态系统服务功能当量；sk′l表示第k′个评价单元第l类生态系统面积；Vk′j,Vk′分别表示第k′个评价单元第j种生态系统服务功能值和功能综合值；M,n,K分别表示生态系统类型数、生态系统服务功能类型数和评价单元个数；l,j,k′分别表示生态系统类型、生态系统服务功能类型和评价单元编号，l＝1,2,…,M；j＝1,2,…,n；k′＝1,2,…,K；
S13：识别生态空间：在长时序生态系统服务综合功能变化幅度分析基础上，结合其功能综合值大小合理划定生态空间，即将生态系统服务功能综合值大于等于某临界值且其平均变化率小于某临界值的单元划为生态空间，其数学表达式为：
Figure FDA0003567359320000033
式中：E表示生态空间；ek表示组成生态空间E的矢量网格；
Figure FDA0003567359320000034
表示t年单元k的生态系统服务功能综合值，计算方法见式(2)；k,m分别表示组成生态空间E的网格编号和个数；t,T分别表示任一时间断面和时间断面个数；δ,ε分别表示划分生态空间与非生态空间的生态系统服务功能综合值、生态系统服务综合功能值年均变化率临界值，采用标准差分级法分析确定，将t年生态系统服务功能综合值不小于δ且年均变化率不大于ε的网格识别为生态空间，否则为非生态空间。
3.根据权利要求1所述的一种生态空间集约利用优化方法，其特征在于，步骤S2具体包括以下子步骤：
S21：根据各单项生态系统服务功能值，识别生态空间主导生态系统服务功能：以矢量网格为识别单元，将网格单元对应的功能值最大的生态系统服务类型作为该矢量网格的主导功能，各网格单元主导功能集合即为生态空间主导功能格局，其数学表达式为：
Figure FDA0003567359320000041
式中：D表示生态空间主导生态系统服务功能格局；
Figure FDA0003567359320000042
表示第k个网格的主导生态系统服务功能，其编号为jk∈[1,n]，编码值1、2、3、4、5、6、7、8分别代表初级产品生产、气体调节、气候调节、环境净化、水文调节、土壤保持、生物多样性、美学景观功能；
Figure FDA0003567359320000043
表示第k个网格主导功能jk的生态系统服务功能值；Vkj表示第k个网格第j类生态系统服务功能值；
S22：进行生态空间主导生态系统服务功能格局模拟：应用Markov-CA模型模拟基期年至目标年主导生态空间生态系统服务功能格局；
S23：进行生态空间主导生态系统服务功能演变分析：基于T个时间断面主导生态空间生态系统服务功能分布图，统计得到t年j项主导生态系统服务功能面积atj，分析得到主导生态空间生态系统服务功能区最大面积max1≤t≤Tatj和最小面积min1≤t≤Tatj。
4.根据权利要求1所述的一种生态空间集约利用优化方法，其特征在于，步骤S3中的所述LP-Markov-CA复合模型的数学表达式为：
U＝(L,S,N,f) (5)
式中：U表示LP-Markov-CA模型框架；L表示元胞空间；S表示元胞状态集合；N表示领域内元胞状态空间向量；f表示元胞转换规则。
5.根据权利要求1所述的一种生态空间集约利用优化方法，其特征在于，步骤S36制作元胞转换适宜性图集的具体步骤为：
将基期年生态空间单项生态系统服务功能值分布栅格图，采用公式(13)标准化后作为该项功能空间分布适宜性图，按主导功能类型编码顺序组合成适宜性图集；并以此为依据对基期年主导生态系统服务功能实际分布进行调整，得到数量结构与面积矩阵x一致的空间分布，如下：
Figure FDA0003567359320000051
式中：Vkj,V′kj表示第j种生态系统服务功能在评价单元k的生态系统服务功能值和标准化值，k,m分别表示组成生态空间E的网格编号和个数。
6.如权利要求1至5中任一所述的一种生态空间集约利用优化方法的集约利用优化系统，其特征在于，该系统包括：
生态空间动态获取模块，用于根据多年生态系统服务功能综合值及其变化特征，动态获取生态空间；所述生态系统服务功能类型包括初级产品生产、气体调节、气候调节、环境净化、水文调节、土壤保持、生物多样性、美学景观功能；
主导生态系统服务功能模拟及演变模块：用于根据各个单项生态系统服务功能值，识别生态空间主导生态系统服务功能；模拟基期年至目标年生态空间主导生态系统服务功能格局，进行主导生态系统服务功能演变分析；
集约利用优化模块，用于采用LP-Markov-CA复合模型，进行目标年主导生态系统服务功能数量结构和空间布局耦合优化，实现生态空间主导生态系统服务功能集约利用优化，得到空间利用效益最大化的生态空间集约利用优化布局；
输出模块，用于输出空间利用效益最大化的生态空间集约利用优化布局。
7.根据权利要求6所述的集约利用优化系统，其特征在于，所述集约利用优化模块中的所述LP-Markov-CA复合模型的数学表达式为：
U＝(L,S,N,f) (5)
式中：U表示LP-Markov-CA模型框架；L表示元胞空间；S表示元胞状态集合；N表示领域内元胞状态空间向量；f表示元胞转换规则。
8.根据权利要求7所述的集约利用优化系统，其特征在于，所述集约利用优化模块的执行过程如下：
将元胞C定义为主导生态系统服务功能分布图的栅格，其大小为栅格图分辨率；
将元胞空间L定义为生态空间E；
将元胞状态S定义为初级产品生产、气体调节、气候调节、环境净化、水文调节、水土保持、生物多样性维持、美学景观8种生态系统服务功能，各功能依次编码为1,2,3,4,5,6,7,8，则S＝{1,2,3,4,5,6,7,8}；
确定元胞领域，采用标准5×5邻近滤波器；
计算元胞转换面积矩阵；
制作元胞转换适宜性图集，将n种经标准化处理的基期年生态空间单项生态系统服务功能值分布栅格图，按主导功能类型编码顺序组合成适宜性图集，作为元胞空间转化规则；
将调整后的基期年主导生态系统服务功能分布图作为元胞空间初始状态，将确定的元胞转换面积矩阵和生成的元胞转换适宜性图集分别作为元胞数量和空间转换规则，利用Markov-CA模型即可模拟得到目标年生态空间集约利用优化布局。
CN202110633616.2A 2021-06-07 2021-06-07 一种生态空间集约利用优化方法及系统 Active CN113361774B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110633616.2A CN113361774B (zh) 2021-06-07 2021-06-07 一种生态空间集约利用优化方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110633616.2A CN113361774B (zh) 2021-06-07 2021-06-07 一种生态空间集约利用优化方法及系统

Publications (2)

Publication Number Publication Date
CN113361774A CN113361774A (zh) 2021-09-07
CN113361774B true CN113361774B (zh) 2022-06-03

Family

ID=77532899

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110633616.2A Active CN113361774B (zh) 2021-06-07 2021-06-07 一种生态空间集约利用优化方法及系统

Country Status (1)

Country Link
CN (1) CN113361774B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115619095A (zh) * 2022-10-17 2023-01-17 四川农业大学 一种基于最优均衡战略的国土空间适宜性评价方法及系统
CN116090707A (zh) * 2023-03-28 2023-05-09 北京世纪农丰土地科技有限公司 一种基于空间生态修复数据的生态问题识别系统

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110069885A (zh) * 2019-05-05 2019-07-30 重庆师范大学 一种基于外部生态功能定位的“三生”空间优化方法
CN110298411A (zh) * 2019-07-04 2019-10-01 中国城市建设研究院有限公司 一种城市群生态空间受损识别评价方法
CN111754132A (zh) * 2020-06-30 2020-10-09 四川省国土整治中心 一种生态用地占补平衡计算与应用方法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
BRPI0909446A2 (pt) * 2008-04-18 2015-12-22 Exxonmobil Upstream Res Co métodos para planejamento do desenvolvimento de reservatório, para suporte à decisão considerando o desenvolvimento de recurso petrolíferos, para otimização do planejamento de desenvolvimento, e para produção de hidrocarbonetos.
US11631022B2 (en) * 2018-03-29 2023-04-18 Daybreak, Llc Forecasting soil and groundwater contamination migration
CN110059855B (zh) * 2019-03-12 2023-04-25 西安理工大学 一种区域土地利用优化配置方法
CN110363356A (zh) * 2019-07-17 2019-10-22 长安大学 一种面向生态的水土资源时空优化配置方法
CN111882245A (zh) * 2020-08-07 2020-11-03 中国环境科学研究院 一种基于生态系统服务视角的生态修复空间识别方法
CN112884226A (zh) * 2021-02-23 2021-06-01 湖南大学 基于多智能体算法的国土空间格局模拟规划方法及系统

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110069885A (zh) * 2019-05-05 2019-07-30 重庆师范大学 一种基于外部生态功能定位的“三生”空间优化方法
CN110298411A (zh) * 2019-07-04 2019-10-01 中国城市建设研究院有限公司 一种城市群生态空间受损识别评价方法
CN111754132A (zh) * 2020-06-30 2020-10-09 四川省国土整治中心 一种生态用地占补平衡计算与应用方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
Delimiting Ecological Space and Simulating its Dominant functions in the Perspective of Dynamic Changes in Ecosystem Service Functions: A Case Study in the Qionglai City of Sichuan Province, China;Nengjun Wu，Yuanxi Li，Dinghua Ou等;《IOP Conference Series: Earth and Environmental Science》;20210512;参见正文第1-14页 *
Nengjun Wu，Yuanxi Li，Dinghua Ou等.Delimiting Ecological Space and Simulating its Dominant functions in the Perspective of Dynamic Changes in Ecosystem Service Functions: A Case Study in the Qionglai City of Sichuan Province, China.《IOP Conference Series: Earth and Environmental Science》.2021, *
城市近郊区景观格局变化特征、潜力与模拟——以成都市龙泉驿区为例;欧定华，夏建国;《地理研究》;20160323;第534-550页 *
基于空间回归模型的城郊区土地景观格局变化驱动力研究——以成都市龙泉驿区为例;欧定华，夏建国，欧晓芳;《江苏农业科学》;20170820;第45卷(第14期);第236-242页 *

Also Published As

Publication number Publication date
CN113361774A (zh) 2021-09-07

Similar Documents

Publication Publication Date Title
Ali et al. Review of urban building energy modeling (UBEM) approaches, methods and tools using qualitative and quantitative analysis
Ren et al. Assessing the reliability, resilience and vulnerability of water supply system under multiple uncertain sources
CN113361774B (zh) 一种生态空间集约利用优化方法及系统
CN101551884B (zh) 面向大规模样本的cvr电力负荷快速预测方法
CN103853106A (zh) 一种建筑能源供应设备的能耗预测参数优化方法
Nichiforov et al. Deep learning techniques for load forecasting in large commercial buildings
Fan et al. Development of PCA-based cluster quantile regression (PCA-CQR) framework for streamflow prediction: Application to the Xiangxi river watershed, China
CN111898258A (zh) 一种水文循环变异驱动下二维干旱灾害评估方法
CN105893669A (zh) 一种基于数据挖掘的全局仿真性能预测方法
CN110457369A (zh) 一种模型的训练方法及相关设备
CN104636834B (zh) 一种改进的联合概率规划模型系统优化方法
CN106599585A (zh) 基于并行蜂群算法的水文模型参数优化方法及装置
CN108764527B (zh) 一种土壤有机碳库时空动态预测最优环境变量筛选方法
CN111311026A (zh) 一种顾及数据特征、模型和校正的径流非线性预测方法
CN103353895A (zh) 一种配电网线损数据的预处理方法
CN114897264A (zh) 一种基于迁移学习的小样本场景下光伏出力区间预测方法
de Chalendar et al. A physics-informed data reconciliation framework for real-time electricity and emissions tracking
Fayer et al. A temporal fusion transformer deep learning model for long-term streamflow forecasting: a case study in the funil reservoir, Southeast Brazil
CN114091776A (zh) 一种基于K-means的多分支AGCNN短期电力负荷预测方法
CN102663562B (zh) 一种基于克隆选择算法的土地资源评价因子级别划分方法
CN115936242B (zh) 空气质量与交通状况的溯源关系数据获取方法及装置
CN117114176A (zh) 基于数据分析和机器学习的土地利用变化预测方法及系统
CN111859787A (zh) 基于气象水文预报因子分类和深度学习模型耦合的日径流预报方法
CN110852496A (zh) 一种基于lstm循环神经网络的天然气负荷预测方法
CN108596781A (zh) 一种电力系统数据挖掘与预测整合方法

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