CN111274545B - 一种栅格尺度基于地形地貌的多模式产流计算方法 - Google Patents

Abstract

本发明公开了一种栅格尺度基于地形地貌的多模式产流计算方法，包括以下步骤：获取流域DEM数据并计算栅格的河面落差、坡度和地形湿度指数值；计算三个指标的累积分布函数，将非水体栅格分为洼地、坡地、高地超渗区和高地下渗区四种下垫面类型，得到四种类型的先验概率；利用最大似然分类法迭代计算各个栅格对每种下垫面类型的隶属度；将每种下垫面类型的产流作为一种产流模式，根据每个栅格的下垫面类型隶属度判别产流的模式和组成，将各产流模式对应的产流量加权叠加得到最终产流量。本发明通过对流域下垫面进行分类并基于栅格对各产流模式的隶属度加权计算得到最终产流量，减少了人为确定分类阈值带来的主观性，促进分布式产流模型的发展。

Description

一种栅格尺度基于地形地貌的多模式产流计算方法

技术领域

本发明属于水文预报领域，具体涉及一种栅格尺度基于地形地貌的产流模式最大似然判别与产流计算方法。

背景技术

产流是流域水文过程的重要组成部分，因此产流计算也是水文预报中的关键环节。降雨和下垫面因子的空间异质性导致了产流的空间不均匀性，而传统的集总式产流模型无法描述这种空间分布的复杂性。因此，分布式产流模型的构建是当前国内外研究的前沿和热点问题。依据地形、地质、土地利用、土壤类型等下垫面因子的水文和水力特征，分布式产流模型将流域划分为异质性空间计算单元，从而更加客观、合理地刻画产流分布。利用主导产流模式的识别划分流域水文单元是其中的方法之一，但人为确定分类阈值存在较强的主观性与经验性，且下垫面类型与产流模式的一对一构建无法刻画产流的复杂性。为此，需要设计一种新的技术方案给予解决。

发明内容

本发明的目的是提供一种栅格尺度基于地形地貌的产流模式最大似然判别与产流计算方法，利用数理统计原理计算流域栅格对各下垫面类型的隶属度，构建不同下垫面类型下的概念性产流模型，将不同模式下的产流量加权叠加得到各栅格的最终产流量，从而解决现有技术中存在的前述问题。

为了实现上述目的，本发明具体采用以下技术方案：

一种栅格尺度基于地形地貌的多模式产流计算方法，其特征在于，包括以下步骤：

步骤1，获取流域DEM数据并在此基础上计算栅格的河面落差(H)、坡度(Slope)和地形湿度指数(TWI)值；

步骤2，计算三个指标的累积分布函数，将非水体栅格分为洼地、坡地、高地超渗区和高地下渗区四类，得到四种类型的先验概率；

步骤3，利用最大似然分类法迭代得到各个栅格对每种下垫面类型的隶属度；

步骤4，在不同类型下垫面上建立不同的产流模式，将各类型对应的产流量根据下垫面类型隶属度加权叠加得到最终产流量。

所述步骤1中河面落差、坡度、地形湿度指数分别为：

H＝d-d₀ (1)

式中，H为河面落差，Slope为坡度，TWI为地形湿度指数，d为栅格所在高程，d₀为沿着经过该栅格的流线方向所到达的第一个水体栅格的高程，

和

分别为栅格在水平和垂直方向的高程变化率，a为栅格单宽累积面积。

所述步骤2中计算三个指标的累积分布函数，将非水体栅格分为洼地、坡地、高地超渗区和高地下渗区四类，得到四种类型的先验概率，具体步骤包括：

步骤21，提取除水体外所有栅格的H、Slope和TWI值，分别从大到小排列得到各自的累积分布函数，从而得到栅格各指标的概率分布：

F_H(x)＝n₁/N_t (4)

F_Slope(x)＝n₂/N_t (5)

F_TWI(x)＝n₃/N_t (6)

