Статистическая обработка рядов гидрологических характеристик (программы FreqShort, CorRows и ConCor)

Задачи, выдвигаемые проектной практикой перед статистической обработкой гидрологической информации, столь своеобразны, что существующие мощные пакеты программ математической статистики не содержат достаточных средств их решения. В известном смысле развитие некоторых направлений и прикладных методов статистики малых рядов ориентировалось именно на решение задач инженерной гидрологи. В других областях науки и практики эти методы применяются редко, обычно в менее разработанных вариантах и остаются мало известными.

Типовая задача применения гидрологической информации в практике проектирования сводится к расчету некоторых важнейших характеристик гидрологического режима рек (озер или др. водных объектов) заданной вероятности превышения. Так, при расчетах пропуска паводков в зависимости от класса сооружений требуется определить максимальные расходы (или уровни) воды с вероятностью их превышения 0.0 1, 0.005, 0.001 и даже 0.0001. При расчетах же водоснабжения, напротив, применяются характеристики минимального меженного стока вероятностью превышения 0.8, 0.9, 0.95 и 0.99. Для решения этих задач необходимо построить вероятностное распределение F(x) данной расчетной гидрологической характеристики x. Вероятности превышения Р=1-F(x) в гидрологии принято называть обеспеченностями и выражать в процентах, либо в продолжительности периодов возврата (10, 50, 100, 500, 1000 лет и т. д.).

Сложность и необычность статистической задачи в данном ее применении определяется несколькими условиями:

• Длина имеющихся вариационных рядов изучаемой характеристики лишь в редких случаях превышает 80-100 значений, обычно составляет 30-50, а зачастую бывает и менее 20. Во всех случаях она недостаточна для расчетов эмпирических частот в интервалах значений переменной. По этой причине для построения вероятностных распределений вычисляется эмпирическая обеспеченность каждого наблюденного значения ряда по одной из приближенных формул. Все эти формулы основаны на допущении Гаусса, согласно которому вероятности превышения значений случайной переменной могут определяться по ее порядковому номеру в ранжированном ряду без учета самого ее значения.
• Выходными результатами расчетов являются значения гидрологической характеристики с обеспеченностями, выходящими за пределы эмпирических обеспеченностей членов ряда при имеющейся продолжительности наблюдений. Возникающая задача экстраполяции эмпирического распределения (или его сглаживания и интерполяции в области крайних членов ряда) требует обоснованного применения некоторого теоретического распределения вероятностей как "математического лекала" и тщательного подбора его параметров.
• Вероятностные распределения большинства гидрологических характеристик по их физической природе существенно асимметричны. Так, максимальные расходы паводков заведомо неотрицательны, тогда как верхний их предел не установлен. Расчет коэффициентов асимметрии распределения при имеющихся рядах наблюдений сопряжен с большими случайными ошибками, которые особенно возрастают при наличии значимой внутрирядной корреляции. В других применениях при такой длине рядов математическая статистика обычно отказывается от учета асимметрии распределений. Однако ответственность задач гидрологических расчетов (высокая стоимость сооружений и требования их безопасной эксплуатации) заставляют искать пути наиболее полного использования имеющейся информации для подбора параметров теоретического распределения с учетом асимметрии, проявившейся за период наблюдений.
• При имеющейся длине рядов применение распределений с четырьмя и более параметрами практически невозможно. Из известных в математической статистике трехпараметрических распределений часто не находится такого, которое бы удовлетворительно соответствовало свойствам эмпирических распределений гидрологических характеристик. В тех случаях, когда переменная существенно положительна, но отдельные ее значения могут быть близки к нулю, известное из статистики (применяемое в гидрологии и реже в других областях) распределение Пирсона III типа становится неудовлетворительным. В этих случаях широко применяется распределение Крицкого-Менкеля, специально разработанное для гидрологических и водохозяйственных расчетов. В американской и английской инженерной гидрологии применяется также специальное распределение Гумбеля.

