Interpretation of 2D magnetic anomalies using wavelet transform

Full Text

1. ВВЕДЕНИЕ

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

Для анализа магнитных данных предложено несколько методов, позволяющих интерпретировать геометрию источников. Среди них можно выделить методы, основанные на вычислении производных поля [Miller and Singh, 1994; Verduzco et al., 2004], [Wijns et al., 2005], [Cooper and Cowan, 2006], [Ferreira et al., 2013]. Эти методы просты в реализации и не требуют больших вычислительных затрат. Их обзор и сравнительный анализ можно найти в работах [Gunn, 1997; Nabighian et al., 2005; Cooper and Cowan, 2011].

Другая не менее важная задача при интерпретации магнитных аномалий состоит в определении глубины источника. Для решения этой задачи наибольшую популярность приобрели такие автоматические методы, как деконволюции Вернера [Werner, 1953] и Эйлера [Thompson, 1982; Reid et al., 1990], в которых проблема оценки глубины трансформируется в решение системы линейных уравнений. Другой подход к получению оценки глубины магнитоактивного слоя состоит в использовании спектрального анализа магнитометрических данных. По сути, этот метод представляет собой преобразование магнитного поля из пространственной области в частотную с использованием одно- или двумерного, в зависимости от вида исходных магнитных данных, преобразования Фурье с последующей его интерпретацией. Широкое распространение получил метод, предложенный в работе [Spector and Grant, 1970], где в качестве модели аномалеобразующих тел рассматривался ансамбль блоков различной глубины, ширины, толщины и намагниченности, и анализировалась форма спектра магнитных данных в зависимости от этих параметров. Возможность определения мощности магнитоактивного слоя по автокорреляционным функциям рассмотрена в работе [Портнова и др., 1987].

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

Для идентификации границ аномалеобразующих тел хорошо известен и широко применяется метод аналитического сигнала [Nabighian, 1972; Nabighian, 1974]. Популярность этого метода обусловлена тем, что при его применении не накладывается никаких ограничений на форму искомого источника. Главное преимущество аналитического сигнала состоит в том, что для двумерных магнитных источников форма и расположение максимумов его амплитуды не зависят от направления намагниченности [MacLeod et al., 1993]. Это свойство обеспечило его широкое применение для интерпретации морских магнитных аномалий в силу независимости от направления древнего и современного поля, что исключало необходимость учета косости аномалий при определении местоположения вертикальных границ между блоками противоположной магнитной полярности, например, [Roest et al., 1992].

Однако было показано, что форма аналитического сигнала магнитного поля дипольного или сферического источника зависит от направления намагниченности, и максимумы его амплитуды не всегда расположены непосредственно над такими источниками, их смещение зависит от расстояния источник-наблюдение и сдвиг может составлять до 30% расстояния [Salem et al., 2002].

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

В последние годы широкое распространение при анализе сигналов получило вейвлет-преобразование. Основы теории вейвлет-преобразования и примеры его использования даны в отечественных и зарубежных обзорах [Астафьева, 1996; Moreau et al., 1977; Kumar, Foufoula-Georgiou, 1997]. Вейвлет-преобразования, как и преобразование Фурье и любое другое интегральное преобразование, наследует свойства анализируемой функции и базисного ядра [Глазнев, 2003]. Но в отличие от преобразования Фурье, в котором анализирующая функция покрывает всю временную ось, двухпараметрическая порождающая функция одномерного вейвлет-преобразования хорошо локализована и во времени, и по частоте. Это свойство позволяет получить больше информации о сигнале, но интерпретация вейвлет-спектра не так очевидна, как анализ преобразования Фурье. Предлагаются различные схемы интерпретации геофизических данных с использованием различных типов вейвлет-преобразования: так, в работах [Оболенский и Булычев, 2011; Кузнецов и Булычев, 2017] предложено использование комплексных вейвлетов Пуассона произвольного порядка. Результаты вейвлет-анализа наземных, аэромагнитных и аэростатных съемок в районе Курской магнитной аномалии описаны в работе [Иванов и др., 2002]. Гауссовы вейвлеты подробно исследуются в работе [Хвастунов, 1998]. В работе [Глазнев и Муравина 2020] рассмотрены свойства вейвлета типа “мексиканская шляпа” и предложено использовать в качестве материнской функции решение прямой задачи для данного типа тел. В работе [Sailhac et al., 2000] показано, как использовать комплексные вейвлеты для интерпретации аэромагнитных данных.

Предварительно нами проведено исследование возникающих ошибок и ограничений на моделях тел простой формы и предложены методы, позволяющие снизить уровень ошибок. Интерпретации магнитных аномалий над телами простой формы рассматривается в обширной геолого-геофизической литературе по разведочной геофизике, см., например [Gay, 1963; Захаров и Логачев, 1979], поскольку форма этих тел составляет значительный процент известных геологических и минеральных проявлений (дайки, жилы, разломы, контакты и слоистые породы). Кроме того, тела простой формы и их комбинации лежат в основе более сложных моделей, используемых при решении задач фундаментальных геофизики, одной из которых является интерпретация морских магнитных аномалий. Разработанные методические подходы были протестированы на спрединговой модели, а затем применены к реальному магнитометрическому профилю, пересекающему хребет Рейкьянес. Результаты вейвлет-анализа на модельных и реальных данных мы сопоставили с результатами применения метода аналитического сигнала.