式中，F_H(x)，F_Slope(x)，F_TWI(x)分别为该栅格H、Slope和TWI值的累积分布函数，n₁，n₂，n₃分别为H、Slope、TWI值小于等于x的非水体栅格数，N_t为非水体栅格总数。

步骤22，计算非水体栅格对洼地、坡地、高地超渗区和高地下渗区四种类型的隶属度,将栅格类型定义为隶属度最高者对应的类型：

M_w＝1-F_H (7)

M_h＝F_H×F_Slope (8)

M_pe＝[F_H×(1-F_Slope)]×(1-F_TWI) (9)

M_pp＝[F_H×(1-F_Slope)]×F_TWI (10)

式中，M_w、M_h、M_pe、M_pp分别为栅格对洼地、坡地、高地超渗区和高地下渗区的隶属度，F_H、F_Slope、F_TWI分别为该栅格H、Slope和TWI值的概率分布；

步骤23，确定洼地、坡地、高地超渗区和高地下渗区四种类型的先验概率并计算每种类型下的栅格各指标的均值和方差(假定呈三维正态分布)：

P(w_i)＝N_i/N_t (11)

式中，P(w_i)为第i种类型下垫面的先验概率，w_i为第i种下垫面类型，N_i为第i种类型的栅格数量。

所述步骤3中利用最大似然分类法迭代得到各个栅格对每种下垫面类型的隶属度，具体步骤包括：

步骤31，将步骤2中的分类结果作为初始场，用最大似然分类法对非水体栅格进行再分类；

步骤32，将前一次的分类结果作为后一次分类的初始场进行迭代；

步骤33，根据迭代终点，计算所有非水体栅格对四种类型的隶属度并将水体栅格定义为100％隶属于水体分类。

最大似然分类的判别函数及判别规则为：

如果d_i(X)＝max_j＝1,2,3,4d_j(X),则X∈w_i (13)

式中，d_i(X)为栅格对第i种类型的隶属度，X为栅格的三维指标向量，D_i为第i种类型中栅格各指标的协方差矩阵，E_i为第i种类型的三维栅格指标均值向量。

最大似然分类法的迭代终点为：

|P_n(w_i)-P_n-1(w_i)|<0.001,i＝1,2,3,4 (14)

式中，n为迭代次数，P_n(w_i)为第n次迭代中第i中类型栅格所占流域的比例。

最终非水体栅格对四种分类的隶属度为：

式中，P(w_i|X)为栅格对第i种类型的隶属度，P(X|w_i)为第i种类型情况下该栅格出现的概率密度。

所述步骤4中在不同类型下垫面上建立不同的产流模式，将各类型对应的产流量根据隶属度加权叠加得到最终产流量，具体内容包括：水体、洼地、坡地、高地超渗区和高地下渗区五种下垫面类型对应的产流模式分别为：直接径流、饱和地面径流、快速壤中流、超渗地面径流和深层下渗(不产流)。各产流模式的计算公式为：

洼地：

坡地：

高地超渗区：

各非水体栅格最终产流量为：

R_s＝P(w₁|X)×R_sof+P(w₃|X)×R_eof (20)

R_i＝P(w₂|X)×R_ssf (21)

R_g＝P(w₁|X)×R_g1+P(w₂|X)×R_g2+P(w₃|X)×R_g3 (22)

式中，R_sof、R_ssf、R_eof、R_g1、R_g2、R_g3分别为饱和地面径流、快速壤中流、超渗地面径流以及洼地、坡地、高地超渗区对应的地下径流，a₁为R_sof与R_g1的分割比例，a₂为R_ssf与R_g2的分割比例，P_e为净雨量，S为蓄水容量，F为下渗能力，R_s、R_i、R_g分别为最终各栅格的地面径流、壤中流、地下径流产流量。