Применение обеспеченных значений гидрологических характеристик основывается на предположении о стационарности природного процесса, отражаемого многолетними рядами наблюдений. К настоящему времени в гидрологии все чаще отмечаются значимые изменения гидрологического режима объектов. Они могут выражаться в непрерывных трендах характеристик, проявляющихся с постоянной или переменной (чаще нарастающей) интенсивностью. Тренды такого рода связывают с нарастающими антропогенными изменениями ландшафтов бассейнов или климата. В других случаях отмечаются резкие разовые изменения гидрологического режима, как правило, под воздействием крупных водохозяйственных мероприятий. Программа должна быть снабжена средствами просмотра и анализа многолетнего хода исследуемой гидрологической характеристики. Необходима также оценка статистической значимости наблюдаемых проявлений неоднородности, т.к. субъективная оценка их ненадежна и часто усматривает значительные тренды там, где они могут быть результатом случайных сочетаний. Такой анализ многолетнего ряда данных позволяет при расчетах обеспеченных характеристик использовать ту его часть, которая наилучшим образом отражает современный гидрологический режим водного объекта.

Средства для решения всех изложенных выше задач сосредоточены в программе FreqShort.

В тех случаях, когда длина имеющегося ряда данных недостаточна для определения параметров распределения изучаемой характеристики, производят его продление с использованием достаточно тесной корреляционной связи данной характеристики расчетного бассейна и бассейна-аналога с более длительным периодом наблюдений.

В некоторых ситуациях могут применяться два или несколько аналогов. Методика продления рядов гидрологических характеристик реализована в программе Corrows, описанной далее отдельно. Здесь же отметим только, что любая, самая совершенная методика продления рядов наблюдений по связям с аналогами по ряду причин способна приводить к смещению оценок параметров распределения и статистических характеристик изучаемого элемента гидрологического режима.

Характеристики корреляции, полученные за период совместных наблюдений на расчетном объекте и объекте аналоге, могут не оставаться постоянными, но изменяться в широких пределах. Чем короче период совместных наблюдений, тем более вероятна неустойчивость полученных параметров связи. В то же время, устойчивость связей проявляется для разных объектов существенно различно. При выборе аналога желательно проверить устойчивость корреляционных связей его характеристик с данными других длиннорядных постов. Для выполнения такого анализа разработана программа ConCor. По всем изложенным причинам продление исходных рядов оправдано при заведомо недостаточной длине исходного ряда. В противном случае прибегать к нему следует с большой осторожностью.

Построение хронологических графиков, оценка однородности исходного ряда, расчет и построение кривых обеспеченностей.

Анализ временного хода исследуемой переменной и проверка однородности ряда должна предшествовать расчету обеспеченных характеристик. Хронологические графики переменной x_i могут быть построены в трех формах:

• собственно график многолетнего хода x_i = f(t);
• интегральная кривая - график нарастания суммы ежегодных значений переменной x_i = f(t);
• разностно-интегральная кривая - график нарастания суммы отклонений ежегодных значений переменной от среднего (x_i - x_ср) = f(t)

Любая неоднородность ряда данных наиболее явно проявляется в форме разностно-интегральных кривых. В период повышенных значений переменной наблюдается рост ординат, а период, когда значения ниже средних обозначается спадом кривой. Графики двух последних форм для обеспечения их сравнимости по разным объектам могут перестраиваться путем приведения к единичному среднему или к единичной дисперсии. Эти операции выполняются выбором соответствующих позиций дополнительного меню, вызываемого правой кнопкой мыши. На графике хронологического хода переменной может быть проведена линия тренда, рассчитываемая одним из двух способов: как линейная регрессии переменной на шкале лет или как прямая, соединяющая средние ординаты двух полупериодов.

