一种缺乏资料区河渠冬季冰情发展过程的模拟方法

摘要

本发明涉及一种缺乏资料区河渠冬季冰情发展过程的模拟方法，包括：根据明渠非恒流理论，模拟无冰状态下水情发展过程，率定河道真实糙率；根据热扩散原理，模拟河道水温变化，率定河道的热交换系数hwa、jwa和kwa；以及临界弗劳德数和最大弗劳德数Frm，确定冰盖发展模式，模拟冰情演变过程。本发明根据对实测流量和水位进行模拟，推演出一种新的沿途流量分配的方法，确保了计算河段流量的平衡，另外本发明针对不同时段、不同冰情状况下太阳反射率的不同，通过模型率定推演出反应太阳反射率的系数，进而计算出时均单位面积太阳辐射值，解决了冰情计算中太阳辐射资料缺乏的难题。从而实现在计算资料不完善情况下，较为准确的计算江、河、渠道的冰情发展演变过程，为北方河渠冬季调控、灾害评估和预防提供科学参考依据。

权利要求

1.一种缺乏资料区河渠冬季冰情发展过程的模拟方法，其特征在于，所述方法包括以下步骤：

步骤1，率定糙率：根据明渠非恒流理论，模拟无冰状态下水情发展过程，率定河道真实糙率，包括如下子步骤：

S11，准备计算过程所需的包括河道断面资料、流量和水位在内的水文资料；

S12，使用无资料条件下支流流量推算法计算各汇入流量和各支流的流出流量，以及在取水口无资料条件下，进行流量分配的计算；

根据明渠非恒定流的连续性方程和运动方程进行计算，确定所述流量分配满足流量平衡和能量守恒；

S13，根据河道条件，在无冰状态下，分段采用达西‑魏斯巴赫公式和谢才公式计算明渠沿程水头损失，率定河道真实糙率n；

步骤2，率定热交换系数：根据热扩散原理，模拟河道水温变化，率定河道的热交换系数hwa、jwa和kwa，包括如下子步骤：S21，准备计算过程需要的太阳辐射、气温、水温、流量和水位在内的现有资料，使用无资料条件下每时太阳辐射量推算法和无资料条件下每时气温推算法进行以小时为时段的太阳辐射和气温的计算；

S22，计算水温，率定热交换系数hwa、jwa和kwa；

根据热扩散方程计算水温：

式中：ρ为水的密度；Cp为水的比热；A为渠道断面面积；Tw为水的断面平均温度；V为水的断面平均流速；Ex为热扩散系数；B水表面宽度； 为水体与周围环境的单位面积热交换率，包括：明流水面与大气的热交换率 水面与飘浮冰块和冰盖的热交换率 河底水体与河床的热交换率

水面与大气的单位面积热交换率采用线性传热方法近似：式中： 为明流水面与大气的热交换率， 为太阳的短波辐射热交换率，hwa、jwa和kwa为热交换系数，Ta为气温；

步骤3，模拟冰情演变：确定冰盖发展模式，模拟冰情演变过程，包括如下子步骤：S31，准备冬季计算所需的包括流量、水位、水温、太阳辐射、气温、冰花浓度、岸冰和锚冰在内的现有资料，其中流量包括冬季支流汇入流量，冬季支流汇入流量使用无资料条件下支流流量推算法进行计算；使用无资料条件下每时太阳辐射量推算法和无资料条件下每时气温推算法进行以小时为时段的太阳辐射计算和气温计算；

S32，根据包括河道过流能力、流量、水位、水温、冰盖厚度和冰盖前沿发展在内的情况，分段率定冰盖平铺上溯模式发展的最大弗劳德数Frc和水力加厚模式发展的最大弗劳德数Frm，确定冰盖的发展模式；

冰盖平铺上溯模式发展的最大弗劳德数Frc表示为：式中：v为冰盖前沿上游水流的平均流速；H为冰盖前沿的水深；ti为冰块厚度； 为系数；li为冰块长度；e为冰块孔隙率；ρ水体的密度；ρi为冰体的密度；g重力加速度；

当弗劳德数超过Frc时，冰块将出现翻转、下潜，冰盖将以水力加厚模式推进，这时冰盖初始厚度h0的计算公式为：

式中，ec＝ep+(1‑ep)e，为整个冰块的孔隙率，ep为冰块间空间间隙的孔隙率；上式存在一个最大弗劳德数Frm，当冰盖前沿弗劳德数Fr超过Frm，Fr＞Frm，冰盖不能向前发展；

S33由水流的热扩散方程、冰花扩散方程、冰盖下水流的输冰能力方程、水面浮冰的输运方程、冰盖和冰块厚度的发展方程，计算结冰前水温下降过程、冰花浓度的增长过程、冰盖前沿的推进过程、岸冰和冰厚的演变过程、开河前水温的升高过程、冰盖消融过程在内的冰情演变过程各要素；

所述的无资料条件下支流流量推算法：各支流分配的流量计算公式：

