Двойственная регуляризация в обратной задаче УНЧ-зондирования земной коры Текст научной статьи по специальности «Математика»

РАДИОФИЗИКА

УДК 621.371+537.87

ДВОЙСТВЕННАЯ РЕГУЛЯРИЗАЦИЯ В ОБРАТНОЙ ЗАДАЧЕ УНЧ-ЗОНДИРОВАНИЯ ЗЕМНОЙ КОРЫ

© 2009 г. К.П. Гайкович 1, Ф.А. Кутерин 3, А.И. Смирнов 2, М.И. Сумин 3

1 Институт физики микроструктур РАН 2 Институт прикладной физики РАН 3 Нижегородский госуниверситет им. Н.И. Лобачевского

gai@ipm. sci-nnov.ru

Поступила в редакцию 11.11.2008

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

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

Введение

Низкочастотное электромагнитное зондирование земной коры (см., например, в [1-4]) осуществляется на частотах от тысячных долей герца до сотен герц с использованием как естественных источников (геомагнитная активность), так и когерентных сигналов, излучаемых антеннами [2, 3], роль которых могут, в частности, выполнять линии электропередачи. Измеренные на поверхности Земли частотные спектры электрического и магнитного полей являются исходной информацией для восстановления профиля геоэлектрического разреза (проводимости), который представляет большой интерес для геофизики и геологии. Соответствующая обратная задача чрезвычайно сложна. Впервые она была сформулирована А.Н. Тихоновым [1] и решалась для многослойной среды методом минимизации функционала невязки. В [4] ее удалось рассмотреть на классе непрерывных функций и свести к нелинейному интегральному уравнению 1-го рода. В [5] представлены результаты численного моделирования этой задачи на основе разработанного итерационного алгоритма с использованием теории Тихонова [6]. Этот метод, однако, с очевидностью ограничен условиями применимости теории возмущений.

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

Общая характеристика метода двойственной регуляризации

При решении данной обратной задачи, как и большинства обратных задач физической диагностики [4], считается естественным переходить от постановки исходной обратной задачи в терминах дифференциального уравнения (обыкновенного или в частных производных) к постановке, в которой основным объектом является интегральное уравнение. Как правило, это линейные уравнения первого рода типа Фред-гольма или Вольтерра, для решения которых имеется широкий арсенал эффективных вычислительных алгоритмов, например, таких как метод Тихонова, метод регуляризации на компактах и ряд других [6, 7].

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

Оказывается, однако, что этап перехода от дифференциальной формы обратной задачи к интегральным уравнениям вовсе не обязателен и может быть исключен из процесса решения обратной задачи, в том числе и в нелинейных случаях. Существуют методы, являющиеся составной частью общей теории методов оптимизации и оптимального управления, которые позволяют решать сложные обратные задачи непосредственно в рамках их естественных дифференциальных постановок. К таким методам можно отнести, в частности, методы, основанные на минимизации нелинейного функционала невязки (см., например, [12]). К тому же классу относится и рассматриваемый в данной работе метод, основанный на теории двойственности, или, другими словами, метод двойственной регуляризации [8-11], устойчивый к ошибкам исходных данных. Его особенностью является то, что, во-первых, он развивает классический в теории оптимизации лагранжев подход, во-вторых, служит удобной универсальной основой для решения различных обратных задач физической диагностики и, в-третьих, позволяет решать такие обратные задачи в их естественных дифференциальных постановках.

