一种喀斯特地区岩石风化成土速率计算的方法与流程

本发明属于土壤学技术领域，尤其涉及一种喀斯特地区岩石风化成土速率计算的方法。

背景技术：

土壤侵蚀预测模型是进行水土流失预测、水土保持规划的有效手段，因此预测值的准确性直接关系到区域生态评估和安全。而喀斯特地区受岩溶环境制约，成土十分缓慢，因此该地区的开展该地区的土壤侵蚀模数的预测时需要结合其土壤允许流失量进行综合考虑。

首先，基于rusle土壤侵蚀模型和gis、rs技术，获取土壤侵蚀各因子，包括降雨侵蚀力因子(r)、土壤可蚀性因子(k)、地形因子(ls)、植被覆盖与作物管理因子(c)以及土壤保持措施因子(p)，并计算土壤侵蚀模数；然而，rusle模型主要是针对国内外缓坡条件开发而成的，并不完全适用于我国地形地貌复杂的喀斯特地区；也就是说通过常规的rusle土壤侵蚀模型来对喀斯特地区土壤水土流失预测、水土保持规划设计等存在误差大，精度低等问题。因此需要基于土壤允许流失量(等同于喀斯特地区岩石风化成土速率)对土壤侵蚀模数进行校正。

技术实现要素：

本发明要解决的技术问题是：提供一种喀斯特地区岩石风化成土速率计算的方法，以解决现有技术通过rusle土壤侵蚀模型对喀斯特地区土壤水土流失预测、水土保持规划设计等存在误差大，精度低等问题。

本发明的技术方案是：

一种喀斯特地区岩石风化成土速率计算的方法，它包括：

步骤1、采集数据及数据预处理；采集包含日降水量数据、日地表最高及最低温度、月总蒸散发数据、土壤质地与有机质含量、数字高程模型、归一化植被指数、土地利用类型以及岩性数据、碳酸盐岩溶蚀速率、酸不溶物含量百分比、碳酸盐岩含量百分比、碳酸盐岩密度和非碳酸盐岩成土速率，并将所有数据投影到同一坐标系中；

步骤2、土壤侵蚀模数的计算：基于修正的通用土壤流失方程(rusle模型)，通过对各因子的计算，得到喀斯特地区的理论侵蚀模数；

步骤3、碳酸盐岩最大溶蚀速率的计算：基于气象水文及地球化学时空数据集，利用碳酸盐岩热力学溶蚀模型，计算碳酸盐岩最大溶蚀速率；

步骤4、喀斯特地区岩石风化成土速率计算：利用碳酸盐岩最大溶蚀速率、酸不溶物含量百分比、碳酸盐岩含量百分比、碳酸盐岩密度和非碳酸盐岩成土速率，计算喀斯特不同岩层组合类型的岩石风化成土速率。

所述的一种喀斯特地区岩石风化成土速率计算的方法，它还包括：步骤5、土壤侵蚀模数检验：将步骤2的结果与步骤4的结果相减，得到一个包含0及正负值的图层，大于0的部分表示该区域理论土壤侵蚀量大于其实际成土量；小于等于0的部分表示该区域理论土壤侵蚀量小于等于实际成土量；

步骤6、土壤侵蚀模数校正：利用岩石风化成土速率图层对步骤5中大于0的部分进行土壤侵蚀模数校正。

步骤2所述土壤侵蚀模数的计算公式为：a＝r×k×l×s×c×p，其中a为土壤侵蚀模数、r为降雨侵蚀力因子、k为土壤可蚀性因子、l为坡长因子、s为坡度因子、c为地表植被覆盖与管理因子、p为土壤保持措施因子。

步骤3所述碳酸盐岩最大溶蚀速率的计算方法为：

步骤3.1、建立碳酸盐岩溶蚀平衡反应式：

步骤3.2、构建碳酸盐岩热力学溶蚀模型：

式中：dmax为在该平衡反应下碳酸盐岩最大的溶蚀速率,p,e分别为降雨及蒸散发总量，ks为方解石溶度积常数，k1为co2水化并离解为hco3^-的平衡常数，k0为co2溶于水的平衡常数，k2为co3^2-形成的平衡常数，为溶液中ca²⁺的活度系数，pco2为土壤或含水层中co2的分压；

所述的ki因子的计算公式为：

log(ks)＝-171.91-0.078tk+2839.32/tk+71.59log(tk)

k0＝1.7×10^-4/k1

式中：tk为开尔文温度；

所述的ca²⁺离子活度系数因子的计算公式为：

式中：a和b取决于摄氏温度t(℃),ai为离子半径(ca²⁺半径为),zi

为离子电荷数，i为离子强度；

所述的a和b因子的计算公式为：

a＝0.4883+8.074×10^-4t

b＝0.3241+1.6×10^-4t

所述的i因子的计算公式为:

式中：ci为离子浓度(mol/l)；

所述的pco2因子的计算公式为:

log(pco2)＝-3.47+2.09×(1-e^-0.00172e)。

步骤4所述喀斯特地区岩石风化成土速率计算公式为：

wi＝v·q·ρ·m+n·(1-m)

式中：wi为岩石风化成土速率，单位t·km^‐2·yr^‐1；v为碳酸盐岩溶蚀速率，单位mm/a，换算为m³·km^‐2·yr^‐1；q为酸不溶物含量，单位％；m为碳酸盐岩含量，单位％；ρ为碳酸盐岩密度，单位t/m³；n为非碳酸盐岩成土速率，单位t·km^‐2·yr^‐1。

