CN113901640B - 一种基于元胞自动机的土壤养分流失仿真模拟方法 - Google Patents
一种基于元胞自动机的土壤养分流失仿真模拟方法 Download PDFInfo
- Publication number
- CN113901640B CN113901640B CN202111059917.5A CN202111059917A CN113901640B CN 113901640 B CN113901640 B CN 113901640B CN 202111059917 A CN202111059917 A CN 202111059917A CN 113901640 B CN113901640 B CN 113901640B
- Authority
- CN
- China
- Prior art keywords
- runoff
- soil
- rainfall
- cell
- infiltration
- 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
Links
- 239000002689 soil Substances 0.000 title claims abstract description 225
- 235000015097 nutrients Nutrition 0.000 title claims abstract description 112
- 238000000034 method Methods 0.000 title claims abstract description 37
- 238000004088 simulation Methods 0.000 title claims abstract description 28
- 230000001413 cellular effect Effects 0.000 title claims abstract description 16
- 238000001764 infiltration Methods 0.000 claims abstract description 104
- 230000008595 infiltration Effects 0.000 claims abstract description 104
- 230000008859 change Effects 0.000 claims abstract description 38
- 230000008569 process Effects 0.000 claims abstract description 25
- 230000002093 peripheral effect Effects 0.000 claims abstract description 5
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 18
- 238000004364 calculation method Methods 0.000 claims description 16
- 238000009825 accumulation Methods 0.000 claims description 12
- 235000021049 nutrient content Nutrition 0.000 claims description 12
- 235000016709 nutrition Nutrition 0.000 claims description 11
- 230000035764 nutrition Effects 0.000 claims description 11
- CIWBSHSKHKDKBQ-JLAZNSOCSA-N Ascorbic acid Chemical compound OC[C@H](O)[C@H]1OC(=O)C(O)=C1O CIWBSHSKHKDKBQ-JLAZNSOCSA-N 0.000 claims description 6
- 239000011159 matrix material Substances 0.000 claims description 6
- 230000000717 retained effect Effects 0.000 claims description 6
- 229920006395 saturated elastomer Polymers 0.000 claims description 6
- 239000002352 surface water Substances 0.000 claims description 6
- 238000003491 array Methods 0.000 claims description 3
- 230000000903 blocking effect Effects 0.000 claims description 3
- 238000010276 construction Methods 0.000 claims description 3
- 238000006243 chemical reaction Methods 0.000 claims description 2
- 230000009471 action Effects 0.000 abstract description 4
- 230000000694 effects Effects 0.000 abstract description 4
- 230000003993 interaction Effects 0.000 abstract description 3
- 230000002776 aggregation Effects 0.000 abstract description 2
- 238000004220 aggregation Methods 0.000 abstract description 2
- 238000009960 carding Methods 0.000 abstract 1
- 238000010586 diagram Methods 0.000 description 4
- 238000011160 research Methods 0.000 description 2
- 238000013528 artificial neural network Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 230000015556 catabolic process Effects 0.000 description 1
- 230000001186 cumulative effect Effects 0.000 description 1
- 238000006731 degradation reaction Methods 0.000 description 1
- 230000003628 erosive effect Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000011010 flushing procedure Methods 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 239000007788 liquid Substances 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 238000007637 random forest analysis Methods 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 238000012876 topography Methods 0.000 description 1
- 238000012795 verification 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
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
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
本发明公开了一种基于元胞自动机的土壤养分流失仿真模拟方法，本发明不仅可对降雨、土壤入渗与径流等水文过程与土壤养分、坡度等地形条件间的相互作用关系进行梳理，提高了土壤养分动态变化模拟的鲁棒性，并且仿真模拟的过程阐明了土壤养分在降雨、土壤入渗与径流作用下的变化过程，实现了对土壤养分动态变化的全过程反映，准确度高；且本发明还可对区域内部每块土壤与其周边区域在径流作用下的养分流失或聚集关系进行确定，可实现对土壤养分空间分布的全区域预测。
Description
技术领域
本发明涉及水文模型的技术领域，尤其涉及一种基于元胞自动机的土壤养分流失仿真模拟方法。
背景技术
土壤在降雨径流冲刷下容易出现养分流失的现象，尤其在植被覆盖度低、土壤质地松散、降雨集中的黄土高原区尤其明显。土壤养分流失将加剧区域内部立地条件的空间差异，从而导致植被格局破碎度的增加，进一步还可能引起部分植物群落的逆向演替。因此，为了避免景观破碎和植被退化，首先需要预测可能出现的土壤养分流失情况，包括其剧烈程度、空间分布、时间变化等特征，以采取针对性的防治措施。然而，降雨、径流、入渗、侵蚀等水文过程的随机性、复杂性和关联性使得定量且直观的反映土壤养分流失的时间变化和空间差异十分困难，而且传统的调查和监测手段也难以预测未来长时间尺度下的土壤养分流失情况。
为了预测大范围区域尺度内土壤养分流失的时空变化，可以采取仿真模拟方法。目前，土壤养分流失的预测主要依靠降雨数据与土壤养分变化之间的经验模型确定。通过野外观测或室内试验，持续记录降雨发生过程中的土壤养分，使用线性/非线性回归、随机森林、神经网络等算法确定土壤养分变化与降雨时间或降雨强度间的关系，建立经验函数模型以预测这一区域的土壤养分流失情况。但这类方法地域特征强、适用性弱、鲁棒性差，一旦超过经验函数的观测年限后精准度将大幅降低。
发明内容
因此，基于以上背景，本发明要解决的技术问题是：开发一种基于元胞自动机的土壤养分流失仿真模拟系统，用于模拟降雨径流作用下土壤养分的流失情况，满足生态环境保护和修复的需要。
本发明通过以下技术方案实现：
一种基于元胞自动机的土壤养分流失仿真模拟方法，其包括如下步骤：
其包括如下步骤：
S1：根据模拟区域实际大小、植被覆盖、降雨情况建立元胞数组，通过随机生成小于降雨概率的元胞数组模拟降雨事件的发生；
S2：根据各个元胞的植被覆盖情况确定土壤入渗速率，结合降雨量计算土壤入渗量和残留在地表的雨量；
S3：通过坡度对地表径流量进行转换，根据高程关系确定每个元胞内径流的流入或流出关系；
S4：根据元胞在径流过程中的径流的流入或流出关系判断土壤养分的流失或汇集，并依据植被覆盖确定土壤养分变化速率，最后计算径流发生后的土壤养分情况；
S5：建立循环函数重复步骤S1至S4，模拟降雨发生后不断发生的降雨、入渗、径流及土壤养分变化过程，当循环函数运行超过设定的降雨时间raintime后停止循环，并输出此时土壤养分流失情况，即对模拟区域的土壤养分流失情况完成了仿真模拟。
进一步地，步骤S1中的元胞数组的构建过程为：在Matlab中建立若干个大小为n*n的矩阵，矩阵中的每个元素代表一个元胞，每个矩阵则代表一个元胞数组。
进一步地，所述元胞数组的数量等于所涉及的元胞状态参数的数量，每个元胞数组内的元胞值代表大小为n*n m2的区域内每块土壤所对应参数的值；涉及到元胞状态参数包括高程t、元胞的坡度s、降雨量rain、入渗量infil、入渗速率infilrate、径流流出量flowout、径流流入量flowin、邻居元胞径流流出量surflowout、径流速率runoffrate、土壤养分含量nutrition以及土壤养分变化量nutrichange。
进一步地，步骤S1的实现方式如下：
1)降雨发生时，使用生成值域区间内随机数的方式，确定此次降雨强度arain，具体公式为
其中sumrain为降雨总量，rainday为降雨天数，raintime为降雨持续时间；
2)使用元胞数组降雨量rain反映当前每个元胞内的降雨情况，对于生成的n*n大小的随机数元胞数组，小于降雨落入概率的元胞则认为该元胞存在降雨，存在降雨的元胞状态为1，其余为0，具体公式为
rain＝arain×(rand(n)<prain) (2)
其中rand(n)为大小为n*n的随机数矩阵，prain为雨水落入概率。
进一步地，步骤S2的实现方式如下：
1)对步骤S1得到的存在降雨的元胞数组进行判断，对于存在降雨的元胞，地表将会滞留等同于降雨量的地表水，随后发生土壤入渗；
2)计算发生土壤入渗的元胞的土壤的入渗速率，根据地表覆盖类型、坡度、入渗时间以及养分含量共同确定，计算公式如下
infilrate＝infilratebareland+infilrateplant (3)
infilratebareland＝cos(s)×(1.91+0.326×e-0.215infil time) (4)
其中infilrate为所有元胞的土壤的入渗速率，infilratebareland与infilrateplant分别为裸土和植被覆盖区域的土壤的入渗速率，infiltime为入渗时间，nutrition为土壤养分含量，inutrition为土壤的初始养分含量。
3)利用元胞数组径流量入渗量infil记录每分钟内的入渗量，并判断当前元胞发生土壤入渗后的地表径流量：
当通过降雨获取的入渗量infil小于元胞此时的入渗速率时，所有降水将转化为土壤水；
当降雨量大于元胞此时的入渗速率时，土壤水的增加量等于单位时间内土壤的入渗量，其余的雨水则会滞留在地表，从而形成地表径流，此外，生成一个n*n大小的元胞数组的累计入渗量suminfil，积累每个元胞累计的入渗量，当累计入渗量超过土壤的饱和含水量时，土壤入渗随即停止。将植被覆盖和裸土区域的土壤饱和含水率分别设定为0.35和0.3，因此，土壤入渗量具体公式为
infil＝infilrate×[(suminfil＜0.3)+(suminfil＜0.35)] (6)
其中，其中，infilrate为入渗速率，suminfil为累计入渗量。
然后，使用降雨量减去土壤入渗量，随机得到土壤入渗后残留在地表的雨量，具体公式为：
rain′＝rain-infil (7)
其中，rain’为土壤入渗后的降雨量，rain为降雨量，infil为入渗量速率。
进一步地，步骤S3的实现方式如下：
1)将降雨获取的径流量进行转化，根据勾股定理可得，坡度为s的元胞的坡面实际单位面积为1/cos(s)，则土壤入渗前的地表径流量runoff的计算公式为：
funoff′＝rain′ (9)
其中runoff’为未修正前的土壤入渗后的地表径流量；
2)根据每个元胞与周边元胞之间的高程关系判断径流的流入或流出情况；
当元胞是相邻元胞中的相对最低点时，相邻元胞的地表径流将会汇聚于此元胞中，则此元胞的径流的流出量flowout等于0，而流入量则等于周边邻居元胞流出量的总和，具体公式为：
flowin(x，y)＝flowout(x-1，y-1)+flowout(x-1，y)+flowout(x-1，y+1)+flowout(x，y-1)+flowout(x，y+1)+flowout(x+1，y-1)+flowout(x+1，y)+flowout(x+1，y+1) (10)
其中，flowin(x，y)为第x行第y列的元胞的径流流入量，flowout(x，y)为第x行第y列的元胞的径流流出量；当元胞不是相对最低点时，此元胞所有的地表径流将会流入其邻居元胞中高程最小的元胞，此元胞的径流的流入量flowin等于0，而其径流的流出量与其径流量以及径流速率有关，而径流的流入量则等于其邻居元胞的流向元胞C(x,y)的径流的流出量总和；
当元胞内的径流量小于元胞的径流速率时，径流流出量等于此时的径流量即
flowout＝runoff (11)
当元胞内的径流量大于元胞的径流速率时，径流流出量则等于元胞的径流速率乘以单位时间，即
flowout＝runoffrate (12)
其中，元胞C(x,y)的径流速率表示了单位时间内通过元胞的最大径流量，所述径流速率通过曼宁公式计算得到；所述曼宁公式如下：
其中k是转换系数，国际标准制为1；n为曼宁粗糙系数；Rh为水力半径，计算方式是过水截面积与湿周长的比值；S为水力坡线或是线性扬程损失的斜率；
3)根据地表覆盖情况确定不同元胞的径流速率，定义裸土的粗糙系数为0.02，植被覆盖土壤的粗糙系数为0.095；斜率S则等于坡度的正切值tan(s)，则地表径流速率的具体公式为
runoffrate＝runoffrateplant+runoffratebareland (14)
runoffrateplant＝(cos(s).*runoff.^(2/3).*(tan(s).^0.5))./0.095 (15)
runoffratebareland＝(cos(s).*runoff.^(2/3).*(tan(s).^0.5))./0.02 (16)
其中runoffrate为所有元胞的径流速率，runoffrateplant为植被覆盖区域的径流速率，runoffratebareland为裸土的径流速率；
最终，当发生流入和流出过程后的径流量的计算公式为
runoff＝runoff″+flowin-flowout (17)
其中runoff”代表土壤入渗后但残余的地表水尚未流动的径流量，flowin是径流流入量，flowout是径流流出量。
进一步地，步骤S4的实现方式如下：
1)根据每个元胞的地表径流和覆盖类型确定土壤养分流失或积累，及其速率；
对于地表径流的流出量大于流入量的元胞，径流此时主要通过携带裹挟带走土壤养分；对于地表径流的流入量大于流出量的元胞，径流主要将携带其他元胞的土壤养分汇集于此元胞；
结合覆盖类型的差异，土壤养分流失和积累的速率随之变化，存在植被覆盖的元胞具有更强的阻拦和保持土壤的能力，因此土壤养分流失的速率较慢、土壤养分积累的速率较快；裸土的元胞则正好相反；参照以往研究中已经确定的土壤养分流失速率与土壤类型之间确定的经验公式参数后，
可得到土壤养分变化速率的计算公式如下：
①当径流流出量>＝径流流入量，根据植被覆盖的元胞土壤养分变化速率为：
裸土的元胞的土壤养分变化速率为：
②当径流流出量<径流流入量，有植被覆盖的元胞土壤养分变化速率为：
裸土的元胞的土壤养分变化速率为：
2)在确定土壤养分变化速率后，使用土壤养分元胞数组减去土壤养分变化元胞数组，从而得到当前时刻的土壤养分情况，公式如下：
nutrition＝nutrition+nutrichange (22)
采取上述技术方案，具有的有益效果如下：
(1)本发明对降雨、土壤入渗与径流等水文过程与土壤养分、坡度等地形条件间的相互作用关系进行梳理，提高了土壤养分动态变化模拟的鲁棒性。
(2)本发明的仿真模拟模型及仿真模拟的过程阐明了土壤养分在降雨、土壤入渗与径流作用下的变化过程，实现了对土壤养分动态变化的全过程反映。
(3)本发明能够对区域内部每块土壤与其周边区域在径流作用下的养分流失或聚集关系进行确定，可实现对土壤养分空间分布的全区域预测。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明仿真模拟方法的流程图；
图2为本发明实施例的相对最低点元胞的地表径流流入示意图；
图3为本发明实施例的非相对最低点元胞的地表径流流出示意图；
图4为本发明实施例的降雨径流模拟过程中的效果图；
图5为本发明实施例的土壤养分流失仿真模拟效果示意图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
下面结合附图1至附图5对本发明做进一步说明。
实施例：本发明选取山西省大同市作为研究区，对其土壤养分流失情况进行模拟仿真，其具体步骤如下：
S1：根据模拟区域实际大小、植被覆盖、降雨情况建立元胞数组，通过随机生成小于降雨概率的元胞数组模拟降雨事件的发生；
首先在Matlab中建立若干个大小为n*n的矩阵，矩阵中的每个元素代表一个元胞，每个矩阵则代表一个元胞数组；所述元胞数组的数量等于所涉及的元胞状态参数的数量，每个元胞数组内的元胞值代表大小为n*n m2的区域内每块土壤所对应参数的值；涉及到元胞状态参数包括高程t、元胞的坡度s、降雨量rain、入渗量infil、入渗速率infilrate、径流流出量flowout、径流流入量flowin、邻居元胞径流流出量surflowout、径流速率runoffrate、土壤养分含量nutrition以及土壤养分变化量nutrichange。
在此步骤的具体实现过程如下：
1)降雨发生时，使用生成值域区间内随机数的方式，确定此次降雨强度arain，具体公式为
其中sumrain为降雨总量，rainday为降雨天数，raintime为降雨持续时间；
其中根据大同市过去50年的降雨统计数据，确定大同市的降雨总量sumrain在220-280mm之间、降雨天数rainday在20-40天之间、降雨时间raintime在60-720min之间、雨水落入概率prain在0.2-0.5之间，根据公式(1)可计算出大同市的降雨雨强。
2)使用元胞数组降雨量rain反映当前每个元胞内的降雨情况，对于生成的n*n大小的随机数元胞数组，小于降雨落入概率的元胞则认为该元胞存在降雨，存在降雨的元胞状态为1，其余为0，具体公式为
rain＝arain×(rand(n)<prain) (2)
其中rand(n)为大小为n*n的随机数矩阵，prain为雨水落入概率。
S2：根据各个元胞的植被覆盖情况确定土壤入渗速率，结合降雨量计算土壤入渗量和残留在地表的雨量；
此步骤的具体实施方式为：
1)对步骤S1得到的存在降雨的元胞数组进行判断，对于存在降雨的元胞，地表将会滞留等同于降雨量的地表水，随后发生土壤入渗；
2)计算发生土壤入渗的元胞的土壤的入渗速率，根据地表覆盖类型、坡度、入渗时间以及养分含量共同确定，计算公式如下：
infilrate＝infilratebareland+infilrateplant (3)
infilratebareland＝cos(s)×(1.91+0.326×e-0.215infiltime) (4)
其中infilrate为所有元胞的土壤的入渗速率，infilratebareland与infilrateplant分别为裸土和植被覆盖区域的土壤的入渗速率，infiltime为入渗时间，nutrition为土壤养分含量，inutrition为土壤初始养分含量。
3)利用元胞数组径流量infil记录每分钟内的入渗量，并判断当前元胞发生土壤入渗后的地表径流量；
当通过降雨获取的径流量infil小于元胞此时的入渗速率时，所有降水将转化为土壤水；
当降雨量大于元胞此时的入渗速率时，土壤水的增加量等于单位时间内土壤的入渗量，其余的雨水则会滞留在地表，从而形成地表径流，此外，生成一个n*n大小的元胞数组suminfil，积累每个元胞累计的入渗量，当累计入渗量超过土壤的饱和含水量时，土壤入渗随即停止，将植被覆盖和裸土区域的土壤饱和含水率分别设定为0.35和0.3，因此，土壤入渗量具体公式为：
infil＝infilrate×[(suminfil＜0.3)+(suminfil＜0.35)] (6)
其中，infilrate为入渗速率，suminfil为累计入渗量。
然后，使用降雨量减去土壤入渗量，随机得到土壤入渗后残留在地表的雨量，具体公式为：
rain′＝rain-infil (7)
其中，rain’为土壤入渗后的降雨量，rain为降雨量，infil为入渗量速率。
S3：通过坡度对地表径流量进行转换，根据高程关系确定每个元胞内径流的流入或流出关系；
此步骤的具体实现方式如下：
1)将降雨获取的径流量进行转化，根据勾股定理可得，坡度为s的元胞的坡面实际单位面积为1/cos(s)，则土壤入渗前地表径流量runoff的计算公式为
funoff′＝rain′ (9)
其中runoff’为未修正前的土壤入渗后的地表径流量；
2)根据每个元胞与周边元胞之间的高程关系判断径流的流入或流出情况；
当元胞是相邻元胞中的相对最低点时，相邻元胞的地表径流将会汇聚于此元胞中，则此元胞的径流的流出量flowout等于0，而流入量则等于周边邻居元胞流出量的总和，具体公式为
flowin(x，y)＝flowout(x-1，y-1)+flowout(x-1，y)+flowout(x-1，y+1)+flowout(x，y-1)+flowout(x，y+1)+flowout(x+1，y-1)+flowout(x+1，y)+flowout(x+1，y+1) (10)
其中，flowin(x，y)为第x行第y列的元胞的径流流入量，flowout(x，y)为第x行第y列的元胞的径流流出量；当元胞不是相对最低点时，此元胞所有的地表径流将会流入其邻居元胞中高程最小的元胞，此元胞的径流的流入量flowin等于0，而其径流的流出量与其径流量以及径流速率有关，而径流的流入量则等于其邻居元胞的流向元胞C(x，y)的径流的流出量总和；
当元胞内的径流量小于元胞的径流速率时，径流流出量等于此时的径流量即
flowout＝runoff (11)
当元胞内的径流量大于元胞的径流速率时，径流流出量则等于元胞的径流速率乘以单位时间，即
flowout＝runoffrate (12)
其中，元胞C(x，y)的径流速率表示了单位时间内通过元胞的最大径流量，所述径流速率通过曼宁公式计算得到；所述曼宁公式如下：
曼宁公式是一个估测液体在开放管道(即明渠流)或非满管流中平均速度的经验公式，其中k是转换系数，国际标准制为1；n为曼宁粗糙系数；Rh为水力半径，计算方式是过水截面积与湿周长的比值；S为水力坡线或是线性扬程损失的斜率；
3)将每个元胞作为径流流动的基本空间单元，则Rh为通过元胞的的水深，为了便于降雨量、径流量之间的水量交互与计算，对径流量与降雨量均统一使用mm作为单位，因此在计算最大径流速率的公式中将径流量直接作为水深计算；经过查询曼宁系数表后，定义裸土的粗糙系数为0.02，植被覆盖土壤的粗糙系数为0.095；斜率S则等于坡度的正切值tan(s)，则地表径流速率的具体公式为
runoffrate＝runoffrateplant+runoffratebareland (14)
runoffrateplant＝(cos(s).*runoff.^(2/3).*(tan(s).^0.5))./0.095 (15)
runoffratebareland＝(cos(s).*runoff.^(2/3).*(tan(s).^0.5))./0.02(16)
其中runoffrate为所有元胞的径流速率，runoffrateplant为植被覆盖区域的径流速率，runoffratebareland为裸土的径流速率；
最终，当发生流入和流出过程后的径流量的计算公式为
runoff＝runoff″+flowin-flowout (17)
其中runoff”代表土壤入渗后但残余的地表水尚未流动的径流量，flowin是径流流入量，flowout是径流流出量。
S4：根据元胞在径流过程中的径流的流入或流出关系判断土壤养分的流失或汇集，并依据植被覆盖确定土壤养分变化速率，最后计算径流发生后的土壤养分情况；
此步骤的具体实现方式如下：
1)根据每个元胞的地表径流和覆盖类型确定土壤养分流失或积累，及其速率；
对于地表径流的流出量大于流入量的元胞，径流此时主要通过携带裹挟带走土壤养分；对于地表径流的流入量大于流出量的元胞，径流主要将携带其他元胞的土壤养分汇集于此元胞；
结合覆盖类型的差异，土壤养分流失和积累的速率随之变化，存在植被覆盖的元胞具有更强的阻拦和保持土壤的能力，因此土壤养分流失的速率较慢、土壤养分积累的速率较快；裸土的元胞则正好相反；参照以往研究中已经确定的土壤养分流失速率与土壤类型之间确定的经验公式参数后，可得到土壤养分变化速率的计算公式如下：
①当径流流出量>＝径流流入量，有植被覆盖的元胞土壤养分变化速率为：
裸土的元胞的土壤养分变化速率为：
②当径流流出量<径流流入量，有植被覆盖的元胞土壤养分变化速率为：
裸土的元胞的土壤养分变化速率为：
2)在确定土壤养分变化速率后，使用土壤养分元胞数组减去土壤养分变化元胞数组，从而得到当前时刻的土壤养分情况，公式如下：
nutrition＝nutrition+nutrichange (22)
S5：建立循环函数重复步骤S1至S4，模拟降雨发生后不断发生的降雨、入渗、径流及土壤养分变化过程，当循环函数运行超过设定的降雨时间raintime后停止循环，并输出此时土壤养分流失情况，即对模拟区域的土壤养分流失情况完成了仿真模拟。
在本实例中，设置循环函数重复步骤S1-S4，模拟降雨发生后不断发生的降雨、土壤入渗、径流及土壤养分变化过程。建立一个初始值为1的计数器j，每完成步骤S1-S4一次该计数器便增加1，当j小于步骤S1中确定的降雨时间raintime时，持续循环S1-S4，如图4所示，图中灰色部分表示地表径流。随着降雨时间的持续发生，土壤入渗量逐渐下雨降雨量，地表径流逐渐增多，相邻元胞内的地表径流根据高程关系，从高处向低处流动。当j等于raintime时，便停止循环在降雨、入渗、径流过程的持续作用下，土壤养分也随之发生变化，当降雨发生8小时后，土壤养分发生变化的区域如图5所示，由此可以获取研究区域内土壤养分的变化情况及空间分布特征。
实施例2：在本实例中，对实施例1的模拟区域在一次降雨前后的土壤养分情况进行了实地调查，随机选择了若干个样点，以土壤有机质作为土壤养分的替代指标，土壤实际与模拟养分变化情况如表1所示。
表1：模拟结果验证
以上对本发明及其实施方式进行了描述，这种描述没有限制性，附图中所示的也只是本发明的实施方式之一，实际的结构并不局限于此。总而言之如果本领域的普通技术人员受其启示，在不脱离本发明创造宗旨的情况下，不经创造性的设计出与该技术方案相似的结构方式及实施例，均应属于本发明的保护范围。
Claims (1)
1.一种基于元胞自动机的土壤养分流失仿真模拟方法，其特征在于，
其包括如下步骤：
S1：根据模拟区域实际大小、植被覆盖、降雨情况建立元胞数组，通过随机生成小于降雨概率的元胞数组模拟降雨事件的发生；
S2：根据各个元胞的植被覆盖情况确定土壤入渗速率，结合降雨量计算土壤入渗量和残留在地表的雨量；
S3：通过坡度对地表径流量进行转换，根据高程关系确定每个元胞内径流的流入或流出关系；
S4：根据元胞在径流过程中的径流的流入或流出关系判断土壤养分的流失或汇集，并依据植被覆盖确定土壤养分变化速率，最后计算径流发生后的土壤养分情况；
S5：建立循环函数重复步骤S1至S4，模拟降雨发生后不断发生的降雨、入渗、径流及土壤养分变化过程，当循环函数运行超过设定的降雨时间raintime后停止循环，并输出此时土壤养分流失情况，即对模拟区域的土壤养分流失情况完成了仿真模拟；
步骤S1中的元胞数组的构建过程为：在Matlab中建立若干个大小为n*n的矩阵，矩阵中的每个元素代表一个元胞，每个矩阵则代表一个元胞数组；
所述元胞数组的数量等于所涉及的元胞状态参数的数量，每个元胞数组内的元胞值代表大小为n*n m2的区域内每块土壤所对应参数的值；涉及到元胞状态参数包括高程t、元胞的坡度s、降雨量rain、入渗量infil、入渗速率infilrate、径流流出量flowout、径流流入量flowin、邻居元胞径流流出量surflowout、径流速率runoffrate、土壤养分含量nutrition以及土壤养分变化量nutrichange；
步骤S1的实现方式如下：
1)降雨发生时，使用生成值域区间内随机数的方式，确定此次降雨强度arain，具体公式为
其中sumrain为降雨总量，rainday为降雨天数，raintime为降雨持续时间；
2)使用元胞数组降雨量rain反映当前每个元胞内的降雨情况，对于生成的n*n大小的随机数元胞数组，小于降雨落入概率的元胞则认为该元胞存在降雨，存在降雨的元胞状态为1，其余为0，具体公式为
rain＝arain×(rand(n)<prain) (2)
其中rand(n)为大小为n*n的随机数矩阵，prain为雨水落入概率；
步骤S2的实现方式如下：
1)对步骤S1得到的存在降雨的元胞数组进行判断，对于存在降雨的元胞，地表将会滞留等同于降雨量的地表水，随后发生土壤入渗；
2)计算发生土壤入渗的元胞的土壤的入渗速率，根据地表覆盖类型、坡度、入渗时间以及养分含量共同确定，计算公式如下
infilrate＝infilratebareland+infilrateplant (3)
infilratebareland＝cos(s)×(1.91+0.326×e-0.215infiltime) (4)
其中infilrate为所有元胞的土壤的入渗速率，infilratebareland与infilrateplant分别为裸土和植被覆盖区域的土壤的入渗速率，infiltime为入渗时间，nutrition为土壤养分含量，inutrition为土壤的初始养分含量；
3)利用元胞数组径流量入渗量infil记录每分钟内的入渗量，并判断当前元胞发生土壤入渗后的地表径流量：
当通过降雨获取的入渗量infil小于元胞此时的入渗速率时，所有降水将转化为土壤水；
当降雨量大于元胞此时的入渗速率时，土壤水的增加量等于单位时间内土壤的入渗量，其余的雨水则会滞留在地表，从而形成地表径流，此外，生成一个n*n大小的元胞数组的累计入渗量suminfil，积累每个元胞累计的入渗量，当累计入渗量超过土壤的饱和含水量时，土壤入渗随即停止；将植被覆盖和裸土区域的土壤饱和含水率分别设定为0.35和0.3，因此，土壤入渗量具体公式为
infil＝infilrate×[(suminfil<0.3)+(suminfil<0.35)] (6)
其中，其中，infilrate为入渗速率，suminfil为累计入渗量；
然后，使用降雨量减去土壤入渗量，随机得到土壤入渗后残留在地表的雨量，具体公式为：
rain′＝rain-infil (7)
其中，rain’为土壤入渗后的降雨量，rain为降雨量，infil为入渗量速率；
步骤S3的实现方式如下：
1)将降雨获取的径流量进行转化，根据勾股定理可得，坡度为s的元胞的坡面实际单位面积为1/cos(s)，则土壤入渗前的地表径流量runoff的计算公式为：
runoff′＝rain′ (9)
其中runoff’为未修正前的土壤入渗后的地表径流量；
2)根据每个元胞与周边元胞之间的高程关系判断径流的流入或流出情况；
当元胞是相邻元胞中的相对最低点时，相邻元胞的地表径流将会汇聚于此元胞中，则此元胞的径流流出量flowout等于0，而流入量则等于周边邻居元胞流出量的总和，具体公式为：
flowin(x，y)＝flowout(x-1，y-1)+flowout(x-1，y)+flowout(x-1，y+1)+flowout(xy-1)+flowout(x，y+1)+flowout(x+1，y-1)+flowout(x+1，y)+flowout(x+1，y+1) (10)
其中，flowin(x，y)为第x行第y列的元胞的径流流入量，flowout(x，y)为第x行第y列的元胞的径流流出量；当元胞不是相对最低点时，此元胞所有的地表径流将会流入其邻居元胞中高程最小的元胞，此元胞的径流的流入量flowin等于0，而其径流的流出量与其径流量以及径流速率有关，而径流的流入量则等于其邻居元胞的流向元胞C(x，y)的径流的流出量总和；
当元胞内的径流量小于元胞的径流速率时，径流流出量等于此时的径流量即
flowout＝runoff (11)
当元胞内的径流量大于元胞的径流速率时，径流流出量则等于元胞的径流速率乘以单位时间，即
flowout＝runoffrate (12)
其中，元胞C(x，y)的径流速率表示了单位时间内通过元胞的最大径流量，所述径流速率通过曼宁公式计算得到；所述曼宁公式如下：
其中k是转换系数，国际标准制为1；n为曼宁粗糙系数；Rh为水力半径，计算方式是过水截面积与湿周长的比值；S为水力坡线或是线性扬程损失的斜率；
3)根据地表覆盖情况确定不同元胞的径流速率，定义裸土的粗糙系数为0.02，植被覆盖土壤的粗糙系数为0.095；斜率S则等于坡度的正切值tan(s)，则地表径流速率的具体公式为
runoffrate＝runoffrateplant+runoffratebareland (14)
runoffrateplant＝(cos(s).*runoff.^(2/3).*(tan(s).^0.5))./0.095 (15)
runoffratebareland＝(cos(s).runoff·^(2/3).*(tan(s).^0.5))./0.02 (16)
其中runoffrate为所有元胞的径流速率，runoffrateplant为植被覆盖区域的径流速率，runoffratebareland为裸土的径流速率；
最终，当发生流入和流出过程后的径流量的计算公式为
runoff＝runoff″+flowin-flowout (17)
其中runoff”代表土壤入渗后但残余的地表水尚未流动的径流量，flowin是径流流入量，flowout是径流流出量；
步骤S4的实现方式如下：
1)根据每个元胞的地表径流和覆盖类型确定土壤养分流失或积累，及其速率；
对于地表径流的流出量大于流入量的元胞，径流此时主要通过携带裹挟带走土壤养分；对于地表径流的流入量大于流出量的元胞，径流主要将携带其他元胞的土壤养分汇集于此元胞；
结合覆盖类型的差异，土壤养分流失和积累的速率随之变化，存在植被覆盖的元胞具有更强的阻拦和保持土壤的能力，因此土壤养分流失的速率较慢、土壤养分积累的速率较快；裸土的元胞则正好相反；参照以往研究中已经确定的土壤养分流失速率与土壤类型之间确定的经验公式参数后，
可得到土壤养分变化速率的计算公式如下：
①当径流流出量>＝径流流入量，根据植被覆盖的元胞土壤养分变化速率为：
裸土的元胞的土壤养分变化速率为：
②当径流流出量<径流流入量，有植被覆盖的元胞土壤养分变化速率为：
裸土的元胞的土壤养分变化速率为：
2)在确定土壤养分变化速率后，使用土壤养分元胞数组减去土壤养分变化元胞数组，从而得到当前时刻的土壤养分情况，公式如下：
nutrtion＝nutrition+nutrichange (22)。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111059917.5A CN113901640B (zh) | 2021-09-10 | 2021-09-10 | 一种基于元胞自动机的土壤养分流失仿真模拟方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111059917.5A CN113901640B (zh) | 2021-09-10 | 2021-09-10 | 一种基于元胞自动机的土壤养分流失仿真模拟方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113901640A CN113901640A (zh) | 2022-01-07 |
| CN113901640B true CN113901640B (zh) | 2024-05-28 |
Family
ID=79027519
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111059917.5A Active CN113901640B (zh) | 2021-09-10 | 2021-09-10 | 一种基于元胞自动机的土壤养分流失仿真模拟方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113901640B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116453583B (zh) * | 2023-03-16 | 2024-05-31 | 宁夏蓝怡生物工程有限公司 | 一种lucc的元胞自动机矫正的蛋白核藻类变化模拟方法 |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103984839A (zh) * | 2014-05-29 | 2014-08-13 | 河南理工大学 | 一种基于元胞自动机的地表下沉仿真方法及系统 |
| KR20160085990A (ko) * | 2015-01-08 | 2016-07-19 | 강릉원주대학교산학협력단 | 토양침식 실험장치 |
| KR20180000619A (ko) * | 2016-06-23 | 2018-01-03 | (주)해동기술개발공사 | Ｇｉｓ 기반 토사유실 평가 방법 |
| CN107860703A (zh) * | 2017-11-02 | 2018-03-30 | 北京师范大学 | 基于室外降雨模拟和室内土壤测试的土壤侵蚀特性分析方法 |
| CN107894498A (zh) * | 2017-11-02 | 2018-04-10 | 北京师范大学 | 基于室外降雨模拟和室内土壤测试的土壤理化演变分析方法 |
| CN108875242A (zh) * | 2018-06-29 | 2018-11-23 | 集美大学 | 一种城市元胞自动机情景模拟方法、终端设备及存储介质 |
| CN110415346A (zh) * | 2019-07-10 | 2019-11-05 | 华中师范大学 | 利用面向对象的三维元胞自动机进行水土流失模拟的方法 |
-
2021
- 2021-09-10 CN CN202111059917.5A patent/CN113901640B/zh active Active
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103984839A (zh) * | 2014-05-29 | 2014-08-13 | 河南理工大学 | 一种基于元胞自动机的地表下沉仿真方法及系统 |
| KR20160085990A (ko) * | 2015-01-08 | 2016-07-19 | 강릉원주대학교산학협력단 | 토양침식 실험장치 |
| KR20180000619A (ko) * | 2016-06-23 | 2018-01-03 | (주)해동기술개발공사 | Ｇｉｓ 기반 토사유실 평가 방법 |
| CN107860703A (zh) * | 2017-11-02 | 2018-03-30 | 北京师范大学 | 基于室外降雨模拟和室内土壤测试的土壤侵蚀特性分析方法 |
| CN107894498A (zh) * | 2017-11-02 | 2018-04-10 | 北京师范大学 | 基于室外降雨模拟和室内土壤测试的土壤理化演变分析方法 |
| CN108875242A (zh) * | 2018-06-29 | 2018-11-23 | 集美大学 | 一种城市元胞自动机情景模拟方法、终端设备及存储介质 |
| CN110415346A (zh) * | 2019-07-10 | 2019-11-05 | 华中师范大学 | 利用面向对象的三维元胞自动机进行水土流失模拟的方法 |
Non-Patent Citations (1)
| Title |
|---|
| 模拟降雨条件下秦岭北麓土壤磷素流失特征;陈曦;王雪松;贺京哲;代允超;刘衡;吕家珑;;水土保持学报;20160415(第02期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113901640A (zh) | 2022-01-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Zhang et al. | A comprehensive assessment framework for quantifying climatic and anthropogenic contributions to streamflow changes: A case study in a typical semi-arid North China basin | |
| Zhu et al. | Sediment flux sensitivity to climate change: A case study in the Longchuanjiang catchment of the upper Yangtze River, China | |
| CN113610264B (zh) | 一种精细化电网台风洪涝灾害预测系统 | |
| Shen et al. | Relating landscape characteristics to non-point source pollution in a typical urbanized watershed in the municipality of Beijing | |
| CN109614655B (zh) | 一种河流径流量的研究方法 | |
| Yang et al. | Long-time series ecological environment quality monitoring and cause analysis in the Dianchi Lake Basin, China | |
| Miao et al. | Spatio-temporal variability of streamflow in the Yellow River: possible causes and implications | |
| CN109409694A (zh) | 基于实测值的地块海绵城市绩效考核指标计算方法 | |
| Gao et al. | Attribution analysis of climatic and multiple anthropogenic causes of runoff change in the Loess Plateau—A case‐study of the Jing River Basin | |
| Zakizadeh et al. | A novel study of SWAT and ANN models for runoff simulation with application on dataset of metrological stations | |
| Yang et al. | Integrating an hourly weather generator with an hourly rainfall SWAT model for climate change impact assessment in the Ru River Basin, China | |
| CN111737853A (zh) | 一种基于swmm模型的低影响开发多目标区间优化配置方法 | |
| Chen et al. | Assimilating multi-source data into a three-dimensional hydro-ecological dynamics model using Ensemble Kalman Filter | |
| CN113901640B (zh) | 一种基于元胞自动机的土壤养分流失仿真模拟方法 | |
| CN110728062A (zh) | 一种基于swmm的农村非点源污染模拟方法 | |
| CN117491585A (zh) | 基于时序网络的水生态污染监测方法及装置、系统 | |
| Tian et al. | Preemptive warning and control strategies for algal blooms in the downstream of Han River, China | |
| Zhang et al. | Quantification of human and climate contributions to multi-dimensional hydrological alterations: A case study in the Upper Minjiang River, China | |
| CN106156529A (zh) | 一种基于数理统计的下垫面径流关系确定方法 | |
| Gillies | Managing the effect of infiltration variability on the performance of surface irrigation | |
| CN112541611B (zh) | 一种雨养农用地面源污染排放量预测方法及系统 | |
| CN116955516A (zh) | Gis和swat结合的流域多要素控制单元划分方法 | |
| Liu et al. | Natural and reservoir-induced channel changes in the Yangtze River Tidal Reach | |
| Li et al. | Variations in water conservation function and attributions in the Three-River Source Region of the Qinghai–Tibet Plateau based on the SWAT model | |
| Kotei et al. | Estimation of flow-duration and low-flow frequency parameters for the Sumanpa Stream at Mampong-Ashanti in Ghana for the 1985-2009 period |
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 |