Идентификация коэффициента фильтрации неоднородного пласта в условиях напорной фильтрации жидкости тема автореферата и диссертации по механике, 01.02.05 ВАК РФ

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

Елесин Андрей Викторович

ИДЕНТИФИКАЦИЯ КОЭФФИЦИЕНТА ФИЛЬТРАЦИИ НЕОДНОРОДНОГО ПЛАСТА В УСЛОВИЯХ НАПОРНОЙ ФИЛЬТРАЦИИ ЖИДКОСТИ

01.02.05 - механика жидкости, газа и плазмы

АВТОРЕФЕРАТ

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

'КАЗАНЬ - 2005

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

Научный руководитель: кандидат физико-математических наук

П.А. Мазуров

Официальные оппоненты: доктор физико-математических наук,

профессор А.В. Лапин

доктор физико-математических наук, профессор Н.Д. Якимов

Ведущая организация: Казанский государственный технический

университет им. А.Н. Туполева

Защита диссертации состоится 24 февраля 2005 г. в 14 часов 30 минут в аудитории Физ.2 на заседании диссертационного Совета Д212.081.11 при Казанском государственном университете по адресу: 420008, г. Казань, ул. Кремлевская, 18.

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

Автореферат разослан "_"_2005 г.

Ученый секретарь диссертационного совета, кандидат физико-математических наук

А.А. Саченков

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

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

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

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

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

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

1. Введено понятие запасов чувствительности для анализа процесса идентификации коэффициента фильтрации.

2. Построены с использованием запасов чувствительности новый квазиградиентный алгоритм минимизации функции невязки и новый критерий прерывания итерационного процесса идентификации.

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

Апробация. Основные результаты диссертационной работы докладывались и обсуждались на: международной конференции "Математические модели и численные методы механики сплошных сред" (Новосибирск, 1996), XI международной конференции по вычислительной механике и современным прикладным программным системам (ВМСППС2001, Москва-Истра, 2001), IV научно-практической конференции молодых ученых и специалистов Республики Татарстан (Казань, 2001), международной конференции 4th International Conference "Calibration and Reliability in

Groundwater Modelling" (Prague, 2002), VIII Четаевской международной конференции "Аналитическая механика, устойчивость и управление движением" (Казань, 2002), Всероссийской конференции "Высокопроизводительные вычисления и технологии (ВВТ-2003)" (Ижевск, 2003), XVII сессии Международной школы по моделям механики сплошной среды (Казань, 2004), итоговых научных конференциях КазНЦ РАН, семинарах Института механики и машиностроения КазНЦ РАН.

Публикации. По теме диссертации опубликовано 13 работ, список которых приведен в конце автореферата.

Структура и объем работы. Диссертация состоит из введения, четырех разделов, заключения и списка литературы. Работа изложена на 138 страницах, содержит 76 рисунков и 12 таблиц. Список литературы насчитывает 91 наименование.

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

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

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

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

д(„ dh\ Э/Л df^dfi) п , . _ „

при граничных условиях 1-го и 2-го рода, где h = h(x,y,z) - напор, Kiy = Kjy(x,y,z), Kz = Kz(x,y,z) - коэффициенты фильтрации, Q - область решения

задачи. Коэффициент фильтрации ищется в классе кусочно-постоянных функций

Kv = ¿Л^Ф^*,**), Kz = f,Kzl^k(x,y,z) (1.2)

*=i *=i

где число зон однородности, и

m

зоны однородности, неизвестные значения коэффициента фильтрации в

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

