CN106372277A - 森林立地指数时空估测中的变异函数模型优化方法 - Google Patents
森林立地指数时空估测中的变异函数模型优化方法 Download PDFInfo
- Publication number
- CN106372277A CN106372277A CN201610694561.5A CN201610694561A CN106372277A CN 106372277 A CN106372277 A CN 106372277A CN 201610694561 A CN201610694561 A CN 201610694561A CN 106372277 A CN106372277 A CN 106372277A
- Authority
- CN
- China
- Prior art keywords
- variation function
- model
- spatial
- computer
- site index
- 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
- 238000000034 method Methods 0.000 title claims abstract description 67
- 238000005457 optimization Methods 0.000 title claims abstract description 18
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 20
- 230000006870 function Effects 0.000 claims description 145
- 230000002123 temporal effect Effects 0.000 claims description 32
- 238000004458 analytical method Methods 0.000 claims description 15
- 238000009826 distribution Methods 0.000 claims description 13
- 238000012360 testing method Methods 0.000 claims description 11
- 239000011159 matrix material Substances 0.000 claims description 9
- 238000006243 chemical reaction Methods 0.000 claims description 8
- 238000005070 sampling Methods 0.000 claims description 8
- 238000002790 cross-validation Methods 0.000 claims description 6
- 238000005516 engineering process Methods 0.000 claims description 6
- 238000002474 experimental method Methods 0.000 claims description 6
- 238000011160 research Methods 0.000 claims description 6
- 230000008569 process Effects 0.000 claims description 5
- 241001269238 Data Species 0.000 claims description 4
- 230000005856 abnormality Effects 0.000 claims description 4
- 230000000694 effects Effects 0.000 claims description 4
- 230000008859 change Effects 0.000 claims description 3
- 238000003745 diagnosis Methods 0.000 claims description 3
- 238000007689 inspection Methods 0.000 claims description 3
- 230000035772 mutation Effects 0.000 claims description 3
- 238000012545 processing Methods 0.000 claims description 3
- 230000008080 stochastic effect Effects 0.000 claims description 3
- 238000005309 stochastic process Methods 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 3
- 238000013519 translation Methods 0.000 claims description 3
- 241000208340 Araliaceae Species 0.000 claims 1
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 claims 1
- 235000003140 Panax quinquefolius Nutrition 0.000 claims 1
- 235000008434 ginseng Nutrition 0.000 claims 1
- 238000011002 quantification Methods 0.000 abstract description 4
- 238000007726 management method Methods 0.000 description 6
- 241000218597 Picea engelmannii Species 0.000 description 4
- PCHJSUWPFVWCPO-UHFFFAOYSA-N gold Chemical compound [Au] PCHJSUWPFVWCPO-UHFFFAOYSA-N 0.000 description 4
- 239000010931 gold Substances 0.000 description 4
- 229910052737 gold Inorganic materials 0.000 description 4
- 239000011435 rock Substances 0.000 description 4
- 239000002689 soil Substances 0.000 description 4
- 238000013459 approach Methods 0.000 description 3
- 238000013461 design Methods 0.000 description 3
- 238000011835 investigation Methods 0.000 description 3
- 241001069925 Orestes Species 0.000 description 2
- 230000008901 benefit Effects 0.000 description 2
- 238000004364 calculation method Methods 0.000 description 2
- 238000005520 cutting process Methods 0.000 description 2
- 238000013506 data mapping Methods 0.000 description 2
- 230000006735 deficit Effects 0.000 description 2
- 238000009795 derivation Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 239000003864 humus Substances 0.000 description 2
- 230000007935 neutral effect Effects 0.000 description 2
- 238000000926 separation method Methods 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- 238000004162 soil erosion Methods 0.000 description 2
- 238000012732 spatial analysis Methods 0.000 description 2
- 238000012876 topography Methods 0.000 description 2
- 238000002604 ultrasonography Methods 0.000 description 2
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 230000004927 fusion Effects 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 238000003786 synthesis reaction Methods 0.000 description 1
- 230000009897 systematic effect Effects 0.000 description 1
- 238000005303 weighing Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开的是森林立地指数时空估测中的变异函数模型优化方法，包括以下步骤：(1)利用计算机进行建立一个可靠的立地指数变异函数模型优选方法，确保变异函数选择的定型化与定量化的有效结合；(2)利用计算机进行构建一个统一的立地指数变异函数多尺度套盒模型表达形式和计算机自动参数最优拟合方法，在主流地统计软件中实现多尺度套合模型功能的扩展算法，确保空间插值算法的有效预测；(3)利用计算机进行构建一个立地指数时空变异函数模型表达形式，进而实现Kriging时空插值算法，提升时空估测精度和可靠性。可有效降低人为主观因素，有效刻画立地指数各向异性和多尺度依赖性，确保空间插值算法有效预测，提升时空估测精度和可靠性。
Description
技术领域
本发明涉及林业技术领域，具体是涉及天山云杉立地指数时空估测中的变异函数模型优化方法。
背景技术
森林生态系统立地质量的评估是森林生态系统经营管理和造林营林的重要理论基础与规划方法,也是研究森林生态系统生产力的重要内容。
立地在生态学上又被称为“生境”，指的是“林地环境和该环境所决定的林地上的植被类型及质量”。更确切地说，立地是森林或其它植被类型生存的空间及相关的自然因子的综合。科学地划分立地类型、评价立地质量，在森林经营中具有重大意义。在制定森林经营方案、确定森林经营的目标和方向；造林调查设计中安排树种规划、初植密度、营林活动中考虑林木培育、间伐、确定主伐年龄、轮伐期，计算材积和生长、预估收货及产量等都涉及到立地质量问题。
森林立地指数是一种已被普遍接受的评价森林立地质量的方法,然而由于技术所限,传统的森林立地指数调查方法因其耗时长、花费大,很难大面积的对森林立地指数进行估测和空间分析,地统计学结合地理信息系统技术为研究森林立地及空间格局提供了一条新的途径。
本发明选择天山中部的天山云杉典型分布区(天山中部的新疆农业大学实习林场)作为实验区,选取合适的变异函数理论模型,对其森林立地指数进行了估测和空间分析,并在精度验证的基础上对结果进行了评价。研究结果表明,天山中部的天山云杉立地指数有中等程度的空间相关，说明地统计学在森林立地指数的估测和空间分析中有一定的应用价值。
发明内容
本发明解决的技术问题是，结合已有变异函数套盒模型表达方法，构建一个统一的立地指数变异函数多尺度套盒模型，再基于地统计学中的变异函数理论及其公式，构建时空变异函数模型表达形式，并实现普通空间Kriging插值的时空扩展。
本发明的技术方案是：
森林立地指数时空估测中的变异函数模型优化方法,其特征在于，包括以下步骤：
(1)选取计算机设备，利用计算机建立一个可靠的立地指数变异函数模型优选方法，确保变异函数选择的定型化与定量化的有效结合；
(2)利用计算机构建一个统一的立地指数变异函数多尺度套盒模型表达形式和计算机自动参数最优拟合方法，在主流地统计软件中实现多尺度套合模型功能的扩展算法，有效刻画立地指数的各向异性和多尺度依赖性，确保空间插值算法的有效预测；
(3)利用计算机构建一个立地指数时空变异函数模型表达形式，进而实现Kriging时空插值算法，提升时空估测精度和可靠性。
进一步地，在上述方案中，步骤(1)所述的建立一个可靠的立地指数变异函数模型优选方法具体为：
①选择实验区，考虑不同海拔、坡向、坡位、郁闭度和立地条件等因素，采用空间平衡抽样方法设计布置临时样地，并确定样地面积，确保每个样地的立地指数均匀分布；
②数据采集和预处理，在每块样地内测定林木胸径、树高、年龄、冠幅和枝下高以及各样地的林分优势木平均高、郁闭度、海拔和公顷株数等测树学因子；用超声波/激光测高仪获取树高或优势高；在条件允许的情况下，在每块样地内测定气候、土壤和地形3类因子；气候因子：温度、降雨量、干旱时间、干燥指数和水分亏损等；土壤因子：土层厚度、裸岩率、PH值、含沙量、容重、质地、腐殖质厚度以及土壤侵蚀程度等；地形因子：坡度、坡位、坡向、岩石类型和坡形等；同时获取实验区1∶1万地形图和DEM数据。对采集的各类原始数据输入到计算机设备中，并进行统计检查，剔除异常和错误数据；
③利用计算机进行样本数据计算，根据实验数据拟合一个立地指数模型，并计算所有样地的立地指数作为样本数据集；应用GIS技术绘制和检查样本数据空间，基于数据值的空间统计分布，找出异常值和错误值以及全局趋势和空间自相关分析的主导方向；
④探索性空间数据分析，利用计算机进行研究样本数据的分布特征，若不符合正态分布，分析数据中的错误并进行处理，再使用对数变换、Box-Cox变换和正态计分变换等数据变换方式对处理后的样本数据进行变换，使之符合正态分布；查找样本数据的全局异常值和局部异常值；计算多种局部统计量(均值、众数、标准差、熵和集群等)，检查立地指数局部变化；分析样本数据中立地指数的空间自相关性和各向异性。
进一步地，在上述方案中，步骤(2)所述的构建一个统一的立地指数变异函数多尺度套盒模型表达形式和计算机自动参数最优拟合方法具体为：
①样本数据随机抽样，利用计算机进行样本数据进行随机抽样，使用70％的数据作为建模样本，30％的数据作为检验样本；
②利用计算机进行变异函数理论模型参数分析，基于球状模型、指数模型、高斯模型等主要变异函数理论模型，分别拟合对立地指数影响最大的几种随机效应因子所对应的实验变异函数模型，并进行块金值、基台值、变程等参数分析；
③利用计算机进行立地指数变异函数多尺度套盒建模，借鉴已有的区域化变量变异函数套盒模型表达方法，通过数学公式推导，选择一种合适的数据结构将多个模型合并，有效刻画立地指数变异函数空间变异特征，形成最优变异函数套盒方法，进而构建一个统一的立地指数变异函数多尺度套盒模型表达形式；
④计算机自动参数最优拟合方法，应用非线性回归理论、进化计算和积分面积法等新兴方法解决立地指数变异函数多尺度套盒模型存在的多参数(且参数数量可变)非线性求解问题，确定最优拟合方法进行计算机自动参数优化估计；
⑤利用计算机进行模型验证与诊断，采用标准平均值、标准均方根预测误差、块金值、基台值、变程和块金比值等检验指标对多尺度套盒模型与传统变异函数理论模型进行比较；
⑥在主流地统计软件中实现多尺度套合模型功能的扩展算法，在ArcGIS10.1中，利用Python语言实现立地指数变异函数多尺度套盒模型扩展模块。
更进一步地，步骤②所述的变异函数理论模型为：
假设变异函数表示分隔距离为h的两点xi和xi+hi的区域化变量Z(xi)和Z(xi+hi)之间的变异，可用其增量〔Z(xi)-Z(xi+hi)〕平方的数学期望表示：
式中，N(h)为距离等于h的点对数，Z(xi)为处于点xi处变量的实测值，Z(xi+hi)为与点xi偏离h处变量的实测值，变异函数中有3个最重要的参数；快金值(Nugget)，基台值(SiII)和变程(Range)。
进一步地，在上述方案中，步骤(3)所述的构建一个立地指数时空变异函数模型表达形式，进而实现Kriging时空插值算法，具体为：
①立地指数时序分解，计算立地指数时间序列中的趋势项T、季节项S和随机项R，获取时间有序间隔，用作后续变异函数建模和Kriging插值；
②时空变异函数建模，在变异函数理论模型公式中引入时间变量，通过数学推导构建时空变异函数模型表达形式；
③时空变异函数曲线拟合，利用建模样本数据集，编程构建协方差函数曲面，对时空变异函数模型的有效性和精度进行分析和评价；
④Kriging时空插值方法，分析普通Kriging空间预测公式，引入时间变量，构建立地指数、时间、空间这一3维空间中的Kriging插值算法；
⑤时空交叉验证，采用k折交叉验证法，利用标准平均值、标准均方根预测误差、块金值、基台值、变程和块金比值等验证指标对时空变异函数模型进行验证。
更进一步地，步骤②所述的空变异函数建模如下：
假设Z(s，t)是定义在Rk×T上的时空随机过程，其中Rk代表k维的欧式空间，T代表时间，(si，ti)，i＝1，2，…，n为时空场中的任意样本点位置，并且h为样本点间的时空间隔距离；设定时空距离h＝(hs，ht)，当Z(s，t)满足二阶平稳时，可定义其协方差函数为：
C(hs，ht)＝Cov(Z(s+hs，t+ht)-Z(s，t)) (1)
同时，变异函数为：
γ(hs，ht)＝1/2E(Z(s+hs，t+ht)-Z(s，t))2＝σ2-C(hs，ht) (2)
式中σ2为Z(s，t)的方差，在满足相应的正定条件下，变异函数是有效的；
计算采样点的样本变异函数并选用时空变异函数模型进行拟合，采用一类积和式变异函数来拟合月均气温的时空变异结构，如下：
Cst(hs，ht)＝k1Cs(hs)Ct(ht)+k2Cs(hs)+k3Ct(ht) (3)
γst(hs，ht)＝(k1Ct(0)+k2)γs(hs)+
(k1Cs(0)+k3)γt(ht)-k1γs(hs)γt(ht) (4)
其中，Cst为时空协方差，Cs为空间协方差，Ct为时间协方差，γst、γs、γt分别是对应的变异函数，而Cst(0，0)、Cs(0)、Ct(0)分别是对应的基台值。模型中的系数k1、k2、k3由下式决定
式(4)的时空变异函数模型，其具体实现过程如下：
a.分别计算纯时间域和纯空间域样本变异函数并进行拟合。变异函数模型拟合通过tem.fit<-fit.variogram(object，vgm(psill，model，range，nugget))，其中object是通过variogram()函数得到的样本变异函数值，model是选用的变异函数模型(如“Sph”，“Gau”，“Exp”等)，psill，range和nugget分别是model所指变异函数的基台值、变程和块金值。然而psill，range和nugget这三个参数是用户预先设定的估计值，并不代表模型实际的拟合参数，实际值被返回至数组tem.fit中。
b.根据步骤a结果分别构建空间变异函数rs(hs)和时间变异函数rt(ht)，并由式(5)计算出系数k1、k2、k3的值。本实验中，rs(hs)采用高斯模型拟合，rt(ht)采用球形模型拟合，分别由以下R程序段得到：
rs<-function(hs){nugget_s+psill_s*(1－exp(－(hs/range_s)2))}rt<-function(ht){if(hs>range_t)return(psill_t+nugget_t)(psill_t+nugget_t)*(1.5*ht/range_t－0.5*(hs/range_t)3)}
c.将rs(hs)和rt(ht)代入式(4)得到时空变异函数模型rst(hs，ht)。
更进一步地，步骤④所述Kriging时空插值采用普通Kriging方法实现数据的时空插值，如下：
式中：Z*(s0，t0)为时空点(s0，t0)处的估计值，λi为邻近观测值Z(si，ti)的加权系数，引入拉格朗日系数μ进行推导可得：
式(7)中等号左侧的系数矩阵定义为coef，右侧的系数矩阵定义为o_coef，上式中的加权系数λ和拉格朗日系数μ的值可通过lamda<-solve(coef，o_coef)得到，继而代入式(6)可得研究区域内任意点的插值估计。
本发明采用野外实验调查、理论研究、数理统计、地统计、时空分析与模拟、算法分析与设计相结合的方法，利用计算机进行系统研究森林立地指数时空估测中变异函数模型优化方法，实现静态模型模拟与动态时空估测的融合。与传统方法相比，本发明的有益效果是：
(1)利用计算机进行建立一个可靠的立地指数变异函数模型优选方法，确保变异函数优选过程的定型化与定量化的结合，有效降低人为主观因素；
(2)利用计算机进行构建一个统一的立地指数变异函数多尺度套盒模型表达形式和计算机自动参数最优拟合方法，有效刻画立地指数各向异性和多尺度依赖性，确保空间插值算法有效预测；
(3)利用计算机进行构建一个立地指数时空变异函数模型表达形式，进而实现Kriging时空插值算法，有效反应立地指数的时空过程与分布，提升时空估测精度和可靠性。
附图说明
图1为本发明实施例各主要研究内容之间的关联关系的研究技术路线框图。
具体实施方式
下面结合具体实施方式来对本发明进行进一步详细的说明：
本发明主要研究内容之间的关联关系的研究技术路线框图如图1所示。
森林立地指数时空估测中的变异函数模型优化方法,包括以下步骤：
(1)选取计算机设备，利用计算机进行建立一个可靠的立地指数变异函数模型优选方法，确保变异函数选择的定型化与定量化的有效结合；
所述立地指数变异函数模型优选方法具体为：
①选择天山中部的天山云杉典型分布区(天山中部的新疆农业大学实习林场)作为实验区，考虑不同海拔、坡向、坡位、郁闭度和立地条件等因素，采用空间平衡抽样方法设计布置临时样地，并确定样地面积为400m2(20m*20m)，确保每个样地的立地指数均匀分布；
②数据采集和预处理，在每块样地内测定林木胸径、树高、年龄、冠幅和枝下高以及各样地的林分优势木平均高、郁闭度、海拔和公顷株数等测树学因子；用超声波/激光测高仪获取树高或优势高；在条件允许的情况下，在每块样地内测定气候、土壤和地形3类因子；气候因子：温度、降雨量、干旱时间、干燥指数和水分亏损等；土壤因子：土层厚度、裸岩率、PH值、含沙量、容重、质地、腐殖质厚度以及土壤侵蚀程度等；地形因子：坡度、坡位、坡向、岩石类型和坡形等；同时获取实验区1∶1万地形图和DEM数据。对采集的各类原始数据输入到计算机设备中，并进行统计检查，剔除异常和错误数据；
③利用计算机进行样本数据计算，根据实验数据拟合一个立地指数模型，并计算所有样地的立地指数作为样本数据集；应用GIS技术绘制和检查样本数据空间，基于数据值的空间统计分布，找出异常值和错误值以及全局趋势和空间自相关分析的主导方向；
④探索性空间数据分析，利用计算机进行研究样本数据的分布特征，若不符合正态分布，分析数据中的错误并进行处理，再使用对数变换、Box-Cox变换和正态计分变换等数据变换方式对处理后的样本数据进行变换，使之符合正态分布；查找样本数据的全局异常值和局部异常值；计算多种局部统计量(均值、众数、标准差、熵和集群等)，检查立地指数局部变化；分析样本数据中立地指数的空间自相关性和各向异性。
(2)利用计算机进行构建一个统一的立地指数变异函数多尺度套盒模型表达形式和计算机自动参数最优拟合方法，在主流地统计软件中实现多尺度套合模型功能的扩展算法，有效刻画立地指数的各向异性和多尺度依赖性，确保空间插值算法的有效预测；
所述的构建一个统一的立地指数变异函数多尺度套盒模型表达形式和计算机自动参数最优拟合方法具体为：
①利用计算机进行样本数据随机抽样，对样本数据进行随机抽样，使用70％的数据作为建模样本，30％的数据作为检验样本；
②变异函数理论模型参数分析，利用计算机进行基于球状模型、指数模型、高斯模型等主要变异函数理论模型，分别拟合对立地指数影响最大的几种随机效应因子所对应的实验变异函数模型，并进行块金值、基台值、变程等参数分析；
所述的变异函数理论模型为：
假设变异函数表示分隔距离为h的两点xi和xi+hi的区域化变量Z(xi)和Z(xi+hi)之间的变异，可用其增量〔Z(xi)-Z(xi+hi)〕平方的数学期望表示：
式中，N(h)为距离等于h的点对数，Z(xi)为处于点xi处变量的实测值，Z(xi+hi)为与点xi偏离h处变量的实测值，变异函数中有3个最重要的参数；快金值(Nugget)，基台值(SiII)和变程(Range)；
③利用计算机进行立地指数变异函数多尺度套盒建模，借鉴已有的区域化变量变异函数套盒模型表达方法，通过数学公式推导，选择一种合适的数据结构将多个模型合并，有效刻画立地指数变异函数空间变异特征，形成最优变异函数套盒方法，进而构建一个统一的立地指数变异函数多尺度套盒模型表达形式；
④计算机自动参数最优拟合方法，利用计算机，应用非线性回归理论、进化计算和积分面积法等新兴方法解决立地指数变异函数多尺度套盒模型存在的多参数(且参数数量可变)非线性求解问题，确定最优拟合方法进行计算机自动参数优化估计；
⑤利用计算机进行模型验证与诊断，采用标准平均值、标准均方根预测误差、块金值、基台值、变程和块金比值等检验指标对多尺度套盒模型与传统变异函数理论模型进行比较；
⑥在主流地统计软件中实现多尺度套合模型功能的扩展算法，在ArcGIS10.1中，利用Python语言实现立地指数变异函数多尺度套盒模型扩展模块。
(3)构建一个立地指数时空变异函数模型表达形式，进而实现Kriging时空插值算法，提升时空估测精度和可靠性，具体为：
①立地指数时序分解，计算立地指数时间序列中的趋势项T、季节项S和随机项R，获取时间有序间隔，用作后续变异函数建模和Kriging插值；
②时空变异函数建模，在变异函数理论模型公式中引入时间变量，通过数学推导构建时空变异函数模型表达形式：
步骤②所述的空变异函数建模如下：
假设Z(s，t)是定义在Rk×T上的时空随机过程，其中Rk代表k维的欧式空间，T代表时间，(si，ti)，i＝1，2，…，n为时空场中的任意样本点位置，并且h为样本点间的时空间隔距离；设定时空距离h＝(hs，ht)，当Z(s，t)满足二阶平稳时，可定义其协方差函数为：
C(hs，ht)＝Cov(Z(s+hs，t+ht)-Z(s，t)) (1)
同时，变异函数为：
γ(hs，ht)＝1/2E(Z(s+hs，t+ht)-Z(s，t))2＝σ2-C(hs，ht) (2)
式中σ2为Z(s，t)的方差，在满足相应的正定条件下，变异函数是有效的；
计算采样点的样本变异函数并选用时空变异函数模型进行拟合，采用一类积和式变异函数来拟合月均气温的时空变异结构，如下：
Cst(hs，ht)＝k1Cs(hs)Ct(ht)+k2Cs(hs)+k3Ct(ht) (3)
γst(hs，ht)＝(k1Ct(0)+k2)γs(hs)+
(k1Cs(0)+k3)γt(ht)-k1γs(hs)γt(ht) (4)
其中，Cst为时空协方差，Cs为空间协方差，Ct为时间协方差，γst、γs、γt分别是对应的变异函数，而Cst(0，0)、Cs(0)、Ct(0)分别是对应的基台值，模型中的系数k1、k2、k3由下式决定，
式(4)的时空变异函数模型，其具体实现过程如下：
a.分别计算纯时间域和纯空间域样本变异函数并进行拟合。变异函数模型拟合通过tem.fit<-fit.variogram(object，vgm(psill，model，range，nugget))，其中object是通过variogram()函数得到的样本变异函数值，model是选用的变异函数模型(如“Sph”，“Gau”，“Exp”等)，psill，range和nugget分别是model所指变异函数的基台值、变程和块金值。然而psill，range和nugget这三个参数是用户预先设定的估计值，并不代表模型实际的拟合参数，实际值被返回至数组tem.fit中。
b.根据步骤a结果分别构建空间变异函数rs(hs)和时间变异函数rt(ht)，并由式(5)计算出系数k1、k2、k3的值。本实验中，rs(hs)采用高斯模型拟合，rt(ht)采用球形模型拟合，分别由以下R程序段得到：
rs<-function(hs){nugget_s+psill_s*(1－exp(－(hs/range_s)2))}rt<-function(ht){if(hs>range_t)return(psill_t+nugget_t)(psill_t+nugget_t)*(1.5*ht/range_t－0.5*(hs/range_t)3)}
c.将rs(hs)和rt(ht)代入式(4)得到时空变异函数模型rst(hs，ht)；
③时空变异函数曲线拟合，利用建模样本数据集，编程构建协方差函数曲面，对时空变异函数模型的有效性和精度进行分析和评价；
④Kriging时空插值方法，分析普通Kriging空间预测公式，引入时间变量，构建立地指数、时间、空间这一3维空间中的Kriging插值算法：
所述Kriging时空插值采用普通Kriging方法实现数据的时空插值，如下：
式中：Z*(s0，t0)为时空点(s0，t0)处的估计值，λi为邻近观测值Z(si，ti)的加权系数，引入拉格朗日系数μ进行推导可得：
式(7)中等号左侧的系数矩阵定义为coef，右侧的系数矩阵定义为o_coef，上式中的加权系数λ和拉格朗日系数μ的值可通过lamda<-solve(coef，o_coef)得到，继而代入式(6)可得研究区域内任意点的插值估计。
⑤时空交叉验证，采用k折交叉验证法，利用标准平均值、标准均方根预测误差、块金值、基台值、变程和块金比值等验证指标对时空变异函数模型进行验证。
Claims (6)
1.森林立地指数时空估测中的变异函数模型优化方法,其特征在于，包括以下步骤：
(1)选取计算机设备，利用计算机建立一个可靠的立地指数变异函数模型优选方法；
(2)利用计算机构建一个统一的立地指数变异函数多尺度套盒模型表达形式和计算机自动参数最优拟合方法；
(3)利用计算机构建一个立地指数时空变异函数模型表达形式，进而实现Kriging时空插值算法。
2.如权利要求1所述的森林立地指数时空估测中的变异函数模型优化方法,其特征在于，步骤(1)所述的建立一个可靠的立地指数变异函数模型优选方法具体为：
①选择实验区，设置样地，确保每个样地的立地指数均匀分布；
②数据采集和预处理，选取计算机设备，对采集的各类原始数据输入到计算机中，并进行统计检查，剔除异常和错误数据，对采集的各类原始数据进行统计检查，剔除异常和错误数据；
③样本数据计算，根据实验数据拟合一个立地指数模型，利用计算机计算所有样地的立地指数作为样本数据集；应用GIS技术绘制和检查样本数据空间，基于数据值的空间统计分布，找出异常值和错误值以及全局趋势和空间自相关分析的主导方向；
④探索性空间数据分析，利用计算机研究样本数据的分布特征，若不符合正态分布，分析数据中的错误并进行处理，再使用对数变换、Box-Cox变换和正态计分变换等数据变换方式对处理后的样本数据进行变换，使之符合正态分布。
3.如权利要求1所述的森林立地指数时空估测中的变异函数模型优化方法,其特征在于，步骤(2)所述的构建一个统一的立地指数变异函数多尺度套盒模型表达形式和计算机自动参数最优拟合方法具体为：
①利用计算机进行样本数据随机抽样，对样本数据进行随机抽样，使用70％的数据作为建模样本，30％的数据作为检验样本；
②利用计算机进行变异函数理论模型参数分析，基于球状模型、指数模型、高斯模型等主要变异函数理论模型，分别拟合对立地指数影响最大的几种随机效应因子所对应的实验变异函数模型，并进行块金值、基台值、变程等参数分析；
③利用计算机立地指数变异函数多尺度套盒建模，借鉴已有的区域化变量变异函数套盒模型表达方法，通过数学公式推导，选择一种合适的数据结构将多个模型合并，有效刻画立地指数变异函数空间变异特征，形成最优变异函数套盒方法，进而构建一个统一的立地指数变异函数多尺度套盒模型表达形式；
④计算机自动参数最优拟合方法，应用非线性回归理论、进化计算和积分面积法等新兴方法解决立地指数变异函数多尺度套盒模型存在的多参数非线性求解问题，确定最优拟合方法进行计算机自动参数优化估计；
⑤利用计算机进行模型验证与诊断，采用标准平均值、标准均方根预测误差、块金值、基台值、变程和块金比值等检验指标对多尺度套盒模型与传统变异函数理论模型进行比较；
⑥在主流地统计软件中实现多尺度套合模型功能的扩展算法，在ArcGIS10.1中，利用Python语言实现立地指数变异函数多尺度套盒模型扩展模块。
4.如权利要求1所述的森林立地指数时空估测中的变异函数模型优化方法,其特征在于，步骤(3)所述的构建一个立地指数时空变异函数模型表达形式，进而实现Kriging时空插值算法，具体为：
①立地指数时序分解，利用计算机计算立地指数时间序列中的趋势项T、季节项S和随机项R，获取时间有序间隔，用作后续变异函数建模和Kriging插值；
②时空变异函数建模，利用计算机在变异函数理论模型公式中引入时间变量，通过数学推导构建时空变异函数模型表达形式；
③时空变异函数曲线拟合，在计算机中利用建模样本数据集，编程构建协方差函数曲面，对时空变异函数模型的有效性和精度进行分析和评价；
④Kriging时空插值方法，利用计算机分析普通Kriging空间预测公式，引入时间变量，构建立地指数、时间、空间这一3维空间中的Kriging插值算法；
⑤利用计算机进行时空交叉验证，采用k折交叉验证法，利用标准平均值、标准均方根预测误差、块金值、基台值、变程和块金比值等验证指标对时空变异函数模型进行验证。
5.如权利要求4所述的森林立地指数时空估测中的变异函数模型优化方法,其特征在于，步骤②所述的空变异函数建模如下：
假设Z(s，t)是定义在Rk×T上的时空随机过程，其中Rk代表k维的欧式空间，T代表时间，(si，ti)，i＝1，2，…，n为时空场中的任意样本点位置，并且h为样本点间的时空间隔距离；设定时空距离h＝(hs，ht)，当Z(s，t)满足二阶平稳时，可定义其协方差函数为：
C(hs，ht)＝Cov(Z(s+hs，t+ht)-Z(s，t)) (1)
同时，变异函数为：
γ(hs，ht)＝1/2E(Z(s+hs，t+ht)-Z(s，t))2＝σ2-C(hs，ht) (2)
式中σ2为Z(s，t)的方差，在满足相应的正定条件下，变异函数是有效的；
计算采样点的样本变异函数并选用时空变异函数模型进行拟合，采用一类积和式变异函数来拟合月均气温的时空变异结构，如下：
Cst(hs，ht)＝k1Cs(hs)Ct(ht)+k2Cs(hs)+k3Ct(ht) (3)
γst(hs，ht)＝(k1Ct(0)+k2)γs(hs)+
(k1Cs(0)+k3)γt(ht)-k1γs(hs)γt(ht) (4)
其中，Cst为时空协方差，Cs为空间协方差，Ct为时间协方差，γst、γs、γt分别是对应的变异函数，而Cst(0，0)、Cs(0)、Ct(0)分别是对应的基台值。模型中的系数k1、k2、k3由下式决定，
式(4)的时空变异函数模型，其具体实现过程如下：
a.分别计算纯时间域和纯空间域样本变异函数并进行拟合。变异函数模型拟合通过tem.fit<-fit.variogram(object，vgm(psill，model，range，nugget))，其中object是通过variogram()函数得到的样本变异函数值，model是选用的变异函数模型(如“Sph”，“Gau”，“Exp”等)，psill，range和nugget分别是model所指变异函数的基台值、变程和块金值。然而psill，range和nugget这三个参数是用户预先设定的估计值，并不代表模型实际的拟合参数，实际值被返回至数组tem.fit中。
b.根据步骤a结果分别构建空间变异函数rs(hs)和时间变异函数rt(ht)，并由式(5)计算出系数k1、k2、k3的值。本实验中，rs(hs)采用高斯模型拟合，rt(ht)采用球形模型拟合，分别由以下R程序段得到：
rs<-function(hs){nugget_s+psill_s*(1－exp(－(hs/range_s)2))}rt<-function(ht){if(hs>range_t)return(psill_t+nugget_t)(psill_t+nugget_t)*(1.5*ht/range_t－0.5*(hs/range_t)3)}
c.将rs(hs)和rt(ht)代入式(4)得到时空变异函数模型rst(hs，ht)。
6.如权利要求4所述的森林立地指数时空估测中的变异函数模型优化方法,其特征在于，步骤④所述Kriging时空插值采用普通Kriging方法实现数据的时空插值，如下：
式中：Z*(s0，t0)为时空点(s0，t0)处的估计值，λi为邻近观测值Z(si，ti)的加权系数，引入拉格朗日系数μ进行推导可得：
式(7)中等号左侧的系数矩阵定义为coef，右侧的系数矩阵定义为o_coef，上式中的加权系数λ和拉格朗日系数μ的值可通过lamda<-solve(coef，o_coef)得到，继而代入式(6)可得研究区域内任意点的插值估计。
Applications Claiming Priority (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610319829 | 2016-05-13 | ||
| CN2016103198297 | 2016-05-13 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN106372277A true CN106372277A (zh) | 2017-02-01 |
| CN106372277B CN106372277B (zh) | 2021-12-28 |
Family
ID=57878550
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610694561.5A Expired - Fee Related CN106372277B (zh) | 2016-05-13 | 2016-08-19 | 森林立地指数时空估测中的变异函数模型优化方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN106372277B (zh) |
Cited By (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107315722A (zh) * | 2017-07-03 | 2017-11-03 | 南京大学 | 基于克里金和信息熵理论相结合的水文站网优化模型 |
| CN107391877A (zh) * | 2017-08-15 | 2017-11-24 | 广西壮族自治区林业科学研究院 | 一种基于r软件的入侵植物频度数量变化调查分析方法 |
| CN107871042A (zh) * | 2017-11-06 | 2018-04-03 | 中国水利水电科学研究院 | 一种田块尺度的土壤入渗参数测算方法 |
| CN108920655A (zh) * | 2018-07-03 | 2018-11-30 | 山东大学 | 一种道路天气信息系统时空覆盖范围量化方法及装置 |
| CN108960517A (zh) * | 2018-07-10 | 2018-12-07 | 湖南城市学院 | 一种Kriging时空统一预测模型及其构建方法和应用 |
| CN109523066A (zh) * | 2018-10-29 | 2019-03-26 | 东华理工大学 | 一种基于克里金插值的pm2.5新增移动站点选址方法 |
| CN110348277A (zh) * | 2018-11-30 | 2019-10-18 | 浙江农林大学 | 一种基于自然背景下的树种图像识别方法 |
| CN110516405A (zh) * | 2019-09-11 | 2019-11-29 | 新疆农业大学 | 硅酸盐水泥基胶凝材料体系水化热无假定预测模型的构建方法 |
| CN110795846A (zh) * | 2019-10-29 | 2020-02-14 | 东北财经大学 | 边界森林模型的构建方法、面向复杂工业过程的多工况软计算模型更新方法及其应用 |
| CN112784493A (zh) * | 2021-01-27 | 2021-05-11 | 武汉轻工大学 | 基于自适应深度q网络的地理空间预测方法及系统 |
| CN113343347A (zh) * | 2021-02-01 | 2021-09-03 | 复旦大学 | 一种翼型前缘cst的垂直补偿修正方法及装置 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101718775A (zh) * | 2009-11-12 | 2010-06-02 | 上海交通大学 | 围垦地土壤中重金属含量的空间变异分布图生成方法 |
| CN103425513A (zh) * | 2013-08-19 | 2013-12-04 | 北京林业大学 | 一种森林经营决策支持模型自动更新方法 |
| CN104573393A (zh) * | 2015-01-28 | 2015-04-29 | 北京师范大学 | 一种基于贝叶斯理论的土壤水分站点数据升尺度方法 |
| CN104764868A (zh) * | 2015-04-02 | 2015-07-08 | 中国科学院南京土壤研究所 | 一种基于地理加权回归的土壤有机碳预测方法 |
-
2016
- 2016-08-19 CN CN201610694561.5A patent/CN106372277B/zh not_active Expired - Fee Related
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101718775A (zh) * | 2009-11-12 | 2010-06-02 | 上海交通大学 | 围垦地土壤中重金属含量的空间变异分布图生成方法 |
| CN103425513A (zh) * | 2013-08-19 | 2013-12-04 | 北京林业大学 | 一种森林经营决策支持模型自动更新方法 |
| CN104573393A (zh) * | 2015-01-28 | 2015-04-29 | 北京师范大学 | 一种基于贝叶斯理论的土壤水分站点数据升尺度方法 |
| CN104764868A (zh) * | 2015-04-02 | 2015-07-08 | 中国科学院南京土壤研究所 | 一种基于地理加权回归的土壤有机碳预测方法 |
Cited By (15)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107315722A (zh) * | 2017-07-03 | 2017-11-03 | 南京大学 | 基于克里金和信息熵理论相结合的水文站网优化模型 |
| CN107315722B (zh) * | 2017-07-03 | 2020-01-21 | 南京大学 | 基于克里金法和信息熵理论耦合的水文站网优化方法 |
| CN107391877A (zh) * | 2017-08-15 | 2017-11-24 | 广西壮族自治区林业科学研究院 | 一种基于r软件的入侵植物频度数量变化调查分析方法 |
| CN107391877B (zh) * | 2017-08-15 | 2020-09-29 | 广西壮族自治区林业科学研究院 | 一种基于r软件的入侵植物频度数量变化调查分析方法 |
| CN107871042B (zh) * | 2017-11-06 | 2019-08-06 | 中国水利水电科学研究院 | 一种田块尺度的土壤入渗参数测算方法 |
| CN107871042A (zh) * | 2017-11-06 | 2018-04-03 | 中国水利水电科学研究院 | 一种田块尺度的土壤入渗参数测算方法 |
| CN108920655A (zh) * | 2018-07-03 | 2018-11-30 | 山东大学 | 一种道路天气信息系统时空覆盖范围量化方法及装置 |
| CN108920655B (zh) * | 2018-07-03 | 2022-07-22 | 山东大学 | 一种道路天气信息系统时空覆盖范围量化方法及装置 |
| CN108960517A (zh) * | 2018-07-10 | 2018-12-07 | 湖南城市学院 | 一种Kriging时空统一预测模型及其构建方法和应用 |
| CN109523066A (zh) * | 2018-10-29 | 2019-03-26 | 东华理工大学 | 一种基于克里金插值的pm2.5新增移动站点选址方法 |
| CN110348277A (zh) * | 2018-11-30 | 2019-10-18 | 浙江农林大学 | 一种基于自然背景下的树种图像识别方法 |
| CN110516405A (zh) * | 2019-09-11 | 2019-11-29 | 新疆农业大学 | 硅酸盐水泥基胶凝材料体系水化热无假定预测模型的构建方法 |
| CN110795846A (zh) * | 2019-10-29 | 2020-02-14 | 东北财经大学 | 边界森林模型的构建方法、面向复杂工业过程的多工况软计算模型更新方法及其应用 |
| CN112784493A (zh) * | 2021-01-27 | 2021-05-11 | 武汉轻工大学 | 基于自适应深度q网络的地理空间预测方法及系统 |
| CN113343347A (zh) * | 2021-02-01 | 2021-09-03 | 复旦大学 | 一种翼型前缘cst的垂直补偿修正方法及装置 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN106372277B (zh) | 2021-12-28 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106372277A (zh) | 森林立地指数时空估测中的变异函数模型优化方法 | |
| Pecchi et al. | Species distribution modelling to support forest management. A literature review | |
| Wang et al. | Large-scale spatial variability of dried soil layers and related factors across the entire Loess Plateau of China | |
| Almeida et al. | Needs and opportunities for using a process-based productivity model as a practical tool in Eucalyptus plantations | |
| Oueslati et al. | Vegetation and topographic control on spatial variability of soil organic carbon | |
| Wang et al. | Mapping and spatial uncertainty analysis of forest vegetation carbon by combining national forest inventory data and satellite images | |
| Ramos et al. | The INFOSOLO database as a first step towards the development of a soil information system in Portugal | |
| Bishop et al. | Digital soil-terrain modeling: the predictive potential and uncertainty | |
| Rezaei et al. | Sensitivity of water stress in a two-layered sandy grassland soil to variations in groundwater depth and soil hydraulic parameters | |
| Baik et al. | Agricultural drought assessment based on multiple soil moisture products | |
| CN106980750B (zh) | 一种基于典型对应分析的土壤氮储量估算方法 | |
| CN113011372B (zh) | 一种盐碱化土地自动监测和识别方法 | |
| Trasobares et al. | A climate-sensitive empirical growth and yield model for forest management planning of even-aged beech stands | |
| CN114169161A (zh) | 一种土壤有机碳时空变异和固碳潜力估计方法和系统 | |
| Ramos et al. | Field-scale assessment of soil water dynamics using distributed modeling and electromagnetic conductivity imaging | |
| Tilse et al. | Mapping the impact of subsoil constraints on soil available water capacity and potential crop yield | |
| Bansod et al. | An application of PCA and fuzzy C-means to delineate management zones and variability analysis of soil | |
| Kumar et al. | Soil properties prediction for agriculture using machine learning techniques | |
| Tian et al. | Developing crown width model for mixed forests using soil, climate and stand factors | |
| Scannavino et al. | Spatial variability on soil pH gradient: A case study in vineyards | |
| Coe et al. | Designing experiments and analysing data. | |
| Yorulmaz et al. | Benchmarking multi-component spatial metrics for hydrologic model calibration using MODIS AET and LAI products | |
| Prusova et al. | Comparison of geostatistical approaches dealing with the distribution of snow | |
| Nurudeen et al. | INDEX AGE FOR ESTIMATING SITE INDEX OF Gmelina arborea (ROXB.) STANDS IN OLUWA FOREST RESERVE, SOUTH WESTERN NIGERIA | |
| Carre et al. | Mapping the CN ratio of the forest litters in Europe-lessons for global digital soil mapping |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | 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 | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20211228 |