автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Численное решение задач идентификации коэффициента фильтрации на основе двухшаговых методов минимизации функции невязки

003494588

На правах рукописи

Кадырова Альфия Шамилевна

ЧИСЛЕННОЕ РЕШЕНИЕ ЗАДАЧ ИДЕНТИФИКАЦИИ КОЭФФИЦИЕНТА ФИЛЬТРАЦИИ НА ОСНОВЕ ДВУХШАГОВЫХ МЕТОДОВ МИНИМИЗАЦИИ ФУНКЦИИ НЕВЯЗКИ

05.13.18 - математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРАТ

диссертации на соискание ученой степени кандидата физико-математических наук

Казань-2010

2 5 МА? 29:9

003494588

Работа выполнена в Учреждении Российской академии наук Институте механики и машиностроения Казанского научного центра РАН

Научный руководитель:

кандидат физико-математических наук, старший научный сотрудник Мазуров Петр Алексеевич

Официальные оппоненты:

доктор физико-математических наук, профессор Лапин Александр Васильевич

доктор физико-математических наук, Шешуков Евгений Геннадьевич

Ведущая организация:

Учреждение Российской академии наук Институт прикладной математики им. М.В. Келдыша РАН (г. Москва)

Защита состоится «08» апреля 2010 г. в 14 часов 30 минут на заседании диссертационного совета Д 212.081.21 в Казанском государственном университете по адресу: 420008, Казань, ул. Кремлёвская, 18, корп. 2, ауд. 218.

С диссертацией можно ознакомиться в Научной библиотеке им. Н.И. Лобачевского Казанского государственного университета.

Автореферат разослан «05» марта 2010 г.

Ученый секретарь диссертационного совета

Д 212.081.21 д.ф.-м.н.

Задворнов О. А.

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

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

Цели диссертационной работы:

- построение эффективных методов минимизации функции невязки;

- разработка на основе предложенных методов численных алгоритмов решения задач идентификации коэффициента фильтрации водоносных пластов.

Научная новизна результатов. 1

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

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

Достоверность полученных результатов обеспечивается корректным применением методов вычислительной математики и математических моделей механики сплошной среды, сравнением полученных решений с известными решениями тестовых и модельных задач.

Практическая ценность. Разработанные двухшаговые методы могут быть использованы для решения систем нелинейных уравнений и в задачах ми-

з

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

Апробация работы. Основные результаты диссертационной работы докладывались и обсуждались на: XI международной конференции по вычислительной механике и современным прикладным программным системам (ВМСППС'2001, Москва-Истра, 2001); IV научно-практической конференции молодых ученых и специалистов Республики Татарстан (Казань, 2001); международной конференции ModelCARE 2002, 4th International Conference "Calibration and Reliability in Groundwater Modelling" (Prague, 2002); VIII Четаевской международной конференции "Аналитическая механика, устойчивость и управление движением" (Казань, 2002); II Международной конференции «Идентификация систем и задачи управления» (Москва, 2003); Всероссийской конференции "Высокопроизводительные вычисления и технологии (ВВТ-2003)" (Ижевск, 2003); XIII Всероссийской конференции-школе «Современные проблемы математического моделирования» (Абрау-Дюрсо, 2009); региональной научно-практической конференции «Актуальные вопросы геолого-гидродинамического моделирования и переоценки нефтяных ресурсов Республики Татарстан» (Казань, 2009); итоговых научных конференциях КазНЦ РАН; семинарах Института механики и машиностроения КазНЦ РАН.

Структура диссертационной работы. Диссертация состоит из введения, четырех разделов, заключения и списка литературы. Общий объем диссертации составляет 150 страниц, включая 26 таблиц и 46 рисунков.

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

СТРУКТУРА И КРАТКОЕ СОДЕРЖАНИЕ РАБОТЫ

Во введении дается обзор литературы, формулируются цели исследования и положения, выносимые на защиту.

В первом разделе дана постановка нелинейной задачи о наименьших квадратах. Приведены результаты поиска минимума тестовых функций стандартными методами безусловной минимизации. Для анализа существующих и построения новых методов минимизации функции невязки введены понятия запасов чувствительности переменной минимизации и функции невязки. На основе методов Ньютона и Гаусса-Ньютона минимизации функции невязки по-

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

В п.1.1 рассмотрена нелинейная задача о наименьших квадратах

где J(x) - функция невязки, х = (хи...,хп)т - вектор переменных минимизации, К = (>],..,,гт)т - вектор невязки, г, = г,(х) - нелинейные функции, п — число переменных минимизации, т - число невязок. Приведены ньютоновские и квазиньютоновские методы минимизации функции невязки: Ньютона и Гаусса-Ньютона без поиска и с поиском шага, Давидона-Флетчера-Пауэлла, Бройдена-Флетчера-Голдфарба-Шэнно, Левенберга-Марквардта.

В методе Левенберга-Марквардта новые значения переменных минимизации на к- ой итерации вычисляются по формуле хк = хкА + с1к, где

Нк -АТА, Л = [8>)Iдх]} - матрица чувствительности, цк >0 - параметр Мар-

квардга, Е - единичная матрица, gk - градиент функции невязки (вектор чувствительности). Различные модификации метода Левенберга-Марквардта отличаются стратегией изменения параметра в процессе минимизации. В работе использовались три варианта метода Левенберга-Марквардта. В первом варианте на первой итерации //, берётся на порядок больше максимального собственного числа матрицы Я1. На каждой итерации проверяется условие убывания функции невязки <У(дгм), и в случае его нарушения параметр Мар-квардта увеличивается в два раза до тех пор, пока это условие не выполнится. При переходе к следующей итерации параметр Марквардта уменьшается в два раза. Во втором варианте в начале каждой итерации = 0. Если условие убывания функции невязки нарушается, то цк меняется по формуле цк=\.5/1к+ 0.001 до тех пор, пока это условие не выполнится. В третьем варианте на каждой итерации цк определяется из условия минимума функции невязки.

Для остановки процесса минимизации в работе использовались два критерия: 1) J(xk~1)~ J(xk)<0.01J(xk~l) в течение трёх итераций (медленная схо-

приближённая матрица вторых производных,

димость итерационного процесса); 2) гтзх = тах[^(д:)| < Ю"6 (достижение заданной точности по невязкам).

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

В п.1.2 введены понятия запасов чувствительности переменной минимизации и функции невязки. Запас чувствительности г'-ой переменной минимизации определяется как величина Pj=-dfNgi, где (¡¡''х' и ^ ■ компоненты вектора спуска Гаусса-Ньютона с1ш и вектора чувствительности g. Вектор спуска Гаусса-Ньютона находится из решения системы НсР}> - . Запас чувствительности функции невязки определяется как сумма запасов чувствительности всех переменных минимизации. При п = т запас чувствительности функции невязки равен квадрату вектора невязки, а при п<т его не превосходит. В пространстве переменных минимизации с помощью сингулярного разложения матрицы Я = УЕУТ введена главная система координат, где V - ортогональная матрица, 2 = ((т,,...,сги) - диагональная матрица, - сингулярные числа,

сг, > <у2 > ...> ап > 0. В главной системе координат запас чувствительности в г'-м направлении РУ[ =(с1у^)2(т, является неотрицательной величиной и характеризует потенциальную возможность минимизации функции невязки в данном направлении - компоненты вектора с1у'4 = Утс1т). Величина запаса чувствительности функции невязки не меняется при переходе к главной системе координат пространства переменных минимизации.

В п.1.3 приведен алгоритм сингулярного разложения симметрической матрицы (БУТ) - разложение). С использованием преобразований Хаусхолдера матрица приводится к трёхдиагональному виду. Затем строится итерационный процесс <311-итераций со сдвигом Уилкинсона для преобразования трёхдиаго-нальной матрицы к диагональному виду.

В п.1.4 с использованием понятия запаса чувствительности проведён анализ метода Гаусса-Ньютона при минимизации функции Розенброка