где г = (/z, -ht ,...,hN - h"N) - вектор невязки, Ir* = , Л™^}^, -

соответственно заданные и вычисленные значения напо_ра в наблюдательных точках, hj = hj(K), К= {К,}ff, = {inЛГ^, InАГ^(, М= 2т - число идентифицируемых параметров, N - число наблюдательных точек. При известной области решения £2 постановка задачи идентификации для напорного пласта сохраняется и в случае напорно-безнапорной фильтрации.

Определены модельные задачи идентификации коэффициента фильтрации пятислойного напорного пласта от-

ражающие особенности задачи идентификации коэффициента фильтрации водоносного пласта речного водосборного бассейна. На кровле пласта задавались граничные условия 2-го рода (в пределах от -9.2x10° м/сут до 2х10"3 м/сут). Подошва и боковая поверхность пласта непроницаемы за исключением участков боковой поверхности нижнего слоя, на которых задавались граничные условия 1-го рода. В каждой зоне однородности задавались значения коэффициента фильтрации (в пределах от 0.1 м/суг до

100м/сут), К^ь (в пределах от 0.0001 м/сут до 0.09 м/сут) для вычисления напоров h* в наблюдательных точках. Далее значения опреде-

лялись по значениям напора ftj = h* без погрешностей и по значениям напора h*=h'J +6j с погрешностяб^иВ процессе идентификации как погрешности, так и значения коэффициента фильтрации К^, К£, считались неизвестными. Начальные значения коэффициента фильтрации К^, определялись в процессе минимизации функции невязки при условии КфЬ =К®, k = 1.....т. Модельные задачи идентификации реша-

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

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

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

относительно параметров) ^ = §га<1У(А")=|^-| и матрицы чувстви-

Гэй/,

тельности А = у^^Г• Для вычисления производных при минимизации

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

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

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

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

(2.1)

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

= (2.2) назовем запасом чувствительности 1-го параметра. Запас чувствительности Р1 характеризует приращение квадрата вектора невязки за счет г-го параметра от текущей точки Къ пространстве параметров до точки К-в предполагаемого минимума функции невязки ./(¿Г-5)= 0 при условии близости указанных точек в пространстве параметров. Показано, что если матрица А квадратная (Ы = М) и невырожденная, то полный запас чувствительности (сумма запасов чувствительности всех параметров) равен хвадрату вектора невязки

м

гтг

независимо от близости точек К и К-5 в пространстве параметров. В случае И> М (число наблюдательных точек больше числа параметров) показано, что сумма запасов чувствительности всех параметров не больше квадрата вектора невязки.

В диссертации используется распределение запасов чувствительности по параметрам в главной системе координат, получаемой вУО - разложением матрицы Н = и£ит, где V - ортогональная матрица, Е = с11од(а1,...|Од/) - диагональная матрица, с^ >ст2 ^ — -ам >0 - сингулярные числа. Система уравнений (2.1) в главной системе координат запишется в виде

где йц. - компоненты вектора отклонений Яц = 1/тх, - чувствительности функции J в направлении осей главной системы координат (компоненты градиента =Uтg). Запас чувствительности в направлении ьой оси преобразованной системы координат является неотрицательной величиной

(2.4)

Полный запас чувствительности является инвариантной величиной относительно данного преобразования матрицы Н

В п. 2.2 рассмотрен алгоритм построения вУО - разложения приближенной матрицы вторых производных функции невязки Н = ААТ для получения распределения запасов чувствительности по главным направлениям в пространстве параметров. Алгоритм включает две части. В первой части матрица Н приводится к трехдиагональному виду с помощью преобразований Хаусхолдера. Во второй части проводится итерационный процесс, основанный на QR - итерациях с одинарным сдвигом. Числа обусловленности матрицы Н для модельных задач при значениях параметров К^ = К? лежат в пределах от »7x10® до «6x10й.

В п. 2.3 распределение запасов чувствительности используется для объяснения медленной сходимости метода наискорейшего спуска. На рис. 2.1-2.4 приведены распределения запасов чувствительности Ри. и абсолютных значений чувствительностей для начального и конечного состояний функции невязки при минимизации методом наискорейшего спуска. Рисунки приведены для модельной задачи идентификации 142-х значений коэффициента фильтрации без погрешностей в замерах напора. Результаты для других модельных задач имеют аналогичный вид.

Рис.2.1. Метод наискорейшего спуска. Чувствительности начального состояния.

Рис.2.2. Метод наискорейшего спуска. Запасы чувствительности начального состояния.

Рис. 2.3. Метод наискорейшего спуска. Чувствительности конечного состояния.

Рис. 2.4. Метод наискорейшего спуска. Запасы чувствительности конечного состояния.

В начальном состоянии (рис. 2.1, 2.2) видно соответствие между значениями чувствительностей и запасов чувствительности. Сравнение чувстви-тельностей и запасов чувствительности конечного состояния (рис. 2.3, 2.4) показывает, что максимальной чувствительности в направлении первой оси соответствует малое значение запаса чувствительно-

сти PVx = 0.025 м2 (максимальное значение = 0.67 м2). Направление

минимизации в методе наискорейшего спуска для конечного состояния данной задачи определяется в основном чувствительностью в направлении первой оси (рис. 2.3). Небольшой запас чувствительности в этом направлении (рис. 2.4) приводит к малому изменению функции невязки и медленной сходимости итерационного процесса.

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

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

=

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

1 Т /<.\ 1 О . . - -

- метод наискорейшего спуска

- метод сопряженных градиентов

- квазиградиентный метод

в п й-! —I

О 50 100 150 215

Рис. 2.5. Сравнение квазиградиентного метода с методами наискорейшего спуска и сопряженных градиентов: (а) - функция невязкй|((р)еднеквадратиче-ское отклонение - номер итерации.

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

Третий раздел посвящен исследованию устойчивости задачи идентификации.

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

методом Левенберга-Марквардта и квазиградиентным алгоритмом при

идентификации 142-х значений коэффициента фильтрации с погрешностями в замерах напора 5 j = 0.01м, 5у =0.1 м.

19 i,,„fO

1.5 -

метод Левенберга-Марквардта квазиградиентный метод

I ~ I-1-~—I п 0.85 -|-Г-1-1 щ

0 20 40 60 87 0 30 60 90 119

Рис. 3.1. Средцеквадратическое отклонение Д1пК^-- (а) погрешности 5^=0.1 м, (б) погрешности = 0.01 м.

Результаты решения модельных задач идентификации при различных погрешностях показали, что методы наискорейшего спуска и сопряженных градиентов имеют устойчивый характер, но значения идентифицируемых параметров, как и в задаче без погрешностей, далеки от истинных. Квазиградиентный метод обладает большей устойчивостью к погрешностям по сравнению с методом Левенберга-Марквардга. В методе Л^венберг<1-Марква{.'дта параметры, начиная с некоторой итерации, значительно удаляются от истинных значений.

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

Обозначим ÜT°° = lim К" как решение задачи идентификации при

Л—»oo

значениях напора как множество всех решений данной задачи

при различных значениях напора А*. Решение К" является предельной точкой итерационного процесса в котором В

пространстве параметров в точке Кт не существует направления, вдоль которого функция невязки убывает. Множество |ä"" } можно разделить на два класса |if" } и {лг^}. Для решен йТ^и з первого класса функция невязки равна нулю Первый класс решений, кроме решения

соответствующего напорам h* —hf, включает решения + дК„., соответствующие напорам h* = hj + , Н а п о р^-о^р еделяются в наблюдательных точках решением прямой задачи при возмущенных параметрах Кц.+ЪК,,., т.е. jiT,,. + 8Jf„. je {лТ^1}. Основным свойством решений К„. + ЬК„. является стремление к нулю отклонений Sy. —* О при К" ->ЛГ„. + 5JST„. в случае невырожденной матрицы Н = ААГ. Для решений Kf из второго класса функция невязки сходится к некоторой конечной величине -И ИГ?) = .1" П Г[ля решений из второго класса при N=Mпока-

па

зано, что ds = J^Sm —>ао в случае maxgy. —* 0 при и—юо. Аналогичные результаты для решений из второго кл <"-•-<' "^.место в случае N>M при условии Рц —> Ру > 0. Величина ds = в окрестности

решения К,г +5К1г характеризует расстояние между текущей точкой К" в пространстве параметров и точкой К„ +5#Г„. предполагаемого минимума 3{К№ +5.ЙГ(Г) = 0. Увеличение величи^ыв окрестности минимума функции невязки при итерационном процессе свидетельствует об удалении параметров от решений {Кw + ЬК# } из первого класса, что лежит в основе предлагаемого критерия прерывания процесса идентификации. Номер итерации прерывания Ияр выбирался с наименьшим значением ds из множества номеров определяющих по некоторым правилам окрестность конечного состояния процесса минимизации. Графики изменения по итерациям значений без прерывания и с прерыванием, полученных методом

Левенберга-Марквардта при идентификации 142-х значений коэффициента фильтрации с погрешностями 5j — -0.1 м, приведены на рис. 3.2. Близкие результаты получены при решении других модельных задач.

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

0 10 и.,»25 40 50 0 10 Яар = 25 40 50 0 10 20 л„ = 25

Рис. 3.2. Метод Левенберга - Марквардтапри погрешности 5у=-0.1 м: (а) - среднеквадратическое отклонение А1пК^ без прерывания, (б) - функция ¿я (в) - среднеквадратическое отклонение Д1п.ЙГ^г «с прерыванием, и5(р = 25 - номер итерации с минимальной величиной «/, среди множества номеров 5= {23-50}.

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

В п. 4.1 показано, что решение задачи идентификации по наблюдательным точкам, расположенным внутри пласта, более неустойчиво по сравнению с решением, полученным по наблюдательным точкам, расположенным на кровле пласта. Во всех задачах, рассмотренных в разделах 1, 2, 3 диссертации, имелись наблюдательные точки внутри пласта. На рис. 4.1 приведены результаты решения задачи идентификации 142-х значений коэффициента фильтрации методом Левенберга-Марквардта при погрешностях 5 , = 0.1 м, 5 , = -0.1 м для различных множеств наблюдательных точек Д|42, о247, -©389 -

• (142 наблюдательные точки внутри пласта, в каждой зоне однородности по две наблюдательные точки) ■ - />247 (247 наблюдательных точек на кровле пласта)

4 Л А - />389 = ¿>,42 О Х>247.

0 20 40 69 0 50 100 176

Рис. 4.1. Среднеквадратические отклонения ДЬК^, полученные методом Левенберга - Марквардта по наблюдательным точкам Х>|42. ^247, Ада : (а) погрешности погрешности

На рис. 4.1 видно, что добавление дополнительных наблюдательных точек к наблюдательным точкам приводит к появлению неустой-

чивости итерационного процесса идентификации коэффициента фильтрации. Для объяснения полученных результатов вводится запас чувствительности /-Й наблюдательной точки P¿, равный запасу чувствительности функции невязки в случае ]-й невязки, равной единице, и остальных невязок, равных нулю. На рис. 4.2 приведены распределения запасов чувствительности наблюдательных точек начального состояния в задаче идентификации коэффициента фильтрации по наблюдательным точкам ■®389- Распределения запасов чувствительности наблюдательных точек конечного состояния имеют аналогичный вид.

1 20 40 60 >0 100 120 142 ' 1 20 40 60 80 100 120 142 1

867 1.33 5 51Е-2 1.66Е-6 867 1.33 5.31Е-2 1.66Е-6 а,

Рис. 4.2. Запасы чувствительности наблюдательных точек Dm и Du7 начального состояния задачи идентификации 142-х значений коэффициента фильтрации по наблюдательным точкам Du?, (а) наблюдательные точки D\n, (б) наблюдательные точки D247.. Ниже порядковых номеров направлений приведены значения сингулярных чисел

Запасы чувствительности наблюдательных точек D142 сосредоточены в направлениях с меньшими сингулярными числами по сравнению с запасами чувствительности наблюдательных точек Дг47- Наблюдательные точки Лиг вносят большую неустойчивость в процесс идентификации по сравнению с наблюдательными точками

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

которого будем определять коэффициент фильтрации. При погрешности S

Пусть в точке х = х* е(0,£) известно значение напора h'(x'), по значению

в замере напора имеем

где К и ЛГ8 - соответственно значения коэффициента фильтрации, полученные по значениям напора и И (х )-6. Из (4.2) видно, что чем дальше точка х' от граничного условия 1-го рода и чем больше расход д, тем меньше ошибка в определении коэффициента фильтрации. Уравнение (4.1) и условие (4.2) отражают сущность неустойчивости итерационного процесса идентификации коэффициента фильтрации по наблюдательным точкам внутри пласта для приведенных модельных задач. При приближении наблюдательной точки к граничному условию 1-го рода задача вырождается, и даже при небольшой погрешности в замере напора решение задачи может не существовать из-за нарушения принципа максимума, что показано в диссертации на примере идентификации коэффициента фильтрации одномерной задачи.

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

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

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

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

2. Построены с использованием запасов чувствительности новый квазиградиентный алгоритм минимизации функции невязки и новый критерий прерывания итерационного процесса идентификации.

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

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

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

1. Елесин А.В. Идентификация параметров пласта и априорная информация / А.В. Елесин, П.А. Мазуров // Тезисы докладов международной конференции "Математические модели и численные методы механики сплошных сред". - Новосибирск, 1996. - С. 268-269.

2. Елесин А.В. Идентификация коэффициента фильтрации (напорный режим, трехмерный случай) / А.В. Елесин // Тезисы докладов международной научно - технической конференции "Механика машиностроения". - Набережные Челны, 1997. - С. 29.

3. Елесин А.В. К идентификации параметров водоносного пласта / А.В. Елесин, О.В. Коляскина, П.А. Мазуров, В.А. Ходина // The VI-th International Symposium on Application of Mathematical Methods and Computers in Mining, Geology and Metallurgy. Proceedings Volume: The Application of Mathematical Methods and Computers in Hydrogeology. -Prague, 1997. - P. 1-6.

4. Елесин А.В. К идентификации параметров водоносного пласта / А.В. Елесин, П.А. Мазуров // Тр. I международной конференции "Модели механики сплошной среды, вычислительные технологии и автоматизированное проектирование в авиа- и машиностроении". - Казань, 1997.-Т.2.-С. 60-64.

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

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

7. Elesin A.V. Use of minimization along the slope for estimation of aquifer parameters / P.A. Mazurov, A.V. Elesin, A.N. Gabidullina, A.Sh. Kadyirova // Proceedings ofthe 4th International Conference on Calibration and Reliability in Groundwater Modelling. - Prague, Czech Republic, 2002. - V. 1. - P. 278-281.

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

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

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

П.Елесин А.В. Квазиградиентный метод минимизации функции невязки при решении обратных коэффициентных задач / П.А. Мазуров. А.Н. Габидуллина, А.В. Елесин, А.Ш. Кадырова // Тезисы Всероссийской конференции "Высокопроизводительные вычисления и технологии (ВВТ-2003)". - Ижевск, 2003. - С. 95-96.

12.Елесин А.В. К минимизации функции невязки квазиградиентным методом при идентификации коэффициента фильтрации трехмерного анизотропного пласта / А.В. Елесин, П.А. Мазуров // Математическое моделирование. - 2004. - Т. 16. - № 8. - С. 99-113.

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

Отпечатано полиграфическим комплексом физического факультета

КГУ

Заказ №01-21-01/05 _бумага офсетная, тираж 12 0 экз._

г. Казань, ул. Кремлевская, дом 16-А, к. 010, тел. (8432) 36-90-16

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

Введение.

Раздел 1. Задача идентификации коэффициента фильтрации трехмерного неоднородного пласта по значениям напора в наблюдательных точках.

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

1.2. Численное решение уравнения стационарной фильтрации жидкости.

1.3. Техника вычисления производных в задачах минимизации функции невязки.

1.4. Идентификация коэффициента фильтрации методом наискорейшего спуска.

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

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

2.2. Сингулярное разложение приближенной матрицы вторых производных для получения распределения запасов чувствительности.

2.3. Использование распределения запасов чувствительности для анализа медленной сходимости метода наискорейшего спуска.

2.4. Построение квазиградиентного метода минимизации функции невязки для определения коэффициента фильтрации.

Раздел 3. Исследование устойчивости задачи идентификации коэффициента фильтрации.

3.1. Сравнение устойчивости квазиградиентного метода с методами наискорейшего спуска, сопряженных градиентов и Левенберга -Марквардта к погрешностям в замерах напора.

3.2. Критерий прерывания итерационного процесса идентификации.

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

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

4.2. Идентификация коэффициента фильтрации по наблюдательным точкам, расположенным на кровле пласта.

При решении различных гидрогеологических задач (управление водными ресурсами, прогнозирование распространения загрязнений и др.) необходимо знание коэффициента фильтрации К пласта. Существуют различные подходы к определению коэффициента фильтрации [10]: (1) по геофизическим данным, (2) по отобранному керну, (3) локальные гидродинамические методы, (4) нелокальные гидродинамические методы. В данной работе будут рассматриваться нелокальные гидродинамические методы, в которых определяется сразу все поле коэффициента фильтрации. Одним из основных способов определения коэффициента фильтрации является минимизация функции невязки J, являющейся суммой квадратов разностей между измеренными значениями напора h* = \h* ], характеризующими состояние гидрогеологического объекта, и значениями напора h = {hj(K))IJl> вычисленными с использованием математической модели,

J = J(K) = yTr, где r = {h\-- h*NJ - вектор невязки. Задачи такого типа относятся к классу некорректно поставленных задач [1-3, 5, 10, И, 42-45, 49, 56, 66, 67, 85]. По определению задача называется некорректно поставленной, если нарушается одно из условий:

1) решение задачи существует,

2) решение задачи единственно,

3) решение непрерывно зависит от исходных данных.

В задачах идентификации гидрогеологических параметров любое из условий 1) - 3) может оказаться не выполненным. Соответствующие примеры можно найти в [56, 77, 85, 88]. Характерным свойством минимизации функции невязки при наличии погрешностей является удаление параметров от своих истинных значений, начиная с некоторой итерации, хотя функция невязки продолжает уменьшаться. В основе методов решения некорректных задач лежит понятие регуляризирующего алгоритма, устойчивого к погрешностям в исходных данных [1-3, 12, 25, 44, 45, 49]. Прерывание итерационного процесса идентификации является одним из регуляризирующих элементов решения таких задач [1-3, 36, 44]. При известных погрешностях задачи критерий прерывания обычно согласуется с погрешностями [1-3, 36, 44].

