Двумерная инверсия магнитотеллурических данных с учётом влияния рельефа поверхности Текст научной статьи по специальности «Математика»

УДК 550.837.211

Вестник СПбГУ. Сер. 4. 2013. Вып. 3

С. А. Вагин, А. В. Козлова, И. Л. Варданянц

ДВУМЕРНАЯ ИНВЕРСИЯ МАГНИТОТЕЛЛУРИЧЕСКИХ ДАННЫХ С УЧЁТОМ ВЛИЯНИЯ РЕЛЬЕФА ПОВЕРХНОСТИ

Введение. Складчатые области, как правило, характеризуются неоднородным по горизонтали строением, и 1D интерпретация при этом неприменима. Кроме того, объём достоверной априорной информации при изучении складчатых областей обычно невелик. В этих условиях интерпретацию начинают с 2D-инверсии, позволяющей получить сглаженное распределение сопротивлений и отражающей основные геоэлектрические структуры [1]. Успех интерпретации зависит не только от выбора стартовой модели, которую желательно согласовать с известными представлениями о регионе, но и от применяемой методики расчета. Для определения влияния рельефа поверхности авторами статьи была разработана и опробована методика 2D-инверсии, основанная на методе минимальных невязок с адаптивной регуляризацией.

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

1.1. Метод минимальных невязок (ММН). Рассмотрим линейное операторное уравнение

f = Lm, (1)

где f, m € M, а L — линейный непрерывный оператор в вещественном гильбертовом пространстве M. Итерационное решение (1) выразим формулой

m„+i = m„ + Дт„ = m„ - knrn, n = 0,1, 2,..., (2)

где mo — начальное приближение, Дmn — итерационный шаг,

Дmn = -knrn, (3)

rn — невязка на n-м шаге,

r„ = Lm„ - f, (4)

kn — итерационный коэффициент.

r„+i = Lmn+i - f = Lmn - knLrn - f = rn - knLrn. (5)

Выбор kn заключается в минимизации нормы невязки rn+i:

KU2 = ||rn+i - knLrn\\2 = ) ^ min. (6)

Станисло,в Александрович Вагин — доктор физико-математических наук, профессор, Санкт-Петербургский государственный университет; e-mail: [email protected]

Александра Вячеславовна Козлова — соискатель, кафедра физики Земли, Санкт-Петербургский государственный университет; e-mail: [email protected]

Изабелла Левоновна Варданянц — кандидат физико-математических наук, научный сотрудник, Санкт-Петербургский государственный университет; e-mail: [email protected] © С. А. Вагин, А.В.Козлова, И. Л. Варданянц, 2013

Функцию ||Ф||2(&„) можно задать, используя операцию скалярного произведения в M: Ф2(кп) = (rn - knbVn, Гп - knLrn) = (rn, r„) - 2кп(тп, Lrn) + кn(Lrn,Lrn) ^ min. (7) Дифференцируя Ф2 по кп, получим условие минимума ||rn+i||2:

= -2(гп, Lrn) + 2kn{Lvn, Lrn) = 0;

dkn

J (rn ,Lrn)

(8)

(Lrrn Lrn)

Соответствующий минимум Ф(kn), согласно (7),

(rn,£rn)2 _ _ (rn, Lrn)2

(Lr„,Lr„) ||rn||2||Lr„

цкп) = iir„+iii = \i(r„,r„)- Ti\ = iir„nji-(9)

Из (9) видно, что на (п + 1)-й итерации всегда

1К+111 < \М. (10)

Равенство возможно, если Ьгп и гп взаимно ортогональны: (Ьгп, гп) = 0. Рекуррентные формулы (2)—(8) составляют метод минимальных невязок (ММН). В работе [2] этот подход был распространен на более общий случай линейного операторного уравнения в комплексном гильбертовом пространстве М. В комплексном случае свойство симметрии скалярного произведения выглядит иначе:

(Г, &) = (& Г)*, (11)

где звёздочка означает комплексное сопряжение. Из этой формулы следует

(Г, к%) = к*(¥, %)*.

Итерационный коэффициент кп в (8) может иметь комплексное значение, так как (тп,Ьтп) — комплексное число. Формула (9) примет вид

