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

РАДИОФИЗИКА

УДК 537.87+621.371

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

К.П. Гайкович 1, А.И. Смирнов 2

© 2007 г.

1 Институт физики микроструктур РАН 2 Институт прикладной физики РАН

gai@ipm.sci-nnov.ru

Поступклв воедвкцкю 9.02.2007

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

Введение

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

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

невязки Тихонова [5] решается интегральное уравнение Фредгольма 1-го рода, и результаты численного моделирования на основе этого алгоритма.

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

Методы решения задач магнитотеллурического зондирования

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

(Е0 = Е0х ■ Х0, Н0 = Н0у ■ у0), котоРая распространяется нормально к земной поверхности г = 0. Другое существенное упрощение связано с тем, что на низких частотах измерений земную кору можно считать проводником с чисто мнимой комплексной диэлектрической прони-

цаемостью є = є -is" x-i-

4na

где проводи-

а

мость о^) зависит от глубины. При этом уравнения для комплексных амплитуд электромагнитного поля в гауссовой системе единиц будут иметь вид:

d2 E.

dz2

, — i 40° E = o,

с2 x

H = i C dEx y a dz

Ho(0) = ^ і і H(a,z")dz"

Ее отклонение от реально измеренной зависимости АН0(а) свидетельствует о неоднородности подповерхностного профиля проводимости. По значениям полей на поверхности из решения прямой задачи (хорошо известного для однородной среды и получаемого численно для неоднородной) определяются и их глубинные распределения Е(с,г), Н(ю^).

В первом приближении, которое по существу является борновским, ядро интегрального уравнения (6) можно найти, считая среду однородной:

(1)

(2)

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

Электрические и магнитные поля на поверхности г = 0 считаются известными (экспериментально измеренными) функциями частоты с :

Е (г = 0,с) = Е0(с), (3)

Н (2 = 0,с) = Н0(с). (4)

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

Н(с, 2) = ]я(с, 2" ^2"

е о, z) = Eo(а) exP^o+i о)

H (0, z) = C(i—1) E (0, z) = ао

= Ho(0)exp[-| +iOJ}

(7)

(8)

где o =

с

толщина скин-слоя.

ст

a(z )dz . (5) =

Если в (5) положить г = 0, получаем интегральное уравнение, связывающее измеренную на поверхности компоненту поля Н0 = Н(с,0) и профиль о(г):

a(z')dz'. (6)

Аналогичные (5) и (6) соотношения можно получить и для электрического поля.

Ядром уравнения (6) служит интеграл от неизвестного распределения поля H (z, 0) и, на первый взгляд, задача определения профиля a(z) кажется неразрешимой. Тем не менее, возможен содержательный подход к её решению. Так, если на поверхности измерена частотная зависимость, например, электрического поля E0(a), то в случае однородной среды a(z) = a0 = = const ей соответствует вполне определенная частотная зависимость магнитного поля H0(a).

В качестве о0 естественно выбрать поверхностное значение проводимости, которое полагается известным. Для однородной среды частотная зависимость одного из полей (в данном случае - электрического) на поверхности навязывается источником поля. При решении задачи на основе реальных данных Е0(с) определяется экспериментально. Вообще говоря, излишне измерять Е0(с) и Н0(а>) по отдельности; достаточно иметь информацию о спектре поверхностного комплексного импеданса среды X = = Ех/Ну^ = 0,с).

Такой подход позволяет свести задачу к двум уравнениям типа Фредгольма 1-го рода (для действительной и мнимой частей комплексной амплитуды магнитного поля И1 '?) и решать их итерационным методом относительно о(г):

0

АН,?’1 (с) =\ К?л (с, г')Аог+1 (г ")йг", (9)

где ядро К пересчитывается на (/+1)-м шаге итерационного процесса по восстановленному профилю о(г). Регуляризированное решение некорректной задачи (9) получается на каждом шаге методом обобщённой невязки [5]. В этом методе решение ищется на классе функций с квадратично интегрируемой производной; невязка минимизируется до уровня погрешности

ад

ад

0

z

—со

—со

ад

z

—ад

данных, мера которой ЪН задается в интегральной (по спектру) метрике (Ь2).

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

Интервал частот измерений выбирался из естественных физических соображений таким образом, чтобы толщина скин-слоя на минимальной частоте была больше слоя, в котором ищутся неоднородности проводимости, а на максимальной частоте зондировался бы тонкий поверхностный слой, в котором проводимость не отличается от поверхностной. В рассматриваемом примере, где восстанавливается неоднородность проводимости на глубине 1 км, при изменении частоты от 3 до 200 Гц толщина скин-слоя изменяется от 3 до 0,35 км. Для достижения достаточно высокой точности расчетов интервал интегрирования в (9) составлял от 0 до 20 км.

На рис. 1 представлены результаты расчёта поля Н0(а>) для неоднородного профиля, показанного жирной линией на рис. 2б и для соответствующего по значению на поверхности од-

Рис. 2. а) Исходные «данные измерений» при среднеквадратичном значении случайной погрешности 8Н0 = 10 в зависимости от частоты f = с/2п; б) жирная сплошная кривая - исходный профиль проводимости, пунктир - результат восстановле-