/(х) = Юо(л:2 -л,2)2 +(1-Х\У- Линии уровня функции Розенброка имеют овражную структуру, дно оврага расположено вдоль параболы х2 =х?. В начальной точке Р0(-1;1), расположенной на дне оврага, запас чувствительности в главной системе координат в направлении первой оси = 0.0064 мал по сравнению с запасом в направлении второй оси РУ1 =3.9936. Направление вдоль первой оси, соответствующей большему сингулярному числу ег, =676.85, является направлением спуска ко дну оврага, а направление вдоль второй оси, соответствующей меньшему сингулярному числу а, =0.15, является смещением вдоль склона оврага. Значение функции Розенброка в начальной точке

б

/(Р0) = 4. Смещение переменных по направлению Гаусса-Ньютона с единичным шагом вследствие изгиба дна оврага привело в точку Р^^-З) на склоне оврага (рис.1, слева). Значение функции Розенброка /(Р,) = 1600 выросло по сравнению со значением в начальной точке, запас чувствительности в направлении первой оси увеличился РУ[ =1597.44, а в направлении второй оси уменьшился Руг = 2.56. По классическому методу Гаусса-Ньютона переход из точки Р0 в точку Р] невозможен, /(Р,)> /(Р0)- Смещение из точки Р1 по направлению Гаусса-Ньютона с единичным шагом привело в точку минимума Рт(1;1) (рис.1, справа), где =Руг= 0 и /(Рт) = 0. Таким образом, подъём на склон

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

РисЛ. Направления Гаусса-Ньютона функции Розенброка (слева- в точке Ро(-1;1), справа - в точке Р1(1;-3)). Линии уровня построены по логарифмам значений функции.

На основе проведённого анализа построены двухшаговые методы Гаусса-Ньютона и Ньютона без поиска и с поиском шага. На каждой к-ой итерации метода Гаусса-Ньютона с поиском шага:

1) Полагается Лк = 1.

2) Проверяется условие

■/(**-'+Л4 «О <У(*Ы),

где д°ык - направление Гаусса-Ньютона, вычисленное в точке хкА. Если это условие выполняется, то

хк^хк-1+Лк<1шк,

и осуществляется переход к следующей итерации, в противном случае - к пункту 3.

3) Проверяется условие убывания функции невязки при смещении переменных минимизации в два шага

+ Лк (с1тк + Зсмк)) < ДхкА), где ¿аык - направление Гаусса-Ньютона, вычисленное в точке хк'} + Лк ¿/СЛ к. Если это условие выполняется, то