Задача определения коэффициента фильтрации водоносного пласта по замерам напора была впервые рассмотрена в работе [54]. Первые попытки использования вычислительной техники для решения данной задачи предприняты в [84]. Замеренные значения напора в [84] интерполировались на все узлы конечно-разностной сетки. Поле коэффициента фильтрации также определялось во всех узлах сетки по полученным значениям напора в этих узлах. Такое решение задачи оказалось неустойчивым. Для устранения неустойчивости предлагалось разбить пласт на зоны с постоянными значениями коэффициента фильтрации (уменьшить число неизвестных параметров), и использовать метод наименьших квадратов для определения значений коэффициента фильтрации в этих зонах.

Существует два подхода к уменьшению числа идентифицируемых параметров [91]: (1) метод разбиения на зоны и (2) интерполяционный метод. В первом подходе область решения задачи разбивается на зоны, каждая из которых характеризуется постоянным значением параметра [59, 62, 84]. Неизвестный параметр аппроксимируется кусочно-постоянной функцией с числом неизвестных значений, равным числу зон. Во втором подходе могут использоваться различные методы: конечноэлементное представление параметра [61, 90], аппроксимация параметра с помощью сплайна [83, 89], полиномиальный метод [65], kriging [58].

Различные методы решения задачи идентификации можно разделить на две группы [80, 85, 91]: (1) явные методы, (2) неявные методы. В явных методах

