Физика Земли, 2020, № 3, стр. 98-111

ВВЕДЕНИЕ

При интерпретации гравиметрических и магнитометрических данных, вычислении поправок за рельеф, построении карт уклонений отвесной линии и т.п. возникает необходимость нахождения аналитических аппроксимаций функций, заданных приближенно на нерегулярной сети точек. В настоящее время предлагаются разнообразные методики решения проблемы выбора аппроксимационных конструкций (см., например, [Костицын и др., 2013; 2017; Мартышко и др., 2003; Кризский, 2011]. Для построения аналитических аппроксимаций рельефа местности и геопотециальных полей можно применить описанный в работах [Страхов, Степанова, 2002а; 2002б] метод S-аппроксимаций, а также улучшенный вариант этого метода – модифицированные S-аппроксимации [Степанова и др., 2017; 2018]. Методика построения S-аппроксимаций основана на применении фундаментальной формулы теории гармонических функций (третьей формулы Грина) [Кошляков и др., 1962]. В настоящей статье будет рассматриваться локальный вариант с декартовой системой координат и вертикальной осью, направленной вверх.

Если известны компоненты магнитного или гравитационного поля (например, первая производная потенциала по x₃ на некотором рельефе), то можно представить потенциал поля в виде суммы потенциалов простого и двойного слоев, распределенных на нескольких горизонтальных плоскостях, которые залегают ниже заданного рельефа. Если систему координат выбрать так, чтобы одна из плоскостей простого и двойного слоев задавалась уравнением х₃ = 0 (однако данное условие не является обязательным), то потенциал представляется в виде:

В формуле (1) и далее использованы следующие обозначения: $M = ({{x}_{1}},{{x}_{2}},{{x}_{3}})$ – точка измерений; $\,{{\xi }_{1}},\;{{\xi }_{2}}$ – точки на соответствующих поверхностях; ${{H}_{r}}$ – глубина залегания r-ой горизонтальной плоскости (расстояние от точки на плоскости до дневной поверхности) в рамках метода S-аппроксимаций; ${{\rho }_{{1,r}}}\left( {\hat {\xi }} \right),$ ${{\rho }_{{2,r}}}\left( {\hat {\xi }} \right) \in {{L}_{2}}\left( { - \infty , + \infty } \right)$ – неизвестные функции, подлежащие определению, $\hat {\xi } = ({{\xi }_{1}},{{\xi }_{2}}).$ Тогда производная по х₃ потенциала V, взятая с обратным знаком, будет иметь вид:

Пусть компоненты поля заданы в конечном множестве точек М_i, M_i = $({{x}_{{1,i}}},{{x}_{{2,i}}},{{x}_{{3,i}}})$, i = 1, 2, …, N; обозначим известную часть подынтегральной функции в первом слагаемом для r-го слоя в (2) в точке М_i через $Q_{{1,r}}^{{(i)}}$, а во втором слагаемом – через $Q_{{2,r}}^{{(i)}}$. Тогда получим:

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

получим, что искомые функции должны иметь вид [Страхов, Степанова, 2002а]:

Таким образом, приходим к следующей системе линейных алгебраических уравнений (СЛАУ):

элементы матрицы которой в нашем случае имеют вид:

Коэффициенты а_ij матрицы А могут быть вычислены явно с помощью интеграла Пуассона:

В формулах (6)–(9) элементы матрицы системы представляют собой сумму выражений, зависящих от R-носителей – либо горизонтальных плоскостей, либо горизонтальных слоев, если в формуле (8) слагаемые проинтегрировать по вертикальной координате (в последнем случае, формула (9), разумеется, несколько изменится, но принцип построения аналитических аппроксимаций останется прежним). На практике требуется разделять поля, создаваемые различными носителями как гравитационных, так и магнитных масс. Это необходимо, в первую очередь, для выявления глубоко расположенных геологических структур. Поля, генерируемые приповерхностными, зачастую хаотически рассеянными источниками быстро затухают при аналитическом продолжении вверх. Однако, если имеются сравнимые по мощности массы, залегающие на разных горизонтах, то вычисление линейных трансформант от всего поля не позволит, в общем случае, идентифицировать источники. Модифицированные S-аппроксимации вместе с усовершенствованной методикой решения СЛАУ являются весьма эффективным инструментом решения описанной проблемы изучения глубинного строения земной коры. Следует отметить, что в упомянутых выше работах авторов речь шла исключительно о нахождении решения системы линейных алгебраических уравнений, общего для всех носителей. В настоящей работе подчеркивается роль структурно-параметрического подхода к решению вариационных задач [Страхов, Степанова, 2002а]: для каждого из носителей задается свой вектор решений; при этом элементы матрицы выглядят совершенно аналогично (9), только без знака суммы по R-носителям. Такой способ решения аппроксимационной задачи позволяет свободнее “маневрировать” при деформации блоков матрицы системы, что приводит к повышению качества решения.

ОПИСАНИЕ МОДИФИЦИРОВАННОГО МЕТОДА S-АППРОКСИМАЦИЙ

При применении “классического” метода S-аппроксимаций теряется существенная часть информации о геометрии носителя масс: о его гладкости, симметрии, приблизительной локализации. Постановка обратной задачи по поиску распределения плотности в классе равномерно непрерывных или равномерно дифференцируемых функций значительно усложнит алгоритм, что приведет к неизбежным ошибкам при его использовании. Поэтому одной из проблем является выбор такой постановки обратной задачи в рамках S-аппроксимации, в которой сочетались бы простота и наиболее полный учет априорной информации об исследуемой области. При построении аппроксимационной конструкции поля, заданного в узлах неравномерной сети точек с резкими перепадами высот, обусловленность матрицы СЛАУ (7) сильно возрастает. При решении системы с плохо обусловленной матрицей S-аппроксимации становятся физически неадекватны и приводят к бессмысленным с геофизической точки зрения результатам.

Одним из способов устранения указанных проблем является введение особого функционала, позволяющего уточнять некоторые пробные решения СЛАУ (7) в зависимости от имеющейся априорной информации. Основная конструкция S-аппроксимации на начальном этапе сохраняется, поэтому назовем такой метод модифицированным методом S-аппроксимаций [Степанова и др., 2018].

Не ограничивая общности, рассмотрим локальный вариант S-аппроксимации, когда сферичностью Земли можно пренебречь. Пусть задано N точек наблюдения ${{M}_{i}} = \left( {{{x}_{{1,i}}},{{x}_{{2,i}}},{{x}_{{3,i}}}} \right)$, i = 1, 2, …, N, в которых известны значения некоторого элемента поля ${{f}_{{i,\delta }}} = {{f}_{i}} + \delta {{f}_{i}}$, i = 1, 2, …, N. Если рассматривается вертикальная производная поля, то аппроксимация поля суммой простого и двойного слоев, расположенных на R-уровнях, задается соотношением (2). Далее ставится следующая вариационная задача:

где: $\rho = \left( {{{\rho }_{{k,r}}}\left( {\hat {\xi }} \right)} \right),$ $k = 1,2$; r = 1, 2, …, R – неизвестные функции; ${{\rho }_{{1,r}}}\left( {\hat {\xi }} \right) \in {{L}_{2}}( - \infty , + \infty )$ – плотность r‑го простого слоя, ${{{\rho }}_{{2,r}}}\left( {{\hat {\xi }}} \right) \in {{L}_{2}}( - \infty , + \infty )$ – плотность r-го двойного слоя, а подынтегральные функции $Q_{{1,r}}^{{(i)}}\left( {\hat {\xi }} \right)$ и $Q_{{2,r}}^{{(i)}}\left( {\hat {\xi }} \right)$ имеют тот же вид, что и в формуле (2). Отличие такой постановки задачи от аналогичной ей (4)–(5) состоит в наличии функционала качества решения Ξ(ρ) и параметра регуляризации μ. При μ = 0 приходим к классическому варианту метода S-аппроксимаций. В зависимости от функционала качества решения (который может характеризовать информацию об искомых геологических объектах: геометрию источников поля, их локализацию, среднюю плотность и т.д.) с помощью параметра регуляризации будет учитываться степень участия априорной информации в подборе качества решения. Граничные условия для квадрата нормы невязки, соответствующие различным значениям параметра регуляризации μ и ограниченных как сверху, так и снизу, запишем в следующем виде:

Для каждого значения параметра μ находится свое приближенное решение СЛАУ ${{\hat {x}}^{\mu }}$ (10), для которого выполняется условие (11). Нижнее и верхнее значения квадрата нормы невязки ${{\left( {\delta _{{\min }}^{\mu }} \right)}^{2}}$ и ${{\left( {\delta _{{\max }}^{\mu }} \right)}^{2}}$ корректируются для различных значений параметра μ. Делается это с той целью, что постановка обратной задачи (10) накладывает определенные ограничения на выбор решения системы. Из теории катастроф известно, что даже незначительное изменение функции ведет к изменению ее особенностей [Постон и др., 1980; Арнольд, 1990]. С физической точки зрения это приводит к тому, что меняется поведение всей системы. Поэтому для каждого последующего значения μ постоянные ${{\left( {\delta _{{\min }}^{\mu }} \right)}^{2}}$ и ${{\left( {\delta _{{\max }}^{\mu }} \right)}^{2}}$ необходимо корректировать.

На начальном этапе решения полагаем параметр μ = 0. В этом случае мы приходим к классическому варианту метода S-аппроксимаций, в котором система описывается соотношениями (6)–(9). Таким образом, находится решение $\hat {x}$, при котором плотности простого и двойного слоев описываются следующими соотношениями:

Из формулы (12) следует, что аппроксимации функций ${{\rho }_{{k,r}}}\left( {\hat {\xi }} \right)$ зависят от найденного решения системы. Для наиболее быстрого решения системы с приемлемой точностью, следует выбрать константы ${{\left( {\delta _{{\min }}^{0}} \right)}^{2}}$ и ${{\left( {\delta _{{\max }}^{0}} \right)}^{2}}$ из следующих соображений:

То есть задавая точность аппроксимации ε, определяем значение констант ${{\left( {\delta _{{\min }}^{0}} \right)}^{2}}$ и ${{\left( {\delta _{{\max }}^{0}} \right)}^{2}}$. Константы c и $\varepsilon $ задаются из условий поставленной задачи и зависят от уровня помех.

Так как выбор функционала качества решения Ξ(ρ) зависит от имеющейся априорной информации, то выписать аналитические выражения для безусловной вариационной задачи (10) не представляется возможным. Однако, если известно решение системы (12) при μ = 0, то, увеличивая параметр регуляризации μ постепенно от 0 до 1, можно определенным образом варьировать вектор решения ${{\hat {x}}^{0}}$ СЛАУ, чтобы аппроксимации искомых функций удовлетворяли поставленным условиям. Не ограничивая общности, будем считать, что значения параметра регуляризации μ меняются по линейному закону:

где s соответствует определенному решению. Так, при s = 0 решением задачи (10) являются выражения (12). Аппроксимации искомых функций $\rho _{{k,r}}^{{({{a}^{\mu }})}}\left( {\hat {\xi }} \right)$ соответствующих разным значениям μ, можно представить в виде однопараметрического семейства функций:

Как было отмечено, функции ${{\tilde {\rho }}_{{k,r}}}\left( {\hat {\xi },{{{\hat {x}}}^{{{{\mu }_{s}}}}},{{\mu }_{s}}} \right)$ при различных значениях μ могут по-разному отображать физико-геологические параметры среды. Запишем выражение для s-й деформации аппроксимированной функции:

где ${{\varepsilon }_{\mu }} = {{\mu }_{s}} - {{\mu }_{{s - 1}}}$, а ${{f}_{{k,r}}}\left( {\hat {\xi },{{\mu }_{s}}} \right) \in {{L}_{2}}\left( { - \infty , + \infty } \right)$ соответствует s-й деформационной функции аппроксимации ${{\tilde {\rho }}_{{k,r}}}\left( {\hat {\xi }{\kern 1pt} ,{{{\hat {x}}}^{{{{\mu }_{{s - 1}}}}}},{{\mu }_{{s - 1}}}} \right)$. Согласно (14) ε_μ = 0.01 = const. Перепишем (16) в следующем виде:

При таком представлении функция ${{\tilde {\rho }}_{{k,r}}}\left( {\hat {\xi }{\kern 1pt} ,{{{\hat {x}}}^{0}},0} \right)$ играет роль нулевого приближения для однопараметрического семейства (15), в котором каждая функция является решением безусловной вариационной задачи (10). Поэтому под s-й деформацией функции подразумевается деформация именно нулевого приближения ${{\tilde {\rho }}_{{k,r}}}\left( {\hat {\xi }{\kern 1pt} ,{{{\hat {x}}}^{0}},0} \right)$, $k = 1,2$; r = 1, 2, …, R.

Целесообразным представляется выбор такого функционала качества решения, который давал бы наиболее содержательные результаты для поставленной задачи. Например, если требуется построить аналитическое продолжение поля вверх для выделения регионального фона [Мартышко и др., 2003] вместе с вычислением высших производных поля, то функционал качества решения для каждой из поставленных задач следует выбирать отдельно.

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

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

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

Кроме того, можно ввести дополнительную переменную ${{\xi }_{4}}$ и перейти в пространство двух комплексных переменных: ${{z}_{1}} = {{\xi }_{1}} + i{{\xi }_{2}},$ ${{z}_{2}} = {{\xi }_{3}} + i{{\xi }_{4}}$. Тогда требование плюригармоничности функции (т.е. решение должно удовлетворять условию $\frac{{{{\partial }^{2}}u}}{{\partial {{z}_{i}}\partial {{{\bar {z}}}_{j}}}} = 0,$ $i = 1,2;$ $j = 1,2.$), представляющей собой значения магнитного поля в заданной области, позволяет ограничить класс возможных решений обратной задачи.

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

Если требуется найти аналитическое продолжение поля в нижнее или верхнее полупространство, то дополнительно целесообразным представляется выделение области $\tilde {D}$, в которой наблюдается наименьшая интенсивность поля. После такого качественного анализа, функционал качества решения можно выбрать следующим образом:

Региональный вариант. Для аналитического продолжения полей вверх или вниз функционал Ξ(ρ), определяемый выражением (20), остается без изменений. В том случае, когда находятся высшие производные поля, выражение (19) принимает следующий вид:

где: ${{D}_{r}}$ – r-я область, представляющая из себя сферу определенного радиуса, на которой распределены простой и двойной слои; $\tilde {\vartheta }$ – ко-широта; $\tilde {\varphi }$ – долгота.

АНАЛИТИЧЕСКАЯ АППРОКСИМАЦИЯ МАГНИТНОГО ПОЛЯ ПО ДАННЫМ ДЕТАЛЬНОЙ ГИДРОМАГНИТНОЙ СЪЕМКИ

Данные гидромагнитной съемки озера Кинерет (рис. 1) были проинтерпретированы нами с помощью усовершенствованного метода блочного контрастирования. Модифицированные S-аппроксимации были построены с использованием структурно-параметрического подхода, как это было описано выше. Шаг сетки $h = 25$ м, общее количество валидных измерений равняется 264 442. В предыдущих работах авторов уже исследовался описываемый регион [Степанова и др., 2018]. Однако ранее размеры блоков в возникающей системе уравнений и глубины залегания плоскостей, на которых залегали носители магнитных масс, не варьировались. В настоящей статье предлагается алгоритм деформации блока, позволяющий выявлять наиболее “значимые” для формирования локальных аномалий геологические структуры.

Исходя из графического анализа аномального магнитного поля можно видеть, что исходная область характеризуется как крупными региональными, так и локальными аномалиями. Исходные данные были получены в результате многопрофильной гидромагнитной съемки, затем были интерполированы на равномерную сетку с шагом $h = 25$ м. Интерполирование данных эксперимента привело к появлению весьма специфических “дефектов” аномального магнитного поля, которые заглушают локальные аномалии исходного поля. Для демонстрации подобного рода “дефектов” на рис. 2 приведена северная часть аномального магнитного поля озера Кинерет.

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

В рамках модифицированной S-аппроксимации поле аппроксимировалось суммой простого и двойного слоев, распределенных на 4 семействах плоскостей, залегающих на глубинах 250 + 5t, 600 + 10t, 1100 + 15t и 2000 + 20t м соответственно ниже высотной отметки, на которой была проведена съемка (для данной гидромагнитной съемки она равняется –220 м); t – параметр, $t \in \left[ {0;5} \right],\;\Delta t = 0.1$ – шаг изменения параметра. Константы, являющиеся границами квадрата нормы невязки снизу и сверху, выбраны следующим образом: $\delta _{{\min }}^{2} = 6000$ нТл², $\delta _{{\max }}^{2} = 30\,000$ нТл². Размеры блоков деформировались по формуле

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

В формуле (21) дискретный индекс i связан с непрерывным параметром t линейной зависимостью:

При построении аналитических аппроксимаций применяется способ трехступенчатого контроля, предложенный В.Н. Страховым, который заключается в следующем. На первом шаге из исходных пунктов исключается около 20% пунктов, имеющих минимальные по модулю значения аномального поля, которые включаются в число контрольных, а оставшиеся 80% пунктов используются для построения аппроксимационной конструкции. На втором шаге из контрольных пунктов выделяется половина, имеющая максимальные по модулю отклонения от исходного поля, и они добавляются к числу пунктов, используемых для построения аппроксимации. Контроль точности производится по оставшимся контрольным пунктам (10%). На третьем шаге для построения аппроксимации используются все исходные пункты.

В табл. 1 приводятся результаты S-аппроксимации. Система линейных алгебраических уравнений (СЛАУ) решалась регуляризованным трехслойным итерационным методом Чебышева с использованием блочного метода контрастирования. Все вычисления проводились на суперкомпьютере “Ломоносов” НИВЦ МГУ. При вычислениях число процессоров выбрано равным 1024.

Приняты следующие обозначения: N – число исходных точек, N₁ – число точек без контрольных точек I типа; N₂ – число точек без контрольных точек II типа; ${{\sigma }_{{\min }}} = \sqrt {\frac{{\delta _{{\min }}^{2}}}{N}} ;$ ${{\sigma }_{{\max }}} = \sqrt {\frac{{\delta _{{\max }}^{2}}}{N}} $ – границы среднеквадратического отклонения, ${{\sigma }_{0}} = \frac{{\left\| {Ax - {{f}_{\delta }}} \right\|_{E}^{{}}}}{{\sqrt N }}$ – среднеквадратическое отклонение; ${{\sigma }_{1}}$ – среднеквадратичное отклонение, полученное с помощью аппроксимации аномального магнитного поля без контрольных точек I типа; ${{\sigma }_{{{\text{конт}}\,1}}}$ – среднеквадратичное отклонение, полученное в контрольных точках I типа; ${{\sigma }_{2}}$ – среднеквадратичное отклонение, полученное с помощью аппроксимации аномального магнитного поля без контрольных точек II типа; ${{\sigma }_{{{\text{конт}}\,2}}}$ – среднеквадратичное отклонение, полученное в контрольных точках II типа; ${\Delta } = \frac{{{{{\left\| {Ax - {{f}_{\delta }}} \right\|}}_{E}}}}{{{{{\left\| {{{f}_{\delta }}} \right\|}}_{E}}}}$ – показатель качества решения; T – время решения СЛАУ (ч : мин : с); роль контрольных точек второго типа играли дефектные зоны. Таким образом, с помощью предлагаемой методики можно выделить поля, создаваемые различными по своему вкладу в результирующее магнитное поле структурами. Основной вклад в поле дефектных зон (68% приблизительно) вносили массы, распределенные на самой верхней плоскости, как показали наши расчеты.

Время решения СЛАУ для данной модели составило всего 5 мин на персональном компьютере. Значение ${{\sigma }_{{\min }}} < {{\sigma }_{0}} < {{\sigma }_{{\max }}}$ и среднеквадратические отклонения ${{\sigma }_{{{\text{конт}}\,1}}}$ и ${{\sigma }_{{{\text{конт}}\,2}}}$, полученные в контрольных точка I и II типа, превышают ${{\sigma }_{0}}$ почти в 3 раза, что связано с наличием вышеуказанных “дефектов” в исходных интерполированных данных. На карте разностного магнитного поля (рис. 3) также выделяются аналогичные кривые, в которых разность между аппроксимацией и реальным полем максимальна. Во второй строке таблицы приводятся значения среднеквадратических отклонений, времени счета и т.д. для деформированного блока, т.е. число элементов исходной матрицы системы было уменьшено в соответствии с формулами (21)–(22). Было выполнено 10 серий расчетов с различными значениями параметра t в формуле (22). При уменьшении размерности блока на 25–35% существенных изменений в погрешности метода не происходит (при редукции блока на 50% среднеквадратическое отклонение ${{\sigma }_{0}}$ возросло на 15–18%), что позволяет сделать вывод об эффективности усовершенствованного метода блочного контрастирования при решении такого рода задач.

РАЗДЕЛЕНИЕ МАГНИТНЫХ ПОЛЕЙ В МОДЕЛЬНОМ ПРИМЕРЕ

В соответствии с изложенным во Введении к настоящей работе для нахождения производных геопотенциального поля в некоторой совокупности точек ${{M}_{\nu }}$ с координатами ${{x}^{{(\nu )}}} = ({{\hat {x}}^{{(\nu )}}},{{x}_{3}}^{{(\nu )}})$, из внешней области СD объема D в рамках метода линейных интегральных представлений необходимо действовать следующим образом.

По заданным плотностям ${{\rho }_{1}}({{\xi }_{1}},{{\xi }_{2}}) \in {{L}_{2}}(D)$ и ${{\rho }_{2}}({{\xi }_{1}},{{\xi }_{2}}) \in {{L}_{2}}(D)$ простого и двойного слоев (обозначим плоскость, на которой распределены простой и двойной слои через D) мы должны найти совокупность ограниченных линейных функционалов вида

где для всех s:

Действительно, значения аномалий силы тяжести $g({{M}_{\nu }};\rho )$ выражаются соотношением (23) при s = ν:

а значения элементов аномального поля $u({{M}_{\nu }};\rho )$ соотношением вида: где $\alpha = ({{\alpha }_{1}},{{\alpha }_{2}},{{\alpha }_{3}})$ – мультииндекс; |α| = α₁ + $ + \,\,{{\alpha }_{2}} + {{\alpha }_{3}} \geqslant 1;$и соотношением (23) при s = ν:

В соотношении (23) ${{\rho }_{1}}({{\xi }_{1}},{{\xi }_{2}}),\,\,{{\rho }_{2}}({{\xi }_{1}},{{\xi }_{2}})$ находятся по заданным значениям ${{g}_{{k,\delta }}}$ аномалий силы тяжести из решения системы уравнений, к которой сводится вариационная задача (10), т.е. системы вида (7)–(9).

Аналогичные выражения можно получить для магнитного потенциала и его производных, если учесть соотношение:

Здесь: $\bar {I}$ – вектор постоянной намагниченности; U – магнитный потенциал; $\gamma ,\;\rho $ – гравитационная постоянная и средняя плотность породы соответственно; V – гравитационный потенциал. Поэтому в дальнейшем, без ограничения общности, производные геопотенциальных полей вычислялись по формуле (28).

Поскольку ${{\rho }_{1}},{{\rho }_{2}}$ представляется в виде (12), то из (28) следует, что

где λ = ${{({{\lambda }_{1}},{{\lambda }_{2}}, \ldots ,{{\lambda }_{N}})}^{T}},$ ${{b}^{k}}(P_{k}^{{(s)}}) = (b_{1}^{k}(P_{k}^{{(s)}}),$ $b_{2}^{k}(P_{k}^{{(s)}}), \ldots ,b_{N}^{k}(P_{k}^{{(s)}}){{)}^{T}};$

Вектора ${{b}^{k}}(P_{k}^{{(s)}})$ вычисляются из соотношений:

Если D-горизонтальная плоскость ${{\varsigma }_{3}} = {\text{const}}$, то

где $\rho _{l}^{2}(\hat {x}) = {{({{\hat {x}}_{{1,l}}} - {{\hat {x}}_{1}})}^{2}} + {{({{\hat {x}}_{{2,l}}} - {{\hat {x}}_{2}})}^{2}}.$

Например, если требуется вычислить первую производную силы тяжести по x₁ в точке М_i, то элементы матрицы ${{a}_{{ij}}}$ нужно продифференцировать по х₁.

Получим:

Здесь $ru = {{({{x}_{{1,i}}} - {{x}_{{1,j}}})}^{2}}$ + ${{({{x}_{{2,i}}} - {{x}_{{2,j}}})}^{2}},$ R = $ = {{[{{({{x}_{{3,i}}} + {{x}_{{3,j}}})}^{2}} + ru]}^{{1/2}}}.$

Затем необходимо матрицу ${{A}_{{{{x}_{1}}}}} = \{ {{a}_{{ij,{{x}_{1}}}}}\} $ умножить на вектор $\lambda $. Получим значения производной ${{g}_{{{{x}_{1}}}}}$ в совокупности точек i = 1, 2, …, S.

Аналогично вычисляются производные ${{g}_{{}}}$ по ${{x}_{2}},{{x}_{3}}$ и вторые производные.

Опишем модельный участок, на котором была апробирована новая методика. Общее количество расчетных гравиметрических пунктов составляет 6000, сеть нерегулярная. Рельеф земной поверхности горный, перепад высот составляет около 3 км. Для данной сети было также рассчитано модельное гравитационное поле, обусловленное 14 прямоугольными призмами, расположенными двумя группами по 7 призм в каждой. Одна группа призм залегала на глубине 7 км, другая –15 км. В целом, данный участок соответствует модельному примеру № 2, подробно описанному в работе [Страхов, Степанова, 2002б].

Для вычисления производных гравитационного поля ${{g}_{{{{x}_{1}}}}},{{g}_{{{{x}_{2}}}}},{{g}_{{{{x}_{3}}}}},$ создаваемого группой призматических тел, находились значения поля в близких по каждой координате точках, и затем разность этих значений делилась на приращение.

С целью проверки эффективности разработанной методики была вычислена сумма первых производных гравитационного поля ${{g}_{{{{x}_{1}}}}} + {{g}_{{{{x}_{2}}}}},$ а затем построена модифицированная S-аппроксимация этой суммы с помощью усовершенствованного метода блочного контрастирования. При этом элементы матрицы системы вычислялись следующим образом: в формуле (33) ${{x}_{1}}$ заменилось на ${{x}_{2}}$, затем находилась сумма значений элементов ${{a}_{{ij,{{x}_{1}}}}} + {{a}_{{ij,{{x}_{2}}}}}$ в каждой точке. Результаты разделения полей (т.е. значений каждой из производных в отдельности) представлены на рис. 4–рис. 8. При вычислении значений конкретной производной вектор решения системы линейных алгебраических уравнений умножался на соответствующую часть ${{a}_{{ij,{{x}_{k}}}}},k = 1,2.$

Относительная ошибка (показатель качества решения ${\Delta } = {{{{{\left\| {Ax - {{f}_{\delta }}} \right\|}}_{E}}} \mathord{\left/ {\vphantom {{{{{\left\| {Ax - {{f}_{\delta }}} \right\|}}_{E}}} {{{{\left\| {{{f}_{\delta }}} \right\|}}_{E}}}}} \right. \kern-0em} {{{{\left\| {{{f}_{\delta }}} \right\|}}_{E}}}}$) не превысила 10% в данном модельном примере.

Для дальнейшего улучшения качества “разделенных” решений была построена аналитическая аппроксимация элемента потенциального поля (т.е. суммы производных ${{g}_{{{{x}_{1}}}}} + {{g}_{{{{x}_{2}}}}}$) полем источников, распределенных по граням двух двугранных углов (Ed.1 и Ed.2):

Выяснилось, что поле первой группы призм с точностью 10–12% восстанавливается с помощью первой конструкции, а поле второй (более глубоко залегающей и расположенной правее первой) – с помощью второго двугранного угла. Время решения системы линейных алгебраических уравнений на Pentium –IV (3 ГГц , ОЗУ -8 Гбайт) составило 17 мин при применении метода усовершенствованного блочного контрастирования.

ВЫВОДЫ

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

2. Усовершенствованный блочный метод контрастирования позволяет найти такое нулевое приближение для СЛАУ, что время решения всей системы существенно сокращается (см. табл. 1), если исследуемое поле не характеризуется малым количеством четко выраженных экстремумов. В противном случае заметного улучшения при использовании блочного метода контрастирования не наблюдается. Однако на практике исследуемое поле всегда характеризуется достаточным количеством геологических объектов, искажающих поле и усложняющих его картину, так что использование блочного метода оправдано при интерпретации данных детальной аэромагнитной съемки (см. табл. 1) и рис. 3–рис. 6.

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

4. На примерах съемки озера Кинерет показано, что при помощи модифицированного метода S-аппроксимаций восстанавливаются первые производные аномального магнитного поля.