qi＝αi·(Qdown‑Qi)式中：Qdown为下游实测流量；Qi为一次洪水传播到支流边界i的流量，i＝1，2，3……f；αi为流量分配系数，

所述的无资料条件下每时太阳辐射量推算法：无云状态下的短波净太阳辐射φcl，由以下公式计算：φcl＝(0.99‑0.17m)φso‑1.253 ‑1

mo＝[sinα+0.15(α+3.885) ]式中：φcl为短波净太阳辐射；φso为每单位面积总的外来太阳辐射；Iso为太阳辐射常数；ω是时角；δ是太阳倾角；为纬度；dn为一年中的天数；m是当地气压为pa时海拔z米的光学气团；p0为水平面气压；m0是水平面处的光学气团；

太阳纬度α＝90‑θz；

其中：

ωi是前一个小时与当前小时时角的均值；E0是地球轨道的偏心校正系数；

在云的覆盖程度没有资料记录的情况下，采用每天日照数同统计时段内最大日照数比值衡量云层覆盖状态，有云状态下的太阳辐射φri计算公式为：式中：φri为有云状态下的太阳辐射；SSH为计算日日照数，SSHMAX为统计时段内最大日照数；

针对不同时段、不同河道冰情状况，采用不同的系数反应太阳反射率的变化，通过模型率定反演获得表达太阳反射变化的系数，对太阳辐射公式再次修正为φR：Rbi＝f(Rti，Rii，Wai，Ici)式中：Rbi为表达太阳反射变化的系数，分为结冰前Rb1、流冰期Rb2、封冻期Rb3、开河前期Rb4模型率定、开河期Rb5、冰盖表面有积雪覆盖Rb6；Rbi同计算河段太阳反射率Rti、河道情况Rii、河道水质Wai、冰情Ici有关；

所述的无资料条件下每时气温推算法：hn＝Y—c

式中：Tk(daytime)是白天第k个小时的温度；Tk(nighttime)是夜晚的第k个小时的温度；Y是白天时长；Z是夜晚时长；Tmax为日最高气温；Tmin是为日最低气温；Ts是日落时的气温；hn是最低温度出现后到日落的小时数；γ是日落后到最低温度时间的小时数；a是最高气温的滞后系数；b是夜间气温系数；c是从日出时开始的最低温度的滞后系数；HR为计算温度对应的小时数，取值为1‑24小时。

2.根据权利要求1所述的模拟方法，其特征在于，步骤1中S11所述的河道断面资料，当实测的河道断面间距大于1km时，则按照最大河道断面间距不超过1km，对断面进行插值加密。

说明书

一种缺乏资料区河渠冬季冰情发展过程的模拟方法

技术领域

[0001] 本发明涉及一种缺乏资料区河渠冬季冰情发展过程的模拟方法，属于水利工程技术领域。

背景技术

[0002] 归纳起来河冰发展演变包括：流凌、冰盖形成、冰盖热增长和封河、冰盖的热衰减和消融开河等过程，是复杂的水文、气象、水力学、热力学和动力学交互作用的结果，当前已
有系统的一维数学模型可以模拟冰演变过程。系统的一维模型可以推溯到美国Clarkson大
学Shen开发的RICE模型，之后改进到一维RICEN模型。在此基础上杨开林等开发的一维冰水
动力学模型实现了浮冰块在冰盖下堆积过程的模拟，郭新蕾等开发了一维树状数学模型模
拟调水工程等人工渠道冰盖发展过程。但当前天然河道冰情模拟存在最大的问题是满足模
拟需要的资料缺乏：水文站主要设立在大流域上，而且观测站之间断面间距长，用于计算的
测量数据空间步长太大；测量的水文和气象资料通常为日均值，在时间分布上，现有数据的
时间步长太大；天然江河支流纵横，但支流通常未进行水文观测，支流流入和流出的水量未
知；大多数水文站跟冰情发生发展相关的水温数据未观测或观测精度不够、误差较大；太阳
辐射、气温等相关气象数据水文站未能系统观测，只能采用当前气象局数据，气象站同水文
观测站的距离较远，数据同步性差。我国黄河是国内水文站系统较为健全、水文观测数据较
为系统的河流，黄河内蒙古河段823公里只有4个水文站，2013年以后又增加了包头水文站，
黑龙江上游900km只有8个基本水位站，主要观测黑龙江水位信息。但是20世纪90年代，黄河
水利委员会先后同芬兰和美国Shen(1994，1996)合作开发了黄河下游冰情预报经验性数学
模型，但是这些模型对实测资料精度要求高，黄河上的观测资料也无能满足计算中对实测
资料的需求，又因为黄河为游荡性河道，河床变化较快，模型建成后基本不能在实际运行中
得以应用。目前，天然河道普遍存在实测资料不全或缺失、河道资料难以测量以及河道上修
建的水工建筑物使得河道流态变化受到人为因素影响较大等条件制约。

