CN114282403A - 一种耦合生境适宜模型的高效高精度栖息地模拟方法 - Google Patents

一种耦合生境适宜模型的高效高精度栖息地模拟方法 Download PDF

Info

Publication number
CN114282403A
CN114282403A CN202111387049.3A CN202111387049A CN114282403A CN 114282403 A CN114282403 A CN 114282403A CN 202111387049 A CN202111387049 A CN 202111387049A CN 114282403 A CN114282403 A CN 114282403A
Authority
CN
China
Prior art keywords
habitat
flux
grid
unit
water depth
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.)
Pending
Application number
CN202111387049.3A
Other languages
English (en)
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.)
Xian University of Technology
Original Assignee
Xian University of Technology
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 Xian University of Technology filed Critical Xian University of Technology
Priority to CN202111387049.3A priority Critical patent/CN114282403A/zh
Publication of CN114282403A publication Critical patent/CN114282403A/zh
Pending legal-status Critical Current

Links

Images

Landscapes

  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种耦合生境适宜模型的高效高精度栖息地模拟方法，采用矩形网格的高分辨率DEM地形，由已知条件赋值相关参数；基于Godunov格式中心有限体积法空间离散二维浅水方程，获取更高空间离散精度；由MUSCL重构和Runge‑Kutta法进一步实现时空的二阶精度；基于IFIM理论，采用适宜指数HSI构建生境适宜模型，以实现水动力过程与生境适宜模型的耦合；并引入GPU加速技术提高模拟效率。至此，耦合生境适宜模型的高效高精度栖息地模拟方法构建完成。本发明提高了生境模型中生境因子模拟精确性，为准确模拟河流水力特性和生物栖息地提供有效技术支撑，使之成为评价生物生境质量的有力工具。

Description