10, 11, 85] значения параметров определяются из решения нелинейной системы уравнений, в которой замеренные значения напора являются известными значениями коэффициентов этой системы. В результате при определении параметров явными методами не требуется решение прямой задачи. В неявных методах [5, 49, 85] оценка неизвестных параметров итерационно улучшается так, чтобы значения напора, полученные при решении прямой задачи, совпадали с известными замерами напора. В этом случае требуется многократное решение прямой задачи с различными значениями идентифицируемых параметров. Один из основных неявных методов заключается в минимизации функции невязки. Алгоритмы, использующиеся при минимизации функции невязки, можно разбить на три группы [86]: методы прямого поиска, градиентные методы, различные модификации метода Гаусса-Ньютона.

В алгоритмах прямого поиска процесс минимизации функции строится только по значениям минимизируемой функции, полученным при различных значениях идентифицируемых параметров. Различные алгоритмы методов прямого поиска можно найти в [4, 6, 26, 27, 43]. Как правило, методы прямого поиска обладают низкой скоростью сходимости и редко используются в задачах идентификации.

При построении градиентных методов на каждой итерации необходимо вычислять производные минимизируемой функции по отношению к искомым параметрам. Широко используемыми в задачах идентификации являются методы наискорейшего спуска и сопряженных градиентов [2-6, 26, 27, 43]. В [56] исследуется зависимость сходимости трех различных методов сопряженных градиентов (метод Флетчера-Ривса, метод Бройдена, метод Флетчера-Пауэлла-Дэвидона) от выбора начальных значений неизвестных параметров. Для улучшения сходимости предлагается использовать комбинированный метод: на первых итерациях минимизация проводится методом Флетчера-Ривса, а на последующих итерациях продолжается методом Бройдена.