[0003] 鉴于此，如何在计算资料不完善的情况下，较为准确的计算江、河、渠道冰盖形成、发展和消融过程，为北方河渠冬季调控、灾害评估和预防提供科学参考依据，是一个需要解
决的问题。

发明内容

[0004] 为了克服现有技术的问题，本发明提出了一种缺乏资料区河渠冬季冰情发展过程的模拟方法。所述的方法可用于江、河、渠道在计算资料不完善情况下较为准确计算冰盖形
成、发展和消融过程，为北方河渠冬季调控、灾害评估和预防提供科学参考依据。

[0005] 本发明的目的是这样实现的：一种缺乏资料区河渠冬季冰情发展过程的模拟方法，其特征在于，所述方法包括以下步骤：

[0006] 步骤1，率定糙率：根据明渠非恒流理论，模拟无冰状态下水情发展过程，率定河道真实糙率，包括如下子步骤：

[0007] S11，准备计算过程所需的包括河道断面资料、流量和水位在内的水文资料；

[0008] S12，使用无资料条件下支流流量推算法计算各汇入流量和各支流的流出流量，以及在取水口无资料条件下，进行流量分配的计算；

[0009] 根据明渠非恒定流的连续性方程和运动方程进行计算，确定所述流量分配满足流量平衡和能量守恒；

[0010] S13，根据河道条件，在春、夏、秋季，分段采用达西‑魏斯巴赫公式和谢才公式计算明渠沿程水头损失，率定河道真实糙率n；

[0011] 步骤2，率定热交换系数：根据热扩散原理，模拟河道水温变化，率定河道的热交换系数hwa、jwa和kwa，包括如下子步骤：

[0012] S21，准备计算过程需要的太阳辐射、气温、水温、流量和水位在内的现有资料，使用无资料条件下每时太阳辐射量推算法和无资料条件下每时气温推算法进行以小时为时
段的太阳辐射和气温的计算；

[0013] S22，计算水温，率定热交换系数hwa、jwa和kwa；

[0014] 根据热扩散方程计算水温：

[0015]

[0016] 式中：ρ为水的密度；Cp为水的比热；A为渠道断面面积；Tw为水的断面平均温度；V为水的断面平均流速；Ex为热扩散系数：B水表面宽度； 为水体与周围环境的单位面积热交
换率，包括：明流水面与大气的热交换率 水面与飘浮冰块和冰盖的热交换率 以及
河底水体与河床的热交换率

[0017] 水面与大气的单位面积热交换率采用线性传热方法近似：

[0018]

[0019] 式中： 为水面与大气的热交换率， 为太阳的短波辐射热交换率，hwa、jwa和kwa为热交换系数，Ta为气温；

[0020] 步骤3，模拟冰情演变：确定冰盖发展模式，模拟冰情演变过程，包括如下子步骤：

[0021] S31，准备冬季计算所需的包括流量、水位、水温、太阳辐射、气温、冰花浓度、岸冰和锚冰在内的现有资料，其中流量包括冬季支流汇入流量，冬季支流汇入流量使用无资料
条件下支流流量推算法进行计算；使用无资料条件下每时太阳辐射量推算法和无资料条件
下每时气温推算法进行以小时为时段的太阳辐射计算和气温计算；

[0022] S32，根据包括河道过流能力、流量、水位、水温、冰盖厚度和冰盖前沿发展在内的情况，分段率定冰盖平铺上溯模式发展的最大弗劳德数Frc和水力加厚模式发展的最大弗
劳德数Frm，确定冰盖的发展模式；

[0023] 冰盖平铺上溯模式发展的最大弗劳德数Frc表示为：

[0024]

[0025] 式中：V为冰盖前沿上游水流的平均流速；H为冰盖前沿的水深；ti为冰块厚度；为系数；li为冰块长度；e为冰块孔隙率；ρ水体的密度；ρi为冰体的密度；g重力加速度；

[0026] 当弗劳德数超过Frc时，冰块将出现翻转、下潜，冰盖将以水力加厚模式推进，这时冰盖初始厚度h0的计算公式为：

[0027]

[0028] 式中:ec＝ep+(1‑ep)e，为整个冰块的孔隙率，ep为冰块间空间间隙的孔隙率；上式存在一个最大弗劳德数Frm，当冰盖前沿弗劳德数Fr超过Frm(Fr＞Frm)，冰盖不能向前发展；

[0029] S33由水流的热扩散方程、冰花扩散方程、冰盖下水流的输冰能力方程、水面浮冰的输运方程、冰盖和冰块厚度的发展方程，计算结冰前水温下降过程、冰花浓度的增长过
程、冰盖前沿的推进过程、岸冰和冰厚的演变过程、开河前水温的升高过程、冰盖消融过程
在内的冰情演变过程各要素；

[0030] 所述的无资料条件下支流流量推算法：

[0031] 各支流分配的流量计算公式：