Обнаружив на любом из этих графиков признаки неоднородности ряда и вызвав стрелку указателя, можно разделить период наблюдений на 2 или 3 части. Проверка гипотезы однородности выделенных частей ряда может быть выполнена в том случае, если для этого достаточна их длина. Проверка однородности ряда в отношении среднего значения переменной выполняется по критериям Стьюдента и Колмогорова-Смирнова. Однородность ряда в отношении дисперсии проверяется по критериям Колмогорова-Смирнова и Фишера. Для применения критерия Колмогорова-Смирнова необходимо найти максимум отклонения эмпирических распределений сравниваемых рядов по шкале обеспеченностей. Выполнение этого расчета отображается графически на экране монитора. При применении этого критерия к оценке однородности рядов по дисперсии исходные ряды приводятся к единичному среднему. Критерии Стьюдента и Фишера относятся к наиболее широко применяемым и не требуют пояснений. В выводных таблицах приводятся хронологические границы выделенных и сравниваемых периодов, длины рядов и число пропусков, основные статистические характеристики переменной, эмпирические значения примененных критериев, их критические значения для данного случая и заключение об однородности периодов по каждому из критериев. Во избежание недоразумений отметим, что данные за год, выбранный в качестве границы между двумя частями периода наблюдений, включаются в оба ряда. По умолчанию проверка гипотезы однородности производится при уровне значимости в 5%, однако расчет может быть повторен при уровне значимости в 1%, 2% или 10% по выбору пользователя. На заключительном этапе этого анализа может быть построен совмещенный график кривых обеспеченностей для двух сравниваемых частей ряда. Их интерполяция в данном случае производится кривой Пирсона III типа.

Построение кривых обеспеченностей P(x) в общем случае производится на специальной клетчатке вероятностей, применяемой в практике гидрологических расчетов. Горизонтальная шкала обеспеченностей P оцифрована в интервале от 0.5% до 0.95%. Для облегчения задачи подбора параметров интерполяционной кривой шкала P трансформирована пропорционально квантилям нормального распределения U(P), вычисляемым для обозначенных на шкале значений P. Обозначенной на шкале обеспеченности 50% соответствует значение U=0, поэтому значения P=0 и P=100% бесконечно удалены от осевой линии графика.

Эмпирическая кривая обеспеченностей обозначается точками с координатами (P_i, x_i), где P_i - эмпирическая обеспеченность, вычисляемая для каждого i-го значения переменной в зависимости от ее порядкового номера m в ранжированном ряду (в порядке убывания) и от общей длины ряда n. Вычисление обеспеченностей P_i производится по одной из формул приведенных в таблице 1.

Формула Крицкого-Менкеля создает некоторый запас надежности в оценке значений переменной малой расчетной обеспеченности (1%). Поэтому эта формула наиболее широко применяется в расчетах максимального стока (и наивысших уровней) при проектировании ответственных сооружений. Формула Гумбеля, напротив, значимо занижает характеристики малых расчетных обеспеченностей и завышает значения больших обеспеченностей (близких к 100%). Значения P_i, рассчитанные по медианной формуле Чегодаева занимают среднее положение между ними и чаще применяются при построении кривых обеспеченностей годового стока. Формула Алексеева дает значения близкие к формуле Чегодаева, но имеет некоторые дополнительные теоретические преимущества. Ее можно рекомендовать для применения при обобщениях характеристик гидрологического режима по множеству объектов и в теоретических расчетах. Формула Гаусса в гидрологии практически не применяется, т.к. обеспеченность минимального наблюденного значения переменной оценивается ею в 100%. По умолчанию в программе установлено применение формулы Алексеева, но опция основного меню "Параметры" позволяет выбрать для применения любую из этих формул. Другими изменяемыми с помощью этой опции параметрами программы являются длина ряда, минимальная для вычисления статистических характеристик (по умолчанию 7) и длина любой из сравниваемых частей ряда, минимальная для проверки его однородности (по умолчанию 5).

Для интерполяции и экстраполяции эмпирических распределений в программе предусмотрена возможность применения кривой Пирсона III типа, кривой Крицкого-Менкеля, кривой Гумбеля, нормального распределения, а также ломаной линии, соединяющей все точки эмпирического распределения. Для каждого анализируемого ряда первоначально производится построение кривой Пирсона III типа. В том случае, когда эта кривая неудовлетворительно соответствует эмпирической в области значений ниже среднего, может быть произведен переход к использованию кривой Крицкого-Менкеля, которая является трансформацией кривой Пирсона III типа, удовлетворяющей условиям P(0)=100% и P(x>0)<100%. Построение кривых этих двух типов производится с использованием в качестве параметров статистических характе6ристик ряда: его среднего x_ср., коэффициента вариации Cv=σ/x_ср. и коэффициента асимметрии Cs=μ₃/σ³, где μ₃ - третий центральный момент. Если кривая, построенная при вычисленных значениях этих характеристик, не вполне удовлетворительно соответствует эмпирическим точкам, можно попытаться подобрать кривую того же типа, но с другим параметром Cs. Для этого необходимо задавать некоторые значения отношения Cs/Cv. Исходное эмпирическое значение этого отношения и примененное на предыдущем этапе подбора, выводятся на экран.