Ф(кп) = ||гп+1|| = у/(гп -кпЬгп,гп -кпЬгп) = \\rjjl - п1(г"2(12)

V \\тп\\ \\Ьгп\\

В комплексном случае оператор Ь должен быть строго положительно определённым:

|(т,£т)| ^ у(т, т) = у\\т\\2, у > 0 (13)

для всех т € М.

1.2. Линеаризация оператора прямой задачи. Магнитотеллурическая задача является нелинейной, т. е. нелинеен оператор А прямой задачи в уравнении

Аш = а. (14)

Для дискретной задачи вектор данных а и вектор модельных параметров ш при помощи операции транспонирования запишем в виде

а = [^1 ,й2,...,йГ1\Т ,

Т (15)

ш = [Ш1,Ш2,...,ШЙ] .

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

к дА1 (щ(°>) 1 (т(°У)

дтл 2\ дттдт!

3 о —1 /— 1 3 1

3=1 3 3 = 1 1=1 г = 1,...,п; ] = I = !,...,к

к

0)\ . '

111

' *—' дт3

3=1

При малых поправках связь между Дт и отклонениями наблюдений Дd приближённо будет линейной:

Л(т) = Л, (т(0)) + ЛДт, (16)

где в случае конечномерного пространства моделей М А — матрица с элементами

Аго(т(°))= , '' I...../<:./ I...../■•• (17)

\ / ^ от^ ) т=т(о)

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

ЛДт(я) = Дd(s-1) = d - Л(т(в_1)), (18)

для которой элементы матрицы А находятся по формуле для модели, полученной на предыдущей итерации:

^'=[тГ~Ч , г = 1,. ..,п; = 1 (19)

^ °т3 ) т=т(="1)

Размер матрицы чувствительности А — п х к, где п — число наблюдений, к — число параметров модели.

Сведение нелинейного оператора (14) прямой задачи к линейному оператору вида (18) называется линеаризацией задачи. Основное преимущество линеаризованной обратной задачи в том, что для линейных задач существуют хорошо разработанные методы и алгоритмы. Кроме того, в тех случаях, для которых оператор решения прямой задачи относительно исходной модели т может быть описан только алгоритмически, можно построить оператор преобразования поправок к модели Дт в поправки к наблюдениям Дd в виде явной функциональной зависимости. 1.3. Регуляризация. В общем случае обратная задача (18)

ЛДт = Дd

является плохо обусловленной. Поэтому линейный оператор Л может не удовлетворять условию (13). В этом случае нужно использовать регуляризацию. Пусть М и Б

или

вещественные или комплексные гильбертовы пространства. Рассмотрим минимизацию параметрического функционала Тихонова

P"(Дт) = ЩДш - Ad||2 + as (Дш) ^ min, (20)

где a > 0, а s^m) — некоторый квадратичный функционал, например,

s (Дш) = (WДш, WДш)м = HWДш||2,

W — положительно определённый линейный непрерывный оператор в M.

Применив вариационный оператор к функционалу Pa и найдя необходимое условие его минимума, получим уравнение Эйлера:

A* Дd = (A*A + aW *W )Дш. (21)

Линейный оператор

La = A*A + aW *W

всегда строго положительно определён, и, следовательно, уравнение Эйлера имеет единственное решение Дш"1, которое можно получить рассмотренным методом минимальных невязок, если положить f = A* Дd и L = La = A*A + aW*W.

Однако можно получить более эффективный регуляризованный ММН, если минимизировать параметрический функционал (20), как сделано в [2]. Окончательно регуляризованный алгоритм ММН для линейной дискретной обратной задачи на n-й итерации можно записать в виде

rn = АДшп - Дd,

Ln = A*rn + "^¥Дшп,

ка = ||Ь«||2 (22)

" || AL" ||2 + ||WL" ||2'

ДшП+1 = ДшП - knLn.

Для линейной дискретной задачи оператору A соответствует матрица A (19), которая должна обновляться на каждой итерации. В программе предусмотрена также возможность вычисления её один раз для нескольких последовательных итераций. Начальная модель поправок Дшо для (22) берётся в виде вектора дш, см. (19). Весовая матрица модели W в простейшем случае равна единичной. Начальное значение параметра адаптивной регуляризации ao берётся равным отношению нормы невязки к норме стабилизирующего функционала (20). Затем параметр регуляризации изменяется на каждой итерации по формуле

an = aoqn-1, n = 1, 2,3,..., q> 0. (23)

В программе взято q = 0,75. Процесс итераций прекращается, когда невязка достигает заданного уровня.

1.4. Блок-схема алгоритма двумерной инверсии МТ данных с учётом рельефа профиля. В качестве магнитотеллурических данных при инверсии были взяты кажущиеся сопротивления р^ (x,t), при этом вектор данных d = log р^ (x,t), вектор модельных параметров ш = log р(х, z), где р(х, z) — вектор удельных сопротивлений геоэлектрического разреза. Согласно п. 1.2, для применения инверсии необходимо MX х MZ х NT раз, исключая точки дискретной сетки модели для воздуха, решать прямую двумерную магнитотеллурическую задачу на каждой итерации, где MX, MZ

