УДК 556.5.048
А.М. Бондаренко
МОДЕЛЬ ФОРМИРОВАНИЯ ПОТЕРЬ ДОЖДЕВЫХ ПАВОДКОВ
Исследована зависимость потерь дождевых паводков от суммарного слоя осадков и предшествующего увлажнения водосбора.
Выполнено математическое моделирование потерь дождевого стока с оптимизацией параметров модели на примере сахалинских малых рек и дальневосточных рек Приморья. Описан алгоритм оптимизации параметров модели методом случайного поиска.
Дождевой сток, слой осадков, оптимизация параметров, случайный поиск
A.M. Bondarenko MODEL FOR RAINFALL RUNOFF LOSS FORMATION
The rainfall runoff loss dependence of the total precipitation and previous watering of watershed area has been studied. A mathematical model of rainfall loss, which parameters can be improved through the optimization methods, has been stated as exemplified by Sakhalin minor rivers and the Far Eastern rivers of the Primorie Territory. The algorithm of model parameters optimization by the means of random search has been described.
Rainfall runoff, precipitation, parameter optimization, random search
Без определения потерь стока на склонах водосборов невозможны прогноз и расчет речного стока. Потерями стока в русле при прохождении паводка на малых водосборах можно пренебречь, поэтому слой стока Y вычисляется как разность слоя водообразования ST на склонах в период дождя и потерь стока после периода водообразования (на спаде склонового стока) Pсп по формуле:
Y = Хт - Рт - Рсп = St - Рсп , (1)
где ХТ и Рт — слой осадков и слой потерь дождевого стока в период водообразования Т. Прогнозирование слоя стока позволяет, в свою очередь, выполнить оценку максимального модуля руслового стока по формулам стока или с применением корреляционных связей.
При расчете потерь стока на водосборе пренебрегают потерями дождевых вод на перехват растительностью и потерями на испарение в период прохождения паводка, как очень малыми величинами. Анализируется только основная составляющая потерь -потери на впитывание. Процесс впитывания на горных склонах состоит из двух фаз [1]. В первую фазу входит заполнение свободной емкости слоя рыхлых отложений и подстилающего относительного водоупора до наименьшей влагоемкости. Вторую фазу можно считать установившейся, и связана она в основном с действием гравитационных сил. Обе фазы впитывания нашли свое отражение в известной формуле интенсивности впитывания, получившей теоретическое и эмпирическое обоснование в работах Г.А. Алексеева, А.Н. Бефани [1,2] и других:
к, = ко + А / гп, (2)
где ко - гравитационная квазиконстанта, не зависящая от времени; А - параметр, зависящий от дефицита влажности почвы, капиллярного напора и коэффициента фильтрации. Показатель редукции впитывания п зависит от типа почв. В [2, 3] показано, что параметр А связан, прежде всего, с индексом предшествующего увлажнения речного водосбора 1ц. Оценивая величины суммарных потерь стока на горных водосборах, прогнозист имеет в своем распоряжении, как правило, только данные наблюдений на дождемерах, регистрирующих суммарные слои осадков. Поэтому имеет смысл рассматривать формулу впитывания в ином виде, разделив ее составляющие на интенсивность осадков, что позволит ввести в расчет непосредственно слой выпавших осадков, как один из основных факторов, определяющих потери дождевого стока. Интенсивность осадков ш можно представить в дифференциалах осадков и времени, как ёх/& и записать отношение
к( _ кДі _ ДР а аДі Дх
(3)
где ДР и Дх - дифференциалы слоя потерь стока и слоя осадков.
Отношение -к- есть мгновенный коэффициент потерь, поскольку — = 1 - —-, где
а а а
к
Ьі - интенсивность водообразования в момент времени і. Введем обозначение кх = — и
а
,„4 , АДі ,, кп Дх тх , ,
запишем (2) в виде кх = ка+-, где ка=— и аі = —. Интегрируя дифференциал
іпДх аг Ді
потерь ДР = кхДх, для фазы водообразования получаем выражения
РТ Хт Т
| ДР = | к'аДх +1Д*, (4)
0 0 0 ^
А Т1~ п
Рт = кО Хт + -----. (5)
1 - п
В течение всего периода водообразования коэффициент потерь кх меняется за счет изменчивости притока дождевых осадков. Он меняется даже тогда, когда интенсивность впитывания постоянна и равна кО. Однако первое слагаемое в формуле (5) получено при условии постоянства мгновенного коэффициента потерь в период установившейся стадии впитывания, когда интенсивность притока осадков незначительно меняется с момента установления к= кО. Такое условие можно принять для районов с преимущественным выпадением дождей обложного характера - таких, как районы дальневосточного Приморья и о. Сахалин с преобладанием дождей муссонного происхождения, для которых характерна незначительная вариация интенсивности осадков. К тому же на горных склонах обычно формируются рыхлые хорошо дренирующие почвы и соответственно подповерхностный (так называемый «контактный») сток на контакте с малопроницаемыми породами относительного водоупора. Это значительно снижает эффект влияния интенсивности дождя [1]. Поэтому при таких условиях формирования стока допустимо принять, что динамика мгновенного коэффициента потерь во времени аналогична динамике интенсивности впитывания и также характеризуется двумя стадиями - неустановившейся и установившейся.
В период спада паводка, начиная с момента окончания дождя, потери стока можно отнести к установившейся фазе впитывания, но при этом затопленная площадь водосбора, участвующая в формировании контактного (подповерхностного) стока, уменьшается за счет поперечной и продольной концентрации в микротальвегах. В связи с этим потери на спаде должны характеризоваться величиной установившегося коэффициента потерь, 138
осредненного не по всей площади, а лишь по затопленной - к'0. Поэтому потери на спаде можно представить формулой Рсп = к'0 Хсп, где Хсп - доля слоя осадков, не успевшая стечь
за период водообразования и соответствующая объему воды на склонах после окончания дождя. Выполнить непосредственный расчет потерь на спаде не представляется возможным из-за трудностей определения величины Хсп, поэтому применяют косвенные методы оценки этой величины. Предлагается определять значение Рсп графически по отрезку на оси ординат, отсекаемому нижней огибающей потока точек на связи общих потерь стока со слоем осадков (рис. 1). Возможна также оптимизация Рсп как параметра математической модели потерь стока.
Общий (суммарный) слой потерь стока Р складывается из потерь неустановившейся фазы АР, соответствующих второму слагаемому формулы (5), потерь установившейся фазы к'аХт и потерь на спаде Рсп:
Р = АР + кО Хт + Рсп . (6)
Связь суммарных потерь стока со слоем осадков хорошо иллюстрирует рис. 1 для р. Холмской (Р = 12,8 км2) на о. Сахалин, на котором показаны составляющие формулы (6). Такая графическая интерпретация общих потерь стока оказалась возможной потому, что для горных водосборов малых рек очевидна зависимость АР от предпаводочного увлажнения водосбора. Потери АР идут на увлажнение почвы от предпаводочной влажности до величины наименьшей влагоемкости в рыхлом слое почвогрунтов и в относительном водоупоре. Среднее по водосбору предпаводочное увлажнение почв непосредственно определить невозможно. Поэтому применяют так называемые индексы предшествующего увлажнения. В данной работе использован индекс, который обоснован
теоретически и нашел широкое применение в гидрологических прогнозах. В общем виде
т
индекс 1Ш = ^ ах, где х, - слой осадков за интервал времени А^-; Т - общая
1
продолжительность времени учета осадков; а, - доля осадков, оставшихся в почве на расчетную дату, удаленную на / суток от начала паводка. Предлагается вести расчет за период Т = 60 суток. Коэффициенты а, определены эмпирически в зависимости от типа почв и температур воздуха и рекомендованы в [2, 3].
В качестве характеристики, косвенно отражающей предпаводочное увлажнение почв, можно использовать и величину предпаводочного расхода воды в реке. На рис. 2 показан пример хорошей корреляции связи 1№ = Д (0п) для р. Холмской.
С целью построения математической модели потерь стока выполнен анализ ряда паводков за годы наблюдений на малых реках Приморья и о.Сахалин. Для каждого паводка вычислялись общие потери стока как разности между слоями осадков за паводок Хт и слоями стока У: Р = ХТ - У. Слои осадков вычислялись по данным дождемеров на исследуемых водосборах, а слои паводочного стока определялись по гидрографам путем отсечения грунтового питания и шлейфа спада предыдущего паводка. Для всех рек строились зависимости Р = ДХт, 1щ), пример которой дан на рис. 1 для р. Холмской. Далее выполнялась оценка величин Рсп и к'о по уравнению нижней огибающей потока точек.
Потери неустановившейся стадии впитывания определялись как разность АР = Р - к'аХт -Рсп. Эти потери идут на увлажнение почвы от предпаводочной влажности до наименьшей влагоемкости, поэтому их зависимость от предпаводочного увлажнения обычно хорошо выражена для данных условий формирования стока.
P, мм 60
АР
ko' Хт
ХТ, мм
Рис. 1. Зависимость слоя общих потерь (Р) от слоя дождевых осадков для р. Холмской (у точек - значения индекса предшествующего увлажнения !„)
Рис. 2. Связь индекса /„ с предпаводочным расходом для р. Холмской (/^ = ЦОп) - уравнение регрессии)
С увеличением слоя осадков способность почвы удерживать влагу восполняется на постепенно увеличивающихся площадях, приближаясь к своему предельному значению. Следовательно, величина АР зависит не только от предпаводочного увлажнения, но и от слоя осадков. Затухающее влияние осадков на потери неустановившейся фазы впитывания хорошо описывается математической моделью, включающей функцию гиперболического тангенса:
АР = j(Iw) Pmax th (Хт/Pmax). (7)
Пример зависимости АР = f(XT, IW) приведен на рис. 3.
Максимально возможные потери Pmax неустановившейся фазы впитывания соответствуют осредненной по водосбору величине наименьшей влагоемкости рыхлых почв и грунтов подстилающего водоупора. Эта величина определяется графическим анализом по верхней огибающей потока точек связи АР = f(XT,IW) - рис. 3 с учетом точек с наименьшим IW. Введем обозначение Рth = Pmax th(XT/Pmax). Тогда множитель, учитывающий в формуле (7) фактор увлажнения, может быть определен отношением j(IW) = АР/Рй. Такие отношения вычислены по паводкам исследуемых водосборов для построения графиков связи АР/Рth = f(IW). Пример подобной зависимости приведен на рис. 4, на котором переменные корреллируются линейной функцией регрессии j(IW) =1,0 - с IW до некоторого предельного для данного водосбора значения IWmax = 55 мм. Такая структура зависимости вполне соответствует природе формирования потерь на горных 140
склонах, поскольку индекс предшествующего увлажнения косвенно отражает предпаводочный дефицит влажности почв и грунтов на водосборе. При достаточно высокой увлажненности водосбора перед паводком (1ц тах) потери стока приобретают преимущественно гравитационный характер и относятся только к установившейся фазе формирования потерь. В этом случае составляющая модели АР принимается равной нулю.
Рис. 3. Зависимость слоя потерь неустановившейся фазы (АР) от слоя ХТ дождевых осадков для р. Холмской (у точек - значения индекса !„)
А Р/Рт
Рис. 4. Зависимость относительных потерь неустановившейся фазы впитывания от индекса предшествующего увлажнения /„ для р. Холмской (АР/Р№ = f (/^) - уравнение регрессии)
Эмпирическая формула общих потерь стока в соответствии со структурой формулы (6), полученная для р. Холмской (Р = 12,8 км2):
Р = (1 - 0,02 1ц) 50 й (Хт/50) + 0,14 Хт +6 , (8)
то есть принято с = 0,02, Ртах =50 мм, к'о = 0,14, Рсп =6 мм.
Для небольших паводков со слоем осадков Хт <30 мм с незначительным начальным увлажнением водосбора 1ц < 20 мм следует учитывать только первое слагаемое формулы (8), принимая Р = АР. При сильно увлажненных почвах с индексом 1ц > 55 мм расчет выполняется без учета первого слагаемого, то есть Р= к'о Хт + Рсп. Потери на спаде достаточно весомы в балансе потерь только для очень малых рек горных водосборов Приморья и о. Сахалин с густотой речной сети менее 3 км/км2. При большей
густоте речной сети потери на спаде паводка не превышают 1-2 мм. При таких условиях формирования стока потерями на спаде можно пренебречь и расчет вести по двучленной формуле Р = АР + к'о Хт.
Графическая оценка параметров модели потерь Ртах, с, к'о отчасти субъективна, поэтому задачу определения этих параметров лучше рассматривать как задачу автоматизированной оптимизации с применением критерия оптимизации (целевой функции).
При решении гидрологических задач необходимо учитывать ограничения на оптимизируемые параметры. Достаточно просто такие ограничения учитываются в алгоритмах случайного поиска [4]. Шаговый алгоритм случайного поиска опробован на оптимизации параметров двучленной формулы суммарных потерь стока на горных водосборах, слагаемой из потерь на почвенное задержание и потерь в относительный подстилающий водоупор:
Р = (1 - с 1щ) Ртах їЬ (Хт / Ртах) + ко Хт . (9)
При этом оптимизировались одновременно три параметра расчетной формулы -Ртах (максимально возможные потери на почвенное задержание), с - коэффициент влияния предпаводочного увлажнения на потери стока, к'о - установившийся коэффициент потерь.
В качестве критерия оптимизации (минимизируемой функции) рекомендуется использовать среднее квадратическое отклонение рассчитанных значений потерь Р' от фактических Р для всех п паводков на водосборе:
F .
n -1
Ограничения на оптимизируемые параметры определены из природы формирования склонового стока и анализа графиков связи Р = f(XT, IW). На рис. 5 приведена блок-схема алгоритма оптимизации методом случайного поиска с переменным шагом в пространстве оптимизируемых параметров рассмотренной модели потерь.
Описание алгоритма. После ввода исходных данных, ограничений на параметры и ограничений на процесс минимизации формируется начальный случайный вектор оптимизируемых параметров с учетом диапазона параметров:
x,= *Г +dx,-x,,
5-0. г\ 1 с max min
где Xi - случайное число в интервале 0-1; ox, =x, - x, - диапазон оптимизации ,-го
параметра. Далее вычисляются фактические значения потерь Р = Xj—Y, модельные значения потерь Р по формуле (9), начальное значение критерия оптимизации F и вектор
Н начальных шагов параметров: h, = (x,max - x,mn)/2. Далее формируется новый вектор
параметров по формуле Xv+1 = XN + H х, где x — новый единичный случайный вектор, а XN — начальный случайный вектор параметров (на предыдущем шаге). Координаты вектора X проверяются на принадлежность их к заданному диапазону параметров, то есть учитываются ограничения. Если условие выполняется (на схеме — первый ромб), то шаг считается неудачным и выполняется адаптация (коррекция) шага Н. Адаптация связана с уменьшением шага по мере приближения к минимуму целевой функции и с его увеличением при отдалении. Алгоритмически это соответствует выражениям: HN+1 = kHN
при F > F и HN+1 =1HN при F <= F.
k
Ввод исходных данных (ХТ, Ут, /„), ограничений на параметры модели (х,), е, _М, ИМ
Генерирование единичного случайного вектора (х) и формирование начального вектора (X) оптимизируемых параметров: х,. = х™” +8х;
Вычисление фактических потерь стока (Р), модельных потерь (Р) начального критерия оптимизации (Р) начального рабочего шага параметров (вектора Н)
Формирование нового вектора х, вычисление приращения АХ = Н • х, нового вектора параметров
*-
«-
Вычисление потерь Р' и критерия
Счетчик удачных шагов _ / ^ р < е \Нет ^ Р = Р Н = Н / к
или ^ ► і > 1_М
Да
Печать: Р, Р, _, _1
Конец
Рис. 5. Блок-схема оптимизации параметров модели
Здесь Им - предыдущее значение вектора рабочего шага параметров; к - коэффициент, меньший 1 (принят равным 0,9); Г и Г - предыдущее и вновь вычисленное значение критерия оптимизации. В алгоритме применена нелинейная тактика случайного поиска [4], при которой случайный вектор параметров Х возвращается в прежнее состояние при неудачном шаге (Р > Г), а затем формируется новый единичный случайный вектор и вычисляется приращение вектора параметров АХ по скорректированному рабочему шагу. Выходом из цикла управляют блоки сравнения текущего числа шагов (Ь и Ь1) с заданным
143
предельным числом удачных или неудачных шагов (LM и L1M). Использованием тактики «возврата с пересчетом» и адаптацией рабочего шага в пространстве параметров удается добиться небольшого числа неудачных шагов, то есть малых затрат времени на оптимизацию. Анализ результатов оптимизации показал, что стабилизация пошаговых значений параметров модели потерь стока и стабилизация критерия оптимизации и наступает достаточно быстро (на 12-15 шаге).
Проверочные прогнозы суммарного слоя потерь дождевого стока по формуле (9) дали хорошие и удовлетворительные результаты. Параметры модели, полученные оптимизацией и графоаналитическим методом, оказались достаточно близкими величинами, что подтверждает целесообразность применения алгоритма случайного поиска с адаптацией шага к оптимизации параметров данной модели.
ЛИТЕРАТУРА
1. Бефани А.Н. Основные виды паводочного стока с горных водосборов и математические модели паводков / А.Н. Бефани // Метеорология, климатология и гидрология. Киев-Одесса, 1980. Вып. 16. С. 37-46.
2. Бефани Н.Ф. Прогнозирование дождевых паводков на основе территориальнообщих зависимостей / А.Н. Бефани. Л., 1977. 182 с.
3. Бефани Н.Ф. Упражнения и методические разработки по гидрологическим прогнозам / А.Н. Бефани, Г.П. Калинин. Л., 1965. 439 с.
4. Растригин Л. А. Случайный поиск в задачах оптимизации многопараметрических систем / Л. А. Растригин. Рига, 1965. 280 с.
Бондаренко Александр Михайлович — Bondarenko Aleksandr Mikhailovich —
кандидат географических наук, доцент, Candidate of Geographical Sciences,
заведующий кафедрой «Инженерные Associate Professor,
изыскания и информационные технологии Head of the Department
в строительстве» of «Engineering Survey and Information
Саратовского государственного Technology in Building»
технического университета of Saratov State Technical University
Статья поступила в редакцию 22.03.10, принята к опубликованию 30.11.10