В некоторых случаях крайние точки ранжированного ряда резко нарушают общий ход эмпирической кривой и далеко отклоняются от проведенной интерполяционной кривой. Если отклонение этих "отскакивающих" точек столь велико, что несомненно обозначает грубую ошибку измерения или физическую неоднородность порождающего его явления остальному ряду, программа предоставляет возможность исключать одну или несколько таких точек и производить расчеты заново без учета этих данных. Исключенные точки ряда с указанием лет, к которым они относятся, указываются в текстовом файле протокола расчетов. Если попытка исключения не дала желаемого результата или признана неправомерной, исключенные точки могут быть восстановлены.

При производстве гидрологических расчетов большое значение имеют характеристики точности расчета статистических параметров распределений изучаемого элемента - его среднего, Cv и Cs. Стандартные ошибки расчета этих параметров существенно зависят от длины ряда, коэффициентов внутрирядной корреляции, от полученных оценок изменчивости и асимметрии распределения. Незначительные отличия отмечаются также в зависимости от соответствия ряда тому или иному распределению. Наиболее точные оценки стандартных ошибок определения этих параметров к настоящему времени получены А. В. Рождественским [1], путем моделирования рядов с наборами названных параметров. Эти оценки выражены в форме таблиц, составленных по результатам моделирования. Для вычисления характеристик точности расчета параметров распределения изучаемого элемента в программе заложены названные таблицы и алгоритм криволинейной интерполяции их данных. В тех случаях, когда какой-либо из учитываемых аргументов расчета оказывается за пределами таблицы, вычисления производятся по эмпирической формуле, наиболее совершенной среди известных. Полный обзор таких формул также приводится в [1]. В последнем случае в соответствующем окне экрана монитора и в текстовом протоколе результатов расчета полученные характеристики точности приводятся в скобках. Учитывая значение этих данных и известное несовершенство изложенного способа их получения, разработчик программы в настоящее время работает над включением в нее процедуры стохастического моделирования рядов с единственным фактическим набором задаваемых параметров. Такой подход, по нашему мнению, позволит получить исчерпывающее решение данной задачи.

Кривая Гумбеля применяется в практике гидрологических расчетов в США, Англии и в ряде других стран. В зависимости от применения ее к рядам характеристик максимального или минимального стока она строится для верхней или нижней части ранжированного ряда. Кривая имеет два параметра и не учитывает асимметрии распределения. В данной программе она строится первоначально на обычной клетчатке P(x), но может быть перестроена в форме графика Гумбеля, применяемой в зарубежной практике. В последнем случае горизонтальная ось графика обозначается двумя параллельными шкалами. На верхней из них нанесены числа Гумбеля (Gumbel reduced variate), на нижней - периоды возврата Т (в годах) значений данной величины (Return period). Последние связаны с обеспеченностями расчетных максимальных значений соотношением T = 1/P, а с обеспеченностями минимальных значений соотношением T = 1/(1-P). Если по примененному кодовому обозначению изучаемого гидрометеорологического элемента известно, что он относится к числу характеристик минимального стока, кривая Гумбеля строится как понижающаяся в противном случае - как возрастающая. Возможность применения кривой Гумбеля в данной программе предусмотрена в связи с возрастающим числом совместных проектов и для целей международного обмена информацией.

Интерполяция ломаной линией, проходящей через все эмпирические точки графика применяется в тех случаях, когда форма эмпирической кривой исключает возможность ее сглаживания примененными простыми интерполяционными кривыми. Ее применение регламентировано действующими руководствами и СНиП-ами при построении кривых обеспеченностей характеристик уровня воды. Программа позволяет также ограничиться построением точек эмпирической кривой без их интерполяции.