2. ПРИМЕНЕНИЕ ВЕЙВЛЕТ-ПРЕОБРАЗОВАНИЯ К МОДЕЛЬНЫМ ПОЛЯМ

Вейвлет-преобразование W исходного сигнала f относительно вейвлет-функции ψ выражается формулой

$W\left(a,b\right)=\frac{1}{{\left|a\right|}^{1/2}}\underset{-\infty }{\overset{\infty }{}}f\left(x\right){\psi }^{*}\left(\frac{x-b}{a}\right)dx,$

где f(x) – анализируемый сигнал, a – масштаб вейвлета, b – сдвиг, ψ – порождающая функция. За счет a изменяется размер отрезка задания порождающей функции, а b задает локализацию порождающей функции вдоль анализируемого сигнала. Одно из наиболее важных свойств вейвлетов состоит в том, что они обеспечивают частотно-временную локализацию структурных особенностей анализируемого сигнала и позволяют проследить их изменение во времени.

К настоящему времени разработано большое количество семейств вейвлетов на основе разных ψ функций, обладающих широким спектром свойств, которые определяют их применение в той или иной задаче. Так, например, те вейвлет-семейства, которые имеют обратное преобразование, часто используются для фильтрации сигналов (см., например, [Воскобойников, 2015]). Для анализа спектрального состава и его изменения во времени полезным является непрерывный вейвлет Морли. Наша задача состоит в локализации особенностей исходного сигнала, связанных со структурными неоднородностями его источника. В работах [Кузнецов и Булычев 2017; Оболенский и Булычев 2011; Кузнецов и др., 2015] для этой цели был использован комплексный вейвлет Пуассона. Мы в данной работе предлагаем использовать семейство вейвлетов, построенных на производных функции Гаусса, с точностью до константы, выражающихся как

$\psi \left(x\right)={\left(-1\right)}^{m-1}\frac{d}{d{x}^m}\left[\text{exp}\left(\frac{-x^2}{2}\right)\right],$

где m – порядок производной. Выбор этого вейвлета обусловлен основой задачей интерпретации морских магнитных аномалий, состоящей в локализации границ блоков океанической коры с разной намагниченностью, что как раз и обеспечивается выбранным семейством вейвлетов. Действительно, подробный анализ свойств вейвлетов, построенных на производных функции Гаусса [Хвастунов, 2002], показал, что для сигнала по форме близкого к функции Гаусса, экстремумы его вейвлет-спектра, вычисленного при разных m, совпадают с экстремумами соответствующих производных сигнала. Мы использовали вейвлеты с m = 1,2,3. По модулю вейвлет-коэффициентов для каждого масштаба определяются координаты локальных максимумов, которые образуют вертикальные линии, параллельные оси масштабов. На первом этапе мы применили такой подход к анализу поля простых телологических образований.

2.1. Тела простой формы

Наиболее важное для интерпретации магнитных аномалий свойство вейвлет-преобразования с использованием производных функции Гаусса состоит в том, что координаты экстремумов этих вейвлетов совпадают с координатами экстремумов соответствующих производных наблюденного поля. Эти экстремумы определяются свойствами геологического объекта, и в том случае, когда для этих точек могут быть получены аналитические выражения, по найденным координатам экстремумов вейвлет-спектров и, предполагая форму тела, можно найти все его параметры. Для набора тел простой формы нами получены соответствующие формулы, которые могут быть использованы для определения их геометрических параметров (все необходимые формулы приведены в Приложении). Мы показываем это на примере нескольких тел простой формы. В работе рассматриваются только линейные структуры, бесконечного простирания по координате y. Сечения рассматриваемых тел показаны на рис. 1.

Рис. 1. Геометрия и параметры 2D геофизических моделей. (а) – квадрант, полубесконечное по вертикали и горизонтали тело; (б) ‒ пласт, полубесконечный по горизонтали слой; (в) ‒ дайка; (г) ‒ блок. 

2.1.1. Квадрант

Конфигурация, которая может рассматриваться как основной элемент в геометрии любого тела с прямоугольными границами, параллельными осям координат, – это тело с сечением в виде квадранта с вершиной в точке (x₀, z₁) (см. рис. 1а). Поле, создаваемое этим объектом, определяется формулой [Talwani, Heirtzler, 1964]

$B_z=\text{arctg}\left((x-x_0)/z_1\right)$.

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

По вейвлет-спектру мы можем определить положение экстремумов производных разного порядка. У квадранта два параметра: x₀ – координата угловой точки сечения и z₁ – глубина до верхней кромки (далее, просто глубина). Согласно (П1) экстремум первой производной достигается в угловой точке x₀. А два экстремума второй производной с координатами x_± =x₀±√3/3z₁ (П2) позволяют оценить глубину z₁. Третья производная имеет три экстремума: центральный совпадает с положением угловой точки, расстояние от двух боковых до центра равно глубине (П3). Пример определения параметров квадранта с заданными параметрами x₀ = 2 км и z₁ = 3 км показан на рис. 2. На рисунке 2 x-координата задана в км, масштабы пересчитаны в км, сплошными линиями показаны изолинии модулей вейвлет-коэффициентов для трех производных. Для каждой производной точками обозначены координаты локальных максимумов. Хорошо видно, что эти точки выстраиваются в прямые, параллельные оси масштабов. По пересечениям этих прямых с осью x были определены параметры квадранта. В результате по 1-й производной (П1) получено значение x₀ = 2.01 км (рис. 2а), по 2-й производной (П2) – x₀ = 2.005 км и z₁ = 3.003 км (рис. 2б), по 3-й производной (П3) – x₀ = 2.0 км и z₁ = 3.01 км. (рис. 2в).