[0032] qi＝αi·(Qdown‑Qi)

[0033] 式中：Qdown为下游实测流量；Qi为一次洪水传播到支流边界i的流量，i＝1，2，3……n；αi为流量分配系数，

[0034] 所述的无资料条件下每时太阳辐射量推算法：

[0035] 无云状态下的短波太阳辐射φcl，由以下公式计算：

[0036] φcl＝(0.99‑0.17m)φso

[0037]

[0038]

[0039] mo＝[sinα+0.15(α+3.885)‑1.253]‑1

[0040]

[0041]

[0042]

[0043] 式中：φcl为短波净太阳辐射；φso为每单位面积总的外来太阳辐射；Iso为太阳辐射常数；ω是时角；δ是太阳倾角；为纬度度；dn为一年中的天数；m是当地气压为pa时海拔z
米的光学气团；p0为水平面气压；m0是水平面处的光学气团；

[0044] 太阳纬度α＝90‑θz；

[0045] 其中：

[0046] ωi是前一个小时与当前小时时角的均值；E0是地球轨道的偏心校正系数；

[0047] 在云的覆盖程度没有资料记录的情况下，采用每天日照数同统计时段内最大日照数比值衡量云层覆盖状态，有云状态下的太阳辐射φri计算公式为：

[0048]

[0049] 式中：φri为有云状态下的太阳辐射；SSH为计算日日照数，SSHMAX为统计时段内最大日照数；

[0050] 针对不同时段、不同河道冰情状况，采用不同的系数反应太阳反射率的变化，通过模型率定反演获得表达太阳反射变化的系数，对太阳辐射公式再次修正为φR：

[0051]

[0052] Rbi＝f(Rti，Rii，Wai，Ici)

[0053] 式中：Rbi为表达太阳反射变化的系数，分为结冰前Rb1、流冰期Rb2、封冻期Rb3、开河前期Rb4模型率定、开河期Rb5、冰盖表面有积雪覆盖Rb6；Rbi同计算河段太阳反射率Rti、河道
情况Rii、河道水质Wai、冰情Ici有关；

[0054] 所述的无资料条件下每时气温推算法：

[0055]

[0056]

[0057] m＝Y‑c

[0058]

[0059] 式中：Tk(daytime)是白天第k个小时的温度；Tk(nighttime)是夜晚的第k个小时的温度；Y是白天时长；Z是夜晚时长；Tmax为日最高气温；Tmin是为日最低气温；Ts是日落时的气
温；m是最低温度出现后到日落的小时数；n是日落后到最低温度时间的小时数；a是最高气
温的滞后系数；b是夜间气温系数；c是从日出时开始的最低温度的滞后系数；HR为计算温度
对应的小时数，取值为1‑24小时。

[0060] 本发明的优点和有益效果是：本发明根据上游流量向下游传播过程的特征，推演出缺少支流流量条件下河道区间支流流量分配的合理比例，找到一种新的区间支流流量进
行动态分配的方法，确保了计算河段流量的平衡；本发明采用每天日照数同统计时段内最
大日照数比值衡量云层覆盖状态，解决了太阳辐射率计算中云覆盖度没有资料记录的问
题；本发明针对不同时段、不同河道冰情太阳反射率不同，通过模型率定得到表达太阳反射
率的系数。在冰情模拟中这个系数分为结冰前Rb1、流冰期Rb2、封冻期Rb3、开河前期Rb4、开河
期Rb5、冰盖表面有积雪覆盖Rb66种反射率，减少了热扩散计算中存在误差。从而实现在计算
资料不完善情况下，较为准确的计算江、河、渠道的冰盖形成、发展和消融过程，为北方河渠
冬季调控、灾害评估和预防的提供科学参考依据。

附图说明

[0061] 下面结合附图和实施例对本发明作进一步说明。

[0062] 图1是本发明实施例一计算河系示意图；

[0063] 图2是本发明实施例一支流流量分配前后流量模拟曲线；

[0064] 图3为本发明实施例一计算河系中气候区域分区示意图；

[0065] 图4为巴彦高勒‑三湖河口支流分水口位置示意图；

[0066] 图5为三湖河口‑头道拐支流分水口位置示意图；

[0067] 图6为所述黄河内蒙古河段的气象数据分区示意图；

[0068] 图7为三湖河口水温演变过程示意图；

[0069] 图8为头道拐水温演变过程示意图；

[0070] 图9为冰盖前沿发展过程示意图；

[0071] 图10为巴彦高勒冰厚发展过程示意图；

[0072] 图11为三湖河口冰厚发展过程示意图。

具体实施方式

[0073] 实施例一：

[0074] 本实施例是一种缺乏资料区河渠冬季冰情发展过程的模拟方法。所述的方法关键在于提出了无资料条件下支流流量推算法、无资料条件下每时太阳辐射量推算法、无资料
条件下每时气温推算法，在缺乏详细资料的情况下，对支流流量、每时的太阳辐射量和每时
气温的演变进行推算，以弥补没有详细资料而造成的冬季冰情无法模拟的缺憾，为防灾减
灾提供了科学的依据。