步骤5和6所述的土壤侵蚀模数检验和土壤侵蚀模数校正的方法为：(1)将步骤2的结果与步骤4的结果相减，得到差值图层；(2)将差值图层里值大于0的部分提取出来作为待校正图层，将差值图层里值小于等于0的部分保留；用待校正图层对喀斯特岩石风化成土速率图层进行掩膜提取；对待校正图层进行赋值，得到校正赋值图层；(3)将掩膜提取的结果与差值图层里值小于等于0的部分进行拼接融合，得到校正后的喀斯特区域土壤侵蚀模数图层。

本发明有益效果：

本发明充分考虑不同区域碳酸盐岩岩石的实际溶蚀特征，并以此对不同区域的碳酸盐岩溶蚀速率进行了空间化；充分考虑到喀斯特地区复杂的地质背景及成土过程，并基于土壤允许流失量(等同于喀斯特地区岩石风化成土速率)对土壤侵蚀模数进行校正，使岩溶生态脆弱区土壤侵蚀的估算更加准确，解决了现有技术通过rusle土壤侵蚀模型对喀斯特地区土壤水土流失预测、水土保持规划设计等存在误差大，精度低等问题；对于土壤侵蚀的防治及石漠化治理和生态恢复与重建具有重要的理论指导意义。

具体实施方式

本发明提出一种喀斯特地区岩石风化成土速率计算的方法，它包括如下步骤：

第一，获取数据与预处理。计算土壤侵蚀模数所需要的数据分别为：日降水量数据，土壤质地与有机质含量，数字高程模型(dem)，归一化植被指数(ndvi)以及土地利用类型。计算碳酸盐岩成土速率需要的数据有岩性矢量数据、碳酸盐岩溶蚀速率，酸不溶物含量百分比，碳酸盐岩含量百分比，碳酸盐岩密度，非碳酸盐岩成土速率。获取相关数据后对数据进行重投影，使其具有统一的投影系统与坐标；对数据进行初运算使其符合模型运算要求。

第二，计算理论土壤侵蚀模数。基于rusle模型，其公式如下：

a＝r×k×l×s×c×p

式中a为土壤侵蚀模数(t·ha^‐1·yr^‐1)；r为降雨侵蚀力因子(mj·mm·ha^‐1·h^‐1·yr^‐1)；k为土壤可蚀性因子(t·h·mj^‐1·mm^‐1)；l为坡长因子；s为坡度因子；c为地表植被覆盖与管理因子；p为水土保持措施因子。通过该公式计算得到喀斯特地区的理论侵蚀模数；

第三，计算碳酸盐岩最大溶蚀速率。基于气象水文及地球化学时空数据集，利用碳酸盐岩热力学溶蚀模型，计算碳酸盐岩最大溶蚀速率图层，基于以下碳酸盐岩溶蚀平衡反应：

构建如下碳酸盐岩热力学溶蚀模型：

其中dmax为在该平衡反应下碳酸盐岩最大的溶蚀速率,p,e分别为降雨及蒸散发总量.ks为方解石溶度积常数,k1为co2水化并离解为hco3^-的平衡常数,k0为co2溶于水的平衡常数,k2为co3^2-形成的平衡常数.为溶液中ca²⁺的活度系数。pco2为土壤或含水层中co2的分压。

其中，ki因子采用以下模型计算：

log(ks)＝-171.91-0.078tk+2839.32/tk+71.59log(tk)

k0＝1.7×10^-4/k1

其中tk为开尔文温度。

ca²⁺离子活度系数因子采用以下模型计算：

其中，a和b取决于摄氏温度t(℃),ai为离子半径(ca²⁺半径为),zi为离子电荷数，i为离子强度.

a和b因子采用以下模型计算：

a＝0.4883+8.074×10^-4t

b＝0.3241+1.6×10^-4t

i因子采用以下模型计算:

其中，ci为离子浓度(mol/l).

pco2因子采用以下模型计算:

log(pco2)＝-3.47+2.09×(1-e^-0.00172e)

第四，计算喀斯特地区岩石风化成土速率。根据喀斯特地区岩性图，通过模型公式计算得到喀斯特地区不同岩层组合类型的岩石风化成土速率图层，其计算公式如下：

wi＝v·q·ρ·m+n·(1-m)

式中：wi为岩石风化成土速率，t·km^‐2·yr^‐1；v为碳酸盐岩溶蚀速率，mm/a，换算为m³·km^‐2·yr^‐1；q为酸不溶物含量，％；m为碳酸盐岩含量，％；ρ为碳酸盐岩密度，t/m³；n为非碳酸盐岩成土速率，t·km^‐2·yr^‐1，以500t·km^‐2·yr^‐1计算(sl190－2007.土壤侵蚀分类分级标准[s])。

第五，喀斯特区域理论土壤侵蚀模数的检验。用喀斯特土壤侵蚀模数图层减去喀斯特成土速率图层，得到一个包含0及正负值的差值图层。差值大于0的部分表示，该区域理论土壤侵蚀量大于其实际成土量。而喀斯特地区的最大土壤侵蚀量不会超过其成土量，因此，需要对差值大于0的部分进行校正处理。差值小于等于0的部分表示，该区域理论土壤侵蚀量小于等于实际成土量，不需要进行校正处理。

第六，理论土壤侵蚀的校正处理。利用岩石风化成土速率图层对喀斯特区域差值大于0的部分进行土壤侵蚀模数校正。将差值图层里值大于0的部分提取出来作为待校正图层，将差值图层里值小于等于0的部分保留；用待校正图层对喀斯特岩石风化成土速率图层进行掩膜提取；对待校正图层进行赋值，得到校正赋值图层；将掩膜提取的结果与差值图层里值小于等于0的部分进行拼接融合，得到校正后的喀斯特区域土壤侵蚀模数图层。