本发明的有益效果：本发明提供一种栅格尺度基于地形地貌的产流模式最大似然判别与产流计算方法，首先获取流域DEM数据并在此基础上计算栅格的H、Slope和TWI值；然后计算三个指标的累积分布函数，将非水体栅格分为洼地、坡地、高地超渗区和高地下渗区四种下垫面类型，得到四种类型的先验概率；之后利用最大似然分类法迭代得到各个栅格对每种类型的隶属度；最后在不同类型下垫面上建立不同的产流模式，将各类型对应的产流量根据隶属度加权叠加得到最终产流量。本发明通过数学方法对流域下垫面进行分类并基于栅格对各产流模式的隶属度加权计算得到最终产流量，减少了人为确定分类阈值带来的主观性，弥补了现阶段一种下垫面对应单一产流模式存在的缺陷，使得产流模型更具有普适性。

附图说明

图1是本发明提供的基于最大似然分类的栅格尺度多模式产流计算方法流程示意图；

图2是具体实施例中流域下垫面的H、Slope、TWI值分布图；

图3是具体实施例中流域下垫面初始分类图；

图4是具体实施例中流域下垫面栅格对每种类型的隶属度分布图。

具体实施方式

下面结合附图和具体实施例对本发明作进一步描述。

应当理解的是，此处所描述的具体实施方式仅仅用以解释本发明，并不用于限定本发明。

如图1所示，本发明实施例提供了一种栅格尺度基于地形地貌的产流模式最大似然判别与产流计算方法，包括以下步骤：

步骤1，获取流域DEM数据并在此基础上计算栅格的河面落差(H)、坡度(Slope)和地形湿度指数(TWI)值，包括：

H＝d-d₀ (1)

式中，d为栅格所在高程，d₀为沿着经过该栅格的流线方向所到达的第一个水体栅格的高程,

和

分别为栅格在水平和垂直方向的高程变化率，a为栅格单宽累积面积。

下载安徽省屯溪流域的DEM数据，利用Arcmap在DEM数据基础上计算H值、Slope值和TWI值分布图，如图2所示。

步骤2，计算三个指标的累积分布函数，将非水体栅格分为洼地、坡地、高地超渗区和高地下渗区四种下垫面类型。对于高地区，地形指数小的位置容易超渗，所以将其划分为高地超渗区，其余位置为高地下渗区。

根据划分的四种下垫面类型，得到四种类型的先验概率，具体步骤包括：

步骤21，提取除水体外所有栅格的H、Slope和TWI值，分别从大到小排列得到各自的累积分布函数，从而得到栅格各指标的概率分布，包括：

F_H(x)＝n₁/N_t (4)

F_Slope(x)＝n₂/N_t (5)

F_TWI(x)＝n₃/N_t (6)

式中，F_H(x)，F_Slope(x)，F_TWI(x)分别为该栅格H、Slope和TWI值的累积分布函数，n₁，n₂，n₃分别为H、Slope、TWI值小于等于x的非水体栅格数，N_t为非水体栅格总数。

步骤22，计算非水体栅格对洼地、坡地、高地超渗区和高地下渗区四种类型的隶属度,将栅格类型定义为隶属度最高者对应的类型，如图3所示：

M_w＝1-F_H (7)

M_h＝F_H×F_Slope (8)

M_pe＝[F_H×(1-F_Slope)]×(1-F_TWI) (9)

M_pp＝[F_H×(1-F_Slope)]×F_TWI (10)

式中，M_w、M_h、M_pe、M_pp分别为栅格对洼地、坡地、高地超渗区和高地下渗区的隶属度，F_H、F_Slope、F_TWI分别为该栅格H、Slope和TWI值的概率分布；

步骤23，确定洼地、坡地、高地超渗区和高地下渗区四种类型的先验概率并计算每种类型下的栅格各指标的均值和方差(假定呈三维正态分布)：

P(w_i)＝N_i/N_t (11)

式中，P(w_i)为第i种类型下垫面的先验概率，w_i为第i种下垫面类型，N_i为第i种类型的栅格数量。

步骤3，利用最大似然分类法迭代得到各个栅格对每种类型的隶属度，具体步骤包括：