[0075] 本实施例包括以下步骤：

[0076] 步骤1，率定糙率：根据明渠非恒流理论，模拟无冰状态下水情发展过程，率定河道真实糙率，包括S11～S13三个子步骤：

[0077] S11，准备计算过程所需的包括河道断面资料、流量和水位在内的水文资料。

[0078] 天然河道冰情模拟应用时常常会遇到一个难题，即缺乏细致系统的河道断面资料、边界条件、水文、气象和水力学基础资料。具体表现一是国内水文站主要设立在较大流
域上，而且站与站之间断面间距长，用于计算的水文站点数据空间步长大；二是测量的水文
和气象资料通常为日均值或者每日最大、最小值，在时间分布上现有数据的时间步长太大；
三是天然江河支流纵横，但支流通常未设水文观测站，支流流入和流出的水量未知；四是水
文站观测的水温数据通常精度不够或数据不完整；五是太阳辐射和气温等相关气象资料采
用气象局发布的信息，气象站同水文观测站的距离较远，数据同步性差。因此需要对现有资
料的断面进行推演，推算出缺少的水文气象资料，以供精细的计算。

[0079] S12，使用无资料条件下支流流量推算法计算各汇入流量和各支流的流出流量，以及在取水口无资料条件下，进行流量分配的计算。

[0080] 根据明渠非恒定流的连续性方程和运动方程进行计算，确定所述流量分配满足流量平衡和能量守恒。

[0081] 所述的无资料条件下支流流量推算法：

[0082] 天然大河支流纵横，大部分支流都未有流量观测数据，导致河道分叉、流入、流出流量数据信息缺乏，因此需要根据对实测流量和水位进行模拟，找到沿途流量进行分配的
方法，满足计算中流量的平衡。

[0083] 明渠非恒定流的连续性方程和运动方程(以流量Q和水位Z表示)是：

[0084]

[0085]

[0086] 式中：x为距离(m)，t为时间(s)，z为水位(m)，Q为流量(m3/s)，A为过流面积，B为水面宽，ql为单位流程上的侧向出流量，负值表示流入； 是断面扩缩引起的，如果渠道
是棱柱体，则方程无此项。

[0087] 在进行数值计算时，需要给出问题的初始条件和上、下游边界和支流的流入和流出边界条件。边界条件通常有三种类型，分别为水位边界条件Z＝Z(t)，流量边界条件Q＝Q
(t)，水位流量边界条件Q＝Q(z)。计算中通常采用上游边界和支流边界给定初始流量、下游
边界给定初始的水位的方法，前后计算断面应满足流量平衡和能量守恒。

[0088] 河系如图1所示，其中边界条件6个：上游边界①，下游边界②，支流边界③‑⑥，在冰情模拟时只有上游边界和下游边界的流量信息，因为支流③‑⑥流量未监测，计算中为了
满足流量平衡，需要对区间汇入和流出的流量进行合理分配。如果仅仅简单采用上、下游边
界流量差进行比例分配的话，忽略了实际流量自上游向下游传递的时间差导致分配方案的
不合理。本实施例采用一种新的流量分配方法，进行支流流量分配，从而确保计算河段流量
的平衡。

[0089] 在计算中只输入上游边界条件①的流量过程和下游边界②的水位条件，模拟上游边界①一次洪水传播到支流边界③‑⑥的流量Qi(i为支流边界，i＝1，2，3，4。下同)和下游
边界②的流量Qdown曲线(如图2所示)，图2中可以看出，下游边界②模拟流量和实测流量存
在明显差异，因为模拟中区间支流流量未加入，导致流量存在不平衡。为此，需给支流边界
③‑⑥分配适当流量，从而使得此段流量基本达到平衡。由于支流边界③‑⑥并无流量观测
数据，采用实测下游边界②的流量Qdown分别减去一次洪水波传递到③‑⑥断面处的模拟流
量，差值按照比例αi分配到支流边界③‑⑥qi，确保计算河段流量平衡，各支流分配的流量
为：

[0090] qi＝αi·(Qdown‑Qi)   (3)

[0091] 式中：Qdown为下游实测流量；Qi为一次洪水传播到支流边界i的流量；αi为流量分配系数，

[0092] 流量分配系数αi的确定可通过调整各支流边界αi值以计算下游边界②流量模拟值，再与实测值对比，以使关键时间点上的流量值最接近率定得到，或可以尝试利用最小二
乘法寻优达到相对准确辨识区间支流流量分配比例的目的，确保计算河段流量的平衡。通
过支流边界③‑⑥流量的合理分配，计算得到支流分配流量后下游边界②的模拟值如图2所
示，由图可知，下游边界②三湖河口水文站模拟值和实测值吻合良好。