Рис. 2. Вейвлет-спектры поля квадранта для трех функций Гаусса. (а) – m = 1, (б) – m = 2 и (в) – m = 3, соответственно. Непрерывные линии – изолинии коэффициентов вейвлет-спектра, точки – положения локальных максимумов при фиксированном масштабе.

Точность определения x координаты вершины квадранта, практически, не зависит от его глубины. Так, в модельном примере глубина менялась от 0.5 до 4 км, при этом координата угла x₀ определялась с ошибкой не большей шага задания поля. Глубины, полученные по экстремумам 2-й и 3-й производных, совпадают с точностью до первых метров.

2.1.2. Полубесконечный по горизонтали пласт

Простейшим производным от квадранта телом можно считать полубесконечный по горизонтали пласт (см. рис. 1б). Такое тело можно считать составленным из двух квадрантов с намагниченностью разного знака, и оно описывается 3-мя параметрами: 2 координаты верхней угловой точки (x₀ и z₁) и толщина h=z₂-z₁. Поле, создаваемое таким телом, выражается формулой [Никитский, Глебовский, 1990]

$B_z=\text{arctg}\left((x-x_0)/z_1\right)-\text{arctg}\left(\right(x-x_0)/z_2)$.

Один из экстремумов 1-й производной совпадает с x₀ (координатой угловой точки). Два других экстремума $x_{\pm }^1$ (здесь нижние индексы – корни, а верхний – порядок производной) определяются нулями 2-й производной (П4).

В численных экспериментах получено, что ошибка определения координат x₀ по всем производным также, как и в случае квадранта, не превышает 10 м независимо от толщины. Что касается глубины и толщины пласта, то в предельных случаях, когда z₂≫z₁ задача сводится к предыдущей, а координаты дополнительных экстремумов стремятся к ∞. Можно найти оценку глубины z₁ по экстремуму 3-й производной, приняв, что в случае пласта применима формула (П3) для квадранта. Тогда, зная положение экстремумов 1-й производной $x_{\pm }^1$, получаем уравнение для определения z₂ (П5). Возможность применения этого приема для определения параметров пласта была исследована эмпирически. Расчеты проводились для трех значений глубины z₁ и восьми значений толщины h. Результаты приведены на рис. 3. На рисунке 3а разными символами показана зависимость ошибки определения глубины от толщины пласта (от 0.5 км ‒ ● до 8 км ‒ □). Получено, что при глубине 1 км эта ошибка мало зависит от толщины пласта и не превышает по абсолютной величине 200 м. Но с ростом глубины толщина пласта оказывает решающее значение на ошибку определения глубины. Для глубины 4 км и толщины 0.5 км ошибка составляет уже 940 м. Ошибка определения толщины пласта показана на рис. 3б. Символами показана зависимость от глубины пласта. Хорошо видно, что для глубины 1км (█) ошибка определения толщины мало зависит от заданной толщины пласта и не превышает 200 м. В то же время для тонкого слоя ошибка быстро растет с глубиной и составляет для z = 4 км (◆) 1200 м. Если перейти к безразмерным единицам, поделив и толщину, и ошибку ее определения на глубину (рис. 3в), то получаем, что относительная ошибка зависит только от относительной глубины и при h/z₁≫1 пренебрежимо мала.

Рис. 3. Ошибки определения глубины и толщины пласта по 3-й производной. (а) ‒ ошибка определения глубины, символами показана зависимость от толщины h в км: h = 0.5 ‒ ●, h = 1 ‒ ◆, h = 2 ‒ ▲, h = 3 ‒ █, h = 4 ‒ □, h = 6 ‒ ○, h = 7 ‒ △, h = 8 ‒ □; (б) ‒ ошибка определения толщины, символами показана зависимость от z₁ в км: z₁ = 1 ‒ ●, z₁ = 2 ‒ ▲, z₁ = 4 ‒◆; (в) ‒ относительная ошибка определения толщины в безразмерных единицах.

2.1.3. Дайка

Дайка или вертикальный пласт заданной ширины характеризуется также тремя параметрами: центр x₀, полуширина d и глубина z₁ до верхней кромки (см. рис. 1в). Сечение S имеет вид S = {x₀–d<x< x₀+d, z < −z₁}. Создаваемое таким телом поле выражается формулой [Никитский, Глебовский, 1990]

$B_z=\text{arctg}\left(\frac{x-x_0-d}{z_{1 }}\right)-\text{arctg}\left(\frac{x-x_0+d}{z_1}\right)$.