В основе различных модификаций метода Гаусса-Ньютона лежит аппрокdh Т симация Н = АА матрицы Гессе функции невязки, где А =

8kj матрица чувствительности. В базовом методе Гаусса-Ньютона неизвестные параметры на каждой итерации изменяются по правилу

Kn+l=K"-s, где s = Н lg - вектор отклонений Гаусса-Ньютона, g - градиент функции невязки (вектор чувствительности, характеризующий чувствительность функции невязки к параметрам), п - номер итерации. Базовый метод Гаусса-Ньютона в задачах идентификации параметров пласта, как правило, не используется по ряду причин [85]:

• значение функции невязки при переходе на новые значения параметров может увеличиться;

• матрица Н может оказаться вырожденной или плохо обусловленной, и соответственно вектор отклонений Гаусса-Ньютона не может быть вычислен или вычислен с большой погрешностью;

• новые значения параметров могут принять недопустимые значения. Различные модификации метода Гаусса-Ньютона позволяют устранить указанные выше трудности [13, 85]. Одной из модификаций метода Гаусса-Ньютона, широко используемой в задачах идентификации, является метод Левенберга -Марквардта. Вектор отклонений в методе Левенберга - Марквардта в зависимости от поведения функции невязки приближается либо к направлению вектора чувствительности функции невязки, либо к вектору отклонений Гаусса-Ньютона. Существуют различные варианты метода Левенберга-Марквардта [13]. В [87] метод Левенберга - Марквардта модифицируется на основе информации о сингулярном разложении приближенной матрицы вторых производных: в системе координат, полученной SVD - разложением, отклонения, соответствующие малым сингулярным числам, зануляются. На сегодняшний день задача оптимального выбора направлений, в которых зануляются отклонения, остается до конца нерешенной [13].