[0093] S13，根据河道条件，在无冰状态下，分段采用达西‑魏斯巴赫公式和谢才公式计算明渠沿程水头损失，率定河道真实糙率n；

[0094] 所述达西‑魏斯巴赫水头损失公式是：

[0095]

[0096] 式中：hf为沿程水头损失，m；f为明渠沿程阻力系数；R为明渠水力半径，m；V为明渠2 2
断面平均流速，m/s；L为明渠长度，m；g为重力加速度，m/s，一般取g＝9.8m/s。

[0097] 所述谢才公式示出明渠沿程阻力系数f与糙率n的关系：

[0098]

[0099] 步骤2，率定热交换系数：根据热扩散原理，模拟河道水温变化，率定河道的热交换系数hwa、jwa和kwa，包括S21、S22两个子步骤：

[0100] S21，准备计算过程需要的太阳辐射、气温、水温、流量和水位在内的现有资料，使用无资料条件下每时太阳辐射量推算法和无资料条件下每时气温推算法进行以小时为时
段的太阳辐射和气温的计算。

[0101] 因为太阳辐射和气温是热交换主要能量来源，对水温和冰情演变快慢有着决定性作用，这两个数据观测资料远远不能满足计算需求，本实施例在太阳辐射和气温无资料条
件下，进行每时太阳辐射和每时气温的计算。

[0102] 较为系统的太阳辐射资料来自气象局，是以日为单位的净太阳辐射，计算中需要精度为每小时为单位太阳辐射值，另外气象站和水文站距离较远，气象资料和水文资料数
据存在位置不同，导致气象数据和真实的计算断面数据不一致。

[0103] 所述的无资料条件下每时太阳辐射量推算法：

[0104] 无云状态下的短波太阳辐射φcl，由以下公式计算：

[0105] φcl＝(0.99‑0.17m)φso   (6)

[0106]

[0107]

[0108] mo＝[sinα+0.15(α+3.885)‑1.253]‑1   (9)

[0109]

[0110]

[0111]

[0112] 式中：φcl为短波净太阳辐射，单位为W/m2；φso为每单位面积总的外来太阳辐射，2 2
单位为J/m ；Iso为太阳辐射常数，冬季约为1380W/m ；ω是时角，午时为0°，每小时变化15°，
早上为正值，下午为负值；δ是太阳倾角，单位为弧度；为纬度度，北纬为正，南纬为负；dn为
一年中的天数，从1月1日算起，如1月20日＝20days；m是当地气压pa时海拔z米的光学气团，
单位为m；p0为水平面气压，单位为pa；m0是水平面处的光学气团；太阳纬度α＝90‑θz；
为弧度；ωi是前一个小时与当前小时时角的
均值；E0是地球轨道的偏心校正系数。

[0113] 通常在冰情演变过程计算中，根据气象站设置的数量和分布位置，需将河段分为尽可能多的区域，尽量确保计算断面使用的太阳辐射值的准确度，每个区域计算一个站点
的太阳辐射和气温，区域划分见图3所示。

[0114] 在云的覆盖程度没有资料记录的情况下，采用每天日照数同统计时段内最大日照数比值衡量云层覆盖状态，有云状态下的太阳辐射φri计算公式为：

[0115]

[0116] 式中：φri为有云状态下的太阳辐射；SSH为计算日日照数，SSHMAX为统计时段内最大日照数。

[0117] 在冰情演变过程中太阳对冰面和水面的反射率变化情况复杂，不仅跟冰的表面状况有关，还和河流中水的含沙量、浑浊度、悬浮物浓度等相关。实际计算中不同的河流不同
冰状态反射率都不相同，一维模型计算中考虑太阳辐射时通常忽略太阳反射率的影响，导
致热扩散计算中存在误差，本实施例针对不同时段、不同河道冰情状况，采用不同的系统反
应太阳反射率的变化，通过模型率定反演得到该系数。在冰情模拟中这个系数分为结冰前
Rb1、流冰期Rb2、封冻期Rb3、开河前期Rb4、开河期Rb5、积雪覆盖Rb66种，不同时段冰水情状态如
表1所示。

[0118] 表1不同冰情状态表达太阳反射率的系数

[0119]

[0120]

[0121] 当太阳辐射到达水面或冰面，太阳辐射公式再次修正为φR：

[0122]

[0123] Rbi＝f(Rti，Rii，Wai，Ici)   (15)

[0124] 式中：Rbi为表达太阳辐射变化的系数，分为结冰前Rb1、流冰期Rb2、封冻期Rb3、开河前期Rb4、开河期Rb5、冰盖表面有积雪覆盖Rb66种反射率，通过模型率定得到；Rbi同计算河段
太阳反射率Rti、河道情况Rii、河道水质Wai，冰情Ici有关。

[0125] 所述的无资料条件下每时气温推算法为：

[0126] 无资料条件下每时气温的计算，采用如下气温计算公式：

[0127]

[0128]

[0129] m＝Y‑c  (18)

[0130]