Этот вид геологического объекта особенно важен для нас при поиске экскурсов в структуре морских магнитных аномалий. Для дайки аналитически могут быть получены только экстремумы 1-й производной (П6). Что дает положение центра. Если рассматривать дайку как комбинацию двух разнесенных квадрантов, то координаты ее угловых точек могут быть определены по положению экстремумов 1-й и 3-й производных, что дает полуширину. Тогда глубина может быть вычислена из (П6). В случае дайки, экстремумы 3-й производной дают более точные значения координат угловых точек. Заметим, что глубина залегания магнитоактивного слоя в случае интерпретации морских магнитометрических данных часто может быть получена из батиметрии.

Эмпирически исследовался важный для интерпретации вопрос: при какой ширине дайки погрешность определения ее параметров тем же способом, что и для квадранта, пренебрежимо мала. Изменяя полуширину и глубину дайки и определяя точки экстремумов 1-й и 3-й производных, получаем оценки значений искомых параметров. Другими словами, мы полагаем, что

а) экстремумы 1-й и 3-й производных соответствуют угловым точкам дайки;

б) экстремумы 3-й производной определяют ее глубину по формуле (П3) для квадранта.

Ширина дайки принимала значения 0.5, 1, 1.5 и 2 км, а глубина изменялась от 1 до 8 км. Ошибки определения параметров дайки для всех возможных комбинаций заданных параметров показаны на рис. 4а и 4б. Важным аспектом является то, что рассматривались только узкие тела. Можно видеть, что глубина, оцененная по экстремумам 3-й производной (рис. 4а), получается несколько завышенной, но, практически, не зависит от ширины дайки. В то же время ошибка определения полуширины (рис. 4б) быстро растет с ростом глубины. Только в том случае, когда полуширина вдвое больше глубины, ошибка становится мала. Это связано с тем, что найденные по вейвлет-спектру координаты экстремумов и 1-й, и 3-й производных для дайки не совпадают с ее угловыми точками, как это было для квадранта, и расхождение растет с увеличением глубины. При большой глубине дайка представляется значительно шире, чем она есть на самом деле. Однако положение центра определяется с малой ошибкой, как и в случаях угловых точек квадранта и пласта.

Рис. 4. Ошибка определения глубины и полуширины дайки по экстремумам 3-й производной. (а) – ошибка определения z₁; (б) – ошибка определения полуширины; (в) – то же, что (а) после коррекции; (г) – тоже, что (б) после коррекции; символы на (а) и (в) соответствуют разным полуширинам, на (б) и (г) – глубинам дайки.

Чтобы уточнить значение полуширины дайки d, можно использовать явное выражение для координат $x_{\pm }^1$экстремумов ее 1-й производной (П6). При этом полагаем z₁ и x₀ известными. После простых преобразований из (П6) получаем формулу для определения полуширины. Обозначив $x_{\pm }^1$ x₀≡q, запишем уравнение для полуширины как

$d=\sqrt{2q\sqrt{q^2+z_1^2}-q^2-z_1^2} .$ (1)

Из этой формулы следует, что расстояние экстремума до центра x₀ больше полуширины. Если глубина меньше полуширины, ε = z₁/d<1, то формулу (1) можно приблизительно записать как d = q(1-ε⁴/8). Если значение z₁ известно из независимых источников, по формуле (1) можно получить значения полуширины с ошибкой, не превышающей сто метров в самом неблагоприятном случае глубокой узкой дайки. Однако, если значение z₁ было определено по экстремумам 3-й производной (рис. 4а), то ошибка резко возрастает, более того, в ряде случаев в качестве полуширины получаются мнимые значения.

Оценки полуширины и глубины для дайки могут быть получены более точно, если усложнить задачу и для их определения привлечь экстремумы 2-й производной. При этом мы получаем систему двух уравнений с двумя неизвестными d и z₁, x₀ считаем известным. Для получения первого уравнения приравнивалась к нулю 2-я производная (П7) при x, равном найденному экстремуму 1-й производной (П6). Для получения второго уравнения приравнивалась к нулю 3-я производная (П8) при x, равном найденному экстремуму 2-й производной (для простоты в формуле центр дайки имеет нулевую координату). Система решалась нами численно. Полученные из решения системы значения глубин и полуширин показаны на рис. 4в и 4г. Можно видеть, что при применении такого подхода погрешность может быть уменьшена в несколько раз. Однако в некоторых случаях решение системы неустойчиво и требует точного определения положения экстремумов, что напрямую зависит от шага задания исходного сигнала. В случае рассматриваемой модели для больших глубин потребовалось знание исходного сигнала с шагом в несколько метров.

Ранее в работе [Иванов и Меркурьев, 2014] была рассмотрена задача определения параметров дайки с использованием метода аналитического сигнала. Это позволяет на примере дайки сравнить интерпретационные возможности двух методов: вейвлет-анализа и метода аналитического сигнала. Для сравнения приведем результаты определения полуширины дайки к безразмерным единицам, поделив и полуширину, и ошибку ее определения на глубину дайки. На рис. 5 показан результат, полученный экспериментально для вышеприведенных значений параметров дайки и при использовании 1-й (треугольники) и 3-й (круги) производных. Для адекватного сравнения была вычислена теоретическая зависимость относительной ошибки вейвлет-анализа с использованием 1-й производной, на рис. 5 она показана сплошной линией. Аналогичная кривая для метода аналитического сигнала на рис. 5 показана пунктиром. Теоретические кривые для 3-й производной не вычислялись, в силу отсутствия для положения экстремумов аналитических выражений. Приведенные результаты наглядно демонстрируют, что в случае дайки вейвлет-анализ дает меньшую ошибку определения параметров по сравнению с методом аналитического сигнала. Кроме того, хорошо видно преимущество использования 3-й производной для определения полуширины дайки по сравнению с 1-й.