一种耦合生境适宜模型的高效高精度栖息地模拟方法
技术领域
本发明属于生物栖息地模拟方法技术领域，具体涉及一种耦合生境适宜模型的高效高精度栖息地模拟方法。
背景技术
一直以来，人们采取一系列水库生态运行和河流生态恢复的措施来缓解水利工程对河流生态系统的压力，为此提出许多模拟方法来解决此类问题，其中生境模型成为世界上应用最广泛的方法，如PHABSIM模型和River2D 模型。而生境因子的准确性对于生境模型更准确地模拟河流生态环境和评价水生生境的质量至关重要。近些年，为提高生境变量的时空分辨率，引入越来越多的水力模型与生境适宜性模型相结合，如Delft3D、Mike21、SSIIM 模型等二维或三维流体动力学模型。但对于高分辨率的地形数据，以及大型且具有复杂流态的河流而言，这些模型仍然存在模拟精度不足，模拟时间长等缺点。因此，亟需发展一种具有能够模拟高精度二维水动力过程的生境模型，来评价单种或多种水生生物的生境质量。基于此，构建一种耦合生境适宜模型的高效高精度栖息地模拟方法，以准确模拟河流水力特性和生物栖息地，使其成为评价生物生境质量的有力工具。
发明内容
本发明的目的是提供一种耦合生境适宜模型的高效高精度栖息地模拟方法，该方法可提高生境变量的时空分辨率；对于高分辨率地形数据及大型且具有复杂流态的河流而言，还可弥补当前很多模型存在的模拟精度不足，模拟时间长等缺点，能够更准确地模拟河流生态环境，评价单种或多种水生生物的生境质量。
本发明所采用的技术方案是，一种耦合生境适宜模型的高效高精度栖息地模拟方法，具体按照以下步骤实施：
步骤1、利用Arcgis软件生成高分辨率DEM地形数据表征真实地形，以ascii格式输出；
步骤2、调研收集研究河段目标物种以及目标物种的生境因子，利用偏好适宜曲线法确定出目标物种的生境因子适宜曲线，包括水深、流速和底质；
步骤3、将高分辨率DEM地形数据离散成矩形网格，根据已知边界条件对相关参数进行赋值；
步骤4、将得出的DEM地形数据及各生境因子适宜条件在GPU中分配内存，并复制到GPU变量中；
步骤5、基于Godunov格式的中心有限体积法空间离散二维浅水方程，求解通量项、底坡源项和摩阻源项；
步骤6、由MUSCL重组法和龙格库塔法进一步实现时间和空间的二阶精度，时间和空间步长计算稳定性由CFL条件限制；
步骤7、基于IFIM理论，由得到的目标物种高精度生境因子及目标物种的生境适宜条件，确定各生境因子的适宜性指数HSI计算出综合适宜性指数 CSI；
步骤8、通过cudaMemcpy函数将GPU计算得到的各生境因子数据及综合适宜指数CSI复制到CPU中，计算目标物种有效栖息地加权可利用面积WUA；
步骤9、根据计算得到的水深、流速和底质数据，以及水深、流速和底质三种生境因子的综合适宜指数CSI数据，利用MATLAB软件输出水深、流速、底质和综合适宜性指数CSI分布图；至此，耦合生境适宜模型的高效高精度栖息地模拟方法构建完成。
本发明的特点还在于，
步骤5中二维浅水方程表示如下：
Figure RE-GDA0003518950780000031
Figure RE-GDA0003518950780000032
Figure RE-GDA0003518950780000033
式中：t为时间步长，单位为s；x、y分别为水平坐标和纵向坐标； q为流量变量矢量；h表示水深，单位为m；u、v分别为x、y方向上的平均流速，单位为m/s；F、G分别为x、y方向上的通量矢量；S为源项矢量；其中，Sb为底坡源项；Sf为摩阻源项；zb为河床高程，m；Cf为床面糙率系数，Cf＝gn2/h1/3，n为曼宁系数，s/m1/3。g为重力加速度，取 9.81，单位为m/s2。
步骤5具体如下：
步骤5.1、第i个控制单元网格内，采用中心有限体积法对SWEs控制方程组积分，基于高斯散度定理通量项面积线积分形式为：
Figure RE-GDA0003518950780000041
通量F(q)·n线积分为：
Figure RE-GDA0003518950780000042
式中：Ω为第i个控制单元网格的体积，单位为m3；Г为第i个控制单元网格的边界；n为边界Г所对应的外法线方向的单位向量；Sb表示底坡源项通量；Sf表示摩阻源项通量；k为控制单元网格边的编号；lk为第i个控制单元网格第k个边的边长；
步骤5.2、在Godunov方法中实现HLLC黎曼求解器步骤为：
左波速SL、中间波速S*及右波速SR计算公式为，
Figure RE-GDA0003518950780000043
Figure RE-GDA0003518950780000044
Figure RE-GDA0003518950780000045
其中，中间变量h*和
Figure RE-GDA0003518950780000049
根据式计算，
Figure RE-GDA0003518950780000046
Figure RE-GDA0003518950780000047
从而计算通量F(q)·n的HLLC黎曼求解器计算公式为，
Figure RE-GDA0003518950780000048
其中，控制体单元网格界面中间通量F*为，
Figure RE-GDA0003518950780000051
控制体单元网格界面通量中间波两侧，F* L和F* R为，
Figure RE-GDA0003518950780000052
式中，FL和FR分别为单元网格界面左、右侧的通量；F*、F* L和F* R分别为控制体单元网格界面通量中间波两侧的左波通量、中间波通量和右波通量；u⊥为沿交界面的切向速度，u⊥＝-uny+vnx，u⊥L和u⊥R分别为控制体单元网格交界面的左侧、右侧切向速度，单位为m3/s；nx和ny分别为沿x 和y的指示方向；u和v分别为沿x和y的单宽流量，单位为m/s；变量 q⊥＝[h,qxnx+qyny]；qx和qy分别为沿x和y的单宽流量，单位为m2/s；SL、 SM、SR分别为左波速、中间波速、右波速，单位为m/s；u⊥为通过网格单元界面垂直的速度，单位为m/s；u⊥＝unx+vny；h*和u* ⊥为中间变量；g 为重力加速度，取9.81，单位为m/s2；
步骤5.3、底坡源项采用底坡通量法求解，将单元内底坡源项转化为单元边界上的通量。
采用底坡通量法得到某一边f处的底坡通量的矢量形式为，
Figure RE-GDA0003518950780000053
式中，Fsf(q)·为某一边f处的底坡通量；nf为面f单位外法向量；nfx和nfy是nf与(面f单位外法向量)的组成部分，hL和zbL分别为网格单元中心处的水深和底部高程，
Figure RE-GDA0003518950780000064
与zML网格边上的水深和底部高程，g为重力加速度，单位为m/s2；
步骤5.4、摩阻源项采用显隐式摩阻处理法求解，将隐式转化为显式，采用显隐式摩阻处理法求解摩阻源项计算公式为，
Figure RE-GDA0003518950780000061
其中，单宽流量更新为时间步长n+1时刻计算公式为，
Figure RE-GDA0003518950780000062
Figure RE-GDA0003518950780000063
式中，Sf为摩阻源项通量，Sfx和Sfy分别为x和y方向的摩阻源项通量；Cf为床面糙率系数，Cf＝gn2/h1/3，n为曼宁系数，s/m1/3；h为水深，单位为m；qx、qy分别为x、y方向上平均流速，单位为m2/s；qx n+1、qy n+1分别为时间步长n+1时刻在x、y方向上的单宽流量，单位为m2/s；qx n、qy n分别为时间步长n时刻在x、y方向上的单宽流量，单位为m2/s；Δt为时间步长，单位为s；Δqx为时间步长n时刻和n+1时刻的x方向上单宽流量之差，即Δqx＝qx n+1-qx n，单位为m2/s；为了避免分母为零，θ取1.0e-12；
步骤5.5、引入静水重构法求解干湿边界问题，即将界面处重构的负水深值强制为零来维持非负水深，并修改相应地水位、底坡高程和界面两侧的流量，并将重构的新值代入近似黎曼求解器中来计算质量和动量的通量，保证干湿界面处的平衡状态和质量守恒。
步骤6具体如下：
步骤6.1、求解出通量项、底坡源项和摩阻源项之后，采用MUSCL重构法在第i个控制单元网格通过同一限制变量梯度插值求解出第i个控制单元网格界面上的变量值；
步骤6.2、由两步Runge-Kutta方法获得二阶时间精度，在新时间步长中第i个控制单元网格的
Figure RE-GDA0003518950780000071
值为：
Figure RE-GDA0003518950780000072
中间变量值
Figure RE-GDA0003518950780000073
为
Figure RE-GDA0003518950780000074
Figure RE-GDA0003518950780000075
式中：i为控制单元网格数；Ω为第i个控制单元网格的体积，单位为 m3；qn+1为时间步长n+1时刻的单宽流量，单位为m2/s，qn为时间步长n 时刻的单宽流量，单位为m2/s；Δt为时间步长，单位为s；S(qn)为时间步长n时刻的源项；
步骤6.3、在空间步长Δx下，时间步长Δt的取值的计算公式如下：
Figure RE-GDA0003518950780000076
式中：Δxmin表示网格单元中心到对应界面的最小距离，单位为m；u为流速，单位为m/s；g为重力加速度，取值9.81m/s2；h为水深，单位为m。
步骤7具体如下：
步骤7.1、根据目标物种高精度生境因子及目标物种适宜条件，采用偏好适宜曲线法，确定出目标物种各生境因子适宜指数HSI，具体为，根据第 i个控制单元网格内水深veli、流速depi和底质subi三种生境因子的模拟数据，采用偏好适宜曲线法，计算出目标物种水深、流速和底质的适宜指数
Figure RE-GDA0003518950780000081
Figure RE-GDA0003518950780000082
和
Figure RE-GDA0003518950780000083
规定生境适宜指数范围为0～1，最适宜生境适宜指数为1，最不适宜生境适宜指数为0；
步骤7.2、第i个控制单元网格内综合适宜指数CSI由下式得出：
Figure RE-GDA0003518950780000084
Figure RE-GDA0003518950780000085
Figure RE-GDA0003518950780000086
式中：i为控制单元网格数；CSI为第i个控制单元网格内的综合适宜指数；HSI为水深、流速和底质等每个生境因子在第i个控制单元网格内的生境适宜指数；depi,veli,subi分别为第i个控制单元网格内水深、流速、底质或覆盖物，
Figure RE-GDA0003518950780000087
分别为第i个控制单元网格内水深、流速、底质或覆盖物的适宜指数；式(4)将所选择的生境因子适宜指数相乘，表示影响因子综合作用结果；式(5)考虑当某一生境因子较为不利时，组成栖息地生境因子之间的补偿影响；式(6)将最不适于鱼种生存的生境因子适宜指数作为组合适宜指数。
步骤8具体如下：
加权可利用面积WUA计算式；
Figure RE-GDA0003518950780000088
式中：WUA为栖息地加权可利用面积，单位为m2；i为控制单元网格数；ai为第i个控制单元网格的面积，单位为m2。
本发明的有益效果是，一种耦合生境适宜模型的高效高精度栖息地模拟方法，是针对现有数值模型在处理高分辨率的地形数据，以及大型且具有复杂流态的河流方面，存在模拟精度不足，模拟时间长等问题，提出一种具有高精度高效的生境适宜模型，以提高生境模型中生境因子精确性，为准确模拟河流水力特性和生物栖息地提供有效的技术支撑，使其成为评价生物生境质量的有力工具。
附图说明
图1是研究河段的高分辨率地形数据；
图2是花斑裸鲤产卵期生境因子适宜性曲线，其中，图2(a)是水深适宜性曲线图，图2(b)是流速适宜性曲线图，图2(c)是底质适宜性曲线图；
图3是基于Godunov格式的二维浅水方程中心有限体积法求解过程；
图4是基于GPU加速的生境适宜模型与高精度水动力过程耦合过程；
图5是GAST和River2D模型下不同流量与WUA关系模拟结果比较图；
图6是流量74m3/s条件下水深分布图；
图7是流量74m3/s条件下流速分布图；
图8是流量74m3/s条件下综合适宜指数CSI分布图。
具体实施方式
下面结合附图和具体实施方式对本发明进行详细说明。
本发明一种耦合生境适宜模型的高效高精度栖息地模拟方法，具体按照以下步骤实施：
步骤1、利用Arcgis软件生成高分辨率DEM地形数据表征真实地形，以ascii格式输出；
步骤2、调研收集研究河段目标物种以及目标物种的生境因子，利用偏好适宜曲线法确定出目标物种的生境因子适宜曲线，包括水深、流速和底质；
步骤3、将高分辨率DEM地形数据离散成矩形网格，根据已知边界条件对相关参数进行赋值；
步骤4、引入GPU加速技术以提高模拟效率，将得出的DEM地形数据及各生境因子适宜条件在GPU(Graphics Processing Unit,即图形处理器,又称为显示核心)中分配内存，并复制到GPU变量中；
步骤5、基于Godunov格式的中心有限体积法空间离散二维浅水方程，求解通量项、底坡源项和摩阻源项，获取更高空间离散精度，以提高激波捕捉的分辨率；
步骤5中二维浅水方程表示如下：
Figure RE-GDA0003518950780000101
Figure RE-GDA0003518950780000102
Figure RE-GDA0003518950780000103
式中：t为时间步长，单位为s；x、y分别为水平坐标和纵向坐标； q为流量变量矢量；h表示水深，单位为m；u、v分别为x、y方向上的平均流速，单位为m/s；F、G分别为x、y方向上的通量矢量；S为源项矢量；其中，Sb为底坡源项；Sf为摩阻源项；zb为河床高程，m；Cf为床面糙率系数，Cf＝gn2/h1/3，n为曼宁系数，s/m1/3。G为重力加速度，取 9.81，单位为m/s2。
步骤5具体如下：
步骤5.1、第i个控制单元网格内，采用中心有限体积法对SWEs控制方程组积分，基于高斯散度定理通量项面积线积分形式为：
Figure RE-GDA0003518950780000111
通量F(q)·n线积分为：
Figure RE-GDA0003518950780000112
式中：Ω为d第i个控制单元网格的体积，单位为m3；Г为第i个控制单元网格的边界；n为边界Г所对应的外法线方向的单位向量；Sb表示底坡源项通量；Sf表示摩阻源项通量；k为控制单元网格边的编号；lk为第i个控制单元网格第k个边的边长；
步骤5.2、通量项计算中，Godunov格式考虑迎风格式(upwind scheme)的信息源特性；对于求解中出现一个任意间断的阶梯函数(黎曼问题)，通过HLLC黎曼求解器，可直接得到控制体间数值通量的近似值，并可自动满足 Godunov格式。在Godunov方法中实现HLLC黎曼求解器步骤为：
考虑到存在网格单元无水的情况，左波速SL、中间波速S*及右波速SR计算公式为，
Figure RE-GDA0003518950780000113
Figure RE-GDA0003518950780000114
Figure RE-GDA0003518950780000115
其中，中间变量h*和
Figure RE-GDA0003518950780000116
根据式计算，
Figure RE-GDA0003518950780000121
Figure RE-GDA0003518950780000122
从而计算通量F(q)·n的HLLC黎曼求解器计算公式为，
Figure RE-GDA0003518950780000123
其中，控制体单元网格界面中间通量F*为，
Figure RE-GDA0003518950780000124
控制体单元网格界面通量中间波两侧，F* L和F* R为，
Figure RE-GDA0003518950780000125
式中，FL和FR分别为单元网格界面左、右侧的通量；F*、F* L和F* R分别为控制体单元网格界面通量中间波两侧的左波通量、中间波通量和右波通量；u||为沿交界面的切向速度，u⊥＝-uny+vnx，u⊥L和u⊥R分别为控制体单元网格交界面的左侧、右侧切向速度，单位为m3/s；nx和ny分别为沿x和y的指示方向；u和v分别为沿x和y的单宽流量，单位为m/s；变量q⊥＝[h,qxnx+qyny]；qx和qy分别为沿x和y的单宽流量，单位为m2/s； SL、SM、SR分别为左波速、中间波速、右波速，单位为m/s；u⊥为通过网格单元界面垂直的速度，单位为m/s；u⊥＝unx+vny；h*和
Figure RE-GDA0003518950780000126
为中间变量；g为重力加速度，取9.81，单位为m/s2；
步骤5.3、底坡源项采用底坡通量法求解，将单元内底坡源项转化为单元边界上的通量，稳定性更佳。
采用底坡通量法得到某一边f处的底坡通量的矢量形式为，
Figure RE-GDA0003518950780000131
式中，Fsf(q)·为某一边f处的底坡通量；nf为面f单位外法向量；nfx和nfy是nf与(面f单位外法向量)的组成部分，hL和zbL分别为网格单元中心处的水深和底部高程，
Figure RE-GDA0003518950780000132
与zML网格边上的水深和底部高程，g为重力加速度，单位为m/s2；
步骤5.4、摩阻源项采用显隐式摩阻处理法求解，将隐式转化为显式，消除常规隐式的冗余迭代，以兼顾计算精度和效率。采用显隐式摩阻处理法求解摩阻源项计算公式为，
Figure RE-GDA0003518950780000133
其中，单宽流量更新为时间步长n+1时刻计算公式为，
Figure RE-GDA0003518950780000134
Figure RE-GDA0003518950780000135
式中，Sf为摩阻源项通量，Sfx和Sfy分别为x和y方向的摩阻源项通量；Cf为床面糙率系数，Cf＝gn2/h1/3，n为曼宁系数，s/m1/3；h为水深，单位为m；qx、qy分别为x、y方向上平均流速，单位为m2/s；qx n+1、qy n+1分别为时间步长n+1时刻在x、y方向上的单宽流量，单位为m2/s；qx n、qy n分别为时间步长n时刻在x、y方向上的单宽流量，单位为m2/s；Δt为时间步长，单位为s；Δqx为时间步长n时刻和n+1时刻的x方向上单宽流量之差，即Δqx＝qx n+1-qx n，单位为m2/s；为了避免分母为零，θ取1.0e-12；
步骤5.5、引入静水重构法求解干湿边界问题，即将界面处重构的负水深值强制为零来维持非负水深，并修改相应地水位、底坡高程和界面两侧的流量，并将重构的新值代入近似黎曼求解器中来计算质量和动量的通量，保证干湿界面处的平衡状态和质量守恒。以保持干湿界面处的平衡状态和质量守恒，达到满足全稳条件。
步骤6、由MUSCL重组法和龙格库塔法(Runge-Kutta)进一步实现时间和空间的二阶精度，时间和空间步长计算稳定性由CFL(Courant Friedrichs Lewy)条件限制；
步骤6具体如下：
步骤6.1、求解出通量项、底坡源项和摩阻源项之后，采用MUSCL重构法在第i个控制单元网格通过同一限制变量梯度插值求解出第i个控制单元网格界面上的变量值，实现空间二阶精度，以平衡数值扩散和振荡问题；
步骤6.2、由两步Runge-Kutta方法获得二阶时间精度，在新时间步长中第i个控制单元网格的
Figure RE-GDA0003518950780000141
值为：
Figure RE-GDA0003518950780000142
中间变量值
Figure RE-GDA0003518950780000143
为
Figure RE-GDA0003518950780000144
Figure RE-GDA0003518950780000145
式中：i为控制单元网格数；Ω为第i个控制单元网格的体积，单位为m3；qn+1为时间步长n+1时刻的单宽流量，单位为m2/s，qn为时间步长n 时刻的单宽流量，单位为m2/s；Δt为时间步长，单位为s；S(qn)为时间步长n时刻的源项；
步骤6.3、在空间步长Δx下，时间步长Δt的取值的计算公式如下：
Figure RE-GDA0003518950780000151
式中：Δxmin表示网格单元中心到对应界面的最小距离，单位为m；u为流速，单位为m/s；g为重力加速度，取值9.81m/s2；h为水深，单位为m。
步骤7具体如下：
步骤7.1、根据目标物种高精度生境因子及目标物种适宜条件，采用偏好适宜曲线法，确定出目标物种各生境因子适宜指数HSI，具体为，根据第 i个控制单元网格内水深depi、流速veli和底质subi三种生境因子的模拟数据，采用偏好适宜曲线法，计算出目标物种水深、流速和底质的适宜指数
Figure RE-GDA0003518950780000152
Figure RE-GDA0003518950780000153
和
Figure RE-GDA0003518950780000154
规定生境适宜指数范围为0～1，最适宜生境适宜指数为1，最不适宜生境适宜指数为0；
步骤7.2、结合实际情况，选择合适的计算方法，计算出综合适宜指数 CSI(Combined Suitability index)；第i个控制单元网格内综合适宜指数CSI 由下式得出：
Figure RE-GDA0003518950780000155
Figure RE-GDA0003518950780000156
Figure RE-GDA0003518950780000157
式中：i为控制单元网格数；CSI为第i个控制单元网格内的综合适宜指数；HSI为水深、流速和底质等每个生境因子在第i个控制单元网格内的生境适宜指数；depi,veli,subi分别为第i个控制单元网格内水深、流速、底质或覆盖物，
Figure RE-GDA0003518950780000161
分别为第i个控制单元网格内水深、流速、底质或覆盖物的适宜指数；式(4)将所选择的生境因子适宜指数相乘，表示影响因子综合作用结果；式(5)考虑当某一生境因子较为不利时，组成栖息地生境因子之间的补偿影响；式(6)将最不适于鱼种生存的生境因子适宜指数作为组合适宜指数。
步骤7、基于IFIM理论，由得到的目标物种高精度生境因子及目标物种的生境适宜条件，确定各生境因子的适宜性指数HSI(Habitat Suitability Index)，结合实际情况，选择合适的计算方法，计算出综合适宜性指数CSI (Combined Suitability index)；
步骤8、通过cudaMemcpy函数将GPU计算得到的各生境因子数据及综合适宜指数CSI复制到CPU中，计算目标物种有效栖息地加权可利用面积 WUA(Weighted UsableArea)；
步骤8具体如下：
根据综合适宜性指数CSI计算出目标物种适宜有效栖息地加权可利用面积WUA(Weighted Usable Area)。与River2D模拟结果进行比较，以验证该方法的准确性；
加权可利用面积WUA计算式：
Figure RE-GDA0003518950780000162
式中：WUA为栖息地加权可利用面积，单位为m2；i为控制单元网格数； ai为第i个控制单元网格的面积，单位为m2。
步骤9、根据计算得到的水深、流速和底质数据，以及水深、流速和底质三种生境因子的综合适宜指数CSI数据，利用MATLAB软件输出水深、流速、底质和综合适宜性指数CSI分布图；至此，耦合生境适宜模型的高效高精度栖息地模拟方法构建完成。
以黄河上游某河段1～7km为研究河段，如图1所示，河段高分辨率地形空间离散成404,544个矩形网格，单元网格边长为5m，河道曼宁糙率系数取值为0.05s/m1/3。选择花斑裸鲤为目标鱼种，水深、流速和底质为花斑裸鲤的生境因子。由水生专家调查研究知，花斑裸鲤产卵期适宜水深为0.1～ 2.0m，最适宜水深为0.4～1.2m；适宜流速为0.25～1.5m/s，最适宜流速为 0.5～1.0m/s；适宜基质条件为砂砾底质。根据花斑裸鲤产卵期水深、流速和底质三种生境因子条件，采用偏好适宜法绘制生境适宜性曲线，生境适宜指数HSI范围为0～1。其中最适宜生境适宜指数为1，最不适宜生境适宜指数为0。从而得水深、流速和底质三种生境因子的生境适宜曲线。
图2(a)是水深适宜性曲线图。该图根据上述调查产卵期适宜水深条件；利用适宜性指数的规定原则，最适宜的产卵水深适宜性指数为1.0，最不适宜产卵水深适宜性指数为0.0，之间为0～1.0.从得到的水深适宜性曲线图。图 2(b)是流速适宜性曲线图。该图根据上述调查产卵期适宜流速条件，利用适宜性指数的规定原则，最适宜的产卵流速适宜性指数为1.0，最不适宜产卵流速适宜性指数为0.0，之间为0～1.0，从而得到流速适宜性曲线图。图2 (c)是底质适宜性曲线图。该图根据河段基质条件多为砂砾底质，又知花斑裸鲤产卵适宜基质条件为砂砾底质，利用适宜性指数的规定原则，底质的适宜性指数为1.0。
通过将构建的生境适宜模型与高精度水动力过程相耦合，得到一种耦合生境适宜模型的高效高精度栖息地模拟方法，如图3和图4所示。为验证耦合生境适宜模型的高效高精度栖息地模拟方法的计算结果准确性，选取流量 26.5m3/s、53m3/s、74m3/s、79.5m3/s、106m3/s、159m3/s、212m3/s、318m3/s、 434m3/s、530m3/s、689m3/s、848m3/s、1060m3/s，分别采用GAST模型和 River2D模型计算花斑裸鲤有效栖息地面积(以加权可用面积表示)。采用两种模型模拟不同流量条件下花斑裸鲤有效栖息地面积，结果如图5所示。由图5可知，在两种模型模拟下，花斑裸鲤产卵的有效栖息地面积均表现为：随着流量增加，有效栖息地面积不断增大，达到某一流量时，有效栖息地面积出现峰值。之后，随着流量增加，有效栖息地面积不断减少。根据对比图两种模型模拟结果几乎相同，规律一致，计算得到有效栖息地面积误差最大仅为5.366％。表明GAST模型结果具有较好的准确性。
此外，选择流量74m3/s，得出该流量下水深和流速分布图，如图6和图7所示。结合生境适宜模型，考虑花斑裸鲤产卵期水深、流速和底质三种生境因子，计算得花斑裸鲤产卵期的综合适宜指数CSI，CSI的分布如图8所示。发现在流量74m3/s时，花斑裸鲤产卵期有效栖息地面积达到最大，适宜产卵繁殖的有效栖息地分布区域向河道中间聚拢。且在流量74m3/s时，采用GAST模型和River2D模型模拟得到的加权可利用面积WUA分别为202492.08m2和198775.00m2。发现GAST模型得到的有效产卵栖息地面积略大，这是由于GAST模型利用高分辨率地形条件模拟得到的水深和速度更为准确。

