ГЕОМАГНЕТИЗМ И АЭРОНОМИЯ, 2014, том 54, № 4, с. 570-576
УДК 550.83.015;550.383;550.384
ИНТЕРПРЕТАЦИЯ МОРСКИХ МАГНИТНЫХ АНОМАЛИЙ. ЧАСТЬ 2. АНАЛИЗ НОВОГО МЕТОДА И АЛГОРИТМА НА ОСНОВЕ МЕТОДА НАИМЕНЬШИХ КВАДРАТОВ © 2014 г. С. А. Иванов, С. А. Меркурьев
СПб филиал Института земного магнетизма, ионосферы и распространения радиоволн им. Н.В. Пушкова РАН, г. Санкт-Петербург e-mail: sergei.a.ivanov@mail.ru
Поступила в редакцию 19.08.2013 г.
В данной работе рассматривается задача интерпретации линейных морских магнитных аномалий в рамках спрединговой модели. На основе метода наименьших квадратов (МНК) предложен новый алгоритм определения положения блоков прямой и обратной полярности. Глобальный минимум невязки ищется путем комбинации мультистарта и метода Монте-Карло. На примере поля трех блоков, возмущенного нормальным шумом, показано, что метод дает погрешность, близкую к минимально возможной. На модели, построенной по шкале инверсий, выяснено, что предложенный алгоритм позволяет полю (при отсутствии шума) найти положения границ блоков с большой точностью, определяемой, по существу, временем счета.
DOI: 10.7868/S0016794014040208
1. ВВЕДЕНИЕ
Ключевым этапом геоисторической и палео-магнитной интерпретации морских магнитных аномалий является задача нахождения границ блоков прямой и обратной полярности, связанных с изменением направления ГМП в прошлом. Обзор методов, применяемых при интерпретации морских магнитных аномалий, изложенный в работе [Гордин и Золотов, 1990] и выполненный в работе [Иванов и Меркурьев, 2014] анализ некоторых алгоритмов, показал, что существующие методы и алгоритмы не в полной мере приспособлены для решения этой задачи.
Главная цель данной работы состоит в разработке и исследовании метода интерпретации морских магнитных аномалий на основе МНК. Рассматривается обратная задача магнитометрии в рамках модели Вайна-Мэйтьюза. В качестве данных обратной задачи берется аномальное магнитное поле на профиле, а также толщина и верхняя кромка магнитоактивного слоя. Решением обратной задачи мы называем модуль намагниченности и координаты блоков прямой и обратной полярности, создающих магнитные аномалии.
На разобранном в работе [Иванов и Меркурьев, 2013] примере, где рассматривается поле от четырех блоков при добавленном шуме, мы проведем исследование эффективности метода МНК и сравним погрешности определения границ блоков методом МНК с минимально возможными погрешностями, полученными на основе неравенства Рао-Крамера.
На модельном поле в отсутствии шума, создаваемом 11 блоками протестируем метод МНК и проверим точность определения границ блоков. Рассмотрим также задачу интерпретации магнитных аномалий в случае медленного спрединга, когда магнитоактивный слой содержит узкие блоки.
2. МНК В ЗАДАЧЕ ОПРЕДЕЛЕНИЯ ГРАНИЦ БЛОКОВ
2.1. Геофизическая модель
В качестве физической модели океанического магнитоактивного слоя рассматривается блоковая инверсионно-спрединговая модель Вайна и Мэтьюза [1963].
Как правило, номерные магнитные аномалии Ламонтской последовательности представляют собой периоды преимущественно одной полярности главного магнитного поля (хроны), разделенные периодами противоположной полярности. Поэтому мы выбрали простую типичную модель фрагмента магнитоактивного слоя, рожденного на оси срединно-океанического хребта в условиях разрастания океанического дна, чтобы показать потенциальные возможности метода и оценить возможные погрешности выделения короткопе-риодных инверсий для данной модели.
Рассмотрим однородно намагниченные прямоугольные блоки, верхние и нижние грани которых залегают соответственно на глубине а и Ь. Чтобы рассчитать значение поля в данной точке от магнитоактивного слоя, представленного по-
следовательностью призм, необходимо рассчитать вклад, вносимый каждой призмой и просуммировать эти значения. Мы предполагаем, что хребет расположен вдоль оси у и намагниченность зависит только от х. Для простоты расчетов считаем, что главное поле ВЕ и намагниченность М направлены вертикально, причем главное поле направлено вниз, в отрицательном направлении. Пусть X = {XI, Х2, Х3,...., Хм} — точки перемены полярности, или контакты.
2.2 МНК для блоковой модели
Пусть в точках А = {а1, а2, а3, ..., ат} задано аномальное поле {у1, у2,у3,..., ут}. Считаем известной глубину верхней кромки базальтового слоя а, и глубину нижней кромки базальтового слоя Ь. Будем минимизировать невязку как функцию координат контактов X1 и модуля намагниченности М. Эта функция строится следующим образом:
Выбирается число Ь контактов с х-координа-
тами блоков X1 = {х1} . Соответствующие блоки имеют вид Б1 = (X), Х]+1). Намагниченность блоков считается вертикальной, равной по амплитуде и меняющей знак. То есть намагниченность у-го блока есть (—1) М, где величина Ммо-жет быть отрицательной.
В точках А считается поле {иь и2, и3,..., ит}, порожденное блоками Б1 и считается невязка
Е(Х>) = ¿(и, - у()2
(1)
400 м. Наблюдения осуществлялись на поверхности г = 0 со среднеквадратичной ошибкой а = 63.7, составляющей 10% от максимального значения.
В основном, мы интересуемся положением блока (Х2, Х3). Для большего приближения к реальной ситуации мы полагаем неизвестными и внешние границы блоков, задаваемые точками Х1 и Х4.
Пусть наблюдения У; проводятся в точках а1, а2,..., а„. Поле в этих точках запишем в виде функций у1 = у1 (Х) от Х, поскольку все остальные величины фиксированы. Пусть £1 есть ошибка в 1-м наблюдении у, = у{ + бБудем считать ошибки независимыми и имеющими нормальное распределение с нулевым средним и одинаковой дисперсией а2. Совместная плотность ошибок есть произведение плотностей ошибок отдельных наблюдений
РИ1,М2,...,«я) =
1
7(2яо) п
-2(и;- - у, )2/2а 2
МНК состоит в поиске глобального минимума этой функции.
2.3. Неравенство Крамера-Рао.
Статистическая модель
В качестве модельного примера для исследования эффективности МНК рассмотрим, как и в работе [Иванов и Меркурьев, 2013], фрагмент магнитоактивного слоя состоящего из двух блоков прямой полярности, разделенных узким блоком обратной полярности. К аномальному магнитному полю, создаваемому этими блоками, добавляется нормально распределенный шум со среднеквадра-тической ошибкой, составляющий 10% от максимальной амплитуды аномального поля.
Координаты границ Х = {12.644, 17.524, 18.134, 25.576} блоков прямой и обратной полярности были получены на основании реального участка шкалы инверсий, содержащего короткий период противоположной полярности. При расчетах намагниченность блоков принята равной 10 А/м, глубина базальтового слоя а = 4 км, а его толщина
Предположим, что мы располагаем методом, позволяющим оценить одну координату £, одного из блоков. Пусть I = 1о§(р) и введем информацию
Фишера I = Е [(5 )2]. Тогда неравенство Крамера-Рао дает нижнюю оценку дисперсии: а2(^) > 1/1. Для оценки нескольких границ блоков составляется информационная матрица Фишера I, элементы которой есть математические ожидания Е[д ¡¡д ¡¡]. Тогда матрица ковариации допускает оценку через матрицу, обратную к информационной матрице: I-1 . В частности, дисперсии оценок оцениваются снизу соответствующими диагональными элементами матрицы I-1. Вычисления [Иванов и Меркурьев, 2013] показывают, что в данном случае неравенства Крамера-Рао [Боровков, 1984] принимают вид
а(Х2) > а™ = 0.7306 км, а(Х3) > а™ = 0.7373 км.
Для среднеквадратичной ошибки длины промежутка Х3-Х2 имеем
ст5 > а^ = 0.1822 км.
Для среднеквадратичной ошибки определения центра хс = (Х2 + Х3)/2 получено:
ас >ас/ = 0.7283 км.
Здесь ст(Х2), а(Х1), а§, ас, среднеквадратичные отклонения оценок соответствующих величин, полученных каким либо методом (точные формулировки см. в [Иванов и Меркурьев, 2013]).
2.4. Неравенство Крамера-Рао и МНК. Численные результаты
Для данного набора наблюдений функция р, есть функция от координат границ блоков, и на-
км
0
20 30 40 50 60 70 80 90 100 110 120 Число точек наблюдения
Рис. 1. Графики среднеквадратичной погрешности определения ширины блока 5 (верхние графики) и центра блока с (нижний график) по Крамеру-Рао (1) и найденные программой (2) в зависимости от числа точек наблюдения. Условные обозначения: 1 — графики погрешности определения ширины блока 5 и центра блока с по Крамеру-Рао; 2 — графики погрешности определения ширины блока и центра блока с, найденные программой.
зывается функцией правдоподобия. Те значения параметров X на которых функция р достигает максимума, называются оценками метода максимального правдоподобия. В случае независимого нормального шума оценки метода максимального правдоподобия, как легко видеть из явного выражения р, совпадает с оценками МНК. Как хорошо известно, см., например, [Боровков, 1984], при большом числе повторных измерений оценки МНК близки к оптимальным (имеют дисперсию, близкую оценкам Крамера-Рао). В рассматриваемой модели повторные измерения отсутствуют, но если измерения взяты в близких точках, то естественно ожидать оценок МНК близких к оптимальным. Строгие утверждения такого рода нам неизвестны и мы не ставили цель исследовать этот вопрос аналитически или численно, но здесь для выбранного выше примера мы численно покажем аналогичный факт. Конкретно, для разных реализаций шума и для разного числа точек наблюдения, мы найдем минимум невязки (МНК) и сравним полученную погрешность с оценками Крамера-Рао. Также покажем зависимость по-
грешности от шага дискретизации наблюденного поля. Сравнение погрешностей было организовано следующим образом:
1. Бралось число наблюдений ЛоЬ5 = 21, 41, 81, 121 в точках {а,,/ = 1, 2, ... Л^}, равномерно распределенных на фиксированном отрезке [Х1 — 1, Х4 + 1]. Строилось (истинное) поле Е(х) соответствующее выбранным блокам.
2. Считались оценки среднеквадратичных погрешностей по Крамеру-Рао для данных точек наблюдений и шума, описанного выше. Именно,
СЯ О ся
находились а § для ширины промежутка о и а с
для центра хс.
3. Брались 1000 реализаций шума, то есть 1000 наборов случайных чисел, {в/, / = 1, 2, ..., Лоья}, распределенных соответственно модели.
Для каждой из 1000 реализаций оп
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.