Для вычисления компонент вектора чувствительности g и матрицы чувствительности А можно использовать три различных метода [85, 91]: (1) метод конечно-разностных соотношений, (2) метод прямого дифференцирования, (3) вариационный метод. В случае метода конечно-разностных соотношений и метода прямого дифференцирования для вычисления g и А необходимо решить М+1 систем уравнений, аналогичных системе уравнений, полученной при дискретизации уравнения фильтрации, где М - число идентифицируемых параметров. В случае вариационного метода вектор чувствительности g может быть вычислен в результате решения всего лишь двух систем уравнений, а для вычисления матрицы чувствительности А необходимо решить ЛН-1 систем уравнений, где N - число замеров напора. Метод конечно-разностных соотношений [52, 53] редко используется в задачах идентификации параметров из-за сложности вычисления производных нужной точности. Метод прямого дифференцирования и вариационный метод обладают более высокой точностью. При вычислении матрицы чувствительности А в целях уменьшения объема вычислений при N<M следует использовать вариационный метод, при M<N - метод прямого дифференцирования. В задачах идентификации параметров вариационный метод впервые использовался в [70], общая схема построения вариационного метода приведена в [46].

При решении задачи идентификации неявными методами на каждой итерации минимизации функции невязки требуется многократное численное решение прямой задачи. При этом уравнение фильтрации аппроксимируется системой линейных алгебраических уравнений, в основном получаемых методами конечных разностей, конечных элементов или конечных объемов [8, 21, 32, 37, 40-42, 47, 50, 60, 63]. Матрица коэффициентов полученной системы уравнений симметрична и сильно разрежена. Решение системы уравнений можно получить прямыми или итерационными методами. На практике при большом количестве неизвестных обычно используются итерационные методы. Одним из наиболее эффективных итерационных методов решения больших, сильно разреженных систем уравнений является метод сопряженных градиентов с предо-буславливающей матрицей [69, 71, 72, 78, 79]. К преимуществам данного метода относится отсутствие задаваемых итерационных параметров. Скорость сходимости метода сопряженных градиентов в значительной степени зависит от выбора предобуславливающей матрицы. Как правило, улучшение эффективности предобуславливающей матрицы ведет к увеличению памяти, необходимой для численной реализации метода [78]. В [79] рассматривается метод сопряженных градиентов с предобуславливающими матрицами, полученными при помощи диагонального масштабирования и полиномиального метода. В [69] при решении линейных и нелинейных задач фильтрации проводится сравнение метода сопряженных градиентов с тремя различными предобуславливающими матрицами, построенными методами неполного разложения Холесского, модифицированного неполного разложения Холесского и полиномиальным методом. В [72] для построения предобуславливающей матрицы используются диагональное масштабирование, неполное разложение Холесского, неполная факторизация, модифицированная неполная факторизация. При построении различных предобуславливающих матриц во многих случаях необходимо, чтобы матрица системы была с диагональным преобладанием [51, 71, 72]. Матрица коэффициентов системы уравнений, полученная при аппроксимации уравнения фильтрации, в общем случае не обладает свойством диагонального преобладания. Для решения таких задач в [72] предлагается строить предобуславливаю-щую матрицу по измененной матрице системы: все положительные внедиаго-нальные элементы добавляются к диагональному элементу и после этого зану-ляются. Другие подходы к построению предобуславливающей матрицы можно также найти в [64, 74, 81]. Различные алгоритмы решения разреженных систем алгебраических уравнений методами неполной факторизации, по своему смыслу близкими к методам с предобуславливающей матрицей, можно найти в [22].