и NT — количество точек по оси x, z и периодов, соответственно. Блок-схема алгоритма двумерной инверсии изображена на рис. 1. Управление алгоритмом осуществляется файлом "Start", в котором задаются название входного файла и основные параметры алгоритма.

Start

Рис. 1. Алгоритм двумерной инверсии магнитотеллурических данных с учётом рельефа профиля

Начальную модель определяет «Файл входных данных». Значения расстояний по осям x и z даны в километрах. Ось z направлена вниз. Массив PROFD несёт информацию о рельефе профиля, начиная с левой точки оси x, значение по оси z принято равным нулю. Начальное распределение сопротивлений разреза в земле либо постоянное во всех ячейках, либо привязывается к псевдоразрезу. Сопротивления в ячейках воздуха в обоих случаях берутся равными 108 Ом • м и в процессе итераций не изменяются. Прямая 2D-задача решается по методу конечных разностей с учётом рельефа местности путём модификации программы, подробно описанной в [3, 4]. В «Блоке инверсии» осуществлён регуляризованный алгоритм метода минимальнах невязок с адаптивной регуляризацией. Невязки контролируются блоком «RMS». В блоке «Файл выходных данных» выдаются основные файлы инверсии на каждой итерации, что очень удобно для анализа результатов и прекращения процесса инверсии.

2. Тестирование алгоритма. Тестирование алгоритма проводилось на модельных данных с учётом рельефа профиля и без его учёта. Были выбраны две принципиальные модели горизонтально-слоистой геоэлектрической среды с сопротивлением 108, 102, 106, 1 Ом-м соответственно, осложнённой рельефом поверхности: в первом случае «горой», во втором «ямой» (рис. 2). Синтетические данные моделей были получены для пятнадцати периодов в интервале 0,0026-900 с для тридцати одной точки зондирования. Важной характеристикой алгоритма является скорость сходимости в зависимости от числа итераций и устойчивость решения. Как видно на рис. 3, сходимость, полученная для данных, учитывающих изменение рельефа вдоль профиля, быстрая и ровная, и, следовательно, решение устойчиво. В то же время для данных МТЗ, смоделированных без учёта дневной поверхности, в случае с моделью «яма» сходимость хуже и медленнее, а в случае с моделью «гора» ещё и искажена на последних итерациях. Для инверсии в качестве начальной модели была использована модель полупространства с удельным сопротивлением 100 Ом-м с такой же сеткой, как в истинной модели. Характерно, что уже на 20-й итерации модельная и экспериментальная кривые для моделей с введённым массивом данных, характеризующих рельеф, различаются мало, геометрия и геоэлектрические свойства среды определяются стабильно и, в целом, верно (рис. 4).

Для модели «гора» (рис. 4, а): X = [—0,2 : 0] — в эту область попадают модельные станции, лежащие внутри области изменения рельефа, X = [—1,8 : 1,2] — станции на границе области, X = [—15,0 : —2,2] — станции, находящиеся вдоль горизонтальной

а б

Рис. 2. Геоэлектрическая модель, осложнённая рельефом земной поверхности — горой (а) и ямой (б)

Рис. 3. Зависимость скорости сходимости от количества итераций для теоретических моделей «яма» и «гора»

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

Для модели «яма» (рис. 4, б): X = [—0,4 : 0] — в эту область попадают модельные станции, лежащие внутри области изменения рельефа, X = [—2,6 : 2,0] — станции на границе области, X = [—15,0 : —3,0] — станции, находящиеся вдоль горизонтальной поверхности. Сплошной толстой линией изображена теоретическая кривая зондирования, характерная для наблюдаемых точек. Как видно на рисунке, наиболее проблемной областью являются точки зондирования, находящиеся в непосредственной близости от внесённых в модель геоморфологических элементов. Кажущееся сопротивление слоя, лежащего непосредственно под «горой» или «ямой», занижено, а форма наблюдённой кривой наиболее искажена относительно модельной.

Аналогичные построения были выполнены для моделей без учёта рельефа.

Модель «гора». Для кривых точек зондирования, попадающих в интервал [—0,2 : 0], характерно завышение ширины слоя сопротивлением 108 Ом • м. Следующий за ним слой с сопротивлением 100 Ом • м выглядит весьма неоднородно с экстремумом в середине, в то же время ширина слоя с сопротивлением 106 Ом • м меньше модельной, а его

-1 км -16 км -500 км

-2 км

-16 км -500 км

--модель

С учётом рельефа: Без учёта рельефа:

---|х| > 1 ---|х| > 1

-------х е -(-1,5, -2) -------х е -(-1,5, -2)

-х е -(-4, -7) -х е -(-4, -7)

модель

С учётом рельефа:

---|х| > 2

------- х е -(-2, -2,5)

- х е -(-7, -15)

Без учёта рельефа:

- |х| > 2

- х е -(-2, -2,5)

- х е -(-7, -15)

Рис. 4- Графики логарифма кажущегося сопротивления для моделей «гора» (а) и «яма» (б) с учётом и без учёта рельефа

0

сопротивление превышает заданное значение. Для интервала станций зондирования X = [—2 : 1,5], находящихся на минимальном удалении от элемента рельефа, характерно максимальное искажение формы кривой зондирования. Максимально совпадающим с теоретической кривой является интервал, соответствующий слою с сопротивлением 108 Ом • м. Интервал, соответствующий слою 100 Ом • м, выглядит незначительным по ширине минимумом кривой сопротивления. Слой с сопротивлением 106 Ом • м в отражении кривой МТЗ смещён от центра и значительно шире, чем модельный. Наибольшего совпадения теоретической и модельной кривых зондирования удалось достичь в интервале X = [—5,0 : —4,0], где расположены точки зондирования, на которые возмущающий геоморфологический элемент производит минимальный эффект.

Модель «яма». Наибольшее искажение модельной кривой приходится на интервал, находящийся непосредственно на участке со структурным элементом «яма», где X = [—0,2 : 0]. В этом случае для модельной кривой зондирования ширина слоя с сопротивлением 108 Ом • м существенно меньше, чем для теоретической кривой, а сам слой оказывается разбитым на два участка с удельным сопротивлением 108 и 106 Ом • м соответственно. Отклик от модельного слоя 106 Ом • м «размазан» по интерпретационной кривой и перекрывает два интервала: собственно от самого слоя и от слоя с сопротивлением в 100 Ом • м. Кривые зондирования точек, расположенных в пределах X = [—15,0 : —7,0], дают неплохой результат. Однако необходимо отметить, что в пределах данной области амплитуда экспериментальной кривой завышена относительно теоретической. В последнем слое модели добавляется участок высокого сопротивления, около 1000 Ом • м.

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

Рассмотрим результирующие разрезы моделей «гора» (рис. 5) и «яма» (рис. 6) с учётом и без учёта рельефа. Минимальная RMS результирующих разрезов 2D-инверсии для модели «гора» с учётом рельефа составила 0,03, а без учёта 0,17. Для модели «яма» — 0,02 и 0,09 соответственно. Результирующие разрезы моделей «гора», «яма» позволяют сделать вывод, что для моделей с учтённым рельефом максимальная разрешающая способность 2D-инверсии снижается в краевой области аномалий дневной поверхности. Скорее всего, это связано с сильным искажением поля на границе диэлек-

5

0

-5

-10

-15

м -20

к

tN -25

-30

-35

-40

-45

-4 -3 -2 -1 0 1 2 3 4 5

X, км

РГн -i» || ¡= п

1 l" "l" "l" 1

-4-3-2-10 1 2 3 4 X, км

ln R от -2 до 1,9 от 1,9 до 2,6 от 2,6 до 6 от 6 до 8,001

к

N"

Рис. 5. Результирующие разрезы модели «гора» с учётом (а) и без учёта (б) рельефа а б

0-5 --10 --15 --20 --25 --30 --35 -

т-т-^-1 ■ I

1012345 X, км

Г^Т" - I

-3-2-10 1 2 3 4 X, км

ln R

■ от -2 до 1,9

от 1,9 до 2,6

от 2,6 до 6

от 6 до 8,001

Рис. 6. Результирующие разрезы модели «яма» с учётом (а) и без учёта (б) рельефа

б

а

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

Выводы.

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

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

3. Разработанный алгоритм двумерной инверсии МТ данных на профиле с изменением рельефа и соответствующая ему программа могут быть использованы при планировании эксперимента и обработке полевых данных в районах с гористым рельефом.

Литература

1. Бердичевский М. Н., Яковлев А. Г., Алексанова Е. Д. и др. Технология и результаты региональных магнитотеллурических исследований // Разведка и охрана недр. 2004- Вып. 5, № 16. С. 37-40.

2. Жданов М. С. Теория обратных задач и регуляризации в геофизике. М.: Научный мир, 2007. 710 с.

3. Варданянц И. Л. Расчёты методом сеток магнитотеллурических полей над двумерно-неоднородными средами (часть I) // Вопросы геофизики. Л., 1978. Вып. 27. С. 36-40.

4. Варданянц И. Л. Расчёты методом сеток магнитотеллурических полей над двумерно-неоднородными средами (часть II) // Вопросы геофизики. Л., 1979. Вып. 28. С. 153-163.

Статья поступила в редакцию 6 марта 2013 г.