Claims (6)

1.一种耦合生境适宜模型的高效高精度栖息地模拟方法，其特征在于，具体按照以下步骤实施：
步骤1、利用Arcgis软件生成高分辨率DEM地形数据表征真实地形，以ascii格式输出；
步骤2、调研收集研究河段目标物种以及目标物种的生境因子，利用偏好适宜曲线法确定出目标物种的生境因子适宜曲线，包括水深、流速和底质；
步骤3、将高分辨率DEM地形数据离散成矩形网格，根据已知边界条件对相关参数进行赋值；
步骤4、将得出的DEM地形数据及各生境因子适宜条件在GPU中分配内存，并复制到GPU变量中；
步骤5、基于Godunov格式的中心有限体积法空间离散二维浅水方程，求解通量项、底坡源项和摩阻源项；
步骤6、由MUSCL重组法和龙格库塔法进一步实现时间和空间的二阶精度，时间和空间步长计算稳定性由CFL条件限制；
步骤7、基于IFIM理论，由得到的目标物种高精度生境因子及目标物种的生境适宜条件，确定各生境因子的适宜性指数HSI计算出综合适宜性指数CSI；
步骤8、通过cudaMemcpy函数将GPU计算得到的各生境因子数据及综合适宜指数CSI复制到CPU中，计算目标物种有效栖息地加权可利用面积WUA；
步骤9、根据计算得到的水深、流速和底质数据，以及水深、流速和底质三种生境因子的综合适宜指数CSI数据，利用MATLAB软件输出水深、流速、底质和综合适宜性指数CSI分布图；至此，耦合生境适宜模型的高效高精度栖息地模拟方法构建完成。
2.根据权利要求1所述的一种耦合生境适宜模型的高效高精度栖息地模拟方法，其特征在于，所述步骤5中二维浅水方程表示如下：
Figure FDA0003367427940000021
Figure FDA0003367427940000022
Figure FDA0003367427940000023
式中：t为时间步长，单位为s；x、y分别为水平坐标和纵向坐标；q为流量变量矢量；h表示水深，单位为m；u、v分别为x、y方向上的平均流速，单位为m/s；F、G分别为x、y方向上的通量矢量；S为源项矢量；其中，Sb为底坡源项；Sf为摩阻源项；zb为河床高程，m；Cf为床面糙率系数，Cf＝gn2/h1/3，n为曼宁系数，s/m1/3。g为重力加速度，取9.81，单位为m/s2。
3.根据权利要求2所述的一种耦合生境适宜模型的高效高精度栖息地模拟方法，其特征在于，所述步骤5具体如下：
步骤5.1、第i个控制单元网格内，采用中心有限体积法对SWEs控制方程组积分，基于高斯散度定理通量项面积线积分形式为：
Figure RE-FDA0003518950770000024
通量F(q)·n线积分为：
Figure RE-FDA0003518950770000031
式中：Ω为d第i个控制单元网格的体积，单位为m3；Г为第i个控制单元网格的边界；n为边界Г所对应的外法线方向的单位向量；Sb表示底坡源项通量；Sf表示摩阻源项通量；k为控制单元网格边的编号；lk为第i个控制单元网格第k个边的边长；
步骤5.2、在Godunov方法中实现HLLC黎曼求解器步骤为：
左波速SL、中间波速S*及右波速SR计算公式为，
Figure RE-FDA0003518950770000032
Figure RE-FDA0003518950770000033
Figure RE-FDA0003518950770000034
其中，中间变量h*和
Figure RE-FDA0003518950770000039
根据式计算，
Figure RE-FDA0003518950770000035
Figure RE-FDA0003518950770000036
从而计算通量F(q)·n的HLLC黎曼求解器计算公式为，
Figure RE-FDA0003518950770000037
其中，控制体单元网格界面中间通量F*为，
Figure RE-FDA0003518950770000038
控制体单元网格界面通量中间波两侧，F* L和F* R为，
Figure RE-FDA0003518950770000041
式中，FL和FR分别为单元网格界面左、右侧的通量；F*、F* L和F* R分别为控制体单元网格界面通量中间波两侧的左波通量、中间波通量和右波通量；u⊥为沿交界面的切向速度，u⊥＝-uny+vnx，u⊥L和u⊥R分别为控制体单元网格交界面的左侧、右侧切向速度，单位为m3/s；nx和ny分别为沿x和y的指示方向；u和v分别为沿x和y的单宽流量，单位为m/s；变量q⊥＝[h,qxnx+qyny]；qx和qy分别为沿x和y的单宽流量，单位为m2/s；SL、SM、SR分别为左波速、中间波速、右波速，单位为m/s；u⊥为通过网格单元界面垂直的速度，单位为m/s；u⊥＝unx+vny；h*和
Figure RE-FDA0003518950770000042
为中间变量；g为重力加速度，取9.81，单位为m/s2；
步骤5.3、底坡源项采用底坡通量法求解，将单元内底坡源项转化为单元边界上的通量。
采用底坡通量法得到某一边f处的底坡通量的矢量形式为，
Figure RE-FDA0003518950770000043
式中，Fsf(q)·为某一边f处的底坡通量；nf为面f单位外法向量；nfx和nfy是nf与(面f单位外法向量)的组成部分，hL和zbL分别为网格单元中心处的水深和底部高程，
Figure RE-FDA0003518950770000044
与zML网格边上的水深和底部高程，g为重力加速度，单位为m/s2；
步骤5.4、摩阻源项采用显隐式摩阻处理法求解，将隐式转化为显式，采用显隐式摩阻处理法求解摩阻源项计算公式为，
Figure RE-FDA0003518950770000051
其中，单宽流量更新为时间步长n+1时刻计算公式为，
Figure RE-FDA0003518950770000052
Figure RE-FDA0003518950770000053
式中，Sf为摩阻源项通量，Sfx和Sfy分别为x和y方向的摩阻源项通量；Cf为床面糙率系数，Cf＝gn2/h1/3，n为曼宁系数，s/m1/3；h为水深，单位为m；qx、qy分别为x、y方向上平均流速，单位为m2/s；qx n+1、qy n+1分别为时间步长n+1时刻在x、y方向上的单宽流量，单位为m2/s；qx n、qy n分别为时间步长n时刻在x、y方向上的单宽流量，单位为m2/s；Δt为时间步长，单位为s；Δqx为时间步长n时刻和n+1时刻的x方向上单宽流量之差，即Δqx＝qx n+1-qx n，单位为m2/s；为了避免分母为零，θ取1.0e-12；
步骤5.5、引入静水重构法求解干湿边界问题，即将界面处重构的负水深值强制为零来维持非负水深，并修改相应地水位、底坡高程和界面两侧的流量，并将重构的新值代入近似黎曼求解器中来计算质量和动量的通量，保证干湿界面处的平衡状态和质量守恒。
4.根据权利要求3所述的一种耦合生境适宜模型的高效高精度栖息地模拟方法，其特征在于，所述步骤6具体如下：
步骤6.1、求解出通量项、底坡源项和摩阻源项之后，采用MUSCL重构法在第i个控制单元网格通过同一限制变量梯度插值求解出第i个控制单元网格界面上的变量值；
步骤6.2、由两步Runge-Kutta方法获得二阶时间精度，在新时间步长中第i个控制单元网格的
Figure FDA0003367427940000061
值为：
Figure FDA0003367427940000062
中间变量值
Figure FDA0003367427940000063
为
Figure FDA0003367427940000064
Figure FDA0003367427940000065
式中：i为控制单元网格数；Ω为第i个控制单元网格的体积，单位为m3；qn+1为时间步长n+1时刻的单宽流量，单位为m2/s，qn为时间步长n时刻的单宽流量，单位为m2/s；Δt为时间步长，单位为s；S(qn)为时间步长n时刻的源项；
步骤6.3、在空间步长Δx下，时间步长Δt的取值的计算公式如下：
Figure FDA0003367427940000066
式中：Δxmin表示网格单元中心到对应界面的最小距离，单位为m；u为流速，单位为m/s；g为重力加速度，取值9.81m/s2；h为水深，单位为m。
5.根据权利要求4所述的一种耦合生境适宜模型的高效高精度栖息地模拟方法，其特征在于，所述步骤7具体如下：
步骤7.1、根据目标物种高精度生境因子及目标物种适宜条件，采用偏好适宜曲线法，确定出目标物种各生境因子适宜指数HSI，具体为，根据第i个控制单元网格内水深depi、流速veli和底质subi三种生境因子的模拟数据，采用偏好适宜曲线法，计算出目标物种水深、流速和底质的适宜指数
Figure FDA0003367427940000067
Figure FDA0003367427940000071
和
Figure FDA0003367427940000072
规定生境适宜指数范围为0～1，最适宜生境适宜指数为1，最不适宜生境适宜指数为0；
步骤7.2、第i个控制单元网格内综合适宜指数CSI由下式得出：
Figure FDA0003367427940000073
Figure FDA0003367427940000074
Figure FDA0003367427940000075
式中：i为控制单元网格数；CSI为第i个控制单元网格内的综合适宜指数；HSI为水深、流速和底质等每个生境因子在第i个控制单元网格内的生境适宜指数；depi,veli,subi分别为第i个控制单元网格内水深、流速、底质或覆盖物，
Figure FDA0003367427940000076
分别为第i个控制单元网格内水深、流速、底质或覆盖物的适宜指数；式(4)将所选择的生境因子适宜指数相乘，表示影响因子综合作用结果；式(5)考虑当某一生境因子较为不利时，组成栖息地生境因子之间的补偿影响；式(6)将最不适于鱼种生存的生境因子适宜指数作为组合适宜指数。
6.根据权利要求5所述的一种耦合生境适宜模型的高效高精度栖息地模拟方法，其特征在于，所述步骤8具体如下：
加权可利用面积WUA计算式；
Figure FDA0003367427940000077
式中：WUA为栖息地加权可利用面积，单位为m2；i为控制单元网格数；ai为第i个控制单元网格的面积，单位为m2。
CN202111387049.3A 2021-11-22 2021-11-22 一种耦合生境适宜模型的高效高精度栖息地模拟方法 Pending CN114282403A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111387049.3A CN114282403A (zh) 2021-11-22 2021-11-22 一种耦合生境适宜模型的高效高精度栖息地模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111387049.3A CN114282403A (zh) 2021-11-22 2021-11-22 一种耦合生境适宜模型的高效高精度栖息地模拟方法