В некоторых случаях, например, при анализе продленного ряда отличного от исходного ряда фактических данных, возникает необходимость записи анализируемого ряда в отдельный файл. Эта операция может быть выполнена выбором соответствующей опции дополнительного меню. При этом, кроме самого анализируемого ряда переменной по выбору пользователя могут быть записаны также ранжированный ряд, ряды рангов переменной, ее обеспеченностей или соответствующих им квантилей нормального распределения. Эти данные могут использоваться в различных менее распространенных методах статистического анализа, таких как ранговая или нормальная корреляция и других. Записываемые значения обеспеченностей могут быть получены по рангам переменной с использованием формул таблицы 1 или путем интерполяции с помощью примененной теоретической кривой.

Продление исходного ряда данных по связи с данными бассейнов аналогов.

Программа CorRows служит для продления рядов ежегодных гидрологических характеристик расчетного объекта по данным бассейнов-аналогов. Программа может принимать и использовать в расчетах данные до 9 аналогов. Если вводятся данные более чем одного аналога, пользователь может по своему усмотрению выполнить продление с использованием данных только одного из аналогов (по своему выбору), или применить вариант множественной корреляции. В последнем случае оптимальный набор аналогов для каждого года подбирается автоматически.

Во всех случаях решения принимаются с учетом показателей точности и надежности корреляционных связей: коэффициентов парной и множественной корреляции, стандартной ошибки их расчета, относительных ошибок расчета коэффициентов регрессии. Даже при наличии наибольшего возможного числа аналогов (девяти) случаи использования связи более чем с двумя аналогами встречаются крайне редко.

Решение по выбору способа продления ряда принимается пользователем на основании данных выводимой на экран матрицы коэффициентов парной корреляции и матрицы продолжительности совместных наблюдений. При использовании нескольких аналогов вместо обычного графика связи данных по двум объектам строится график связи наблюденных и вычисленных (по правилам применяемой методики) значений исследуемой характеристики для расчетного объекта.

При наличии нескольких аналогов с рядами наблюдений, достаточными для продления ряда данных расчетного объекта, и при достаточной корреляции, выбор единственного аналога может оказаться целесообразным в том случае, когда в результате применения программы ConCor обнаружено, что только этот аналог имеет устойчивые связи с длинными рядами данных по другим бассейнам.

Если на экран выведен график связи данных расчетного объекта и аналога, появляется возможность построить взаимную интегральную кривую как зависимость сумм членов обоих рядов, последовательно нарастающих в хронологическом порядке. На графиках связи данных разных постов и на взаимно-интегральных кривых пользователь может указать курсором мыши интересующую его точку графика. При нажатии кнопки мыши против этой точки появится надпись, обозначающая год, к которому относятся данные этой точки.

В текстовых файлах протокола счета отражаются результаты всех опробованных пользователем вариантов расчета. Продленный ряд записывается в файл в соответствии с условиями подготовки вводных данных для программы FreqShort. Сколько бы раз не выполнялись различные варианты продления исходного ряда, для последующего использования сохраняются только данные последнего расчета.

Программа ConCor, получая данные многолетних рядов ежегодных характеристик по нескольким постам, выполняет расчеты коэффициентов корреляции для любой указанной пары постов за скользящие интервалы продолжительностью 5, 7, 9, 11, 15 и 25 лет (по умолчанию) и строит совмещенный график хронологического хода полученных значений коэффициентов корреляции. Каждое из вычисленных значений при данном положении шаблона относится к середине примененного интервала. Длина скользящих интервалов может изменяться пользователем в широких пределах. Проявляющиеся в результате расчетов эффекты высокой устойчивости коэффициентов корреляции, определяемых даже за наиболее короткие интервалы времени, и, напротив, резких изменений уровня корреляции в пределах от -0.8 до 0.95 не имеют в настоящее время убедительного объяснения, но они, несомненно, должны учитываться при выборе аналогов для продления коротких рядов гидрологических характеристик.

Пример расчета хронологического хода коэффициентов корреляции максимальных расходов весеннего половодья для двух постов Верхне-Волжского УГМС по скользящим интервалам времени дан на рис.1.

Рис. 1. Типичный пример графика временного хода коэффициентов корреляции по скользящим коротким интервалам времени,
построенного при помощи программы ConCor.

Вопросы подготовки данных для описанных программ и организация вывода результатов расчета подробно изложены в их технических описаниях и здесь не рассматриваются. Все графики, создаваемые программами, в любой момент могут быть сохранены в файлах формата BMP.