В данной работе решены модельные задачи идентификации коэффициента фильтрации неоднородного напорного пласта, отражающие особенности задачи идентификации коэффициента фильтрации водоносного пласта речного водосборного бассейна. Пласт трехмерный реальной конфигурации («40 км х 30 км х 200 м), пятислойный, слои зональнонеоднородные. Предполагается, что зоны имеют слоистую структуру с горизонтальным характером напластования прослоек различной проницаемости. Поэтому каждая зона однородности Qk характеризуется двумя значениями коэффициента фильтрации К^, KZk, постоянными в пределах зоны однородности. Значение К^ обычно на несколько порядков больше значения Kzk. Значения коэффициента фильтрации определялись по результатам минимизации функции невязки. На кровле пласта заданы граничные условия 2-го рода, в реальном пласте определяемые инфильтрацион-ным питанием, расходными характеристиками реки и родников. Подошва и боковая поверхность пласта непроницаемы, за исключением участков боковой поверхности нижнего слоя с заданными граничными условиями первого рода.

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

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

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

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

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

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

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

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

Диссертация состоит из введения, четырех разделов, заключения и списка литературы. Работа содержит: страниц -138, рисунков - 76, таблиц - 12, список литературы - 91 наименование.

Выводы.

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

Заключение.

1. В диссертации проведено исследование возможностей гидродинамического метода определения коэффициента фильтрации и поля напоров водоносного пласта речного водосборного бассейна по замерам напора в наблюдательных точках пласта. Поставлены и решены модельные задачи идентификации коэффициента фильтрации пятислойного напорного пласта, отражающие особенности задачи идентификации коэффициента фильтрации водоносного пласта речного водосборного бассейна. С использованием метода Левенберга-Марквардта для минимизации функции невязки проведена идентификация коэффициента фильтрации (от 20-ти до 142-х значений) в модельных задачах с различным числом и расположением зон однородности по наблюдательным точкам, часть которых расположена внутри пласта в зонах однородности. Задача идентификации 20-ти параметров без погрешностей в замерах напора имела несколько решений (площадь зоны однородности да 440 км ). В задачах идентификации 38-ми и 142-х параметров получена единственность решений. При погрешностях в замерах напора задача идентификации 142-х параметров (площадь зоны однородности да 70 км ) имела более устойчивый характер по сравнению с задачей идентификации 38-ми параметров (площадь зоны однородности да 270 км). Для всех решений, полученных методом Левенберга-Марквардта, характерным является неустойчивость к погрешностям замеров напора в наблюдательных точках. На первых итерациях параметры приближаются к истинным значениям, но, начиная с некоторой итерации, начинают от них удаляться. Показана недостаточность минимизирующих свойств методов наискорейшего спуска и сопряженных градиентов для идентификации коэффициента фильтрации. Методы наискорейшего спуска и сопряженных градиентов имеют устойчивый характер, но значения коэффициента фильтрации, полученные этими методами, далеки от истинных значений.

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

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

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

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

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

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

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

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

6. Васильев Ф.П. Численные методы решения экстремальных задач / Ф.П. Васильев. -М.: Наука, 1988. 552 с.

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

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

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

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

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

12. Денисов A.M. Введение в теорию обратных задач / A.M. Денисов. М.: Изд-во МГУ, 1994. - 208 с.

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

14. Елесин А.В. Идентификация коэффициента фильтрации (напорный режим, трехмерный случай) / А.В. Елесин // Тезисы докладов международной научно технической конференции "Механика машиностроения". - Набережные Челны, 1997. - С 29.

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

16. Елесин А.В. Идентификация параметров пласта и априорная информация / А.В. Елесин, П.А. Мазуров // Тезисы докладов международной конференции "Математические модели и численные методы механики сплошных сред". -Новосибирск, 1996. С. 268-269.

17. Елесин А.В. К минимизации функции невязки квазиградиентным методом при идентификации коэффициента фильтрации трехмерного анизотропного пласта / А.В. Елесин, П.А. Мазуров // Математическое моделирование. -2004. Т.16. - № 8. - С. 99-113.

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

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

20. Калиткин Н.Н. Численные методы / Н.Н. Калиткин. М.: Наука, 1978. -512 с.

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

22. Лаврентьев М.М Некорректные задачи математической физики и анализа / М.М. Лаврентьев, В.Г. Романов, С.П. Шишатский. М.: Наука, 1989. - 286 с.

23. Лесин В.В. Основы методов оптимизации / В.В. Лесин, Ю.П. Лисовец. М.: Изд-во МАИ, 1998. - 344 с.

24. Летова Т.А. Экстремум функций в примерах и задачах / Т.А. Летова, А.В. Пантелеев. М.: Изд-во МАИ, 1988. - 376 с.

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

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