Рис. 5. Сравнение интерпретационных возможностей вейвлет-метода и метода аналитического сигнала. Сплошной линией показана теоретическая зависимость для вейвлет-метода, пунктиром – для аналитического сигнала, треугольники – значения, полученные по экстремумам 1-й производной, круги – то же по экстремумам 3-й производной.

2.2. Модель трех блоков

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

Рассмотрим в качестве модели три последовательных блока, находящихся на одной глубине z₁ = 3 км, имеющих толщину h = 0.4 км, отличающихся полушириной и намагниченностью. Положим полуширины d равными 4.5, 0.5 и 5 км. Такую модель можно рассматривать как элемент океанического магнитоактивного слоя. Поле, создаваемое этой моделью, и его производные приведены на рис. 6.

Рис. 6. Поле, создаваемое структурой из 3-х блоков, его производные и результат интерпретации. Жирной черной линией показано поле, длинный пунктир – 1-я производная, короткий пунктир – 2-я производная, сплошная серая линия – 3-я производная; (а) – заданная структура, (б) – первоначальная оценка, (в) – после коррекции.

На рис. 6 в блоке (а) показана заданная структура магнитоактивного слоя. Анализ вейвлет-спектров показал, оценки глубины до верхней кромки для всех блоков оказываются заниженными и колеблются от 2.7 до 2.85 км. Этот результат определяется малой толщиной слоя и соответствует оценкам, полученным ранее для модели слоя. Наилучшее совпадение с границами можно извлечь из положения экстремумов 3-й производной. Результат представлен на рис. 6 в блоке (б), при этом хорошо видно, что узкий блок представляется уширенным. Этот результат соответствует полученному выше для узкой дайки. Там же было показано, что уточнить полуширину можно на основе формулы (1). Параметром в этой формуле является положение центральной точки x₀. Оценить положение центров блоков можно 3 способами: по положению экстремумов поля и по расстоянию между границами, определенными по 1-й и 3-й производным. Результаты представлены в табл. 1.

Таблица 1. Координаты центров блоков до и после коррекции

Можно видеть, что положение центра маленького блока определяется с малой ошибкой любым способом. Для крупных блоков определение центра по экстремумам поля и 3-й производной дают схожие результаты. Далее мы использовали положения центров, полученные по экстремумам 3-й производной, и по ним оценивали полуширины блоков. Значения заданных полуширин и полученных из вейвлет-спектра приведены в табл. 2 в соответствующих строках. Очевидно, что требуется коррекция не только для узкого блока. Для решения этой задачи применим методику, описанную выше для дайки. Будем считать известной глубину до верхней кромки, а за положения центров примем координаты, полученные по 3-й производной. Тогда, применив формулу (1) и решив полученные уравнения, вычислим новые значения полуширин блоков (шаг 1). Результаты, полученные на шаге 1, вынесены в табл. 2.

Таблица 2. Полуширины блоков до и после коррекции

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

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

3. ПРИМЕНЕНИЕ ВЕЙВЛЕТ-АНАЛИЗА ДЛЯ ИНТЕРПРЕТАЦИИ ОКЕАНИЧЕСКИХ МАГНИТНЫХ АНОМАЛИЙ

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

Теоретические вопросы возможности разрешения океанических магнитных аномалий с помощью метода аналитического сигнала были рассмотрены нами ранее в ряде статей [Иванов и Меркурьев, 2014; Иванов и Меркурьев, 2016], где на простейшем примере вертикальной дайки было аналитически и численно показано, при каких соотношениях глубины залегания и ширины источника возникают трудности определения его положения.

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

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

3.1. Аналитический сигнал и его свойства

Напомним определение и основные свойства аналитического сигнала.

Аналитический сигнал определяется как квадратный корень из суммы квадратов вертикальной и двух горизонтальных производных аномалий магнитного поля ΔT [Roest et all., 1992].

$\left|AS\right|=\sqrt{{\left(\frac{\Delta T}{\partial x}\right)}^2+{\left(\frac{\Delta T}{\partial y}\right)}^2+{\left(\frac{\Delta T}{\partial z}\right)}^2}.$

Здесь производная по z есть преобразование Гильберта горизонтальных производных [Nabighian, 1972], что позволяет строить модуль аналитического сигнала |AS(x,y)| по заданной магнитной аномалии ΔT(x, y).

Метод аналитического сигнала полезен для определения местоположения вертикальных границ между блоками противоположной магнитной полярности океанической коры, например, [Roest, et all., 1992]. По мнению некоторых авторов, этот метод гораздо более надежен и объективен, чем традиционный метод выбора путем визуального сравнения наблюдаемых и синтетических моделей аномалий [Chaubey, 2002]. Согласно опыта других исследователей, метод не улучшил способность идентифицировать хроны разной полярности, в качестве возможной причины приводятся сложная геометрия и вращение намагниченных тел в осевой зоне медленно-спредингового хребта [Parnell-Turner et al., 2016].