Далее излагаются результаты исследования возможностей метода двойственной регуляризации [11] для решения нелинейной обратной задачи геомагнитного зондирования в ее естественной дифференциальной постановке. В качестве общей характеристики рассматриваемого подхода следует отметить, что двойственные алгоритмы являются одними из самых эффективных при решении оптимизационных задач с ограничениями [13, 14], в которых непосредственным образом используется классическая идея «снятия» ограничений, заложенная в принципе Лагранжа. Отметим также, что характерной особенностью любого алгоритма, основанного на теории двойственности, является параллельное решение сразу двух задач - исходной и двойственной к ней. При этом процесс решения двойственной задачи, которая всегда является задачей выпуклой оптимизации (с точ-

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

Обратная задача УНЧ геомагнитного зондирования

Пусть электромагнитная монохроматическая волна с векторами электрического и магнитного поля, направленными вдоль осей х и у соответственно (Е = Ехх0, И = Иуу0), распространяется нормально к земной поверхности г = 0 и проводимость ст(г) зависит от глубины (в дальнейшем у проекций Ех и Иу опускаем индексы х и у). Для интересующих нас низкочастотных полей земную кору можно считать проводником с чисто мнимой комплексной диэлектриче-

, ■ „ ■4па

ской проницаемостью є = є — і є ~ —і---------.

а

Уравнения для комплексных амплитуд электрического и магнитного поля в гауссовой системе единиц имеют вид:

d2Е 4па(г)а „ _ ТТ с dE

—т — і-----2^— Е = 0, И = і-------,

dz с а dz

го < г < 0 , (1)

Е = Е1 + іЕ2 , И = И1 + іИ 2 , где с - скорость света, а - циклическая частота, г0 < 0 - нижняя граница интервала анализа, в

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

Е(г = 0, а) = Е0(а),

И(г = 0, а) = И0(а).

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

Подчеркнем и другую особенность постановки обратной задачи УНЧ-зондирования в данной работе. Мы считаем величину г0 конечной и задаем поля Е (г = г0, а), И (г = г0, а) на этой нижней границе, считая их дополнительными вспомогательными неизвестными, подлежащими определению (нахождению) наряду с проводимостью о(г). Данный прием введения вспомогательных дополнительных неизвестных оказывается удобным с математической точки зре-

ния, так как он основывается на использовании классов функций, задаваемых лишь на конечном интервале [z0,0], что облегчает обоснование сходимости метода регуляризации. Отметим также, что рассматриваемая задача линейна по этим вспомогательным переменным. Поэтому введение дополнительных регуляризирую-щих слагаемых, обусловленных такими переменными, при формулировке оптимизационной задачи, которая соответствует рассматриваемой обратной задаче, лишь усиливает, на наш взгляд, общий регуляризирующий эффект. Отметим здесь же, что изучение исходной обратной задачи УНЧ-зондирования на бесконечном полуинтервале [-»,0] приводит формально к иной постановке оптимизационной задачи, в которой все равно приходится определять некий конечный интервал, за пределами которого измерения уже становятся неинформативными. В силу названных причин такая постановка в данной работе не рассматривается.

Введем обозначения:

Е1 = Хр

оЁ1 аЁ2

■ = х, Е = х, —2 = х..

, -„-2= Х3, ^ = Х4. (3)

йч йч

Тогда исходная обратная задача УНЧ-зондирования для восстановления подповерхностного профиля проводимости земной коры ст(ч) в новых переменных примет вид эквивалентной задачи минимизации

!00, Хт ) = ||ст|1 2+|\ХЛ ^ min,

/1(CT, Хт ) = Х0(®),^е D (4)

D = е ^(Ч0,0) :0 ^ СТ(4) ^ СТ0}

D с L2(zo,0), 40 <0, Хт е R4

где сг0 - достаточно большое положительное

число, а = 4п/с2. Существенно отметить, что в рамках данного рассмотрения компоненты электромагнитного поля хт на нижнем уровне

интервала, на котором восстанавливается профиль проводимости, также рассматриваются как неизвестные величины, которые определяются из решения задачи.

Можно доказать, что эта задача имеет решение. Значение минимизируемого функционала

10 на этом решении обозначим через /0 . В то же время нельзя гарантировать единственность оптимального элемента из-за нелинейности в ограничении этой оптимизационной задачи (нелинейность типа произведения управления на траекторию).

В случае линейных обратных задач, которые сводятся к эквивалентным линейно выпуклым оптимизационным задачам, для их решения на основе двойственной регуляризации достаточно конструкции классической функции Лагранжа [8-10]. В нелинейных случаях необходимо использовать конструкции модифицированных функций Лагранжа, т.е. обычных функций Лагранжа с добавленными к ним штрафными слагаемыми того или иного вида. Рассмотрим здесь один из возможных случаев [11] модифицированной функции Лагранжа задачи (4)

