Определение параметров гидрологической модели Текст научной статьи по специальности «Математика»

Определение параметров гидрологической

модели

Е. С. Засухина, С. В. Засухин

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

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

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

I. Введение

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

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

Статья получена 20 декабря 2016 г.

Работа выполнена при поддержке Российского фонда фундаментальных исследований № 15-07-08952.

Е. С. Засухина, научный сотрудник Вычислительного центра им. А. А. Дородницына Федерального исследовательского центра «Информатика и управление» Российской академии наук. [email protected].

С. В. Засухин, аспирант Московского физико-технического института (государственного университета). [email protected].

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

При моделировании процессов формирования стока на водосборе важная роль отводится моделям передвижения воды в почве. Присутствующие в эти моделях гидрофизические характеристики почвы — коэффициент диффузии и гидравлическая проводимость, как правило, рассчитываемые по формулам ван Генухтена [1], содержат трудно определяемые параметры. Задание точных значений этих параметров имеет критическое значение при моделировании и прогнозировании потока воды и переноса растворенных веществ в зоне аэрации. Этому вопросу посвящено огромное число работ. В некоторых из них при нахождении параметров модели используются методы типа градиентных [2-6]. Разработаны компьютерные программы для определения гидрологических параметров. Здесь следует отметить компьютерную программу КЕТС [6], позволяющую, в том числе, определять эти параметры по измеренной функции водоудерживания, а также программу Rosetta [7], которая позволяет, в частности, находить гидрологические параметры с помощью педотрансферных функций, получаемых методом нейронных сетей. В ряде работ для нахождения гидрофизических параметров в процессе численной оптимизации был использован алгоритм имитации отжига, например в [8]. В работах [9, 10] искомые параметры находились с применением генетического алгоритма. В [11] применяется алгоритм глобальной оптимизации с покоординатным дроблением области определения целевой функции. На протяжении последних десятилетий появилось много работ, в которых для поиска искомых гидрофизических параметров применяются оптимизационные алгоритмы, имитирующие поведение биологических популяций в условиях недостатка жизненных ресурсов и мигрирующих с целью найти место с благоприятными условиями проживания, алгоритмы, имитирующие социальное поведение [12-13]. В ряде работ при нахождении искомых гидрофизических параметров применяют стохастические методы оптимизации [14-17].

В предлагаемой статье рассматривается задача

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

II. Постановка задачи

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

Рассмотрим следующую начально-краевую задачу:

(О е Q,

f ( - ' И

0(z,O) = cp(z ),

в(Ь, t ) = w(t ),

_Г D(0)^ - K (0)

z=0

= R(t)-E (t*

z e(0, L), t e(0,T ),

t e(0, T ),

(1)

K (0) = K0S 0 D(0) = K0

-il - S Vm

1 - m

S

0.5-7 m

(2)

am(0s -0mm )

+ 11 -s'm| -2

1/ \-m / 1/ \m

1 -s /m) +11 - s lm

где ^ = (в-вттУ(в, -втт), и К0, а, т , -

некоторые параметры, причем 0 < в, - втах << в,, Назовем описанную задачу (1) прямой задачей.

Сформулируем задачу идентификации параметров а и т . Предположим, что на некотором множестве

Qo с Q задана функция в(х,(). Назовем эту функцию "экспериментальными данными." Поставим задачу подобрать параметры а и т таким образом, чтобы соответствующее решение прямой задачи было как можно ближе к функции в(р,() на множестве Qo с Q . Или, более точно, найти а°> и т°> и соответствующее решение в0°>> (г, t) прямой задачи (1)

такие, чтобы функционал J = — |(в°> -в}

2 Qo

достигал минимума.

Перейдем к дискретному аналогу задачи (1). Разобьем интервалы (0,ь) и (0,Т) на I и N равных подынтервалов с концевыми точками = Ы , 0 < 1 < I,

и ^ = тп, 0 < п < N, соответственно, где т = Т^ и Ы = Ь/1. Аппроксимируем прямую задачу (1) с помощью следующей конечно-разностной схемы:

0n+1 -0n

D_

0n+1 - 0n+1 }

n+1 °i+1 - Kn+1

i+1/2 h Ki+1/2

0П+1 - 0П+1 ^

-.n+1 i-1 n+1

Л-1/2 h V1/2

втШ <в(0,t)<втах, t е (0,Т),

где 2 - пространственная координата по оси,

направленной вниз от поверхности почвы; t - время;

в(х, t) - искомая влажность в точке (г,t), так называемая «объемная» влажность почвы, выражаемая в единицах объема воды в единичном объеме почвы (безразмерная величина); Р = (0, ь)х(0,Т); ф(т.) и ц/({) - заданные функции; в(в) и К (в) - коэффициент диффузии и гидравлическая проводимость -гидрофизические характеристики почвы; я() и Е() -интенсивности осадков и испарения - линейные потоки влаги; 0 < Е(t) < О, t е (0, Т), О — некоторая константа.

Входящие в уравнение коэффициент диффузии в(в) и гидравлическая проводимость К (в) вычисляются по широко применяемым формулам ван Генухтена [1]

1 < 1 < 1, 0 < п < N,

в0 =щ, 1 < 1 < I, в'п =уп, 1 < п < N,

где вв , ^1+12, Кгп_12 есть значения функций в(х, t),

в(в^,, К(в(г, t)) в точках ^, ^ ), (( + 1/2)ы,тп),

(( - 1/2)ы, тп) соответственно. Краевое условие слева имеет вид:

nn+1 nn 0 0о - 0о _ 2

D

,n+1 - 0n+1 n+1 01 0о - kn+1 + rn+1 - ßn+1

h 12

1/2

0 < n < N,

где Яп+1, Еп+1 есть значения функций Я^) и Е^) в

точках = т(п +1).

Дискретный аналог прямой задачи выглядит следующий образом:

ф п =-| 1+А в" 1в0п +Л-в" + 1в,п-1 +

г h2 D1/2 DV2°1

21

K" + Rn - En)= 0, h V V2 > '

0mm <0q" <0max, 1 < n < N,

Ф" =-L D^/ 0P-1 -

h

o«-1

2Di-1/2 i-1 "

1 + d;

Г h

2 Гг+1/2 т г-12^

0in d"„0+1 +

(3)

h

2 Di+1/2"i+1 "

1 (ki-

Г h Ki-V2 Ki+12

= 0, 1 < i < 1, 1 < n < N,

1 < n < N,

ф п =^п-вч = 0,

вв0 0 < 1 <I.

При этом фигурирующие в формулах (3) коэффициент диффузии в и гидравлическую проводимость К в промежуточных точках будем вычислять по следующим формулам:

г

h

г

г

2

m

1

+

X

n—\. пИ-1

D

Di + D

i+1

n-1 . T^n-1

K

K"-L + K

i+1

,+l/2 - 2 > Ki+1/2 - 3 ' (4)

1 < n < N, 0 < i < I. Положим Qo - {(z, t): z - ih, t - It, (i, l) e a} , где

A - {(i, n) : i - kj, n - ql, j - 0,..., [l/к], l - 0,..., [N/q]}, к > 1, q > 1 — некоторые фиксированные натуральные числа. Зададим целевой функционал в виде

W(e,u)-2 Yten4nfhT. (5)

2 (i,n)eA

Дискретная задача оптимального управления формулируется следующим образом: найти оптимальное

управление uopt -{aopt,mopt} и соответствующее оптимальное решение задачи (3) такие, что функционал W (u) (1.5) достигал бы минимального значения.

Поиск численного решения полученной задачи конечномерной оптимизации предлагается проводить методом наискорейшего спуска, при этом градиент будет вычисляться с применением метода быстрого автоматического дифференцирования (БАД) [18-21]. Изложим кратко суть этого метода. Приведем содержащийся в [20] способ получения формул БАД для вычисления производных сложной функции, основанный на теореме о неявной функции.

III. Формулы БАД

Предположим, что для векторов z e Rn и u e Rr дифференцируемые функции W (z, u) и Ф^, u)

определяют отображения W : Rn x Rr ^ R1 и

Ф: Rn x Rr ^ Rn . Пусть z и u удовлетворяют системе из n скалярных алгебраических уравнений

Ф^, u)- 0n , (6)

где 0n есть нулевой n -мерный вектор. Предположим

также, что матрица Ф^ (z, u) неособенная всюду в интересующей нас области. Тогда по теореме о неявной функции система связей (6) определяет непрерывно-дифференцируемую функцию z - z(u). Согласно методу БАД градиент функции W(z(u), u) вычисляется по формуле

dW(z(u),uVdu - Wu (z(u), u) + Ф^ (z(u),u)p. (7)

Входящий в эту формулу вектор p e Rn является множителем Лагранжа и определяется в результате решения системы линейных уравнений

Wz(z(u),u)+Ф^(z(u),u)p - 0n . (8)

Линейная система (8) является сопряженной к исходной системе связей (6).

IV. Численные результаты

Полученные формулы были применены при нахождении численного решения описанной выше дискретной задачи оптимального управления градиентным методом. Использовалась сетка с параметрами: I - 100, N - 96.

Расчеты производились со следующими значениями входных параметров:

L - 100 (см), T - 1 (сут), K0 - 100 (см/сут),

0mln - 0.05 (см3/см3) 0тах - 0.5 (см3/см3)

<p(z)- 0.3, z e(0,L), y/(t)- 0.3, t e (0,T).

Расчеты проводились следующим образом. Сначала

решалась прямая задача (3) со значениями параметров

агие - 0.01 и mtrue - 0.2. Из вида уравнений (3) и с учетом формул для вычисления коэффициентов в промежуточных точках (4) следует, что система уравнений (3) расщепляется на N подсистем, каждая из которых соответствует определенному временному слою, и может решаться отдельно от других. Каждая система решается методом прогонки, так как ее основная матрица является трехдиагональной. Полученное решение d(z, t) задачи (3) принималось в качестве предписанной функции d(z, t).

Далее было проведено несколько серий численных экспериментов. Каждая серия численных экспериментов состояла из решения описанной дискретной задачи оптимального управления с целевой функцией, вычисляемой по формуле (5), при этом множество A имело вид: A - {(i,n): i - kj, j - 0,.,[l/к],n - 1,.,d}, где d - 1,2,3,...,26 , а к оставалось постоянным. Было проведено четыре серии таких численных экспериментов с к - 1,2,5,10. В качестве начального приближения были приняты следующие значения

параметров: aimt - 0.03 и mimt - 0.11. Поиск численного решения каждой задачи производился методом наискорейшего спуска с применением точного градиента, вычисляемого по формулам БАД (7)-(8). Величина шага вдоль выбранного направления определялась в результате выполнения процедуры одномерной минимизации вдоль этого направления функции, интерполирующей целевую функцию с помощью сплайнов по 40 точкам. Процесс численной оптимизации прекращался, когда чебышевская норма

градиента становилась меньше 10 7 .

Результаты численных расчетов представлены в четырех таблицах. Каждая таблица содержит численные результаты отдельной серии экспериментов: таблица I содержит результаты расчетов серии экспериментов, в которой к - 1; таблица II содержит результаты серии экспериментов с к - 2; таблица III содержит численные результаты серии с к - 5; таблица IV содержит результаты численных расчетов серии экспериментов с к -10.

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

возможности позволяют проводить указанные измерения влажности почвы через каждые 15 минут и с шагом в 1, 2, 5 и 10 см по глубине. Выделим четыре группы измерений, в каждой из которых измерения влажности почвы проводятся с постоянным шагом по глубине в 1, 2, 5 и 10 см и через каждые 15 минут. В каждой группе произведем по 26 измерений профилей влажности, следующих одно за другим через каждые 15 минут, и рассмотрим 26 наборов данных, состоящих из 1, 2, ..., 26 последовательно измеренных профилей влажности. И теперь рассмотрим каждый из этих наборов данных в качестве предписанных значений, среднеквадратическое отклонение которых от полученных в результате решения соответствующей прямой задачи при искомых параметрах должно быть минимальным. Зададимся вопросом, какое количество профилей влажности необходимо измерить для того, чтобы в результате решения соответствующей задачи оптимального управления найти искомые параметры с заданной точностью? В предположении о том, что ошибки измерения пренебрежимо малы, можно сказать, что анализ проведенных численных экспериментов позволяет ответить на этот вопрос.

Таблица I содержит результаты решения дискретной задачи оптимального управления с целевой функцией (5), в которой множество А имеет вид: А = {(/',п): I = к],] = 0,...,[¡/к],п = 1,...,ё}, где к = 1, ё = 1,2,3,. ,26, т. е. вычисленные и предписанные значения влажности почвы сравниваются в точках, следующих друг за другом через 1 см.

opt л „2 а г ■ 10

1.16155 1.01417 1.00772 1.00593 1.00517 1.00476 1.00451 1.00434 1.00421

1.00412

1.00405

1.00399

1.00394

1.00390

1.00387

m0pt ■ 10

Таблица I

Число итераций

1.95558 1.99477 1.99689 1.99750 1.99776 1.99790 1.99799 1.99805 1.99809

1.99812

1.99815

1.99817

1.99819

1.99820

1.99821

47813

18320

12151

10056

9093

8564

8237

8017

7860

7742

7651

7578

7519

7470

7428

Значение функции

2.7 1.6 1.3 1.2 1.1 1.1 1.0 1.0

9.9 ■ 10

9.8 ■ 10 9.7 ■ 10 9.6 ■ 10 9.5 ■ 10 9.5 ■ 10

10- -10

10- 11

10- 11

10- 11

10- 11

10- 11

10- 11

10- 11

10- 11

-12

12

12

12

12

12

1.00384 1.00381 1.00379 1.00377 1.00375 1.00406 1.00430 1.00424 1.00417

1.00411

1.00412

1.99822

1.99823

1.99824

1.99824

1.99825 1.99811 1.99799 1.99802 1.99805 1.99808 1.99807

7393 7362 7336 7312 7291 7185 7198 7349 7491 7655 7865

9.4 ■ 10

9.4 ■ 10

-12

9.3 ■ 10

-12

9.3 ■ 10

-12

9.3 ■ 10

-12

1.1 ■ 10

-1

1.2 ■ 10

-1

1.2 ■ 10

-1

1.2 ■ 10

-1

1.1 ■ 10

-1

1.1 ■ 10

-1

Как видно увеличивается

из таблицы I с увеличением n точность нахождения искомых параметров: точность в определении а изменяется от 16% до 0.4% при изменении n от 1 до 26, а точность в определении m — от 2.3% до 0.1%. При этом уменьшается количество итераций, за которое процесс численной оптимизации приводит к решению.

opt opt „ true

Отклонение а г и m г от истинных значений а

true ,

и m резко уменьшается при переходе от n = 1 к

n = 2. Далее с увеличением n от 2 до 10 происходит

дальнейшее уменьшение отклонений вычисленных

значений параметров а и m от их истинных значений.

С последующим увеличением n от 10 до 26 эти

отклонения уменьшаются очень незначительно. Таким

образом, в свете описанной интерпретации при

составлении плана измерений профилей влажности

через каждые 15 минут с шагом по глубине в 1 см для

получения экспериментальных данных, которые будут

использованы в рамках предлагаемого подхода для

нахождения а с точностью 0.4% и m с точностью

0.1%, можно было ограничиться измерением 8-10 таких

профилей.

Таблица II содержит результаты численных расчетов дискретной задачи оптимального управления с целевой функцией (5), в которой множество A имеет вид: A = {(i,n): i = kj, j = 0,.,[l/k],n = 1,.,d}, где к = 2, d = 1,2,3,. ,26 , т. е. сравнение полученных в результате решения прямой задачи (3) значений влажности почвы с предписанными значениями происходит в точках, следующих друг за другом через каждые 2 см.

Таблица II

opt 2 а г ■ 10

m0pt ■ 10

Число итераций

Значение функции

n

n

1.27157 1.02775 1.01516 1.01167 1.01017 1.00936 1.00887 1.00853 1.00829

1.00811

1.00797

1.00785

1.00776

1.00768

1.00761

1.00755

1.00751

1.00746

1.00742

1.00739

1.00736

1.00733

1.00730

1.00728

1.00726

1.00724

1.93120 1.98988 1.99393 1.99510 1.99561 1.99589 1.99606 1.99617 1.99626

1.99632

1.99637

1.99641

1.99644

1.99647

1.99649

1.99651

1.99653

1.99654

1.99655

1.99657

1.99658

1.99659

1.99660

1.99660

1.99661

1.99662

58793 29503 19986 16621 15053 14189 13654 13295 13039

12847

12699

12582

12486

12406

12340

12283

12233

12191

12153

12120

12090

12063

12039

12017

11997

11978

3.3 5.2

3.2 2.5

2.3 2.2 2.1 2.0 2.0

10 1010101010101010-

1.9 • 10 1.9 • 101.9 • 101.9 • 101.9 • 101.9 • 101.8 • 101.8 • 101.8 • 101.8 • 101.8 • 101.8 • 101.8 • 101.8 • 101.8 • 101.8 • 101.8 • 10-

true л

и m резко уменьшается при переходе от n = 1 к

n = 2. Далее с увеличением n от 2 до 15 происходит

дальнейшее уменьшение отклонений вычисленных

значений параметров а и m от их истинных значений.

С последующим увеличением n от 16 до 26 эта разница

снижается гораздо медленнее,. Таким образом, при

составлении плана измерений профилей влажности

через каждые 15 минут с шагом по глубине в 2 см для

получения экспериментальных данных, которые будут

использованы в рамках предлагаемого подхода для

нахождения а с точностью 0.8% и m с точностью

0.2%, можно было ограничиться измерением 10-15

таких профилей.

В таблице III приводятся результаты серии численных экспериментов, когда решалась дискретная задача оптимального управления с целевой функцией (5), в которой множество A имеет вид: A = {(i,n): i = kj, j = 0,...,[//к],n = 1,...,d}, где к = 5, d = 1,2,3,. ,26 , т. е. сравнение полученных в результате решения прямой задачи (3) значений влажности почвы с предписанными значениями происходит в точках, следующих друг за другом через каждые 5 см.

Таблица III

Результаты, содержащиеся в таблице II, аналогичны результатам из таблицы I. С увеличением n уменьшается разница между найденными значениями параметров и их истинными значениями, а именно: разница между найденными значениями а и его истинным значением изменяется от 27% до 0.7% при изменении n от 1 до 26, а в отношении m эта разница изменяется от 3.5% до 0.2%. С увеличением n количество итераций, за которое процесс численной оптимизации приводит к решению, уменьшается.

opt opt „ true

Отклонение а г и m г от истинных значений а

opt • 102 m opt • 10 Число итераций Значение функции

1.50588 1 .88993 72516 5.3 ю-10

1.06549 1 .97679 51740 1.2 10-10

1.03599 1 98581 36323 7.3 10-11

1.02775 1 .98849 30424 6.0 10-11

1.02420 1 .98966 27610 5.4 10-11

1.02230 1 .99029 26046 5.1 10-11

1.02113 1 .99069 25078 4.9 10-11

1.02034 1 .99095 24429 4.7 10-11

1.01977 1 .99114 23968 4.6 10-11

1.01934 1 99129 23625 4.6 10-11

1.01900 1 .99140 23361 4.5 10-11

1.01873 1 .99149 23152 4.5 10-11

1.01851 1 .99157 22982 4.4 10-11

1.01832 1 .99163 22842 4.4 10-11

1.01816 1 .99168 22725 4.4 10-11

1.01803 1 .99173 22625 4.3 10-11

1.01791 1 99177 22539 4.3 10-11

1.01781 1 .99180 22464 4.3 10-11

0

n

1.01771 1.01763 1.01756 1.01749 1.01744 1.01738 1.01733 1.01729

1.99184 1.99186 1.99189 1.99191 1.99193 1.99195

1.99197

1.99198

22399 22341 22289 22243 22201 22163 22129 22097

4.3 • 10" 4.3 • 104.2 • 104.2 • 104.2 • 104.2 • 104.2 • 104.2 • 10-

Как видно из таблицы III, результаты численных расчетов, приведенные в ней, подчиняются тем же закономерностям, что и результаты, содержащиеся в таблицах I и II. С увеличением n увеличивается точность нахождения искомых параметров: точность в определении а изменяется от 50,6% до 1.7% при изменении n от 1 до 26, а точность в определении m — от 5.5% до 0.4%. При этом уменьшается количество итераций, за которое процесс численной оптимизации

opt

приводит к решению. Отклонение

opt m г от

true true

истинных значений а и m резко уменьшается при переходе от n = 1 к n = 2. Далее с увеличением n от 2 до 12 происходит дальнейшее уменьшение отклонений вычисленных значений параметров а и m от их истинных значений, при этом это уменьшение все более и более замедляется. С последующим увеличением n от 12 до 26 такие отклонения практически остаются неизменными. Таким образом, при составлении плана измерений профилей влажности через каждые 15 минут с шагом по глубине в 5 см для получения экспериментальных данных, которые в дальнейшем будут использованы в рамках предлагаемого подхода для нахождения а с точностью 2% и m с точностью 0.5%, можно было ограничиться измерением 10-12 таких профилей.

Таблица IV содержит результаты численных расчетов дискретной задачи оптимального управления с целевой функцией (5), в которой множество A имеет вид: A = {(i,n): i = kj, j = 0,.,[¡/к],n = 1,.,d}, где к = 10, d = 1,2,3,. ,26 , т. е. сравнение полученных в результате решения прямой задачи (3) значений влажности почвы с предписанными значениями происходит в точках, следующих друг за другом через каждые 10 см.

Таблица IV

n a0pt • 102 m0pt • 10 Число итераций Значение функции

1 1.79506 1.85193 79097 -10 7.8 • 10 10

2 1.12074 1.95891 74318 -10 2.1 • 10 10

06660 1.97431 53916 1.3 10

05144 1.97901 45382 1.1 10

04492 1.98109 41202 9.7 10

04142 1.98221 38859 9.2 10

03926 1.98291 37407 8.8 10

03780 1.98339 36436 8.6 10

03675 1.98373 35749 8.4 10

03596 1.98399 35241 8.3 10

03534 1.98419 34851 8.2 10

03484 1.98435 34544 8.1 10

03444 1.98448 34296 8.0 10

03410 1.98459 34092 8.0 10

03381 1.98469 33922 7.9 10

03356 1.98477 33777 7.9 10

03334 1.98484 33654 7.9 10

03315 1.98490 33546 7.8 10

03299 1.98496 33452 7.8 10

03283 1.98501 33370 7.8 10

03270 1.98505 33296 7.8 10

03258 1.98509 33230 7.7 10

03247 1.98513 33171 7.7 10

03237 1.98516 33118 7.7 10

03228 1.98519 33069 7.7 10

03220 1.98522 33025 7.7 10

Анализ

результатов

численных

расчетов,

приведенных в таблице IV, приводит к выводам, аналогичным полученным в результате анализа таблиц I-III. С увеличением n увеличивается точность нахождения искомых параметров: точность в определении а изменяется от 79,5% до 3% при изменении n от 1 до 26, а точность в определении m — от 15% до 0.8%. При этом уменьшается количество итераций, за которое процесс численной оптимизации

opt opt

приводит к решению. Отклонение а г и m г от

true true

истинных значений а и m резко уменьшается при переходе от n = 1 к n = 2 . Далее с увеличением n от 2 до 15 происходит дальнейшее уменьшение отклонений вычисленных значений параметров а и m от их истинных значений, при этом это уменьшение все

более и более замедляется. С последующим увеличением п от 16 до 26 такие отклонения практически остаются неизменными. Таким образом, при составлении плана измерений профилей влажности через каждые 15 минут с шагом по глубине в 10 см для получения экспериментальных данных, которые в дальнейшем будут использованы в рамках предлагаемого подхода для нахождения а с точностью 3% и т с точностью 0.8%, можно было бы ограничиться измерением 12-15 таких профилей.

Таким образом, мы можем ограничиться измерением 10-15 профилей влажности почвы и в зависимости от шага по глубине при проведении измерений рассчитывать на различную точность нахождения параметров а и т в результате применения описанного подхода. Так, при шаге в 1 см эта точность составляет 0.4% для а и 0.1% для т , при шаге в 2 см — 0.8% для а и 0.2% для т , при шаге в 5 см — 2% для а и 0.5% для т , при шаге в 10 см — 3% для а и 0.8% для т .

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

V. Заключение

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

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

Следует отметить недостаток описанного подхода: применяемый в процессе численной оптимизации метод градиентного спуска относится к числу локальных методов. В связи с этим возникает вопрос о возможности применения в данной ситуации метода глобальной оптимизации, например, известного метода неравномерных покрытий [24].

Библиография

[1] M. Th. van Genuchten, "A closed form equation for predicting the hydraulic conductivity of unsaturated soils," Soil. Sci. Soc. Am. J. vol. 44, pp. 892-898, Sept. 1980.

[2] S. O. Eching and J.W. Hopman, "Optimization of hydraulic functions from transient outflow and soil water pressure head data," Soil. Sci. Soc. Am. J. vol. 57, pp. 1167-1175, Sept. 1993.

[3] J. B. Kool and J. C. Parker, "Analysis of the inverse problem for transient unsaturated flow," Water Resour. Res. vol. 24(6), pp. 817830, Jun. 1999.

[4] N. Romano and A. Santini, "Determining soil hydraulic functions from evaporation experiments by a parameter estimation approach: experimental verifications and numerical studies," Water Resour. Res. vol. 35(11), pp. 3342-3359, Nov. 1999.

[5] J. Simunek and M. Th. van Genuchten, "Estimating unsaturated soil parameters from multiple tension disc infiltrometer data," Soil Sci. vol. 162 (6), pp. 383-398, Jun. 1997.

[6] M. Th. van Genuchten, F. J. Leij and S. R. Yates, "The RETC code for quantifying the hydraulic functions of unsaturated soils," U.S. Salinity Lab., USDA, ARS, Riverside, California. EPA Report 600/291/065, 1991.

[7] M. G.Scaap., F .J. Leij and M. Th. van Genuchten, "Rosetta: a computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions," J. Hydrol. vol. 251, pp. 163176, Oct. 2001.

[8] L. H.. Pan and L. S. Wu, "A hybrid global optimization method for inverse estimation of hydraulic parameters: annealing-simplex method," Water Resour. Res. vol. 34, pp. 2261-2269, Nov. 1998.

[9] Y. Takeshita, "Parameter estimation of unsaturated soil hydraulic properties from transient outfl ow experiments using genetic algorithms," in Proc. Int. Worksh.,Riverside, CA. 22-24 Oct. 1997, U.S. Salinity Lab., Riverside, CA, New York, 1999, pp. 761-768.

[10] J. A. Vrugt, A. H. Weerts, and W. Bouten, "Information content of data for identifying soil hydraulic parameters from outfl ow experiments," Soil. Sci. Soc. Am. J. vol. 65, pp. 19-27, Jan. 2001.

[11] S. Lambot, M. Javaux, F. Hupet, and M. Vanclooster, "A global multilevel coordinate search procedure for estimating the unsaturated soil hydraulic properties," Water Resour. Res. vol. 38(11), pp. 1224?, Nov. 2002.

[12] K. C. Abbaspour, R. Schulin, and M.Th . van Genuchten, "Estimating unsaturated soil hydraulic parameters using ant colony optimization," Adv. Water Resour. vol. 24, pp. 827-841, Aug. 2001.

[13] X. Yang and X. You, "Estimating parameters of van Genuchten model for soil water retention curve by intelligent algorithms," Appl. Math. Inf. Sci. vol. 7(5), pp. 1977- 1983, Sept. 2013.

[14] J. A. Vrugt and W. Bouten, "Validity of first-order approximations to describe parameter uncertainty in soil hydrologic models," Soil. Sci. Soc. Am. J. vol. 66(6), pp. 1740-1752, Nov. 2002.

[15] J. A. Vrugt, H. V. Gupta, W. Bouten, L. A. Bastidas, and S. Sorooshian, "Effective and efficient algorithm for multi-objective optimization of hydrologic models," Water Resour. Res. vol. 39(8), pp. 1214-1232, Aug. 2003.

[16] J. Mertens, H. Madsen, M. Kristensen, D. Jacques, and J. Feyen, "Sensitivity of soil parameters in unsaturated zone modeling and the relation between effective, laboratory and in situ estimates," Hydrol. Processes. vol. 19(8), pp. 1611-1633, May. 2005.

[17] J. Mertens, R. Stenger, and G. F. Barkle, "Multiobjective inverse modeling for soil parameter estimation and model verification," Vadose Zone J. vol. 5(3), pp. 917-933, Aug. 2006.

[18] К. Р. 2. Айда-Заде, Ю. Г. Евтушенко, "Быстрое автоматическое дифференцирование на ЭВМ," Математическое моделирование. Т. 1, № 1, С. 121-139, 1989.

[19] Automatic Differentiation of Algorithms. Theory, Implementation and Application, A. Griewank and G. F. Corliss, Ed. Philadelphia: SIAM, 1991.

[20] Yu. Evtushenko, "Automatic differentiation viewed from optimal control theory," in Automatic Differentiation of Algorithms. Theory, Implementation and Application, A. Griewank and

G. F. Corliss, Ed. Philadelphia: SIAM, 1991, pp. 25-30.

[21] Yu. Evtushenko, "Computation of exact gradients in distributed dynamic systems," Optimization methods and software. vol. 9, pp. 45-75, Sept. 1998.

[22] С. А. Лупин, М. А. Посыпкин, Технологии параллельного программирования. М: Форум Инфра-М, 2008. 208 с.

[23] Y. Evtushenko, M. Posypkin and I. Sigal, "A framework for parallel large-scale global optimization," Computer Science - Research and Development. vol. 23, № 3, pp. 211-215, Jun. 2009.

[24] Y. Evtushenko and M. Posypkin, "A deterministic approach to global box-constrained optimization," Optimization Letters. vol. 7, № 4, pp. 819-829, Apr. 2013.

Determining parameters of hydrological model

E. S. Zasukhina, S. V. Zasukhin

Abstract—The study considers a problem of determining parameters of hydro-physical characteristics included in many hydrological models of runoff formation in the catchment area. The problem of parameters determination is formulated as an optimal control problem. The object function is standard deviation of modeled soil moisture profiles from some prescribed values, and the control is defined parameters. The discretized optimal control problem is solved numerically by gradient method and exact values of gradient of the objective function are calculated by fast automatic differentiation method.

Keywords—optimization, optimal control, gradient method, fast automatic differentiation.