Алгоритм и программа, с помощью которой проводились расчеты спрединговой модели и аналитического сигнала, использовали быстрое преобразование Фурье [Трошков и Шалаев, 1961; Schouten and McCamy, 1972], что способствовало быстроте вычислений.

3.2. Тестовый пример для спрединговой модели

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

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

В данной работе в качестве тестового примера был рассчитан теоретический магнитный профиль от спрединговой модели, которая была нами ранее использована при геохронологическом анализе магнитных аномалий, образовавшихся на границе Северо-Американской и Евразийской плит за последние 20 млн лет [Merkouriev and DeMets, 2008]. При расчете мы использовали временную шкалу инверсий [Gradstein et al., 2020], полная скорость спрединга 20 мм/год, толщина слоя 0.4 км, глубина залегания верхней кромки слоя 2 км. На рис. 7а показана структура магнитоактивного слоя в виде блоков прямой и обратной полярности. Блоки прямой полярности закрашены, отмечены границы хронов, которые необходимо определить по магнитным аномалиям с помощью вейвлет-анализа. Сплошной линией показан график теоретических магнитных аномалий, рассчитанных от этого магнитоактивного слоя.

Рис. 7. Сравнение двух методов определения границ источников аномалий для спрединговой модели. (а) –магнитные аномалии, рассчитанные от спрединговой модели, темными прямоугольниками показаны блоки прямой полярности, буквенно-цифровыми символами – границы хронов; (б) – аналитический сигнал; (в) – сопоставление положения границ хронов: отрезки с ▲ – искомые границы тел, отрезки с ▼ – найденные с помощью вейвлет-анализа, отрезки с ◆ – найденные с помощью аналитического сигнала; (г) – вейвлет-преобразование, изолиниями показан модуль вейвлет-коэффициентов, прерывистая линия – локальные экстремумы.

К данному профилю был применен вейвлет-анализ с порождающей функцией, являющейся 3-й производной от функции Гаусса. Одновременно для сравнения проводился анализ методом аналитического сигнала. На рис. 7б показан график модуля аналитического сигнала, рассчитанного от магнитных аномалий. Локальные максимумы аналитического сигнала, согласно его свойствам, расположены над или вблизи границ аномалеобразующих тел. Результат вейвлет-анализа приведен на рис. 7г. Сплошными линиями показаны изолинии модуля вейвлет-коэффициентов, точками ‒ локальные максимумы на каждом масштабе, масштабы пересчитаны в км. Сравнение полученных результатов вынесено в центральный блок (в). Здесь повторена модельная шкала, длинными вертикальными отрезками показаны положения экстремумов 3-й производной вейвлет-преобразования, перевернутыми треугольниками ▼ отмечены экстремумы, наиболее близко расположенные к искомым границам хронов, треугольниками ▲ отмечены сами границы, а короткими отрезками с ромбом ◆ ‒ показаны границы, определенные с помощью метода аналитического сигнала.

Как видно из рис. 7в, наблюдается довольно неплохое соответствие искомых и найденных границ обоими методами. Численные оценки отклонения найденных координат от искомых показывают, что модуль среднего отклонения, которое дает вейвлет-анализ, составляет 0.33 км, тогда как эта величина для метода аналитического сигнала равна 1.2 км. Кроме того, с помощью аналитического сигнала было выделено значительно меньшее количество границ блоков и, в частности, не удалось определить границу хрона 4n.1, тогда как с помощью вейвлет-анализа границы всех искомых хронов были определены.

3.3. Пример использования вейвлет-анализа на реальном магнитометрическом профиле, пересекающем хребет Рейкьянес

Описанную выше методологию мы применили к реальным данным для определения границ блоков аномалеобразующих тел по магнитометрическому профилю, пересекающему срединно-океанический хребет (хребет Рейкьянес). Такая задача является типичной при проведении геоисторического и кинематического анализа аномального магнитного поля, когда по магнитным аномалиям необходимо определить границы одновозрастных хронов, расположенных по обеим сторонам от оси хребта для последующего их совмещения с целью определения кинематических параметров вращения литосферных плит. В качестве такого профиля был выбран профиль KNOR24, пересекающий хребет Рейкьянес, расположенный в Северной Атлантике южнее Исландии (рис. 8). Для того, чтобы представить геологическую обстановку этого региона, напомним, что хребет Рейкьянес является дивергентной границей Северо-Американской и Евразийской плит, которые расходятся с полной скоростью спрединга около 20 мм/год [Merkouriev and DeMets, 2008].

Рис. 8. Магнитометрический профиль, пересекающий хребет Рейкъянес в Северной Атлантике. (а) – рельеф дна; (б) – магнитные аномалии, измеренные на профиле KNOR24, символы ‒ хроны; цифры – номера основных аномалий; (в) – теневая карта рельефа дна с вынесенными положениями изохрон и наблюденных магнитных аномалий на магнитометрическом профиле KNOR24, символы те же, что на (б).