нородного распределения проводимости о = = 0,01 См/м. Можно видеть, что существует область частот, где эти частотные зависимости заметно различаются. Их разность используется в уравнении (9) для решения задачи восстановления профиля проводимости при среднеквадратичном уровне моделируемой погрешности измерения мнимой части комплексной амплитуды ЪН1 = 10 (рис. 2а) и ЪН$ = 1 (рис. 3а).

На рис. 2, 3 показаны результаты численного эксперимента. Видно, что качество решения существенно зависит от точности измерений.

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

0.01 0.011 0.012 0.013 0.014 0.015

ния

Рис. 3. а) Исходные «данные измерений» при среднеквадратичном значении случайной погрешности ЪН0 = 1 в зависимости от частоты f = ®/2л; б) жирная сплошная кривая - исходный профиль проводимости, пунктир - результат восстановления

0.01 0.011 0.012 0.013 0.014 0.015

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

2. Сканирующая низкочастотная томография трёхмерной структуры проводимости земной коры. Следуя [6], рассмотрим область рассеивающей среды с комплексной диэлектрической проницаемостью е(г, с) = е0(с) +

+ £^( г , с) погружена в полупространство г <0, где е = е0(с). Для зависимости полей от времени вида ехр(/ю?) невозмущённое поле без рассеивающей области Е0(г) определяется свёрткой источника с функцией Г рина:

Е (г ) = Е0(—) +—-—| Я (г, —' )Е0(—' )&'. (11)

4п—0 V

Резольвента Я определяется распределением — и функцией Грина в виде известного ряда Неймана. Уравнения (10)-( 11) решают прямую задачу электродинамики, но, к сожалению, нельзя использовать (11) для решения обратной задачи в рамках рассматриваемого подхода, поскольку это выражение не является свёрткой. Поэтому приходится ограничиваться рассмотрением обратной задачи в рамках борновского приближения ( — (—) <<— 0), где Е(г) « Е0(г) + Е1(г) и рассеянная компонента поля определяется как

-* 1 Ґ Т —