LЛa, Хт,Л) НМ12 + 11 Хт ||2 +

®2

+ Г (Л(ю), х[а, хт ](ч = 0,®) - х0(®))й® +

+ мЦ І I х[ст,хт](г = 0,а) — Х0(а)| dа + (6)

Ох

— = Л0( г)а) х, х( г0) = хт,

аг

¡1(0, хт )(а) = х[а, хт ](г = 0) =

=х[ста,хт](г = 0) = х0(а), а є [а1;а2].

( 0 1 0 0^

Л(ст( г)а) =

0

0

0 —аа(г)а 0

0

0

аа(г)а 0 0

х0(а) =

( Е Е01

—И02(а)а/ с Е

Е02

V И01(а)а/ с

1

(5) + а I х[0,хт](г = 0,а) — х0(а)|2 аа) =

НИР + 11 Хт ||2 + (Л ¡1(0 Хт ) — х0) +

+ М(|1 Л0, х т ) — х0|1 + 11 Л0, хт ) — х0||2).

В соответствии с результатами [11] запишем регуляризованную по Тихонову с параметром регуляризации а модифицированную двойственную задачу, представляющую собой задачу максимизации вогнутого функционала на гильбертовом пространстве Х^(а1,а2):

Ум(Л) — а || Л ||2^ тах,

ЛєЛм ={Лє ¿2(а1,а2): 11 Л|1 <м}

^(Л) = тІП ^м0, хт ,Л),

0єD ,fєR4

Л є Еі2(а1,а2), м > 0.

0

L42(o1,o2) =

=L2(o1,o2) х L2(®j,®2) х х L2(®j,®2) х L2(®j,®2).

