Двухшаговые методы Левенберга-Марквардта в задаче идентификации коэффициента фильтрации Текст научной статьи по специальности «Математика»

УДК: 519.6

А.В. Елесин, А.Ш. Кадырова, П.А. Мазуров

Институт механики и машиностроения Казанского научного центра РАН, Казань

mazurov@mail.knc.ru

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

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

Ключевые слова: математическое моделирование, минимизация функции невязки, идентификация коэффициента фильтрации.

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

д_

дх

К

дh) д

— н------

дх J ду

К

** ду J dz

К.

Эй

dz

= 0

(1)

где Ку, Кг - коэффициенты фильтрации, И = И(х,у,2) - напор. Пласт Щпятислойный («40 кмх30 кмх200м), слои зо-

71

нально-неоднородные, □ = 1М . Каждая зона однородно-

к=1

сти Щ, характеризуется двумя значениями коэффициента фильтрации Кф, К*к. Значения коэффициента фильтрации Кф брались в пределах от 0.1 м-сут-1 до 100 м-сут-1, К^ в пределах от 0.0001 м-сут-1 до 0.02 м-сут-1. На кровле пласта заданы граничные условия 2-го рода (от -9.2-10-3 м-сут-1 до 2-10-3 м-сут-1), которые моделируют инфильтрацию, расходные характеристики родников и реки. Подошва и боковая поверхность непроницаемы, за исключением участка боковой поверхности 5-го слоя, на котором заданы граничные условия 1-го рода И = 80 м. Для дискретизации уравнения (1) использовался метод конечных элементов. Полученная в результате дискретизации система линейных алгебраических уравнений решалась методом

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

^тах ДІпХ"0 xyz ^шах дь.*:" xyz it nc

ЛМ 9.27 1.7 ь 1-Н X <Г) 0.71 222** 31934

ДЛМ 9x10 7 0.17 108 15645

ДЛММ 8х10“7 0.16 57 8437

Тaбл. 1. EM - мemoд Лeвeнбepгa-Mapквapдma; ДЛM - deyxrna-гoвыlй мemoд Лeвeнбepгa-Mapквapдma, ДЛMM - мoдuфuцu-poeaurnrn двyxшaгoвыlй мemoд Лeвeнбepгa-Mapквapдma; hmax

г» A I'It

мaкcuмaльнaя нeвязкa,

' xyz ’ xyz

nananbHbie u

umoгoвыle знaчeнuя cpeднeквaдpamuчecкoгo omклoнeнuя, it -чucлo umepaijm, nc - чucлo, xapaкmepuзyющee выlчucлumeль-Hbie 3ampambi, ** - npoijecc мuнuмuзaцuu ocmaнoвлeн no щ>u-mepuю мeдлeннoй cxoдuмocmu.

сопряженных градиентов с предобусловливающей матрицей в виде неполного разложения Холесского (Hill, 1990).

Из решения уравнения (1) при заданных значениях коэффициента фильтрации Кф, определялись значения

напора в наблюдательных точках hj (/'=1,_,247). Наблюда-

тельные точки располагались на кровле пласта.

По значениям h* = hj + S, где 8 - задаваемая погрешность, восстанавливались логарифмы значений коэффициента фильтрации К = {К$ы =§аКхукМКл^=х нимума функции невязки

из ми-

J{K) = -RTR =—V(ft - h* )2 2 2£V ; J’’

(2)

где R = (hl- -й^7) - вектор невязок, h. = h.(K) -

вычисленные значения напора в наблюдательных точках. В процессе идентификации параметры К^к, считались неизвестными.

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

Kn = Kn-1 - (H + MnE)-1g,

где E - единичная матрица, H = ATA - приближённая мат-

Г dhi 1

рица вторых производных, А = г - матрица чувстви-

тельности, g - градиент функции невязки, ц - параметр Марквардта, n - номер итерации. При больших значениях m направление спуска метода Левенберга-Марквардта близко к направлению метода наискорейшего спуска, при малых значениях ц близко к направлению метода Гаусса-Ньютона. Различные варианты метода Левенберга-Марк-вардта отличаются стратегией выбора на каждой итерации параметра Марквардта (Дэннис, 1988; Hill, 1998; Sun, 1994). В статье (Мазуров и др., 2009) построен двухшаговый метод на основе метода Левенберга-Марквардта, в

котором параметр Марквардта определялся методом золотого сечения из минимума функции невязки

тшДГ^-СЯ+^Г^). В предлагаемом двухша-

Мп

говом методе Левенберга-Марквардта (ДЛМ) используется следующая стратегия выбора параметра Марквардта (Пантелеев, Летова, 2005). Начальное значение параметра Марквардта выбирается на порядок больше максимального сингулярного числа матрицы Н. В случае уменьшения функции невязки на текущей итерации 3(Кп) < 3(Кп'г) параметр Марквардта уменьшается в два раза, в случае нарушения условия убывания параметр Марквардта увеличивается в два раза до тех пор, пока это условие не выполнится. При построении двухшагового алгоритма используется главная система координат, полученная с помощью сингулярного разложения И=У£ V т, где V - ортогональная матрица, £ =(о'1,...,о'142) - диагональная матрица, а. - сингулярные числа, а1 > с2...> о142> 0. Направления вдоль осей с большими сингулярными числами соответствуют спуску ко дну оврага, а с малыми сингулярными числами - смещению вдоль дна оврага. Каждая итерация двухшагового метода Левенберга-Марквардта проводится в два шага. На первом шаге допускается возрастание функции невязки за счет её роста в направлениях, соответствующих большим сингулярным числам. На втором шаге проводится уменьшение функции невязки вдоль этих направлений (спуск ко дну оврага). В итоге по итерациям строится убывающая последовательность значений функции невязки. Также рассмотрена модификация двухшагового метода Левенберга-Марквардта (ДЛММ), в которой на втором шаге проводятся дополнительные спуски в направлениях, соответствующих большим сингулярным числам. Использование второго шага в двухшаговых методах ДЛМ и ДЛММ позволило сократить число итераций и вычислительные затраты процесса минимизации функции невязки за счет большего уменьшения функции невязки на каждой итерации по сравнению с классическим методом Левенберга-Марквардта (ЛМ). Приведем алгоритм метода ДЛММ. Каждая итерация состоит из последовательности следующих операций:

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

3^.) = 3(Кп1),

где Кп1 = Кп-1 + 8 1, 8 1 = - (Н+ цп£)-1я. При выполнении условия уменьшения функции невязки Зх(рп) < 3(КпЛ)новые значения переменных определяются, как Кп = Кп-1+81, параметр Марквардта уменьшается в два раза, и итерация заканчивается, иначе - переходим к пункту 2.

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

3/^ = 3(К^), г = 2,3,...

где Кп-‘ = Кп’‘-‘ + Й?, й% =У7у, Яу. =§у./(ст1 +м„), г = !,...,#, Зу. =0, г = д+1,...,142, £у. - компоненты вектора §у =УТ§, £ = А7/?, /? - вектор невязок в точке Кп- '-1 , # - число направлений в главной системе координат, вдоль которых проводится смещение параметров, # выбирается из условия о > цп > <г+1. Значения Jt(^t) вычисляются

до тех пор, пока выполняется условие

-3/^ > °.°Ч/мЛ * ^ 2. (3)

Определим

^п = JtXmn)},