Е1(—) =-------1—(—' 'Ре (— - —' )Е0 (г-Уг' • (12)

4п— 0 V

Если представить (12) как интеграл от эквива-

— т —

лентного источника тока ] =---------—1Е0

4п

Е (—) = —1— | ](—')&е (— - —', (13)

т— 0 V

то электрическую и магнитную компоненты рассеянного поля можно представить в виде разложения по плоским волнам:

Е0(г ) =-

1

гт—

0 V

- г')йг'. (10) ехр {-ікх (х - х')- гку (у - у’)± г^2 -к\ (7 - 7 ')}

При наличии рассеивающей области поле Е(г ) в среде определяется как

і^2 - к2

х{[Е0х (k2 -кх2) - Е0уКуКх + Е07 (±Kx^/kГ-K2)] Х0 +

+[-Е0хкхку + Е0у (02 - КУ ) + Е07 (±кут]к2 - кЛ )]у0 + + [Е0х (±K^^/k7^' - Е0у ^уЛ^ ' +

+ Е07кЛ]70} • ІЇКх(ЇКу(№' , (14)

Н1(х у^7) = -1 \тІг-3 • А —!(х ' ^ у ' ^7 ' '• 2с V' 16п

ехр {-ікх (х - х ') - іку (у - у') ± і^к^-кї(7 - 2 ')}

' ф2 - к х

{[Е0у (х '^ у '^ 7 ' )(+Цк 2 -к2 ' - Е07 (х '^ у '^ 7 ')х х (-гКу )]х0 +' [Е0і (х , у , 7 ' • (±г00 - К ' +

+ Е07 (х '^ у '^ 7 ') • (ІКх )]у0 + [Е0х (х '^ у '^ 7 Г) X

х (-іку' + Е0 у(х' ^ у '^ 7 ') • (ікх )]70 } х

"0 у

х с(кг(ку^'

х у

л

Н0(г' = Н0у ехР| Т + *Т Iу0 =

5 5

= Е

с(і -1) (7 ,7

0 х

где 5 =

т5 ПГ і 5 1 у0,

(17)

л/2

- толщина скин-слоя. Можно

пт о о

видеть, что волна в среде не зависит от поперечных координат и, следовательно, войдет в

интегральные уравнения (14), (15) лишь как множитель при — и тогда эти уравнения будут уравнениями типа свёртки по поперечным координатам. Это позволяет выполнить по этим координатам двумерное фурье-преобразование, которое приводит к одномерным интегральным уравнениям для каждой пары поперечных волновых чисел полученного пространственного спектра рассеянных компонент электрического и магнитного полей:

Е1(кх ,ку т) = -П— 8п —

0

X ЄХр(^ + і К) • ехр [±^02 -к22 (7 - 7 ')] Ъ о

(02 - °1 'х0 - КхКуу0 + кх4°

2 2 — -кЛ 70

і\[к~2

к/02 - кЛ

(7 ' , (18)

(1 + і)т5

0 у 4с2

где k =1 — I е . Таким! образом, если имеется

возможность измерять двумерное распределение компонент рассеянного поля, задача томографии сводится к решению трехмерных интегральных уравнений типа (14)—(15). Такие обратные задачи до сих пор практически не решались из-за исключительно высоких требований к памяти и быстродействию компьютера. Однако в данном случае можно использовать тот факт, что (14)-(15) представляют собой уравнения типа свёртки по поперечным координатам.

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

Е0, Н 0 имеют вид волны, затухающей с ростом глубины в поглощающем полупространстве:

Е0(Г) = Е0х ехР(/Ь)х0 = Е0х ехр( 4 +г' 4 Iх0,(16)

7 . 7

X ехр(— + і—) • ехр|-^02 - кЛ (7 - 7')], 5 5

у0

40

2 2 —

-°л +^2'

іу[к~2

і\102 - кЛ

(19)

где О] - вариации проводимости, определяющие е^. Если, например, методом двумерного сканирования измеряются вариации у-ком-поненты магнитного поля на поверхности грунта (г = 0) в зависимости от частоты падающей волны, то из (19) следует еще более простое уравнение:

Н1 у (кх ,ку ,т) = Н

0 у

(1 + і)т5 4с2

х Ї —1(кх ,ку ^ 7 ')ехр[5 + і 5 |х

:ехр

•і2 2 2 г

-и\0 -к -к7

(20)

(7'

Уравнение (20) является уравнением Фредголь-ма 1-го рода относительно глубинного распределения поперечных спектральных компонент магнитного поля. Можно видеть, что вклад неоднородностей среды в возмущения полей экспоненциально уменьшается на масштабе скин-слоя 4 величина которого зависит от частоты, возрастая с уменьшением последней, что и лежит в основе метода глубинного ОНЧ зондирования. Вместе с тем, существует и другой вид затухания, определяемый второй экспонентой, который для ближнепольных компонент попе-

)

х

х

—со

х

х

с

речного спектра таких, что у2 = К + к2у > |^| ,

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

Измеряемые в сканирующей низкочастотной томографии двумерные распределения полей в зависимости от частоты о связаны с объёмным распределением проводимости среды трёхмерными интегральными уравнениями 1-го рода (18), (19), которые можно записать в общем виде как:

Нх( х, у, о) = К (х - х ', у - у' , г О) ;

V'

х|ст1 (хуг')dx' йу' йг',

Н1Кх Ку,о) =

Далее уравнение (22) решается как одномерное интегральное уравнение Фредгольма 1-го рода относительно глубинного профиля поперечного спектра 71(ух ,у, г') для каждой пары к, К в

интервалах дискретного спектрального разложения.

Наиболее последовательный подход к решению таких одномерных уравнений основан на теории некорректных задач Тихонова [5], в котором <71(кх ,к , г ') может быть найдено методом

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

. < К, (23)

Н - нЛ г <8, ||К - К„

(21)

где ядро уравнения К - аппаратная функция, Н\(х, у, о) - измеренный сигнал, 7^(х', у', г') -искомое распределение вариаций проводимости в полупространстве г < 0. Зависимость ядра уравнения К от частоты позволяет в принципе искать решение (21) относительно трёхмерного распределения 71.

В данной работе используется тот факт, что двумерное фурье-преобразование по поперечным координатам уравнения (21) приводит к одномерным интегральному уравнению типа (20) относительно глубинных распределений компонент поперечного спектра вида:

Н1(кх Ку=

0

= 4^2 | К(кх Ку,г ',°)71(кх Ку,г ')0г ', (22)

где

где Н1 - левая часть (22), соответствующая

ТТ 8 «

точному решению; Н1 - данные измерений с

ошибкой 8 Согласно методу Тихонова приближённое решение (22) находится из условия минимума функционала обобщённой невязки:

II о 812 и м2

МЛ71] = ||4п Кк71 - Н1 ||4 +«N1 ^ (24)

при дополнительном условии, которое выражает принцип обобщенной невязки:

||4п2Кй71 - Н18|[ = (8 + К|7Л 4)2. (25)

При (8, К)—>0 восстановленное распределение равномерно сходится к истинному. Параметр погрешности ядра К определяется обычно ио численного решения прямой задачи, и его стремятся сделать меньше меры погрешности измерений 8, которая согласно (25) с учетом теоремы Планшереля может быть оценена как

82 =--------1----- |Г йкхйку х

Akx Aky Ьо}} х у

-П-Н1(х,уо)ехр(-Кхх - Куу)йхйу , х | [8Н1(}х} ,о)]2 4п о

а под интегралом - аналогичные фурье-транс-форманты ядра уравнения и искомого распределения проводимости.

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

1

4п АхАуАо

|| йхйу | [8НХ (х, у, о)]2

х, у

о

х йо =■

1

(26)

2 82,

4п Е

где 8Н1(К Ку}° =

- Н^ К х, у у ,о) , а интегрирование выполняется по частотной области анализа. Искомое решение получается путём обратного преобразования Фурье полученного глубинного профиля

компонент вариаций латерального спектра проводимости:

71(X У, z) = jj71(yx ,yy . z) x exp(iyxx + iyyy )dyxdyy

exx

(27)

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

def (а) =

| z|| K(x, y, z, а)dxdy dz

(28)

Точность решения быстро убывает на глубинах г < -йед-.

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

Результаты и обсуждение

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

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

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

Представленные результаты получены по программам фундаментальных исследований Отделения физических наук РАН (проекты «Когерентная томография земной коры с использованием ультранизкочастотных электромагнитных полей» и «Ближнепольное электромагнитное зондирование диэлектрических и проводящих сред») и при поддержке РФФИ (грант № 07-02-00191).

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

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

2. Велихов Е.П., Жамалетдинов А.А., Собчаков Л.А. и др. // Доклады АН РАН. 1994. Т. 338. № 1. - С. 106-109.

3. Поляков С.В., Ермакова Е.Н., Поляков А.С., Якунин М.Н. // Геомагнетизм и аэрономия. 2003. Т. 42. № 2. С. 240-248.

4. Gaikovich K.P. Inverse Problems in Physical Diagnostics. - Nova Science Inc.: N.Y., 2004. - 372 p.

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

6. Gaikovich K.P. // Proceedings of 2006 8th International Conference on Transparent Optical Networks, ICTON 2006, Nottingham, United Kingdom, June 1822, 2006. V. 1. P. 250-255.

INVERSE PROBLEMS OF COHERENT LOW-FREQUENCY SENSING OF THE EARTH CRUST

K.P. Gaikovich, A.I. Smirnov

We consider inverse problems of reconstructing the conductivity structure of the Earth’s crust by using the data of spectral VLF measurements. In the case of one-dimentional depth distribution of the conductivity, the Maxwell equations yield nonlinear integral equations relating the measured electric-field components and the conductivity profile. An iterative algorithm for solving such an inverse problem is developed on the basis of the Tikhonov theory of ill-posed problems. A tomographic technique is proposed for probing the three-dimensional conductivity structure of the Earth’s crust. The technique is based on near-ground measurements of the frequency dependences of two-dimensional distributions of the tangential electric- or magnetic-field variations. It is shown that in the Born approximation these variations can be represented by a triple integral in the form of convolution over the transverse coordinates. The Fourier transform over the transverse coordinates reduces the problem to a one-dimensional integral equation for the depth profile of the lateral spectral components of the conductivity. The latter equation can be solved by the Tikhonov method.