Publications (1)

Publication Number Publication Date
CN114282403A true CN114282403A (zh) 2022-04-05

Family

ID=80869567

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111387049.3A Pending CN114282403A (zh) 2021-11-22 2021-11-22 一种耦合生境适宜模型的高效高精度栖息地模拟方法

Country Status (1)

Country Link
CN (1) CN114282403A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116029469A (zh) * 2023-03-30 2023-04-28 长江水资源保护科学研究所 基于水鸟生境适宜性的闸控湖泊湿地水位确定方法及装置

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116029469A (zh) * 2023-03-30 2023-04-28 长江水资源保护科学研究所 基于水鸟生境适宜性的闸控湖泊湿地水位确定方法及装置
CN116029469B (zh) * 2023-03-30 2023-06-27 长江水资源保护科学研究所 基于水鸟生境适宜性的闸控湖泊湿地水位确定方法及装置

Similar Documents

Publication Publication Date Title
CN111768502A (zh) 一种基于gpu加速技术的非结构网格二维洪水模拟系统
Martins et al. 3D modelling in the Sado estuary using a new generic vertical discretization approach
Losasso et al. Two-way coupled SPH and particle level set fluid simulation
van Maanen et al. Modelling the effects of tidal range and initial bathymetry on the morphological evolution of tidal embayments
CN108021780B (zh) 一种基于无规则非结构网格模型的山洪动态仿真方法
Hunke et al. Sea-ice models for climate study: retrospective and new directions
Legrand et al. High-resolution, unstructured meshes for hydrodynamic models of the Great Barrier Reef, Australia
Nakahashi et al. Building-cube method for large-scale, high resolution flow computations
CN103699714B (zh) 一种基于有限元和无网格耦合的柔性物体实时切割仿真方法
CN108563837B (zh) 一种冲积河流水沙模型的模型参数实时校正方法和系统
CN107451372B (zh) 一种运动波与动力波相结合的山区洪水过程数值模拟方法
CN112270115B (zh) 一种基于元胞自动机的复杂地形洪水淹没进程模拟方法
CN111881596B (zh) 一种基于拉格朗日插值的溢油污染源逆时追踪模拟方法
Siadatmousavi et al. On the importance of high frequency tail in third generation wave models
Ibrayev Model of enclosed and semi-enclosed sea hydrodynamics
CN115994496B (zh) 城市公园高分辨率大气co2浓度三维场的数值模拟方法
CN114282403A (zh) 一种耦合生境适宜模型的高效高精度栖息地模拟方法
CN106569270B (zh) 规则网格速度模型的自适应非结构三角网格化方法
Yang et al. Application of habitat suitability model coupling with high-precision hydrodynamic processes
Cordonnier et al. Forming terrains by glacial erosion
CN114580283A (zh) 山溪性强潮河口分汊段中长期动力地貌演变数值模拟方法
CN114925624A (zh) 一种天然河道三维水流数值模拟方法
CN114997006A (zh) 基于多gpu并行框架的河流泥沙输移过程数值模拟方法
Tanaka et al. A consistent finite difference local inertial model for shallow water simulation
CN112364559A (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