где г* - значение г, при котором условие (3) нарушается. Если / < /(К”'1), то новые значения параметров берутся для 1, соответствующего / , параметр Марквардта уменьшается в два раза, и итерация заканчивается. В противном случае параметр Марквардта увеличивается в два раза, и возвращаемся к пункту 1.

Элементы матрицы чувствительности вычисляются методом прямого дифференцирования уравнения фильтрации (1) и соответствующих граничных условий, градиент функции невязки вычисляется по формуле g = ЛТЯ.

Для остановки процесса минимизации использовались два критерия:

1) медленная сходимость итерационного процесса /(К”'1) - /(К”) < 0.01/(Ки'1) в течение трёх итераций,

2) достижение заданной точности по напору в наблюдательных точках = тах^. (Кп)-к* | < КТ6

Численные результаты. Значения максимальной невязки И , среднеквадратического отклонения коэффициента фильтрации от истинных значений

г (Г ^ хук

+

.*=1

0пК‘-1п*1)2)/142

1/2

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

5 3 метод АЫК1 дьа:" хуг и пс

0.1 ЛМ 1.7 1.16 58** 8321

ДЛМ 1.08 64** 9239

ДЛММ 1.06 46** 6735

0.01 ЛМ 1.7 0.76 207** 29774

ДЛМ 0.52 98** 14193

ДЛММ 0.52 63** 9424

-0.1 ЛМ 1.7 0.98 146** 20992

ДЛМ 0.91 74** 10701

ДЛММ 0.92 51** 7487

-0.01 ЛМ 1.7 0.77 202** 29054

ДЛМ 0.56 100** 11481

ДЛММ 0.55 60** 8855

Табл. 2. Уел. обозн. см. табл. 1.

для определения элементов матрицы чувствительности.

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

Задача идентификации коэффициента фильтрации является обратной задачей и относится к классу некорректно поставленных задач (Тихонов, 1986; Бип, 1994). При наличии погрешностей в замерах напора значения идентифицируемых параметров, начиная с некоторой итерации, обычно удаляются от своих истинных значений, при этом функция невязки продолжает уменьшаться. Для выбора номера итерации с итоговыми значениями коэффициента фильтрации применяются специальные правила останова процесса минимизации или прерывания полученной последовательности значений функции невязки. Эти правила являются одним из регуляризирующих элементов решения обратных задач.

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

1) определялся номер итерации к, с которого начинается медленная сходимость процесса минимизации;

2) определялся максимальный номер г = 1,2,..., при котором выполняется условие / к' < 1.5/ к;

3) итоговые значения коэффициента фильтрации брались с итерации с номером к - г.

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

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

Литература

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

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

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

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

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

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

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

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

A.V. Elesin, A.Sh. Kadyrova, P.A. Mazurov. The two-step Levenberg-Marquardt methods in hydraulic conductivity identification task.

The two-step Levenberg-Marquardt methods are proposed for the minimization of a residual function. The proposed methods and the classical Levenberg-Marquardt method are compared on solving model problem of hydraulic conductivity identification for a threedimensional anisotropic confined aquifer. The numerical results are shown the efficiency of the two-step Levenberg-Marquardt methods.

Keywords: mathematical modeling, minimization of residual function, inverse problem.

Казань: «Изд-во ПЛУТОН». 2007. 124 с.

Гидродинамические исследования и моделирование многоствольных горизонтальных скважин

Иктисанов В.А.

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

область применения многоствольных горизонтальных скважин. Рассмотрены вопросы управления выработкой

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

К.ф.-м.н., старший научный сотрудник.

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

Научный сотрудник.

Петр Алексеевич Мазуров

К.ф.-м.н., ведущий научный сотрудник, и.о. зав. лабораторией.

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

420111, Россия, Казань, ул. Лобачевского, 2/31. Тел.: (843) 292-74-90.