Минимум модифицированного функционала Лагранжа Lu(a,xm,Л) в силу его специфики достигается на некотором элементе (ац[Х\,х [Л\), причем такой элемент, вообще

говоря, не единственный. Обозначим множество всех таких элементов Ии[Л\.

Функционал Vи является вогнутым на

L4,(®j,®2), но из-за нелинейности задачи (4) он не является дифференцируемым в обычном смысле. Однако замечательным его свойством является то, что для него может быть выписана явная формула для супердифференциала в каждой точке Ле L4,(®j,®2). При этом под супердифференциалом вогнутого функционала V понимается взятый с обратным знаком субдифференциал в смысле выпуклого анализа функционала — V. Именно явное выражение для супердифференциала и позволяет организовать вычислительную устойчивую к ошибкам исходных данных процедуру решения обратной задачи.

Можно показать (см., например, [8, 10]), что супердифференциал функционала V задается

формулой

dVuW= CO{x[CTu(Л), Хти (Л)\( Z = 0) —

-Хо :[°и(Л),Хти(Л)\ е ЕДЛ = co{Gu[Л]}, (8)

где coA - выпуклая оболочка множества А. Для вычисления (8) нужно решить задачу минимизации модифицированного функционала Лагранжа

L (ст, х ,Л) ^ min, <7 е D,x е R4.

U У ? m5 ' 5 5 т

Приведем выражения для его градиентов по ст и по хт, которые служат основой для нахождения такого решения (см., например, [14\)

daLu0, Хт,Л) =

= 2ст + J (Л(о),Т[ст\(о) :

О

х А(1 • о)х[ст\(о))0 +

(

+ и

1

— + 2

11 ¡1(ст, Хт ) — Х0 у

Т[ст\(z = 0,о)А(1о) х х х[ст, хт\( z = 0,o))do, д L (ст, х ,Л) =

Хт и 4 т * '

: 2хт + f Т [ст\(z = 0, a>a(o)do +

J о

( 1

+ ui

1 11(CT, Хт ) — Х0

-+ 2

/• От

х f Т [ст\(z = 0,o)edo,

J О

(^1[ст\( z,o) ^ ^2[СТ\( z,0) ^з[ст\( z,o) ^4[СТ\( z,0)

¥< + А ((у( ч))ю)^ = 0,

(0) = е,, i = 1,2,3,4, где е., 7 = 1, 2, 3, 4 - единичные орты, * -

знак транспонирования.

Пусть исходные данные Е0®), #0®) заданы с ошибкой. Будем поэтому считать, что

I х^ — х° II2 =

1 х0 х0 11 _

От

= f (Хд(о) — х0(о))2 do<S2.

Ja

(10)

f (х[ст, Хт\(z = 0,о) — Х0О),

J О

На основе приведенных выше явных выражений для градиентов дУ , д L , д L можно

Л Хт Л у /Л

предложить несколько вариантов регуляризи-рующих алгоритмов для решения задачи (4). Остановимся в данной работе на одном из них, представляющем собой итерационный процесс, сходимость которого доказана в [11]. Этот процесс, в соответствии с [11], протекает двояким образом в зависимости от того, существует или нет решение классической двойственной к (4) задачи:

А. Решение X классической (при /л = 0) двойственной задачи (7)

У0(Л) ^ тах, Ле L2(ю1,ю2) (11)

существует, причем У0 (Л* ) = /0 ;

Б. Решение классической двойственной задачи не существует, или в случае его существования К0(Л*) < /0*.

Указанный итерационный процесс, который мы рассмотрим сначала для случая неограниченно убывающей ошибки исходных данных 8 = 8к, к = 1,2, к, состоит в следующем. Берем

х

e = e1 + e2 + e3 + e4

х

х

начальную точку Л0 є Л . В случае А штрафной коэффициент м - произвольная достаточно большая положительная величина. В случае Б его надо стремить к бесконечности согласованным с ошибкой исходных данных образом. Пусть последовательность Л, к = 1,2, к, конструируется по правилу

Лк+1 = РгЛ (Лк + рЪ* — 2вкакЛк),

М

к = 0,1,2,...,

Ък = 1^к (0 м ’к, х м ’к) — х0к,

м 1 ^ ’т' 0’

(12)

р(&м, ъм [Л ]) ,

где РгЛ - оператор проектирования на множество А, р (а, А) - расстояние от точки а до множества А,

(ум ,к, хЛ,к )-

^ (У М,8к,ак,7(М,к) ХМ,8к,ак,7( л,к))

- такой элемент минимизирующей модифицированный функционал Лагранжа £8 (у, хт ,Лк)

/ _и8к ,ак ,7 л 8к ,ак ,7 ч последовательности (У ■ ’ ’ , хт ), г =

= 1, 2, к, что

^¡1 (у, Хт ,Лк )-

Нт0к +а +вк ) = 0,^+г < С0,

к ^да а

| ак+1 —ак |_ вк

Ііт . кч 2 п к--Ііт —=

к(ак)2вк к ^да (ак)

^к тк _

= 1іт аз=0, 1іт а = 0,

к ^ да (а ) к ^ да а

ак = к 1/3, /Зк = к 1/2. В случае А, т.е. при произвольном фиксированном достаточно большом , последовательность (уЛ,к,хЛ,к) ,к = 1,2,...,

должна сходиться, в соответствии с [11], к одному из решений исходной обратной задачи. В случае же Б для того чтобы имела место та же сходимость, необходимо еще стремить к + » и штрафной коэффициент = к так, что выполняется предельное соотношение 8к Л 0, к ^ да.

Обсудим в заключение важнейший для решения исходной обратной задачи случай, когда ошибка задания исходных данных не стремится

к нулю и остается конечной и равной 8 <81. Именно такая ситуация и реализуется в конкретных практических обратных задачах. Пусть в этом случае итерационный процесс задается в случае А рекуррентным соотношением

Лк+1 = РгЛ (Лк + вкОк — 2вкакЛк),

М

к = 0,1,2,...,

(15)

2>к ___ т5 м,к .у.м,к ) х^

т ' 0 5

- 4к (0мк,х^к,Лк)> (13)

м,5к ,ак ,і( м, к )

—є ' Уст є Б, х є Я4,

5 5 т 5

м,5к ,ак ,і(мк) . п г . є^ > > ^ 0, к ^ да.

Здесь последовательности 8к, ак Тк,

к = 0,1,2,..., удовлетворяют условиям согласования

5к > 0,ак > 0,вк >0,

(14)

^ ак вк = +да.

к=1

Легко проверить, что такие последовательности ак, рк, к = 1,2, к, существуют. В качестве одного из возможных примеров таких последовательностей можно взять, в частности,

р(&ц,G¡[Лk]) <тк, а в случае Б - аналогичным соотношением, но с заменой л на Лк • Пусть также итерации в случае А продолжаются до такого наибольшего номера к = к(8), при котором выполняются

неравенства 8к >8, к = 1,2,..., к(8).Тогда в соответствии с [11] оказывается, что любая предельная точка последовательности (ул,к(8(*:»,хЛ ,к(8(5))), 5 = 1,2,..., для любой неотрицательной последовательности 8(5) ^ 0, 5 ^ да является решением исходной обратной задачи. Аналогичное регуляризирующее правило останова может быть сформулировано и в случае Б (см. [11]).

Заключение

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

к

ветствует таким дополнительным данным. Такой способ введения дополнительной априорной информации зачастую приводит к радикальному улучшению качества решения некорректных задач (см. примеры в [4]).

Представленные результаты получены по программе фундаментальных исследований Отделения физических наук РАН (проекты «Когерентная томография земной коры с использованием ультранизкочастотных электромагнитных полей» и «Ближнепольное электромагнитное зондирование диэлектрических и проводящих сред»). Работа выполнена также при финансовой поддержке РФФИ (коды проектов 0802-00117, 07-01-00495) и аналитической целевой ведомственной программы «Развитие научного потенциала высшей школы (2009-2010 гг.)» Минобрнауки РФ (регистрационный номер 2.1.1/3927).

Список литеретуры

1. Тихонов А.Н. // Доклады АН СССР. Нов. сер. 1950. Т.73, № 2. С. 295-297.

2. Велихов Е.П. // Доклады РАН. 1994. Т. 338, № 1. С. 106-109.

3. Поляков С.В. // Геомагнетизм и аэрономия.

2003. Т. 42, № 2. С. 240-248.

6. Тихонов А.Н., Гончарский А.В., Степанов В.В., Ягола А.Г. Регуляризирующие алгоритмы и априорная информация. М.: Наука, 1983. 200 с.

7. Бакушинский А.Б., Гончарский А.В. Некорректные задачи. Численные методы и приложения. М.: Изд-во Моск. ун-та, 1989. 632 с.

8. Сумин М.И. // Журн. вычисл. матем. и матем. физ. 2004. Т. 44, № 11. С. 2001-2019.

9. Сумин М.И. // Вестник Нижегородского университета. Серия Математика. 2004. Вып. 1(2). С. 192-208.

10. Сумин М.И. // Журн. вычисл. матем. и матем. физ. 2007. Т. 47, № 4. С. 602-625.

11. Сумин М.И. // Журн. вычисл. матем. и матем. физ. 2007. Т. 47, № 5. С. 794-814.

12. Короткий А.И., Цепелев И.А. // Труды Института математики и механики УрО РАН. Т. 9. 2. Динамика жидкости и газа. Екатеринбург: Изд-во УрО РАН, 2003. С. 73-86.

13. Мину М. Математическое программирование. Теория и алгоритмы. М.: Наука, 1990. 488 с.

14. Васильев Ф.П. Методы оптимизации. М.: Факториал Пресс, 2002. 824 с.

DUAL REGULARIZATION IN THE INVERSE PROBLEM OF EARTHS CRUST ULF SOUNDING

K.P. Gaikovich, F.A. Kuterin, A.I. Smirnov, M.I. Sumin

A dual regularization method to retrieve subsurface Earth’s crust conductivity profiles by ultra-low frequency (ULF) spectral measurement data has been proposed. In this method, the corresponding inverse problem is solved in its initial differential statement on the basis of the theory of optimization and optimal control without the reduction of Maxwell’s equations to non-linear integral equations, as it was done before.

4. Gaikovich K.P. Inverse Problems in Physical Diagnostics. New York: Nova Science Publishers Inc.,

2004. 372 p.

5. Гайкович К.П., Смирнов А.И. // Вестник Нижегородского университета им. Н.И. Лобачевского. Н. Новгород: Изд-во ННГУ, 2007. № 2. С. 73-80.