步骤31，将步骤2中的分类结果作为初始场，用最大似然分类法对非水体栅格进行再分类；

步骤32，将前一次的分类结果作为后一次分类的初始场进行迭代；

步骤33，根据迭代终点，计算所有非水体栅格对四种下垫面类型的隶属度并将水体栅格定义为100％隶属于水体分类，如图4所示。

最大似然分类的判别函数及判别规则为：

i＝1,2,3,4

如果d_i(X)＝max_j＝1,2,3,4d_j(X),则X∈w_i (13)

式中，d_i(X)为栅格对第i种类型的隶属度，X为栅格的三维指标向量，D_i为第i种类型中栅格各指标的协方差矩阵，E_i为第i种类型的三维栅格指标均值向量。

最大似然分类法的迭代终点为：

|P_n(w_i)-P_n-1(w_i)|<0.001,i＝1,2,3,4 (14)

式中，n为迭代次数，P_n(w_i)为第n次迭代中第i中类型栅格所占流域的比例。

最终非水体栅格对四种分类的隶属度为：

式中，P(w_i|X)为栅格对第i种类型的隶属度，P(X|w_i)为第i种类型情况下该栅格出现的概率密度。

步骤4，在不同类型下垫面上建立不同的产流模式，将各类型对应的产流量根据隶属度加权叠加得到最终产流量。

所述步骤4中在不同类型下垫面上建立不同的产流模式，将各类型对应的产流量根据隶属度加权叠加得到最终产流量，具体内容包括：水体、洼地、坡地、高地超渗区和高地下渗区五种下垫面类型对应的产流模式分别为：直接径流、饱和地面径流、快速壤中流、超渗地面径流和深层下渗(不产流)。各产流模式的计算公式为：

洼地：

坡地：

高地超渗区：

各非水体栅格最终产流量为：

R_s＝P(w₁|X)×R_sof+P(w₃|X)×R_eof (20)

R_i＝P(w₂|X)×R_ssf (21)

R_g＝P(w₁|X)×R_g1+P(w₂|X)×R_g2+P(w₃|X)×R_g3 (22)

式中，R_sof、R_ssf、R_eof、R_g1、R_g2、R_g3分别为饱和地面径流、快速壤中流、超渗地面径流以及洼地、坡地、高地超渗区对应的地下径流，a₁为R_sof与R_g1的分割比例，a₂为R_ssf与R_g2的分割比例，P_e为净雨量，S为蓄水容量(与各点H值成正比)，F为下渗能力，R_s、R_i、R_g分别为最终各栅格的地面径流、壤中流、地下径流产流量。

以上所述仅是本发明的优选实施方式，并非用以限定本发明。本领域技术人员在不脱离本发明构思的前提下做出若干改进和优化，都应当视为本发明保护的范围。

Claims (7)

1.一种栅格尺度基于地形地貌的多模式产流计算方法，其特征在于，包括以下步骤：

步骤1，获取流域DEM数据并在此基础上计算栅格的河面落差、坡度和地形湿度指数值；

步骤2，计算河面落差、坡度和地形湿度指数值三个指标的累积分布函数，将非水体栅格分为洼地、坡地、高地超渗区和高地下渗区四种下垫面类型，得到四种类型的先验概率；

步骤3，利用最大似然分类法迭代得到各个栅格对每种下垫面类型的隶属度；

步骤4，将每种下垫面类型的产流作为一种产流模式并建立相应计算方法，根据每个栅格下垫面类型的隶属度判别产流的模式和组成，将各产流模式对应的产流量根据下垫面类型隶属度加权叠加得到最终产流量；

所述步骤2包括：

步骤21，提取除水体外所有栅格的河面落差、坡度和地形湿度指数值，分别从大到小排列得到各自的累积分布函数，从而得到栅格各指标的概率分布：

F_H(x)＝n₁/N_t (4)

F_Slope(x)＝n₂/N_t (5)

F_TWI(x)＝n₃/N_t (6)