[0131] 式中：Ti(daytime)是白天第k个小时的温度，Ti(nighttime)是夜晚的第k个小时的温度；Y是白天时长，小时；Z是夜晚时长，小时；Tmax和Tmin是日最高和最低气温；Ts是日落时
的气温；m是最低温度出现后到日落的小时数；n是日落后到最低温度时间的小时数；a是最
高气温的滞后系数；b是夜间气温系数；c是从日出时开始的最低温度的滞后系数；HR为计算
温度对应的小时数，取值为1‑24小时；系数a、b、c根据逐日最大气温和最小气温实测值率定
得到。

[0132] S22，计算水温，率定热交换系数hwa、jwa和kwa。

[0133] 根据热扩散方程计算水温：

[0134]

[0135] 式中：ρ为水的密度；Cp＝4.2kJ/为水的比热，kg·℃；A为渠道断面面积；Tw为水的断面平均温度；V为水的断面平均流速；Ex为热扩散系数：B水表面宽度； 为水体与周围环
境的单位面积热交换率，包括：明流水面与大气的热交换率 水面与飘浮冰块和冰盖的
热交换率 大气与冰盖的热交换率 以及河底水体与河床的热交换率

[0136] 水面与大气的单位面积热交换率采用线性传热方法近似：

[0137]

[0138] 式中： 为水面与大气的热交换率， 为太阳的短波辐射热交换率，hwa、jwa和kwa为热交换系数，Ta为气温。

[0139] 步骤3，模拟冰情演变：确定冰盖发展模式，模拟冰情演变过程，包括S31～S33三个子步骤：

[0140] S31，准备冬季计算所需的包括流量、水位、水温、太阳辐射、气温、冰花浓度、岸冰和锚冰在内的现有资料，其中流量包括冬季支流汇入流量，冬季支流汇入流量使用无资料
条件下支流流量推算法进行计算；使用无资料条件下每时太阳辐射量推算法和无资料条件
下每时气温推算法进行以小时为时段的太阳辐射计算和气温计算。

[0141] S32，根据包括河道过流能力、流量、水位在内的情况，分段率定冰盖平铺上溯模式发展的最大弗劳德数Frc和水力加厚模式发展的最大弗劳德数Frm，确定冰盖的发展模式。

[0142] 冰盖平铺上溯模式发展的最大弗劳德数Frc表示为：

[0143]

[0144] 式中：V为冰盖前沿上游水流的平均流速；H为冰盖前沿的水深；ti为冰块厚度；为系数，取值在0.66～1.3之间变化；li为冰块长度；e为冰块孔隙率；ρ水体的密度；ρi
为冰体的密度；g重力加速度。

[0145] 当弗劳德数超过Frc时，冰块将出现翻转、下潜，冰盖将以水力加厚模式推进，这时冰盖初始厚度h0的计算公式为：

[0146]

[0147] 式中，ec＝ep+(1‑ep)e，为整个冰块的孔隙率，ep为冰块间空间间隙的孔隙率；上式存在一个最大弗劳德数Frm，当冰盖前沿弗劳德数Fr超过Frm(Fr＞Frm)，冰盖不能向前发展。

[0148] S33，在上述步骤1、步骤2的数据准备和参数率定的基础上，由水流的热扩散方程、冰花扩散方程、冰盖下水流的冰情能力方程、水面浮冰的输运方程、冰盖和冰块厚度的发展
方程，计算结冰前水温下降过程、冰花浓度的增长过程、冰盖前沿的推进过程、岸冰和冰厚
的演变过程、开河前水温的升高过程、冰盖消融过程等冰情演变过程。

[0149] 模拟举例：

[0150] “黄河宁，天下平”，治理黄河，自古就是安民兴邦的大事。黄河是中国凌汛出现最为频繁的河流之一，凌汛问题一直受到中央和水利部的高度重视，主要凌汛灾害集中在最
北端的内蒙古河段。黄河内蒙古河段位居倒“U”字型河道大拐弯最底部，全长823km多，冬季
最低气温可达‑35℃，特殊的地理位置导致内蒙古河段封河自下游向上游、开河自上游向下
游发展。因此，无论从地理位置、河道走势还是水文气象条件看内蒙古河段都是黄河冰情最
为严重的区域。

[0151] 近年来，受气候变化、人类活动等影响，内蒙古河段的冰情特点产生了很大的变化，出现了许多以前未有过的新情况，主要表现在下述几个方面：1)冬季平均气温偏高，气
温变幅大；2)河道淤积严重；3)槽蓄水增量增大。这些新特点的出现说明内蒙古河段的冰情
规律发生了明显变化，亟待建立能反映冰凌物理特性变化影响因素和规律的数学模型，并
对冰凌生消变化的物理过程进行模拟，才能更准确地把握其变化规律，作出较准确的凌情
演变计算和冰情发展的预报。

[0152] 一、黄河宁蒙河段断面和河道情况