хк=хк~х+Лк{с1апк+2ст), и осуществляется переход к следующей итерации, в противном случае полагается Лк =Лк/2 и возврат к пункту 2. Деление шага Лк прекращается, если два последних вычисленных значения функции отличаются менее чем на один процент от значения J(xk'l).

В двухшаговых методах Ньютона вместо отклонения Гаусса-Ньютона используется отклонение Ньютона.

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

Таблица 1. Численные результаты минимизации тестовых функций.

функция Н НПШ ДФП БФГШ гн гнпш ЛМ1 ЛМ2 ДН днпш ДГН дгнпш

квадратичная простой структуры 1\2 1\2 4\35 4\35 1\2 1\2 14\15 142 142 142 142 142

Розенброка * 21\29 23\333 201204 * 10433 26V31 6464 447 447 143 143

Розенброка, модификация I 7\8 7\8 6\53 6\53 * 3\5 15\16 3414 748 748 143 143

Розенброка, модификация II 4\5 4\5 5\41 5\40 2\3 2V3 16\17 243 445 445 243 243

Розенброка, модификация III * 26Y37 214317 21\225 * 4\10 25\27 3415 446 446 143 143

двумерная экспоненциальная * * 10432 10Ш1 * 7\9 21\22 7410 * * 648 648

Билл 7\8 7\8 8Y73 9Y79 * 5\8 14415 8426 748 748 * 5410

Зангвилла 142 га 3\20 3\29 1\2 142 12413 142 142 142 142 142

Пауэлла * ** 224323 21V244 12\13 12\13 30431 12413 * ** 12413 12413

* - нарушается условие убывания функции,

** - процесс минимизации прерван по критерию медленной сходимости, а\Ь : а - число итераций, Ъ - число вычислений целевой функции.

Методы: Н - Ньютона, НПШ - Ньютона с поиском шага, ДФП - Давидона-Флетчера-Пауэлла, БФГШ - Бройдена-Флетчера-Годдфарба-Шэнно, ГН - Гаусса-Ньютона, ГНПШ - Гаусса-Ньютона с поиском шага, ЛМ1 - Левенберга-Марквардта, вариант 1, ЛМ2 - Левенберга-Марквардга, вариант 2, ДН - двухшаговый Ньютона, ДНПШ - двухшаговый Ньютона с поиском шага, ДГН - двухшаговый Гаусса-Ньютона, ДГНПШ - двухшаговый Гаусса-Ньютона с поиском шага.

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

В п.2.1 описывается численное решение задачи определения поля напора при заданных граничных условиях в трехмерном анизотропном пласте в усло-

виях однофазной стационарной фильтрации жидкости, подчиняющейся линейному закону Дарси,

(31у(А^гаё/г) = 0.

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

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

1

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

В п.2.3 определены шесть модельных задач идентификации коэффициента фильтрации. Модельные задачи различаются числом и расположением наблюдательных точек, а также граничными условиями первого рода. Область решения всех модельных задач представляет собой пятислойный пласт. Каждый слой пласта разбит на зоны однородности, которые характеризуются двумя значениями коэффициента фильтрации К'гху< и К'[г Всего задана 71 зона однородности. При построении модельных задач из решения уравнения фильтрации определялись значения напора Щ в наблюдательных точках. Далее значения

К"ху1 и К" считались неизвестными, и требовалось определить их по значениям напора в наблюдательных точках Н' = Щ + ¿ь, где 3] - задаваемая погрешность. Начальные значения коэффициента фильтрации определялись при условии однородности пласта (два неизвестных значения коэффициента фильтрации).

В п.2.4 приведены результаты решения модельных задач идентификации коэффициента фильтрации без погрешности в замерах напора, полученные при минимизации функции невязки квазиньютоновскими методами (методы Гаусса-Ньютона, Давидона-Флетчера-Пауэлла, Бройдена-Флетчера-Голдфарба-Шэнно, Левенберга-Марквардта). Сложность этих задач заключается в большом числе переменных минимизации и большом числе обусловленности к

приближенной матрицы вторых производных функции невязки (109 <*"<1016). Использование методов Гаусса-Ньютона, Давидона-Флетчера-Пауэлла, Брой-дена-Флетчера-Голдфарба-Шэнно и второго варианта метода Левенберга-Марквардга не позволило получить решение с заданной точностью ни в одной модельной задаче. Заданная точность по напору была достигнута при использовании первого и третьего вариантов метода Левенберга-Марквардта при решении всех модельных задач, кроме третьей.

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

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

В п.3.1 на примере решения первой модельной задачи идентификации коэффициента фильтрации проведён анализ третьего варианта метода Левенберга-Марквардта с использованием понятия запаса чувствительности. Заданная точность по напору в наблюдательных точках в этой задаче достигается после выполнения 332-х итераций. Из рис.2 видно, что, начиная с некоторой итерации, процесс сходимости замедляется. Так, значения функции невязки на 50-ой и 51-ой итерациях соответственно равны 1.88х10"5 и 1.83х10'5 (значение параметра Марквардта на 51-ой итерации /¿5, = 4.11 х 10~5). Распределение запаса

чувствительности по направлениям в главной системе координат на 50-ой и 51-ой итерациях показано на рис.3.

549 <а « 7 1Е-3

ру,

-1-г-П I

Ру,

1Е-7

О 100 200 332 1 50 100 142 1 50 100 142

Рис.2. Значения по итерациям Рис.3. Распределение запаса чувствительности

удвоенной функции невязки. на 50-ой (слева) и 51-ой (справа) итерациях.

Из приведенных результатов видно, что распределение запаса чувствительности после проведения 51-ой итерации практически не меняется, причём весь запас чувствительности в основном сосредоточен на последних осях, соответст-

ю

вующих малым сингулярным числам. Направления минимизации в главной системе координат условно делились на две группы: направления, соответствующие большим сингулярным числам (смещение ко дну оврага), и направления, соответствующие маленьким сингулярным числам (смещение вдоль дна оврага). Далее проводился анализ изменения распределения запаса чувствительности при смещении переменных минимизации на 51-ой итерации с различными значениями параметра /¿51 в два шага. На первом шаге переменные

минимизации смещались на величину с1[ =-(# + g и принимали значе-

ния К50 +с11. На втором шаге переменные минимизации дополнительно смещались на величину с12 = Ус1у, где с!у =(¿/^,...,5^) - вектор смещения переменных в главной системе координат с ненулевыми компонентами в направлениях, соответствующих большим сингулярным числам матрицы Я, --gv /(сг,. + /751),= 1,...,<7, =0, г = 9+1,...,142, номер оси д выбирался из условия ст? > ¡.I,л > ст?+|. Распределения запаса чувствительности на 51-ой итерации после первого и второго шага при различных значениях параметра ¿и5[ показаны на рис.4,5.

1Е-3 ■

1Е-3 -1 1Е-4 1Е-5 1Е-6 -1Е-7

'У,

1Е-3 1Е-4 1Е-5 1Е-6 1Е-7

1Е-4 -1Е-5 1Е-6 I 1Е-7

ГУ,

50 100 142

1

//51 =4.11x10"

50 100 142 -5

50 100 142

50 100 142

(Хц =1x10"

/¿51 = 0.5 х 10 //51 =1x10"

Рис.4. Распределение запаса чувствительности при различных //5| после первого шага.

16-3 1Е-4 1Е-5 1Е-6 -1Е-7

1

—1-1—"1 I

50 100 142

цьх= 4.11x10"

50 100 142 -5

50 100 142

1 50 100 142

-6

цьх =1x10 =0.5x10"' /¿л =1x10"

Рис.5. Распределение запаса чувствительности при различных /У51 после второго шага.

При уменьшении параметра /¿51 запас чувствительности после первого шага растёт на осях, соответствующих большим сингулярным числам, при этом также растёт функция невязки (табл.2), а запас чувствительности на осях, соответствующих малым сингулярным числам, уменьшается. После второго шага запас чувствительности на осях, соответствующих большим сингулярным числам, и значение функции невязки уменьшаются (табл.2). В результате выполнения

двух шагов при /х51 =10 получено значение функции невязки меньшее, чем при /и51 = 4.11х10~5.

Таблица 2. Первая модельная задача идентификации коэффициента фильтрации.

>"51 начальное значение функции невязки функция невязки после 1-го шага функция невязки после 2-го шага

4.11x10° 1.88хЮ'5 1.84x10° 1.83x10°

1.00x10° 2.73x10° 1.75x10°

0.50x10° 4.99x10° 2.28x10°

1.00x10'" б.згхю"1 2.39х10"4

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

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

1) Вычисляется значение

где с1х = -(# + /JnEy1g - вектор спуска, аналогичный вектору спуска метода Левенберга-Марквардта.

2) Вычисляется значение

где О1 = Ус1у характеризует спуск ко дну оврага; ¿у - вектор смещения переменных с компонентами с!у. =-,§>, /(сг, +/'„), г=1,...,д, с1у. =0, г = ^ + 1,...,и; ¡¡у - компоненты вектора -¥т§; # = АтЯ, Я - вектор невязок в точке К"'1 + й?1 ; номер оси д выбирается из условия сг? > ¡л п > сг?+1.

3)Определяется У

Методом золотого сечения по параметру вычисляется = гшп./^. При

Ал

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

(к"-'+с1', если У1, <32, К" = \ ^

1 К"'1 + + с!2, если У1. >

1. Мп Ип

где //* значение , при котором вычислено ./*.

Для определения вектора используются значения невязок в точке К"'1, а для определения вектора с!2 - значения невязок в точке + В обоих случаях

используется матрица чувствительности А, вычисленная в точке К"'1. Аналогичным образом построены двухшаговые методы Левенберга-Марквардта с другими процедурами поиска параметра Марквардта.

В п.3.2 приведены численные результаты решения модельных задач идентификации коэффициента фильтрации и минимизации тестовых функций двухшаговыми методами Левенберга-Марквардта. Проведено сравнение результатов, полученных методами Левенберга-Марквардта и двухшаговыми методами Левенберга-Марквардта. Основные вычислительные затраты при минимизации функции невязки при решении задачи идентификации коэффициента фильтрации приходятся на вычисление функции невязки и на вычисление элементов матрицы чувствительности. Для оценки этих затрат введено число пс = пс\ +пс2, где пс\ - число решений уравнения фильтрации, пс1 — число решений уравнений, полученных прямым дифференцированием уравнения фильтрации (используются для определения элементов матрицы чувствительности). Результаты решения первой модельной задачи методами Левенберга-Марквардта и двухшаговыми методами Левенберга-Марквардта приведены в табл.3.

Таблица 3. Первая модельная задача идентификации коэффициента фильтрации.

методы начальное состояние конечное состояние

г шах Д1п Кт г шах Д1п Кт и ПС

ЛМ1 7.18 1.7 9х10"7 0.24 248 35670

ДЛМ1 9x10'7 0.27 110 15946

ЛМЗ 1х10'6 0.24 332 50700

ДЛМЗ 8х10"7 0.26 82 13368

Д 1п К^ = ((1п К^- 1п Кф / + (1п К % - 1п К!к )2)/1421 - среднеквадратическое отклонение, гтах- максимальная невязка, и - число итерадий, ЛМ1 -метод Левенберга-Марквардта, вариант 1; ЛМЗ - метод Левенберга-Марквардта, вариант 3; ДЛМ1 - двухшаго-вый метод Левенберга-Марквардта, вариант 1; ДЛМЗ - двухшаговый метод Левенберга-Марквардта, вариант 3.

Аналогичные результаты получены и при решении других модельных задач. Двухшаговые методы Левенберга-Марквардта показали более высокую скорость минимизации функции невязки и потребовали меньших вычислительных затрат в задачах без погрешности и с погрешностями в замерах напора.

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

i

1 50 100 142

Рис.6. Распределение запаса чувствительности на 15-ой итерации при {¡15 = 1.34 х 10 5: в начале (слева), после первого шага (в центре) и после второго шага (справа).

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

/«15 начальное значение функции невязки функция невязки после 1-го шага функция невязки после 2-го шага функция невязки после дополнительных шагов

1.34х10"5 7.12x10"5 1.26x1o"4 4.73xl0"i 4.57х10"5

1.0x10"5 2.16x10"4 4.75х10'5 4.1х10"5

0.5х10"5 9.55x1o"4 1.05Х10"4 3.19x10°

В результате выполнения дополнительных шагов при /¿l5 = 0.5x10 5 получено значение функции невязки меньшее, чем при //15 =1.34х10"5 (значение параметра Марквардта, полученное по двухшаговому методу Левенберга-Марквардта). Распределение запаса чувствительности по отдельным шагам при //15 = 0.5 х 10~5 показано на рис.7.

1Е"31 Р

1Е-4 -

1Е-5 - .

1Е-6 - Б

1Е-7 -1-г-Ч i

1 50 1С0 142

Рис.7. Распределение запаса чувствительности на 15-ой итерации при /л[5 - 0.5 х 10~5: после первого (слева), второго (в центре) и восьмого (справа) шагов для.

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

so 100 142

S0 100 142

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

Таблица 5. Первая модельная задача идентификации коэффициента фильтрации.

методы начальное состояние конечное состояние

г шах Д1п г тах Д1п Кху1 и пс

ДЛМ1 9x10"7 0.27 110 15946

ДЛМ1М 7.18 1.7 7x10"7 0.27 59 8666

ДЛМЗ 8х10"7 0.26 82 13368

ДЛМЗМ 4х10"7 0.25 17 4534

ДЛМ1 - двухшаговый метод Левенберга-Марквардта, вариант 1; ДЛМЗ - двухшаговый метод Левенберга-Марквардта, вариант 3; ДЛМ1М - модифицированный двухшаговый метод Левенберга-Марквардта, вариант 1; ДЛМЗМ - модифицированный двухшаговый метод Левенберга-Марквардта, вариант 3; /"тах - максимальная невязка, и - число итераций, Д]п А^

- среднеквадратическое отклонение.

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

Кху\ ~Кху2-■■■-Кхут> К:\ - ■ ■ ■ - К тт■

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

В п.4.1 на основе метода Левенберга-Марквардта с поиском параметра Марквардта по методу золотого сечения построен метод, сохраняющий упорядоченность идентифицируемых парамегров. В этом методе для каждого значения //„ формируются переменные минимизации х-. Переменная минимизации

х^ объединяет все параметры с одинаковыми значениями Кд =... = Кп. Значение переменной х} совпадает со значением параметров в соответствующей ей группе После формирования переменных последовательно

перебираются разбиения каждого набора на две части [к

и где у, <у'с < j2. Каждому такому разбиению соответствует

множество переменных х}, число которых на единицу больше числа переменных х]. Разбиение считается недопустимым, если любое смещение переменных минимизации Xj по вычисленному для них направлению спуска

={а?,...,<1%х)т =-{H + ^nE)~lg ведёт к нарушению упорядоченности параметров (их - число переменных минимизации х^ ). Из всех допустимых разбие-

ний выбирается разбиение с максимальным запасом чувствительности функции невязки /*тах. Если Ртах больше запаса чувствительности функции невязки для переменных х}, то делается переход к соответствующим переменным минимизации ЗЕу (число переменных минимизации увеличивается на единицу). Процедура увеличения числа переменных минимизации повторяется до тех пор, пока возрастает запас чувствительности функции невязки.

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

1) Суммарный шаг ^ полагается равным нулю.

2) Вычисляется р = тт((х,+] -ху- )) для ¿1 >0, р - максимальное значение шага р, при котором выполняется условие сохранения упорядоченности параметров х^ + р < х]А + р . Определяется шаг р'=тт{р,\-р,}.

3) Определяются новые значения переменных х} - + р" ^, вычисляется

суммарный шаг pJ = р5 + р* , и проверяется условие р8 = 1. При его нарушении переходим к пункту 4. В противном случае смещение переменных заканчивается, и вычисляется значение функции невязки.

4) Переменные, удовлетворяющие условию

Х] =■*/+! и ¿1 (1) объединяются, и проводится пересчет отклонений. Объединение переменных и пересчет отклонений для новых переменных повторяется до тех пор, пока остаются переменные, удовлетворяющие условию (1). После этого процедура повторяется, начиная с шага 2.

Аналогичным образом построен метод Гаусса-Ньютона с учетом априорной сравнительной информации.

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

В п.4.3 с учётом априорной сравнительной информации о значениях идентифицируемых параметров построен двухшаговый метод минимизации функции невязки. Результаты, полученные при решении первой модельной задачи, приведены в табл.6. Аналогичные результаты получены при решении остальных модельных задач без погрешности и с погрешностями в замерах напора. Двухшаговый метод Левенберга-Марквардта, сохраняющий упорядоченность идентифицируемых параметров, показал себя наиболее эффективным по вычислительным затратам по сравнению со всеми рассмотренными методами.

Таблица 6. Первая модельная задача идентификации коэффициента фильтрации.

метод начальное состояние конечное состояние

У 'тах Д1п К^ г шах и пс

ДЛМЗ 7.18 1.7 8х10'7 0.26 82 13368

ЛМЗ А 9x10"7 0.01 41 6251

ДЛМЗА 8х10"7 0.01 21 3420

JIM3A - метод ЛМЗ с учётом априорной сравнительной информации; ДЛМЗ - двухшаговый метод с учётом априорной сравнительной информации; ДЛМЗА - двухшаговый метод Левенберга-Марквардта, сохраняющий упорядоченность идентифицируемых параметров; rmax-максимальная невязка, Д In - среднеквадратическое отклонение, /I - число итераций.

В заключении приводятся основные результаты работы.

ОСНОВНЫЕ РЕЗУЛЬТАТЫ

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

2. На основе двухшаговых методов Левенберга-Марквардта разработаны вычислительные алгоритмы для решения задач идентификации коэффициента фильтрации неоднородных пластов по замерам напора в наблюдательных точках. Алгоритмы протестированы на модельных задачах идентификации коэффициента фильтрации трёхмерного пласта с учётом и без учёта априорной сравнительной информации о значениях идентифицируемых параметров. Показана более высокая эффективность по вычислительным затратам разработанных двухшаговых алгоритмов по сравнению с классическими.

СПИСОК ОПУБЛИКОВАННЫХ РАБОТ ПО ТЕМЕ ДИССЕРТАЦИИ

1. Елесин A.B. К решению обратной задачи по определению коэффициента фильтрации трехмерного напорного пласта / A.B. Елесин, А.Н. Габидуллина, А.Ш. Кадырова // Труды математического центра им. Н.И.Лобачевского. - Казань, "Унипресс", 1998.-С. 103-105.

2. Елесин A.B. Исследование минимизационных свойств метода Марквардта при идентификации коэффициента фильтрации напорного анизотропного пласта / A.B. Елесин, А.Ш. Кадырова // Тезисы докладов IV-й научно-практической конференции молодых ученых и специалистов Республики Татарстан. - Казань, 2001. - С. 65.

3. Габидуллина А.Н. К идентификации коэффициента фильтрации трёхмерного напорного анизотропного пласта / А.Н. Габидуллина, A.B. Елесин, А.Ш. Кады-

рова, П.А. Мазуров // Математическое моделирование. -2002. - Т. 14. - №9. - С. 97-102.

4. Мазуров П. А. Определение параметров водоносных пластов с использованием анализа чувствительности / ПЛ. Мазуров, A.B. Елесин, А.Н. Габидуллина, А.Ш. Кадырова // Современные проблемы гидрогеологии и гидромеханики. Сб. докл. конференции. - СПб., 2002. - С. 462-471.

5. Mazurov P.A. Use of minimization along the slope for estimation of aquifer parameters / P.A. Mazurov, А.У. Elesin, A.N. Gabidullina, A.Sh. Kadyirova // Proceedings of the 4th International Conference on Calibration and reliability in groundwater modelling. - Prague, Czech Republic, 17-20 June 2002. -V.l. - P. 278-281.

6. Мазуров П.А. Новый метод минимизации функции невязки при идентификации параметров водоносных слоев / П.А. Мазуров, А.Н. Габидуллина, A.B. Елесин, А.Ш. Кадырова // Труды II Международной конференции "Идентификация систем и задачи управления". - Москва, 29-31 января 2003. - С.714-727.

7. Мазуров П.А. Запасы чувствительности в задачах идентификации коэффициента фильтрации трехмерных пластов / П.А. Мазуров, А.Н. Габидуллина, A.B. Елесин, А.Ш. Кадырова // Вычислительные методы и программирование. -2004. - Т.5. - № 1,-С. 50-61.

8. Габидуллина А.Н. Идентификация коэффициента фильтрации с учётом сравнительной информации о значениях коэффициента фильтрации / А.Н. Габидуллина, A.B. Елесин, А.Ш. Кадырова, П.А. Мазуров // Актуальные проблемы механики сплошной среды. К 15-летию ИММ КазНЦ РАН. - Казань, КГУ, 2006. С. 172-178.

9. Елесин A.B. Учёт априорной сравнительной информации в задачах идентификации коэффициента фильтрации / A.B. Елесин, А.Ш. Кадырова // Вычислительные методы и программирование. - 2008. - Т.9. - № 1. - С. 14-19.

10. Елесин A.B. Двухшаговые методы Левенберга-Марквардта в задаче идентификации коэффициента фильтрации / A.B. Елесин, А.Ш. Кадырова, П.А. Мазуров // Георесурсы. - 2009. - 4(32). - С.40-42.

11. Мазуров П.А. Квазиньютоновский двухшаговый метод минмизации функции невязки / П.А. Мазуров, A.B. Елесин, А.Ш. Кадырова // Вычислительные методы и программирование. - 2009. - Т. 10. - № 1. - С. 64-71.

12. Елесин A.B. Построение двухшаговых методов минимизации функции невязки / A.B. Елесин, А.Ш. Кадырова, П.А. Мазуров // Современные проблемы математического моделирования. Математическое моделирование, численные методы и комплексы программ. Сб.трудов научных молодежных школ. - Ростов-на-Дону, 2009. - с.240-247.

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

в типографии издательства Казанского государственного университета Тираж 120 экз. Заказ 104/2

420008, ул. Профессора Нужина, 1/37 тел.: (843)233-73-59, 292-65-60

Основные обозначения.

Введение.

Раздел 1. Двухшаговые методы Ньютона и Гаусса-Ньютона минимизации функции невязки.

1.1. Задача минимизации функции невязки.

1.2. Запасы чувствительности.

1.3. Сингулярное разложение симметрической матрицы.

1.4. Построение двухшаговых методов Ньютона и Гаусса-Ньютона минимизации функции невязки.

Раздел 2. Идентификация коэффициента фильтрации неоднородного пласта.

2.1. Численное решение уравнения фильтрации.

2.2. Постановка задачи идентификации коэффициента фильтрации.

2.3. Модельные задачи.78*

2.4. Результаты решения модельных задач без погрешностей в замерах напора.

2.5. Решение модельных задач с погрешностями в замерах напора.

Раздел 3. Двухшаговые методы Левенберга-Марквардта минимизации функции невязки.

3.1. Построение двухшаговых методов Левенберга-Марквардта.

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

3.3. Модификации двухшаговых методов5 Левенберга-Марквардта.

Раздел 4. Учет априорной сравнительной информации в задачах идентификации коэффициента фильтрации.

4.1. Метод минимизации функции невязки, учитывающий сравнительную информацию.

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

4.3. Двухшаговые методы минимизации функции невязки, учитывающие сравнительную информацию.

При решении многих важных практических задач управления, идентификации и адаптации возникает нелинейная задача о наименьших квадратах где J(x) - функция невязки (целевая функция), х = (x,,.,xn)r - вектор переменных минимизации, R = (rx,.,rm)T - вектор невязки, rt = ri(х) - нелинейные функции, п — число переменных минимизации, т - число невязок. При численном решении нелинейной задачи о наименьших квадратах, как правило, строразуют убывающую последовательность, т.е. значения целевой функции в точках последовательности удовлетворяют неравенству (условие убывания функции невязки) минимизации различаются способами построения релаксационной последовательности и в зависимости от порядка используемых производных делятся на три группы [3,4,7,44,47,51,56]: методы нулевого порядка (или методы прямого поиска), методы первого (градиентные методы) и второго порядка. Методы нулевого порядка используют только значения целевой функции. К методам нулевого порядка относятся метод конфигураций, метод деформируемого многогранника, метод Розенброка, методы случайного поиска и др. В методах первого порядка используются значения целевой функции и значения её частных производных первого порядка. Наиболее часто используемые методы этого типа: градиентный метод, метод наискорейшего спуска и метод сопряженных градиентов. В методах второго порядка (метод Ньютона, метод Ньютона с по

14,19,90]: ится последовательность точек хА в которых значения целевой функции об

J{xk)< J(xk~l).

Такая последовательность точек |хЧ называется релаксационной [3]. Методы иском шага) используются значения целевой функции и значения её частных производных до второго порядка включительно. При решении практических задач построение матрицы вторых производных (матрица Гессе) требует больших вычислительных затрат. Поэтому на практике обычно применяются методы квазиньютоновского типа. В этих методах используется аппроксимация матрицы вторых производных или её обратной с помощью частных производных первого порядка. При приближении к точке минимума эти методы по своим свойствам приближаются к методам второго порядка. К методам, в которых аппроксимируется матрица обратная к матрице Гессе, относятся методы Дави-дона-Флетчера-Пауэлла, Бройдена-Флетчера-Голдфарба-Шэнно, Пауэлла и др., в методе Гаусса-Ньютона аппроксимируется матрица Гессе.

Метод Гаусса-Ньютона лежит в основе многих методов квазиньютоновского типа и создан специально для минимизации функций, имеющих вид суммы квадратов. В этом методе матрица Гессе аппроксимируется'с помощью матдг

- > - матрица чувствительности, и переход на новые cbcfJ значения переменных осуществляется, по формуле x"ew -xo!d + dGN, dGN = -H~lg - направление спуска Гаусса-Ньютона, g - градиент функции невязки. В процессе-минимизации методом Гаусса-Ньютона матрица Н может оказаться плохо обусловленной или вырожденной, что вносит большие погрешности в определение направления спуска, и значение функции невязки при переходе на новые значения переменных может вырасти. Различные модификации метода Гаусса-Ньютона направлены на устранение этих недостатков (метод Гаусса-Ньютона с поиском шага, метод Левенберга-Марквардта).

В методе Гаусса-Ньютона с поиском шага используются различные процедуры поиска шага вдоль направления спуска. Длина шага определяется из условия минимума функции вдоль направления спуска с использованием методов одномерной минимизации или с помощью процедуры дробления шага. рицы Н = А А, где А =

Другой модификацией метода Гаусса-Ньютона является метод Левенбер-га-Марквардта. В этом методе последовательные приближения определяются формулой x"cw = xold - (Н + juE)~] g, где ju> 0 - параметр Марквардта. Эта формула была независимо предложена Левенбергом [78] и Марквардтом [81]. При больших значениях параметра Марквардта направление минимизации метода Левенберга-Марквардта близко к направлению метода наискорейшего спуска, при малых значениях - к направлению метода Гаусса-Ньютона. В начале процесса минимизации, когда текущая точка расположена вдали от точки минимума, выбираются большие значения параметра Марквардта. При приближении к точке минимума значения параметра Марквардта уменьшаются. Таким образом, метод Левенберга-Марквардта, в отличие от метода Гаусса-Ньютона, не требует хорошего начального приближения. Направление спуска метода Левенберга-Марквардта d ~-{Н + juE)~l g является решением задачи на условный минимум [14,19] mmi|| Ad + R\2 при ограничении \d\2 < Д, где А - матрица чувствительности, R - вектор невязок, А - параметр, связанный с параметром jli , определяет размер доверительной области. Поэтому метод Левенберга-Марквардта можно отнести к классу методов доверительной области [19]. Существуют различные модификации метода Левенберга-Марквардта, отличающиеся стратегиями выбора параметра Марквардта. В [44,90] начальное значение параметра Марквардта выбирается достаточно большим. В случае выполнения условия убывания функции невязки на итерации при переходе на следующую итерацию параметр Марквардта уменьшается, а при нарушении этого условия параметр Марквардта увеличивается до тех пор, пока оно не выполнится. В [44] коэффициент изменения параметра Марквардта берётся равным двум, а в [90] - десяти. В' [75] в начале каждой итерации параметр Марквардта равен нулю. В случае нарушения условия убывания функции невязки параметр Марквардта увеличивается по формуле /л = 1.5// + 0.001 до тех пор, пока это условие не выполнится. Этот вариант метода Левенберга-Марквардта используется в пакетах UCODE и MODFLOWP, предназначенных для решения обратных задач. Вариант метода Левенберга-Марквардта, в котором подбирается А, а // вычисляется по значению А приводится в [86]. Этот вариант включен в пакет программ MINPACK. Процедура подбора значения параметра Марквардта в зависимости от отношения между фактическим изменением функции невязки в результате пробных шагов и ожидаемыми величинами этих изменений построена в [68]. В [34] параметр Марквардта определяется методом золотого сечения из условия минимума функции невязки.

Для проверки эффективности методов минимизации обычно используются специальные тестовые функции. Различные примеры тестовых функций можно найти в [7,59,60,87].

К задаче минимизации функции невязки сводится задача идентификации коэффициента фильтрации водоносных и нефтяных пластов при фильтрации в них жидкости по замерам напора в наблюдательных точках [75,90,92]. В этом случае функция невязки имеет вид суммы квадратов разностей между вычисленными и измеренными значениями напора в наблюдательных точках. Коэффициент фильтрации ищется в виде кусочно-постоянной функции. Область решения задачи разбивается на зоны однородности, каждая из которых характеризуется своими значениями коэффициента фильтрации [64,66,89]. Задача определения поля напора при известном коэффициенте фильтрации является прямой задачей, а задача идентификации коэффициента фильтрации по замерам напора в наблюдательных точках - обратной коэффициентной задачей.

Задача идентификации коэффициента фильтрации относится к классу некорректно поставленных задач [1,5,6,10,17,18,41,50,52,53,58,90]. Задача называется некорректно поставленной, если нарушается одно из условий: решение существует, решение единственно, решение устойчиво относительно изменения исходных данных. В задачах определения параметров водоносных пластов может нарушаться любое из условий корректности. Примеры таких задач приведены в [62,83,90,91]. При решении обратных коэффициентных задач из-за наличия погрешностей обычно на первых итерациях значения параметров приближаются, а, начиная с некоторых итераций, удаляются от своих истинных значений, при этом функция невязки продолжает уменьшаться. Для выбора номера итерации с более достоверными значениями коэффициента фильтрации применяются специальные правила останова процесса минимизации или прерывания полученной последовательности значений функции невязки. Эти правила являются одним из регуляризирующих элементов решения обратных задач [1,5,6,41,52].

При решении задач идентификации коэффициента фильтрации с большим числом идентифицируемых параметров обычно используются квазиньютоновские методы минимизации функции невязки. Для вычисления .градиента функции невязки и матрицы чувствительности используется один из методов: метод конечно-разностных соотношений, вариационный метод, метод прямого дифференцирования [90,92]. Недостатком метода конечно-разностных соотношений является, сложность вычисления производных с требуемой точностью. Более точные значения производных получаются при - использовании метода-прямого дифференцирования и вариационного метода. В [54] приведена общая схема построения вариационного метода, в [76] вариационный метод впервые использован при решении задач идентификации параметров. Сравнение трёх перечисленных выше методов вычисления производных при решении задачи идентификации коэффициента фильтрации приводится в [22]. Для нахождения^ матрицы чувствительности методом конечно разностных соотношений или методом прямого дифференцирования необходимо решить п+1 уравнений, аналогичных уравнению фильтрации (п - число идентифицируемых параметров), а в случае использования вариационного метода т+1 уравнений (т — число замеров напора).

В процессе решения задачи идентификации коэффициента фильтрации на каждой итерации требуется неоднократно решать прямую задачу. Решение прямой задачи проводится численными методами. Система дифференциальных уравнений с частными производными (уравнение фильтрации и соответствующие начальные и граничные условия) сводится к системе линейных алгебраических уравнений. Для этого используются метод конечных разностей, метод конечных элементов, метод контрольных объёмов и др. [8,15,26,36,40,42,4850,55,65,69,71,72,79]. В практических задачах получаемая система алгебраических уравнений имеет большую размерность и обычно решается одним из итерационных методов, учитывающих симметричность и сильную разреженность матрицы системы. Одним из наиболее часто используемых методов является метод сопряженных градиентов с предобусловливанием [74,77,84,85]. Скорость сходимости и эффективность этого метода сильно зависят от используемой предобусловливающей матрицы. Сравнение результатов, полученных методом сопряженных градиентов с использованием трёх различных предобусловли-вающих матриц (неполное разложение Холесского, модифицированное неполное разложение Холесского, полиномиальный метод), проведено в [74]. В [77] приведены результаты решения девяти трёхмерных задач по определению поля напора в условиях однофазной стационарной фильтрации жидкости, полученные методом сопряженных градиентов без предобусловливания и с использованием для построения предобусловливающей' матрицы диагонального масштабирования, неполного разложения Холесского, неполной факторизации, модифицированной неполной факторизации. Другие подходы к построению предобусловливающей матрицы используются в [70,80,88].

Для тестирования методов решения задачи идентификации коэффициента фильтрации используются модельные задачи [63]. В модельных задачах по известным значениям коэффициента фильтрации (истинные значения) определяются значения напора в наблюдательных точках. Далее значения коэффициента фильтрации считаются неизвестными, и требуется их определить по значениям напора в наблюдательных точках. Число и границы зон однородности могут считаться как известными, так и неизвестными. В последнем случае их также требуется определить в процессе идентификации. Методы решения отрабатываются на модельных задачах, как без погрешностей, так и с погрешностями, вносимыми в задачу. В модельных задачах без погрешностей решение всегда существует. Число и положение наблюдательных точек должно быть достаточным для выполнения условия единственности решения в задаче без погрешностей, что проверяется численными экспериментами. Полученные в- результате решения модельной задачи идентификации значения коэффициента фильтрации всегда можно сравнить с истинными значениями, и по результатам сравнения оценить эффективность метода решения задачи. При решении практической задачи идентификации коэффициента фильтрации приходится многократно строить и решать различные модельные задачи. При решении модельных задач приходится определять достаточность исходных данных для идентификации того или иного набора параметров, влияние погрешностей в исходных данных и используемых регуляризирующих процедур. Всё это требует больших вычислительных затрат.

Основной целью данной работы является разработка методов*минимизации функции невязки, сокращающих вычислительных затраты. В работе предложены двухшаговые методы минимизации функции невязки: двухшаговые методы Ньютона, двухшаговые методы Гаусса-Ньютона и двухшаговые методы Левенберга-Марквардта. В двухшаговых методах первый <, шаг каждой итерации проводится по алгоритмам классических методов, но допускается-увеличение функции невязки. Итоговые же значения функции невязки на итерациях, как и в классических методах, образуют убывающую последовательность. Для анализа известных методов минимизации функции невязки и построения двухшаговых методов введено понятие запаса чувствительности. При построении двухшаговых методов используется главная система координат в пространстве переменных минимизации. Главная система координат вводится с помощью сингулярного разложения симметрической приближённой матрицы вторых производных [16,19]. На первом этапе сингулярного разложения с помощью преобразований Хаусхолдера строится подобная трёхдиагональная симметрическая матрица, которая на втором этапе итерационно приводится к диагональному виду с помощью QR-итераций с одинарным сдвигом. В главной системе координат запас чувствительности характеризует потенциальную возможность к минимизации функции невязки вдоль каждого из направлений.

Построены двухшаговые методы Ньютона и Гаусса-Ньютона минимизации функции без поиска и с поиском шага, обобщающие классические методы Ньютона и Гаусса-Ньютона. В двухшаговых методах Ньютона и Гаусса-Ньютона при минимизации овражных функций, когда текущая точка расположена вблизи дна оврага, первый шаг соответствует движению вдоль дна оврага. При изменении направления дна оврага первый шаг может привести к подъёму на его склон (значение функции невязки увеличивается), а второй шаг к спуску ко дну оврага (значение функции невязки становится меньше начального значения на итерации). Проведение таких двух шагов на одной итерации позволяет обойти изгибы дна оврага и тем самым ускорить процесс минимизации. Эти методы протестированы на наиболее часто встречающихся тестовых функциях (квадратичная функция простой структуры, функция Розенброка, двумерная экспоненциальная функция, функции Биля-, Пауэлла и Зангвилла). Результаты тестирования показали, что эти методы наиболее эффективны при минимизации функций с сильным изгибом дна оврага.

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

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

Предложенные двухшаговые методы минимизации функции невязки относятся к монотонным методам. В работах [67,93] рассмотрены немонотонные методы, в которых допускается увеличение функции невязки на отдельных итерациях. В немонотонных методах при переходе на новую итерацию значение функции невязки сравнивается не с предыдущим значением, а с максимальным значением из М предыдущих итераций, где М - задаваемый параметр.

Численные алгоритмы, основанные на двухшаговых методах Левенберга-Марквардта, протестированы на шести модельных задачах идентификации коэффициента фильтрации трехмерного анизотропного пласта по замерам напора в наблюдательных точках в условиях однофазной стационарной фильтрации жидкости, подчиняющейся закону Дарси. Модельные задачи представляют собой водоносный пласт, состоящий из пяти слоев. Каждый слой разбит на зоны однородности коэффициента фильтрации (71 зона однородности). Каждая зона однородности характеризуется двумя значениями коэффициента фильтрации. В задачах требовалось восстановить 142 значения коэффициента фильтрации. На кровле пласта заданы граничные условия второго рода. Подошва и боковая поверхность пласта считались непроницаемыми, за исключением участков боковой поверхности пятого слоя. Модельные задачи отличаются числом и положением наблюдательных точек, граничными условиями. Для определения поля напора уравнение фильтрации с соответствующими граничными условиями методом конечных элементов Галёркина аппроксимировалось системой линейных алгебраических уравнений. Полученная система уравнений решалась методом сопряжённых градиентов с предобусловливанием в виде неполного разложения Холесского. Значения коэффициента фильтрации определялись из минимума функции невязки, имеющей вид суммы квадратов разности между вычисленными и измеренными значениями напора в наблюдательных точках. Сравнение результатов решения модельных задач идентификации коэффициента фильтрации классическими и двухшаговыми методами показало^ что; по вычислительным затратам наиболее эффективными являются двухшаговые методы Левенберга-Марквардта с дополнительными спусками ко дну оврага. При решении задач с погрешностями в замерах напора значения» идентифицируемых параметров, как правило, начиная с некоторой итерации, удаляются от своих истинных значений, хотя функция невязки продолжает уменьшаться. Для получения итоговых значений коэффициента фильтрации более близких к истинным, в работе использовался критерий выбора номера итерации с итоговыми значениями идентифицируемых параметров.

На основе классических методов Гаусса-Ньютона и Левенберга-Марквардта, а также двухшагового метода Левенберга-Марквардта, построены методы, сохраняющие априорную упорядоченность значений идентифицируемых параметров. Учёт различного рода априорной информации о значениях идентифицируемых параметров является одним из регуляризирующих элементов, повышающим устойчивость решения задачи. Результаты решения модельных задач показали, что учёт априорной сравнительной информации о значениях идентифицируемых параметров позволяет сократить вычислительные затраты, при этом итоговые значения идентифицируемых параметров получаются более близкими к своим истинным значениям. Из методов, учитывающих сравнительную априорную информацию, наилучшим показал себя двухшаговый метод.

Цели диссертационной работы:

- построение эффективных методов минимизации функции невязки;

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

Структура и краткое содержание работы.

Диссертация состоит из введения, четырёх разделов, заключения и списка литературы. Работа содержит: страниц - 150, рисунков - 46, таблиц - 26, список литературы - 93 наименование.

Выводы.

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

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

Заключение.

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

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

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

На примерах минимизации тестовых функций (квадратичная функция простой структуры, функция Розенброка, двумерная экспоненциальная функция, функции Биля, Пауэлла и Зангвилла) проведено сравнение известных методов: Ньютона, Гаусса-Ньютона, Ньютона и Гаусса-Ньютона с поиском шага, Давидона-Флетчера-Пауэлла, Бройдена-Флетчера-Голдфарба-Шэнно, Левенберга-Марквардта с двухшаговыми методами Ньютона и Гаусса-Ньютона без поиска и с поиском шага. Результаты тестирования показали, что двухшаговый метод Гаусса-Ньютона с поиском шага позволил достичь точки минимума для всех рассмотренных тестовых функций с наименьшими вычислительными затратами.

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

На основе методов Левенберга-Марквардта, Гаусса-Ньютона и двухшагового метода Левенберга-Марквардта построены алгоритмы минимизации функции невязки, учитывающие априорную сравнительную информацию о значениях идентифицируемых параметров. Учёт априорной сравнительной информации о значениях идентифицируемых параметров позволил сократить вычислительные затраты и получить итоговые значения коэффициента фильтрации более близкие к истинным. Из этих методов наиболее эффективным по вычислительным затратам показал себя двухшаговый метод.

1. Алифанов О.М. Экстремальные методы решения некорректных задач / О.М. Алифанов, Е.А. Артюхин, С.В. Румянцев - М.: Наука, 1988. - 288 с.

2. Андреева Е.А. Вариационное исчисление и методы оптимизации / Е.А. Андреева, В.М. Цирулева М.: Высш.шк., 2006.- 584 с.

3. Аттетков А.В. Методы оптимизации / А.В. Аттетков, С.В. Галкин, B.C. Зарубин М.: Изд-во МГТУ им. Н.Э.Баумана, 2003. - 440 с.

4. Аттетков А.В. Введение в методы оптимизации / А.В. Аттетков, B.C. Зарубин, А.Н. Канатников- М.: Финансы и статистика; ИНФРА-М, 2008. 272 с.

5. Бакушинский А.Б. Итеративные методы решения некорректных задач / А.Б. Бакушинский, А.В. Гончарский М.: Наука, 1989. - 128 с.

6. Бакушинский А.Б. Итерационные методы решения некорректных операторных уравнений с гладкими операторами / А.Б. Бакушинский, М.Ю. Кокурин М.: Едиториал УРСС, 2002. - 192 с.

7. Банди Б. Методы оптимизации. Вводный курс / Б. Банди М.: Радио и связь, 1988.- 128 с.

8. Булыгин В.Я. Гидромеханика нефтяного пласта / В.Я. Булыгин М.: Недра, 1974.-232 с.

9. Васильев Ф.П. Лекции по методам решения экстремальных задач / Ф.П. Васильев М.: Изд-во МГУ, 1974. - 376 с.

10. Ю.Васин В.В. Некорректные задачи с априорной информацией / В.В. Васин, А.Л. Агеев Екатеринбург: УИФ Наука, 1993. 262 с.

11. П.Вержбицкий В.М. Численные методы (линейная алгебра и нелинейные уравнения) / В.М. Вержбицкий М.: Оникс 21 век, 2005. - 432 с.

12. Габидуллина А.Н. К идентификации коэффициента фильтрации трёхмерного напорного анизотропного пласта / А.Н. Габидуллина, А.В. Елесин, А.Ш. Кадырова, П.А. Мазуров // Математическое моделирование. -2002. Т.14. -№9.-С. 97-102.

13. Гилл Ф. Практическая оптимизация / Ф. Гилл, У. Мюррей, М. Райт М.: Мир, 1985.-509 с.

14. Годунов С.К. Разностные схемы / С.К. Годунов, B.C. Рябенький М.: Наука, 1977.440 с.

15. Голуб Дж. Матричные вычисления / Дж. Голуб, Ч. Ван Лоун М.: Мир, 1999. - 548 с.

16. Голубев Г.В. Определение гидропроводности неоднородных нефтяных пластов нелокальными методами / Г.В. Голубев, П.Г. Данилаев, Г.Г. Тумашев -Казань: Изд-во Казанского университета, 1978. 168 с.

17. Данилаев П.Г. Коэффициентные обратные задачи для уравнений параболического типа и их приложения / П.Г. Данилаев. Казань: Изд-во Казанского математического общества, Изд-во УНИПРЕСС, 1998. - 127 с.

18. Дэннис Дж. Численные методы безусловной оптимизации и решения нелинейных уравнений / Дж. Дэннис, Р. Шнабель М.: Мир, 1988. - 440 с.

19. Елесин А.В. К решению обратной задачи по определению коэффициента фильтрации трехмерного напорного пласта / А.В. Елесин, А.Н. Габидуллина, А.Ш. Кадырова // Труды математического центра им. Н.И.Лобачевского. Казань, "Унипресс", 1998. - С. 103-105.

20. Елесин А.В. Идентификация коэффициента фильтрации неоднородного пласта в условиях напорной фильтрации жидкости: Дис. канд. физико-математических наук: 01.02.05., 2005. 138 с.

21. Елесин А.В. Учёт априорной сравнительной информации в задачах идентификации коэффициента фильтрации / А.В. Елесин, А.Ш. Кадырова // Вычислительные методы и программирование. 2008. - Т.9. - № 1. - С. 14-19.

22. Елесин А.В. Двухшаговые методы Левенберга-Марквардта в задаче идентификации коэффициента фильтрации / А.В. Елесин, А.Ш. Кадырова, П.А. Мазуров // Георесурсы. 2009. - 4(32). - С.40-42.

23. Зенкевич О. Конечные элементы и аппроксимация / О. Зенкевич, К. Морган -М.: Мир, 1986.- 318 с.

24. Измаилов А.Ф. Численные методы оптимизации / А.Ф. Измаилов, М.В. Со-лодов М.: ФИЗМАТЛИТ, 2008. - 320 с.

25. Ильин В.П. Методы неполной факторизации для решения алгебраических систем / В.П. Ильин. М.: Физматлит, 1995. - 288 с.

26. Коллинз Р. Течения жидкостей через пористые материалы / Р. Коллинз М.: Мир, 1964.-350 с.

27. Корн Г. Справочник по математике для научных работников и инженеров / Г. Корн, Т. Корн М.: Наука. Гл. ред. физ.-мат. лит., 1973. - 832 с.

28. Мазуров П.А. Запасы чувствительности в задачах идентификации коэффициента фильтрации трехмерных пластов / П.А. Мазуров, А.Н. Габидуллина, А.В. Елесин, А.Ш. Кадырова // Вычислительные методы и программирование.-2004.-Т.5. № 1.-С. 50-61.

29. Мазуров П.А. Квазиньютоновский двухшаговый метод минмизации функции невязки / П.А. Мазуров, А.В. Елесин, А.Ш. Кадырова // Вычислительные методы и программирование. 2009. - Т. 10. - № 1. — С. 64-71.

30. Мальцев А.И. Основы линейной алгебры / А.И. Мальцев М.: Наука. Гл. ред. физ.-мат. лит., 1970. - 400 с.

31. Марчук Г.И. Методы вычислительной математики / Г.И. Марчук М.: Наука, 1989. - 608 с.

32. Маскет М. Течение однородных жидкостей в пористой среде / М. Маскет -Гостоптехиздат, 1949. 628 с.

33. Мину М. Математическое программирование / М. Мину М.: Наука. Гл. ред. физ.-мат. лит., 1990.-488 с.

34. Мироненко В.А. Динамика подземных вод / В.А. Мироненко М. Изд-во МГГУ, 1996.-520 с.

35. Митчелл Э. Метод конечных элементов для уравнений с частными производными / Э. Митчелл, Р. Уэйт М.: Мир, 1981. - 216 с.

36. Морозов В.А. Алгоритмические основы методов решения некорректно поставленных задач / В.А. Морозов // Вычислительные методы и программирование.-2003.-Т.4. № 1.-с. 134-145.

37. Норри Д. Введение в метод конечных элементов / Д. Норри, Ж. де Фриз -М.: Мир, 1981.-304 с.43.0ртега Дж. Итерационные методы решения нелинейных систем уравнений со многими неизвестными/ Дж. Ортега, В. Рейнболдт -М.: Мир, 1975. -560 с.

38. Пантелеев А.В. Методы оптимизации в примерах и задачах / А.В. Пантелеев, Т.А. Летова М.: Высш. шк., 2005. - 544 с.

39. Парлетт Б. Симметричная проблема собственных значений. Численные методы / Б. Парлетт М.: Мир, 1983. - 384 с.

40. Полубаринова-Кочина П.Я. Теория движения грунтовых вод / П.Я. Полубаринова-Кочина. М.: Гос. Изд-во Технико-Теоретической литературы, 1952. - 676 с.

41. Рыжиков Ю.И. Вычислительные методы / Ю.И. Рыжиков СПб.: БХВ-Петербург, 2007. - 400 с.

42. Самарский А.А. Разностные методы для эллиптических уравнений / А.А. Самарский, В.Б. Андреев М.: Наука, 1976. - 352 с.

43. Самарский А.А. Численные методы / А.А. Самарский, А.В. Гулин М.: Наука, 1989.-432 с.

44. Самарский А.А. Численные методы решения обратных задач математической физики / А.А. Самарский, П.Н. Вабищевич. М.: Едиториал УРСС, 2004.-480 с.

45. Сухарев А.Г. Курс методов оптимизации / А.Г. Сухарев А.Г., А.В. Тимохов, В.В. Федоров-М.: Наука, 1986. 328 с.

46. Тихонов А.Н. Методы решения некорректных задач / А.Н. Тихонов, В.Я. Арсенин М.: Наука, 1986. - 288 с.

47. Тихонов А.Н. Численные методы решения некорректных задач / А.Н. Тихонов, А.В. Гончарский, В.В. Степанов, А.Г. Ягола. М.: Наука, 1990. - 232 с.

48. Федоренко Р.П. Введение в вычислительную физику / Р.П. Федоренко. М.: Изд-во Моск. физ.-техн. ин-та, 1994. - 528 с.

49. Флетчер К. Численные методы на основе метода Галёркина / К. Флетчер -М.: Мир, 1988.-352 с.

50. Формалев В.Ф. Численные методы / В.Ф. Формалев, Д.Л. Ревизников М.: ФИЗМАТЛИТ, 2006. - 400 с.

51. Форсайт Дж. Численное решение систем линейных алгебраических уравнений / Дж. Форсайт, К. Молер М.: Мир, 1969.

52. Хайруллин М.Х. О решении обратных коэффициентных задач фильтрации многослойных пластов методом регуляризации / М.Х. Хайруллин // ДАН РАН. 1996. Т.347. - №1. С.103-105.

53. Химмельблау Д. Прикладное нелинейное программирование / Д. Химмель-блау-М.: Мир, 1975. 534 с.

54. Черноруцкий И.Г. Методы оптимизации и принятия решений / И.Г. Черно-руцкий СПб.: Изд-во «Лань», 2001. - 384 с.

55. Carrera J. Estimation of aquifer parameters under transient and steady state conditions: Maximum likelihood method incorporating prior information / J. Carrera, S.P. Neuman // Water Resour. Res. 1986. - Vol. 22. - No. 2. - P. 199-210.

56. Carrera J. Estimation of aquifer parameters under transient and steady state conditions: 2. Uniqueness, Stability, and Solution Algorithms / J. Carrera, S.P. Neuman // Water Resour. Res. 1986. - Vol.22. - No.2. - P. 211-227.

57. Carrera J. Estimation of aquifer parameters under transient and steady state conditions: 3. Application to Synthetic and Field Data / J. Carrera, S.P. Neuman // Water Resour. Res. 1986. - Vol.22. - No.2. - P. 228-242.

58. Coats K.H. A new technique for determining reservoir description from field performance data / K.H. Coats, J.R. Dempsey, J.H. Henderson // Soc. Pet. Eng, J. -1970.- 10(1)-P. 66-74.

59. Emsellem Y. An automatic solution for the inverse problem / Y. Emsellem, G. de Marsily // Water Resour. Res. 1971. - Vol. 7. - No. 5. - P. 1264-1283.

60. Fasano G. A Truncated Nonmonotone Gauss-Newton Method for Large-Scale Nonlinear Least-Squares Problems / G. Fasano, F. Lampariello, M. Sciandrone // Computational Optimization and Applications. 2006. - 34. - P.343-358.

61. Fletcher R. A Modified Marquardt subroutine for nonlinear least-squares / R. Fletcher // Report R6799, Atomic Energy Research Establishment, England. -1971.

62. Gambolati G. A 3-D finite element conjugate gradient model of subsurface flow with automatic mesh generation / G. Gambolati, G. Pini, T. Tucciarelli // Adv. Water Resour. 1986. - 9. - P. 34-41.

63. Gambolati G. Numerical comparison of preconditionings for large sparse finite element problems / G. Gambolati, G. Pini, G. Zilli // Numerical Methods for Partial Differential Equations. John Wiley, New York, 1988. - P. 139-157.

64. Hadamard J. Le problem de Cauchy et les equations aux derivees partielles lin-eares hyperboliques / J. Hadamard. Paris, Hermann, 1932.

65. Hadamard J. Sur les problems aux derivees partielles et leur signification phisique / J. Hadamard. Bull. Univ. Princeton., 1902.

66. Hestenes M.R. Methods of conjugate gradients for solving linear systems / M.R. Hestenes, E. Stiefel // J. Res. Nat. Bur. Stand. 1952. - V. 49. - P. 409-436.

67. Hill M.C. Solving groundwater flow problems by conjugate-gradient methods and the strongly implicit procedure / M.S. Hill // Water Resour. Res. 1990. - Vol.26. - No.9. - P. 1961-1969.

68. Hill M.C. Methods and guidelines for effective model calibration. U.S Geological survey water-resources investigations report 98-4005. Denver, Colorado, 1998.

69. Jacquard P. Permeability distribution from field pressure data / P. Jacquard, C. Jain // Soc. Pet. Eng. J. 1965. - 5(4). - P. 281-294.

70. Larabi A. Solving three-dimensional hexahedral finite element groundwater models by preconditioned conjugate gradient methods / A. Larabi, F. De Smedt // Water Resour. Res. 1994. - Vol. 30. - No. 2. - P. 509-521.

71. Levenberg К. A method for solution of certain problems in least squares / K. Levenberg // Quart. Appl. Math. 1944. - 2. - P. 164-168.

72. Li B. Control volume function approximation methods and their applications to modeling porous media flow / B. Li, Z. Chen, G. Huan // Adv. Water Resour. 2003.-26.-P. 435-444.

73. Manteuffel T.A. An incomplete factorization technique for positive definite linear systems / T.A. Manteuffel // Math. Comput. 1980. - 34(150). - p. 473-497.

74. Marquardt D. An algorithm for least-squares estimation of nonlinear parameters / D. Marquardt // SIAM J. Appl. Math. 1963. - 11. - p. 431-441.

75. McLaughlin D. A reassessment of the groundwater inverse problem / D. McLaughlin, L.R. Townley // Water Resour. Res. 1996. - Vol. 32. - No. 5. - P. 1131-1161.

76. Meijerink J.A. An iterative solution method for linear systems of which the coeff-cient matrix is a symmetric M-matrix / J.A. Meijerink, H.A. Van der Vorst // Math. Comput.- 1977.-31.-P. 148-162.

77. More J.J. The Levenberg-Marquardt algorithm: implementation and theory / J.J. More // Lecture Notes in Math. 1977. -630. - p. 105-116.

78. More J.J. Testing Unconstrained Optimization Software / J.J. More, B.S. Garbow, K.E. Hillstrom // ACM transactions on Mathematical Software. 1981. - Vol. 7. -No.l. - P. 17-41.

79. Ortega J.M. Introduction to Parallel and Vector Solution of Linear Systems / J.M. Ortega. Plenum, New York, 1988. - 503 p.

80. Stallman R.W. Numerical analysis of regional water levels to define aquifer hydrology / R.W. Stallman // Eos. Trans. AGU. 1956. - 37(4). - P. 451-460.

81. Sun N.-Z. Inverse Problems in Groundwater Modeling / N.-Z. Sun. Kluwer Acad., Norwell, Mass., 1994. - 337 p.

82. Yakowitz S. Instability in aquifer identification: Theory and case study / S. Ya-kowitz, L. Duckstein // Water Resour. Res. 1980. - Vol. 16. - No. 6. - P. 10451064.

83. Yeh W. W-G. Review of parameter identification procedures in groundwater hydrology: The inverse problem / W.W-G. Yeh // Water Resour. Res. 1986. - Vol. 22.-No. 2.-P. 95-108.

84. Zhang J.Z. Nonmonotone Levenberg-Marquardt Algorithms and Their Convergence Analysis // J.Z. Zhang, L.H. Chen // Journal of Optimization Theory and Applications. 1997. - Vol. 92. -No. 2. - P. 393-418.