式中，F_H(x)，F_Slope(x)，F_TWI(x)分别为该栅格河面落差、坡度和地形湿度指数值的累积分布函数，n₁，n₂，n₃分别为河面落差、坡度、地形湿度指数值小于等于x的非水体栅格数，N_t为非水体栅格总数；

步骤22，计算非水体栅格对洼地、坡地、高地超渗区和高地下渗区四种类型的隶属度,将栅格类型定义为隶属度最高者对应的类型：

M_w＝1-F_H (7)

M_h＝F_H×F_Slope (8)

M_pe＝[F_H×(1-F_Slope)]×(1-F_TWI) (9)

M_pp＝[F_H×(1-F_Slope)]×F_TWI (10)

式中，M_w、M_h、M_pe、M_pp分别为栅格对洼地、坡地、高地超渗区和高地下渗区的隶属度，F_H、F_Slope、F_TWI分别为栅格的河面落差、坡度和地形湿度指数值的概率分布；

步骤23，确定洼地、坡地、高地超渗区和高地下渗区四种类型的先验概率并计算每种类型下的栅格各指标的均值和方差：

P(w_i)＝N_i/N_t (11)

式中，P(w_i)为第i种类型下垫面的先验概率，w_i为第i种下垫面类型，N_i为第i种类型的栅格数量。

2.根据权利要求1所述的多模式产流计算方法，其特征在于，所述步骤3包括：

步骤31，将步骤2中的分类结果作为初始场，用最大似然分类法对非水体栅格进行再分类；

步骤32，将前一次的分类结果作为后一次分类的初始场进行迭代；

步骤33，根据迭代终点，计算最终非水体栅格对四种下垫面类型的隶属度并将水体栅格定义为100％隶属于水体分类。

3.根据权利要求1所述的多模式产流计算方法，其特征在于，最大似然分类的判别函数及判别规则为：

如果d_i(X)＝max_j＝1,2,3,4d_j(X),则X∈w_i (13)

式中，d_i(X)为栅格对第i种类型的隶属度，X为栅格的三维指标向量，D_i为第i种类型中栅格各指标的协方差矩阵，E_i为第i种类型的三维栅格指标均值向量。

4.根据权利要求3所述的多模式产流计算方法，其特征在于，最大似然分类法的迭代终点为：

|P_n(w_i)-P_n-1(w_i)|<0.001, (14)

i＝1,2,3,4

式中，n为迭代次数，P_n(w_i)为第n次迭代中第i种类型栅格所占流域的比例。

5.根据权利要求4所述的多模式产流计算方法，其特征在于，最终非水体栅格对四种下垫面类型的隶属度为：

式中，P(w_i|X)为栅格对第i种类型的隶属度，P(X|w_i)为第i种类型情况下该栅格出现的概率密度。

6.根据权利要求5所述的多模式产流计算方法，其特征在于，所述步骤4包括：水体、洼地、坡地、高地超渗区和高地下渗区五种下垫面类型对应的产流模式分别为：直接径流、饱和地面径流、快速壤中流、超渗地面径流和深层下渗。

7.根据权利要求6所述的多模式产流计算方法，其特征在于，各产流模式的计算公式为：

洼地：

坡地：

各非水体栅格最终产流量为：

R_s＝P(w₁|X)×R_sof+P(w₃|X)×R_eof (20)

R_i＝P(w₂|X)×R_ssf (21)

R_g＝P(w₁|X)×R_g1+P(w₂|X)×R_g2+P(w₃|X)×R_g3 (22)

式中，R_sof、R_ssf、R_eof、R_g1、R_g2、R_g3分别为饱和地面径流、快速壤中流、超渗地面径流以及洼地、坡地、高地超渗区对应的地下径流，a₁为R_sof与R_g1的分割比例，a₂为R_ssf与R_g2的分割比例，P_e为净雨量，S为蓄水容量，F为下渗能力，R_s、R_i、R_g分别为最终各栅格的地面径流、壤中流、地下径流产流量。

Images