[0153] 黄河内蒙古河段最新河道断面数据为2012年10月的大断面资料，从巴彦高勒到头道拐共测量断面166个，因此本实施例研究冰情计算、验证、校核及冰情模拟以2012‑2013凌
汛期作为研究对象。

[0154] 计算河段巴彦高勒至头道拐，在实测166个断面基础上增加桥梁断面18个，浮桥断面3个，总计断面187个，模拟河段只有巴彦高勒、三湖河口和头道拐三个水文站。通过调整
桥梁和浮桥处断面的糙率，考虑桥梁和浮桥处的阻力作用。

[0155] 实测断面间距较大会影响结冰期模型的稳定性及模拟结果，因此按照最大断面间距不超过1km，对断面进行线性插值，最终得到河道断面580个。计算得到的是两个断面之间
的直线距离，没有充分考虑河道弯曲情况，河道长度约为486km，落差68m，平均坡度0.14‰。

[0156] 二、无资料情况下支流流量分配

[0157] 因模拟河段汇入支流众多，经过国家重要的畜牧业产区和重工业基地，需要取水量大，缺少沿途支流取水或者退水的流量资料，因此需要通过对巴彦高勒、三湖河口和头道
拐三个水文站实测真实流量和水位进行模拟，对沿途支流的流量进行合理分配。为了确保
巴彦高勒和三湖河口区间、三湖河口和头道拐之区间流量的平衡，使用无资料条件下支流
流量推算法将流量进行分配。计算中选取巴彦高勒‑三湖河口之间4个支流边界，分别为二
闸、三闸、四闸和乌梁素海分水口(图4)，按照0.1:0.1:0.3:0.5的比例进行流量分配。三湖
河口和头道拐之间主要分水支流边界6个，分别为打不素扬水站、包钢水源地取水口、画匠
银子岸边取水泵房、磴口净水厂、团结渠电力扬水站、民利渠扬水站分水口(图5)，按照0.1:
0.1:0.1:0.1:0.3:0.3的比例进行分配，通过以上流量的合理分配确保了计算河段流量的
平衡。

[0158] 三、无资料情况下太阳辐射和气温的模拟

[0159] 由于研究河段较长，流经区域气候条件有明显差异，因此需要对计算流域划分不同的区域，不同区域输入对应气象站的气温和太阳辐射数据，根据气象站的分布将计算区
域划分为8个如图6所示。由于目前缺乏实测的逐时气温和太阳辐射数据，采用无资料条件
下每时太阳辐射量推算法和无资料条件下每时气温推算法推算出所需的太阳辐射和气温
值资料。

[0160] 四、黄河内蒙古河段冰情模拟

[0161] 冰情模拟除了需要河段断面资料，要需要完备系统的水文、气象资料和率定的关键参数。

[0162] 所需水文、气象资料：每时气温、每时水温、每时太阳辐射、每时流量、每时水位、流凌密度和岸冰、锚冰和冰盖的发展情况。

[0163] 率定参数：河床糙率，热交换系数hwa、jwa、kwa，冰盖平铺上溯模式发展的最大弗劳德数Frc和水力加厚模式发展的最大弗劳德数Frm等。

[0164] 冰情演变过程计算步骤见上述步骤1‑3，模拟冰情发展过程，包括：水温、冰花密度、岸冰和锚冰的发展、冰厚(冰盖产生、增长和消融过程)、流量、水位、冰盖前沿的发展过
程等。

[0165] 三湖河口和头道拐水温计算见图7、8所示，整体上看，计算值和实测值吻合较好，只有1‑2天出现相对大点误差，三湖河口气温的实测值和模拟值最大误差0.29℃，头道拐气
温的实测值和模拟值最大误差0.56℃。

[0166] 冰盖于2012年11月30日在昭君坟上游附近首封，首封断面距离巴彦高勒287km，随着气温下降，冰盖前沿不断向前推进，到12月23日封冻至巴彦高勒，冰盖前沿发展过程如图
9所示。从计算图上看出，冰盖前沿发展过程模拟值和实测值吻合非常好。

[0167] 水文站测量冰厚仍采用冰面上打孔测量，因此测量冰厚要在冰盖有足够强度、能够支撑人在上面走动时进行，所以冰盖厚度20‑30cm之后才能测量，冰厚测量为每周一次，
巴彦高勒在冰盖增长到稳定期冰厚实测值仅6个，三湖河口冰盖实测值10个，冰厚实测值和
模拟值的比较见图10‑11，巴彦高勒和三湖河口冰盖厚度测量值和模拟值误差如表2所示，
冰厚实测值和模拟值最大误差未超0.065m，模拟效果好。

[0168] 表2冰厚实测值和模拟值

[0169]

[0170] 最后应说明的是，以上仅用以说明本发明的技术方案而非限制，尽管参照较佳布置方案对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术
方案(各种公式的运用、步骤的先后顺序等)进行修改或者等同替换，而不脱离本发明技术
方案的精神和范围。