27. Маскет М. Течение однородных жидкостей в пористой среде / М. Маскет. -М.: ГОСТОПТЕХИЗДАТ, 1949. 627 с.

28. Мироненко В.А. Динамика подземных вод / В.А. Мироненко. М.: Изд-во Московского государственного горного университета, 1996. - 519 с.

29. Мишина А. П. Высшая алгебра / А. П. Мишина, И. В. Проскуряков. М.: Наука, 1965. - 300 с.

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

31. Норри Д. Введение в метод конечных элементов / Д. Норри, Ж. де Фриз. -М.: Мир, 1981.-304 с.

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

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

34. Самарский А.А. Разностные методы для эллиптических уравнений /

35. A.А. Самарский, В.Б. Андреев. М.: Наука, 1976. - 352 с.

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

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

38. Сухарев А.Г. Курс методов оптимизации / А.Г. Сухарев, А.В. Тимохов,

39. B.В. Федоров. -М.: Наука, Гл. ред. физ.-мат. лит., 1986. 328 с.

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

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

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

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

44. Фридман А. Уравнения с частными производными параболического типа /

45. A. Фридман. М.: Мир, 1968. - 427 с.

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

47. Шайдуров В.В. Многосеточные методы конечных элементов /

48. B.В. Шайдуров. -М.: Наука, 1989. 288 с.

49. Axelsson О. Finite Element Solution of Boundary Value Problems / O. Axelsson, V.A. Barker. Theory and Computation, Academic, San Diego, Calif., 1984.432 p.

50. Bard Y. Nonlinear parameter estimation / Y. Bard. John Wiley, New York, 1974.

51. Becker L. Identification of parameters in unsteady open-channel flows / L. Becker, W. W-G.Yeh // Water Resour. Res. 1972.- Vol. 8. - No.4. - P. 956965.

52. Bennett R.R. Geology and groundwater resources of the Baltimore area / R.R. Bennett, R.R. Meyer. Mines and Water Resour. Bull. 4, Maryland Dep. of Geol. - 1952.

53. 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.

54. 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.

55. Chavent G. Identification of distributed parameter system: About the output least square method, its implementation, and identifiability / G. Chavent // Identification and System Parameter Estimation. Pergamon, New York, 1979. -Vol. 1. - P. 85-97.

56. Clifton P.M. Effects of kriging and inverse modeling on conditional simulation of the Avra Valley aquifer in Southern Arizona / P.M. Clifton, S.P. Newman // Water Resour. Res. 1982. - Vol. 18. - No. 4. - P. 1215-1234.

57. 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.

58. DiStefano N. An identification approach to subsurface hydrological systems / N. DiStefano, A. Rath // Water Resour. Res. 1975. - Vol. 11. - No. 6. - P. 10051012.

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. 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.

61. 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.

62. Garay H.L. Distributed parameter identification of groundwater system by nonlinear estimation / H.L. Garay, Y.Y. Haimes, P. Das // J. Hydrol. 1976. - 30. -P. 47-61.

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

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

65. 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.

66. Hill M.S. 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.

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

68. Kershaw D.S. The incomplete Cholesky conjugate gradient method for the iterative solution of systems of linear equations / D.S. Kershaw // J. Compute. Phys. -1978.-No. 26.-P. 43-65.

69. 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.

70. Levenberg K. A method for the solution of a certain nonlinear problems in least squares / K. Levenberg // Q. Appl. Math. 1944. - No. 2. - P 164-168.

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

72. Marquardt D.W. A algorithm for least squares estimation of nonlinear parameters / D.W. Marquardt // SIAM, J. 1963. - 11. - P. 431-441.

73. 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.

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

75. Neuman S.P. Calibration of distributed parameter groundwater flow models viewed as a multiple-objective decision process under uncertainty / S.P. Neuman // Water Resour. Res. 1973. - Vol. 9. - No. 4. - P. 1006-1021.

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

77. Sagar B. A direct method for the identification of the parameters of dynamic non-homogeneous aquifers / B. Sagar, S. Yakowitz, L. Duckstein // Water Resour. Res. 1975. - Vol. 11. - No. 4. - P. 563-570.

78. 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.

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

80. Sun N.-Z. Coupled inverse problem in groundwater modeling, 1, Sensitiviti analysis and parameter identification / N.-Z. Sun, W.W-G. Yeh // Water Resour. Res. 1990. - Vol. 26. - No. 10. - P. 2507-2525.

81. Weiss R. Parameter space methods in joint estimation for groundwater flow models / R. Weiss, L. Smith // Water Resour. Res. 1998. - Vol. 34. - No. 4. - P. 647661.

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

83. Yakowitz S. On the identification of inhomogeneous parameters in dynamic linear partial differential equations / S. Yakowitz, P. Noren // J. Math. Anal. Appl. -1976.-53.-P. 521-538.

84. Yeh W.W.-G. Parameter identification with optimum dimension in parameterization / W.W.-G. Yeh, Y.S. Yoon // Water Resour. Res. 1981. - Vol.17. - No. 3. -P. 664-672.

85. 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.