Для хребта Рейкьянес характерны высокоинтенсивные аномалии с амплитудой 300‒800 нТл и высокоинтенсивной осевой аномалией с амплитудой до 2500 нТл, приуроченной к оси хребта (рис. 8б). Рельеф океанического дна довольно изрезанный, особенно в осевой зоне хребта, которая представляет собой грабен с глубиной около 1000 м и шириной около 70 км, который переходит в плавно спадающий рельеф, характерный для склонов срединно-океанических хребтов (рис. 8а).

На рис. 8в показана теневая карта рельефа дна с вынесенными положениями изохрон и графиком наблюденных магнитных аномалий на магнитометрическом профиле KNOR24, полученном из базы данных NGDC [Hey, 2013]. Изохроны, показанные на рис. 8в различными маркерами, были получены нами в ходе исследования магнитометрических данных на границе Северо-Американской и Евразийской плит с целью построения кинематической модели высокого разрешения [Merkouriev and DeMets, 2008]. Всего были оцифрованы границы 21 инверсии в диапазоне возрастов 0.72‒20 млн лет (аномалии 1n – 6n), средний интервал между хронами составил 1 млн лет. Оцифровка аномалий осуществлялась в ручном режиме с помощью программы, позволяющей путем визуального сопоставления модельного и наблюденного профиля выбирать точку на графике, соответствующую положению (границе) каждого хрона. Отметим, что точка, соответствующая границе хрона, у каждой аномалии не всегда приходится на участок графика, пересекающего горизонтальную нулевую ось или в месте максимального градиента графика поля, поэтому оцифровку аномалий не всегда удается свести к формальной процедуре, которую можно запрограммировать.

С помощью вейвлет-анализа и метода аналитического сигнала по описанной выше методике были определены границы аномалеобразующих тел. Результаты анализа реальных данных приведены на рис. 9, который построен так же, как и рис. 7. График анализируемого магнитного профиля показан на рис. 9а, на который вынесены границы хронов, оцифрованные нами в ручном режиме [Merkouriev and DeMets, 2008; DeMets and Merkouriev, 2020]. Положение этих хронов мы рассматриваем как искомые границы источников, которые надо определить с помощью вейвлет-анализа и метода аналитического сигнала. В этом же блоке ниже показан магнитоактивный слой спрединговой модели, которую мы использовали ранее для тестирования. Напомним, что спрединговая модель была нами рассчитана при постоянной скорости спрединга 20 мм/год и изображенная здесь шкала служит для того, чтобы иметь возможность сопоставлять хроны, оцифрованные на наблюденном профиле, с соответствующими хронами на шкале.

Рис. 9. Сравнение двух методов определения границ источников магнитных аномалий для профиля KNOR24. (а) – наблюденные на профиле магнитные аномалии, темными прямоугольниками показаны блоки прямой полярности, буквенно-цифровыми символами – границы хронов; (б) – аналитический сигнал; (в) – сопоставление положения границ хронов; (г) – вейвлет-преобразование. Условные обозначения те же, словные обозначения те же, что для рис. 7.

На рис. 9б показан график амплитуды аналитического сигнала, вычисленный по магнитным аномалиям, наблюденным на профиле KNOR24. По локальным максимумам аналитического сигнала были определены границы источников, которые вынесены в виде черных отрезков с ромбами ◆ на рис. 9в. Границы, определенные по локальным максимумам, найденным с помощью вейвлет-анализа с порождающей функцией ‒ 3-я производная от функции Гаусса (рис. 9г) показана на рис. 9в перевернутыми треугольниками ▼. На этом же рисунке маркером ▲ показаны искомые границы хронов.

Как видно из рис. 9в, также как и для спрединговой модели, наблюдается довольно неплохое соответствие искомых и найденных границ, полученных с помощью вейвлет-анализа. Полученные для реального профиля численные оценки модуля отклонения найденных координат от искомых составили для вейвлет-анализ 0.25 ± 0.2 км, тогда как для метода аналитического сигнала – 1.5 ± 1.1 км.

4. ЗАКЛЮЧЕНИЕ

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

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

1. Даны рабочие формулы, связывающие положение экстремумов производных магнитных полей этих линейных структур с их параметрами.
2. Показано, что, используя положения экстремумов вейвлет-преобразования с производными функции Гауса разного порядка, можно определить геометрические параметры перечисленных выше тел простой формы.
3. Проведен анализ возникающих ошибок и даны оценки точности получаемых результатов.
4. Численно показано, что вейвлет-преобразование позволяет определять координаты контактов тем точнее, чем ближе поля в окрестности контакта к полю тела с сечением в виде квадранта.
5. Предложены методы уточнения параметров тел при использовании комбинации вейвлет-преобразований с производными разного порядка от функции Гаусса. Показано, что задача сводится к решению системы уравнений, степень которых зависит от порядка производных. При этом для обеспечения устойчивости получаемого решения требуется большая точность определения координат экстремумов вейвлет-спектра.
6. На примере определения параметров дайки аналитически и численно показано, что вейвлет-преобразования дает меньшую ошибку определения параметров по сравнению с методом аналитического сигнала.
7. В качестве примера структуры из нескольких тел детально рассмотрена модель с сечением в виде 3-х блоков. Получено, что в случае присутствия в модели узкого блока, вейвлет-преобразование позволяет установить его наличие в структуре. Уточнение координат границ и ширин блоков требует дополнительных методических разработок.
8. Показана работоспособность метода на примере модели Вайна-Метьюза, охватывающей последние 20 миллионов лет. Численные оценки показали, что модуль среднего отклонения найденных координат границ хронов от искомых по вейвлет-анализу составляет 0.33 км, тогда как эта величина для метода аналитического сигнала равна 1.2 км. Кроме того, с помощью метода аналитического сигнала было выделено значительно меньшее количество границ блоков, тогда как с помощью вейвлет-анализа были определены границы всех искомых хронов.
9. Метод применен к данным реального профиля, пересекающего хребет Рейкьянес. Получено, что численные оценки средних отклонений координат хронов по вейвлет-анализу составляет 0.25±0.2 км, тогда как эта величина для метода аналитического сигнала равна 1.5±1.1 км. Это показало преимущество метода вейвлет-преобразования перед методом аналитического сигнала в задаче определения контактов для реальных данных.

