CN114282403A - 一种耦合生境适宜模型的高效高精度栖息地模拟方法 - Google Patents
一种耦合生境适宜模型的高效高精度栖息地模拟方法 Download PDFInfo
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 74
- 238000004088 simulation Methods 0.000 title claims abstract description 36
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims abstract description 81
- 230000001133 acceleration Effects 0.000 claims abstract description 15
- 238000010168 coupling process Methods 0.000 claims abstract description 12
- 230000008878 coupling Effects 0.000 claims abstract description 11
- 238000005859 coupling reaction Methods 0.000 claims abstract description 11
- 210000003205 muscle Anatomy 0.000 claims abstract description 7
- 238000010276 construction Methods 0.000 claims abstract description 4
- 230000004907 flux Effects 0.000 claims description 62
- 241000894007 species Species 0.000 claims description 31
- 238000004364 calculation method Methods 0.000 claims description 26
- 239000000758 substrate Substances 0.000 claims description 25
- 230000005484 gravity Effects 0.000 claims description 12
- 238000007716 flux method Methods 0.000 claims description 6
- 239000000463 material Substances 0.000 claims description 6
- 241000251468 Actinopterygii Species 0.000 claims description 4
- 230000006870 function Effects 0.000 claims description 4
- 230000009471 action Effects 0.000 claims description 3
- 230000002146 bilateral effect Effects 0.000 claims description 3
- 230000001276 controlling effect Effects 0.000 claims description 3
- 238000005215 recombination Methods 0.000 claims description 3
- 230000006798 recombination Effects 0.000 claims description 3
- 230000001105 regulatory effect Effects 0.000 claims description 3
- 230000004083 survival effect Effects 0.000 claims description 3
- 230000008569 process Effects 0.000 abstract description 5
- 238000005516 engineering process Methods 0.000 abstract description 2
- 241001275921 Gymnocypris przewalskii Species 0.000 description 8
- 230000017448 oviposition Effects 0.000 description 7
- 238000010586 diagram Methods 0.000 description 3
- 241001275923 Gymnocypris Species 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 230000003044 adaptive effect Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000009395 breeding Methods 0.000 description 1
- 230000001488 breeding effect Effects 0.000 description 1
- 238000009792 diffusion process Methods 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 230000010355 oscillation Effects 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 230000035939 shock Effects 0.000 description 1
- 238000011144 upstream manufacturing Methods 0.000 description 1
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中二维浅水方程表示如下:
式中: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控制方程组积分,基于高斯散度定理通量项面积线积分形式为:
通量F(q)·n线积分为:
式中:Ω为第i个控制单元网格的体积,单位为m3;Г为第i个控制单元网格的边界;n为边界Г所对应的外法线方向的单位向量;Sb表示底坡源项通量;Sf表示摩阻源项通量;k为控制单元网格边的编号;lk为第i个控制单元网格第k个边的边长;
步骤5.2、在Godunov方法中实现HLLC黎曼求解器步骤为:
左波速SL、中间波速S*及右波速SR计算公式为,
从而计算通量F(q)·n的HLLC黎曼求解器计算公式为,
其中,控制体单元网格界面中间通量F*为,
控制体单元网格界面通量中间波两侧,F* L和F* R为,
式中,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处的底坡通量的矢量形式为,
式中,Fsf(q)·为某一边f处的底坡通量;nf为面f单位外法向量;nfx和nfy是nf与(面f单位外法向量)的组成部分,hL和zbL分别为网格单元中心处的水深和底部高程,与zML网格边上的水深和底部高程,g为重力加速度,单位为m/s2;
步骤5.4、摩阻源项采用显隐式摩阻处理法求解,将隐式转化为显式,采用显隐式摩阻处理法求解摩阻源项计算公式为,
其中,单宽流量更新为时间步长n+1时刻计算公式为,
式中,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个控制单元网格界面上的变量值;
式中:i为控制单元网格数;Ω为第i个控制单元网格的体积,单位为 m3;qn+1为时间步长n+1时刻的单宽流量,单位为m2/s,qn为时间步长n 时刻的单宽流量,单位为m2/s;Δt为时间步长,单位为s;S(qn)为时间步长n时刻的源项;
步骤6.3、在空间步长Δx下,时间步长Δt的取值的计算公式如下:
式中:Δxmin表示网格单元中心到对应界面的最小距离,单位为m;u为流速,单位为m/s;g为重力加速度,取值9.81m/s2;h为水深,单位为m。
步骤7具体如下:
步骤7.1、根据目标物种高精度生境因子及目标物种适宜条件,采用偏好适宜曲线法,确定出目标物种各生境因子适宜指数HSI,具体为,根据第 i个控制单元网格内水深veli、流速depi和底质subi三种生境因子的模拟数据,采用偏好适宜曲线法,计算出目标物种水深、流速和底质的适宜指数 和规定生境适宜指数范围为0~1,最适宜生境适宜指数为1,最不适宜生境适宜指数为0;
步骤7.2、第i个控制单元网格内综合适宜指数CSI由下式得出:
式中:i为控制单元网格数;CSI为第i个控制单元网格内的综合适宜指数;HSI为水深、流速和底质等每个生境因子在第i个控制单元网格内的生境适宜指数;depi,veli,subi分别为第i个控制单元网格内水深、流速、底质或覆盖物,分别为第i个控制单元网格内水深、流速、底质或覆盖物的适宜指数;式(4)将所选择的生境因子适宜指数相乘,表示影响因子综合作用结果;式(5)考虑当某一生境因子较为不利时,组成栖息地生境因子之间的补偿影响;式(6)将最不适于鱼种生存的生境因子适宜指数作为组合适宜指数。
步骤8具体如下:
加权可利用面积WUA计算式;
式中: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中二维浅水方程表示如下:
式中: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控制方程组积分,基于高斯散度定理通量项面积线积分形式为:
通量F(q)·n线积分为:
式中:Ω为d第i个控制单元网格的体积,单位为m3;Г为第i个控制单元网格的边界;n为边界Г所对应的外法线方向的单位向量;Sb表示底坡源项通量;Sf表示摩阻源项通量;k为控制单元网格边的编号;lk为第i个控制单元网格第k个边的边长;
步骤5.2、通量项计算中,Godunov格式考虑迎风格式(upwind scheme)的信息源特性;对于求解中出现一个任意间断的阶梯函数(黎曼问题),通过HLLC黎曼求解器,可直接得到控制体间数值通量的近似值,并可自动满足 Godunov格式。在Godunov方法中实现HLLC黎曼求解器步骤为:
考虑到存在网格单元无水的情况,左波速SL、中间波速S*及右波速SR计算公式为,
从而计算通量F(q)·n的HLLC黎曼求解器计算公式为,
其中,控制体单元网格界面中间通量F*为,
控制体单元网格界面通量中间波两侧,F* L和F* R为,
式中,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*和为中间变量;g为重力加速度,取9.81,单位为m/s2;
步骤5.3、底坡源项采用底坡通量法求解,将单元内底坡源项转化为单元边界上的通量,稳定性更佳。
采用底坡通量法得到某一边f处的底坡通量的矢量形式为,
式中,Fsf(q)·为某一边f处的底坡通量;nf为面f单位外法向量;nfx和nfy是nf与(面f单位外法向量)的组成部分,hL和zbL分别为网格单元中心处的水深和底部高程,与zML网格边上的水深和底部高程,g为重力加速度,单位为m/s2;
步骤5.4、摩阻源项采用显隐式摩阻处理法求解,将隐式转化为显式,消除常规隐式的冗余迭代,以兼顾计算精度和效率。采用显隐式摩阻处理法求解摩阻源项计算公式为,
其中,单宽流量更新为时间步长n+1时刻计算公式为,
式中,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个控制单元网格界面上的变量值,实现空间二阶精度,以平衡数值扩散和振荡问题;
式中:i为控制单元网格数;Ω为第i个控制单元网格的体积,单位为m3;qn+1为时间步长n+1时刻的单宽流量,单位为m2/s,qn为时间步长n 时刻的单宽流量,单位为m2/s;Δt为时间步长,单位为s;S(qn)为时间步长n时刻的源项;
步骤6.3、在空间步长Δx下,时间步长Δt的取值的计算公式如下:
式中:Δxmin表示网格单元中心到对应界面的最小距离,单位为m;u为流速,单位为m/s;g为重力加速度,取值9.81m/s2;h为水深,单位为m。
步骤7具体如下:
步骤7.1、根据目标物种高精度生境因子及目标物种适宜条件,采用偏好适宜曲线法,确定出目标物种各生境因子适宜指数HSI,具体为,根据第 i个控制单元网格内水深depi、流速veli和底质subi三种生境因子的模拟数据,采用偏好适宜曲线法,计算出目标物种水深、流速和底质的适宜指数 和规定生境适宜指数范围为0~1,最适宜生境适宜指数为1,最不适宜生境适宜指数为0;
步骤7.2、结合实际情况,选择合适的计算方法,计算出综合适宜指数 CSI(Combined Suitability index);第i个控制单元网格内综合适宜指数CSI 由下式得出:
式中:i为控制单元网格数;CSI为第i个控制单元网格内的综合适宜指数;HSI为水深、流速和底质等每个生境因子在第i个控制单元网格内的生境适宜指数;depi,veli,subi分别为第i个控制单元网格内水深、流速、底质或覆盖物,分别为第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计算式:
式中: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分布图;至此,耦合生境适宜模型的高效高精度栖息地模拟方法构建完成。
3.根据权利要求2所述的一种耦合生境适宜模型的高效高精度栖息地模拟方法,其特征在于,所述步骤5具体如下:
步骤5.1、第i个控制单元网格内,采用中心有限体积法对SWEs控制方程组积分,基于高斯散度定理通量项面积线积分形式为:
通量F(q)·n线积分为:
式中:Ω为d第i个控制单元网格的体积,单位为m3;Г为第i个控制单元网格的边界;n为边界Г所对应的外法线方向的单位向量;Sb表示底坡源项通量;Sf表示摩阻源项通量;k为控制单元网格边的编号;lk为第i个控制单元网格第k个边的边长;
步骤5.2、在Godunov方法中实现HLLC黎曼求解器步骤为:
左波速SL、中间波速S*及右波速SR计算公式为,
从而计算通量F(q)·n的HLLC黎曼求解器计算公式为,
其中,控制体单元网格界面中间通量F*为,
控制体单元网格界面通量中间波两侧,F* L和F* R为,
式中,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*和为中间变量;g为重力加速度,取9.81,单位为m/s2;
步骤5.3、底坡源项采用底坡通量法求解,将单元内底坡源项转化为单元边界上的通量。
采用底坡通量法得到某一边f处的底坡通量的矢量形式为,
式中,Fsf(q)·为某一边f处的底坡通量;nf为面f单位外法向量;nfx和nfy是nf与(面f单位外法向量)的组成部分,hL和zbL分别为网格单元中心处的水深和底部高程,与zML网格边上的水深和底部高程,g为重力加速度,单位为m/s2;
步骤5.4、摩阻源项采用显隐式摩阻处理法求解,将隐式转化为显式,采用显隐式摩阻处理法求解摩阻源项计算公式为,
其中,单宽流量更新为时间步长n+1时刻计算公式为,
式中,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个控制单元网格界面上的变量值;
式中:i为控制单元网格数;Ω为第i个控制单元网格的体积,单位为m3;qn+1为时间步长n+1时刻的单宽流量,单位为m2/s,qn为时间步长n时刻的单宽流量,单位为m2/s;Δt为时间步长,单位为s;S(qn)为时间步长n时刻的源项;
步骤6.3、在空间步长Δx下,时间步长Δt的取值的计算公式如下:
式中:Δxmin表示网格单元中心到对应界面的最小距离,单位为m;u为流速,单位为m/s;g为重力加速度,取值9.81m/s2;h为水深,单位为m。
5.根据权利要求4所述的一种耦合生境适宜模型的高效高精度栖息地模拟方法,其特征在于,所述步骤7具体如下:
步骤7.1、根据目标物种高精度生境因子及目标物种适宜条件,采用偏好适宜曲线法,确定出目标物种各生境因子适宜指数HSI,具体为,根据第i个控制单元网格内水深depi、流速veli和底质subi三种生境因子的模拟数据,采用偏好适宜曲线法,计算出目标物种水深、流速和底质的适宜指数 和规定生境适宜指数范围为0~1,最适宜生境适宜指数为1,最不适宜生境适宜指数为0;
步骤7.2、第i个控制单元网格内综合适宜指数CSI由下式得出:
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202111387049.3A CN114282403B (zh) | 2021-11-22 | 2021-11-22 | 一种耦合生境适宜模型的高效高精度栖息地模拟方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202111387049.3A CN114282403B (zh) | 2021-11-22 | 2021-11-22 | 一种耦合生境适宜模型的高效高精度栖息地模拟方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN114282403A true CN114282403A (zh) | 2022-04-05 |
CN114282403B CN114282403B (zh) | 2024-07-23 |
Family
ID=80869567
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202111387049.3A Active CN114282403B (zh) | 2021-11-22 | 2021-11-22 | 一种耦合生境适宜模型的高效高精度栖息地模拟方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN114282403B (zh) |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN116029469A (zh) * | 2023-03-30 | 2023-04-28 | 长江水资源保护科学研究所 | 基于水鸟生境适宜性的闸控湖泊湿地水位确定方法及装置 |
Citations (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20110046927A1 (en) * | 2007-10-30 | 2011-02-24 | University Of Utah Research Foundation | Fast iterative method for processing hamilton-jacobi equations |
US20150242545A1 (en) * | 2014-02-21 | 2015-08-27 | Junghyun Cho | Method of Simulation of Moving Interfaces using Geometry-Aware Volume of Fluid Method |
CN110414088A (zh) * | 2019-07-10 | 2019-11-05 | 武汉大学 | 结合水动力模型的涉禽栖息地适宜性空间模糊评价方法 |
US20190362035A1 (en) * | 2018-05-23 | 2019-11-28 | Nvidia Corporation | Systems and methods for computer simulation of detailed waves for large-scale water simulation |
US20200043230A1 (en) * | 2017-09-13 | 2020-02-06 | Tencent Technology (Shenzhen) Company Limited | Liquid simulation method, liquid interaction method and apparatuses |
CN111767684A (zh) * | 2020-06-30 | 2020-10-13 | 西安理工大学 | 一种优化的摩阻源项隐式格式二维浅水方程建模方法 |
CN111768502A (zh) * | 2020-07-08 | 2020-10-13 | 西安理工大学 | 一种基于gpu加速技术的非结构网格二维洪水模拟系统 |
CN112257313A (zh) * | 2020-10-21 | 2021-01-22 | 西安理工大学 | 一种基于gpu加速的污染物输移高分辨率数值模拟方法 |
CN113191054A (zh) * | 2021-04-30 | 2021-07-30 | 西安理工大学 | 一种基于显卡加速耦合管网的高精度城市雨洪模拟方法 |
-
2021
- 2021-11-22 CN CN202111387049.3A patent/CN114282403B/zh active Active
Patent Citations (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20110046927A1 (en) * | 2007-10-30 | 2011-02-24 | University Of Utah Research Foundation | Fast iterative method for processing hamilton-jacobi equations |
US20150242545A1 (en) * | 2014-02-21 | 2015-08-27 | Junghyun Cho | Method of Simulation of Moving Interfaces using Geometry-Aware Volume of Fluid Method |
US20200043230A1 (en) * | 2017-09-13 | 2020-02-06 | Tencent Technology (Shenzhen) Company Limited | Liquid simulation method, liquid interaction method and apparatuses |
US20190362035A1 (en) * | 2018-05-23 | 2019-11-28 | Nvidia Corporation | Systems and methods for computer simulation of detailed waves for large-scale water simulation |
CN110414088A (zh) * | 2019-07-10 | 2019-11-05 | 武汉大学 | 结合水动力模型的涉禽栖息地适宜性空间模糊评价方法 |
CN111767684A (zh) * | 2020-06-30 | 2020-10-13 | 西安理工大学 | 一种优化的摩阻源项隐式格式二维浅水方程建模方法 |
CN111768502A (zh) * | 2020-07-08 | 2020-10-13 | 西安理工大学 | 一种基于gpu加速技术的非结构网格二维洪水模拟系统 |
CN112257313A (zh) * | 2020-10-21 | 2021-01-22 | 西安理工大学 | 一种基于gpu加速的污染物输移高分辨率数值模拟方法 |
CN113191054A (zh) * | 2021-04-30 | 2021-07-30 | 西安理工大学 | 一种基于显卡加速耦合管网的高精度城市雨洪模拟方法 |
Non-Patent Citations (3)
Title |
---|
LU YANG等: "Application of Habitat Suitability Model Coupling with High-precision Hydrodynamic Process", 《ECOLOGICAL MODELLING》, 15 December 2021 (2021-12-15), pages 1 - 15 * |
侯精明;李桂伊;李国栋;LIANG QIUHUA;支再兴;: "高效高精度水动力模型在洪水演进中的应用研究", 水力发电学报, no. 02, 25 February 2018 (2018-02-25), pages 98 - 109 * |
郭敏鹏;侯精明;付德宇;康永德;石宝山;李昌镐;: "面源污染输移过程高性能数值模拟方法", 环境工程, no. 08, 15 August 2020 (2020-08-15), pages 163 - 169 * |
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN116029469A (zh) * | 2023-03-30 | 2023-04-28 | 长江水资源保护科学研究所 | 基于水鸟生境适宜性的闸控湖泊湿地水位确定方法及装置 |
CN116029469B (zh) * | 2023-03-30 | 2023-06-27 | 长江水资源保护科学研究所 | 基于水鸟生境适宜性的闸控湖泊湿地水位确定方法及装置 |
Also Published As
Publication number | Publication date |
---|---|
CN114282403B (zh) | 2024-07-23 |
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 | |
Hunke et al. | Sea-ice models for climate study: retrospective and new directions | |
CN108021780B (zh) | 一种基于无规则非结构网格模型的山洪动态仿真方法 | |
Nakahashi et al. | Building-cube method for large-scale, high resolution flow computations | |
Legrand et al. | High-resolution, unstructured meshes for hydrodynamic models of the Great Barrier Reef, Australia | |
CN103699714B (zh) | 一种基于有限元和无网格耦合的柔性物体实时切割仿真方法 | |
CN108563837B (zh) | 一种冲积河流水沙模型的模型参数实时校正方法和系统 | |
CN107451372B (zh) | 一种运动波与动力波相结合的山区洪水过程数值模拟方法 | |
CN111881596B (zh) | 一种基于拉格朗日插值的溢油污染源逆时追踪模拟方法 | |
CN112270115B (zh) | 一种基于元胞自动机的复杂地形洪水淹没进程模拟方法 | |
Siadatmousavi et al. | On the importance of high frequency tail in third generation wave models | |
CN115994496B (zh) | 城市公园高分辨率大气co2浓度三维场的数值模拟方法 | |
CN106934192A (zh) | 一种参数优化的浅水方程模型水体建模方法 | |
CN106569270B (zh) | 规则网格速度模型的自适应非结构三角网格化方法 | |
CN114282403A (zh) | 一种耦合生境适宜模型的高效高精度栖息地模拟方法 | |
CN114580283A (zh) | 山溪性强潮河口分汊段中长期动力地貌演变数值模拟方法 | |
Yang et al. | Application of habitat suitability model coupling with high-precision hydrodynamic processes | |
Cordonnier et al. | Forming terrains by glacial erosion | |
CN114997006A (zh) | 基于多gpu并行框架的河流泥沙输移过程数值模拟方法 | |
Tanaka et al. | A consistent finite difference local inertial model for shallow water simulation | |
CN117829028B (zh) | 浮式风机动力响应全耦合数值模拟方法、装置及相关设备 | |
McCann | Long-term morphological modelling of tidal basins |
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 |