ФИНАНСИРОВАНИЕ РАБОТЫ

Работа выполнена в рамках государственного бюджета.

КОНФЛИКТ ИНТЕРЕСОВ

Авторы заявляют об отсутствии конфликта интересов.

ПРИЛОЖЕНИЕ

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

1.Квадрант

2. Пласт

3. Дайка

4. Блок

About the authors

S. А. Merkuriev

Author for correspondence.
Email: sam_hg@hotmail.com

St. Petersburg Branch

S. А. Ivanov

Email: sergei.a.ivanov@mail.ru

St. Petersburg Branch

I. М. Demina

Email: dim@izmiran.spb.ru

St. Petersburg Branch

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Fig. 1. Geometry and parameters of 2D geophysical models. (a) - quadrant, vertically and horizontally semi-infinite body; (b) - layer, horizontally semi-infinite layer; (c) - dike; (d) - block.

Download (70KB)

Indexing metadata

3. Fig. 2. Wavelet spectra of the quadrant field for three Gaussian functions. (a) - m = 1, (b) - m = 2 and (c) - m = 3, respectively. Continuous lines are isolines of wavelet spectrum coefficients, dots are positions of local maxima at fixed scale.

Download (166KB)

Indexing metadata

4. Fig. 3. Errors of determining the depth and thickness of the reservoir by the 3rd derivative. (a) - depth determination error, symbols show the dependence on thickness h in km: h = 0. 5 - ●, h = 1 - ◆, h = 2 - ▲, h = 3 - ”, h = 4 - □, h = 6 - ○, h = 7 - △, h = 8 - □; (b) - thickness determination error, symbols show dependence on z1 in km: z1 = 1 - ●, z1 = 2 - ▲, z1 = 4 -◆; (c) - relative thickness determination error in dimensionless units.

Download (112KB)

Indexing metadata

5. Fig. 4. Error in determining the depth and half-width of the dike from the extrema of the 3rd derivative. (a) - error of z1 determination; (b) - error of half-width determination; (c) - the same as (a) after correction; (d) - the same as (b) after correction; symbols on (a) and (c) correspond to different half-widths, on (b) and (d) - dike depths.

Download (189KB)

Indexing metadata

6. Fig. 5.Comparison of interpretation capabilities of the wavelet method and the analytical signal method.The solid line shows the theoretical dependence for the wavelet method, the dotted line - for the analytical signal, triangles - values obtained by the extrema of the 1st derivative, circles - the same by the extrema of the 3rd derivative.

Download (70KB)

Indexing metadata

7. Fig. 6.The field created by the structure of 3 blocks, its derivatives and the result of interpretation.Bold black line shows the field, long dashed line - 1st derivative, short dashed line - 2nd derivative, solid gray line - 3rd derivative; (a) - given structure, (b) - initial estimation, (c) - after correction.

Download (209KB)

Indexing metadata

8. Fig. 7.Comparison of two methods for determining the anomaly source boundaries for the spreading model.(a) -magnetic anomalies calculated from the spreading model, dark rectangles show the blocks of direct polarity, alphanumeric symbols - chron boundaries; (b) - analytical signal; (c) - comparison of the position of chron boundaries: segments with ▲ are the sought boundaries of the bodies, segments with ▼ are those found by wavelet analysis, segments with ◆ are those found by analytical signal; (d) - wavelet transform, isolines show the modulus of wavelet coefficients, dashed line - local extrema.

Download (480KB)

Indexing metadata

9. Fig. 8.Magnetometer profile crossing the Reykjanes Ridge in the North Atlantic. (a) - bottom relief; (b) - magnetic anomalies measured on the KNOR24 profile, symbols - chrons; numbers - numbers of the main anomalies; (c) - shadow map of the bottom relief with the positions of isochrons and observed magnetic anomalies on the KNOR24 magnetometer profile, symbols are the same as in (b).

Download (507KB)

Indexing metadata

10. Fig. 9. Comparison of two methods of determining the boundaries of magnetic anomaly sources for the KNOR24 profile. (a) - magnetic anomalies observed on the profile, dark rectangles show blocks of direct polarity, alphanumeric symbols - chron boundaries; (b) - analytical signal; (c) - comparison of the chron boundary positions; (d) - wavelet transformation.

Download (428KB)

Indexing metadata

In the print version, the article was published under the DOI: 10.31857/S0016794025020097