Данная утилита решает задачу поиска эквивалентного перераспределения физического параметра (плотности, намагниченности, магнитной восприимчивости) по данным геопотенциального поля. Исходными данными является сеточная 2D модель поля на горизонтальной плоскости, результатом - 3D конечноэлементная модель под плоскостью, задаваемая 3D сеткой. Используемый метод аналогичен методу Приезжева, но содержит ряд отличий, главное из которых состоит в учёте дискретной природы моделей и возможности параметризации метода, позволяющей подбирать результат (модель) под априорную информацию (всю другую геолого-геофизическую информацию о моделируемом объекте). Утилита является обратной для утилиты Прямая задача (геопотенциальные поля) (ДПФ) в том смысле, что решение прямой задачи от результата решения обратной даёт исходное поле с машинной точностью (но с некоторыми оговорками, касающимися экстраполяции, и описанными ниже).

Простое описание

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

Множество эквивалентных решений само по себе очень большое - например, если исходная сетка имеет размеры ${\displaystyle n_x\times n_y}$, а строящаяся 3Д сеть - ${\displaystyle n_x\times n_y\times n_z}$, то 3Д сеть будет иметь очень много степеней свободы - ${\displaystyle n_x\times n_y\times (n_z-1)}$. Эти степени свободы "фиксируются" следующим способом: для заданной 3D сети в каждой её точке задаётся значение функции ${\displaystyle \psi }$, и в соответствии с изложенным ниже алгоритмом значения функции трансформируются так, что решение прямой задачи от результата этой трансформации равняется изначальному полю ${\displaystyle u}$. Подробнее метод изложен в разделе #Математическая формулировка.

В модуле реализовано два способа задания ${\displaystyle \psi }$: явно и аналитически. Аналитический способ задания - это математическая функция ${\displaystyle \psi =\psi (x,y,z)}$; обычно это полосовой фильтр. Важно отметтиь, что решение не будет полем, пропущенным через этот фильтр: алгоритм в общем случае изменит ${\displaystyle \psi }$ перед "фильтрацией".

Явный способ задания заключается в том, что пользователь самостоятельно задаёт функцию ${\displaystyle \psi }$ как данные сетки, обычно с помощью внешних по отношению к утилите инструментов, а утилита "исправляет" данные так, чтобы "исправленная" версия ${\displaystyle g_{\psi }(u)}$ являлась эквивалентной по заданному полю ${\displaystyle A(g_{\psi }(u))=u}$.

Функцию ${\displaystyle \psi }$ можно интерпретировать как первое приближение модели физического параметра. Если функция ${\displaystyle \psi }$ и правда является эквивалентным распределением физического параметра, то алгоритм не изменит её; если же ${\displaystyle \psi }$ близко к решению (в некотором смысле), то алгоритм изменит её очень мало (то есть, ${\displaystyle g_{\psi }(u)-\psi }$ будет мало отличаться от ${\displaystyle \psi }$, или ${\displaystyle |g_{\psi }(u)-\psi |\to 0}$ при ${\displaystyle |A\psi -u|\to 0}$).

Решение ищется в классе моделей эффективного физического параметра, с точностью до константы. То есть, считается, что к каждому горизонтальному срезу модели с номером ${\displaystyle i}$ можно прибавить константу ${\displaystyle C_i}$, и решение фактически не изменится. Абсолютные значения физического параметра следует назначать вне этой утилиты - например, суммированием с градиентной моделью (с нулевым эффектом).

Математическая формулировка

Математическая формулировка метода опирается на определения, данные на странице Прямая задача (геопотенциальные поля) (ДПФ). Для дискретных моделей поля ${\displaystyle \hat{u}}$ и распределения физического параметра ${\displaystyle \hat{g}}$ (причём ${\displaystyle {\hat{g}}_i}$ - горизонтальный ${\displaystyle i}$-ый срез, ${\displaystyle i\in 1..n_z}$), и горизонтального (по ${\displaystyle x}$ и ${\displaystyle y}$ координатам) Фурье-преобразования ${\displaystyle {ℱ}}$ было выведено соотношение:

${\displaystyle \hat{u}=\sum _{i=1..n_z}{{ℱ}}^{-1}Diag\left[{ℱ}{\hat{k}}_i\right]({ℱ}{\hat{g}}_i)}$. (1)

Обратная задача заключается в нахождении таких срезов ${\displaystyle {\hat{g}}_i}$ для данного поля ${\displaystyle \hat{u}}$, что соотношение (1) - верно. Заметим, что если ${\displaystyle n_z=1}$, то решение находится однозначно:

${\displaystyle {\hat{g}}_i={{ℱ}}^{-1}Dia{g}^{-1}\left[{ℱ}{\hat{k}}_i\right]({ℱ}{\hat{u}}_i)}$. (2)

При этом, естественно, ни один коэффициент вектора ${\displaystyle {ℱ}{\hat{k}}_i}$ не равняется нулю (исходя из аналитического вида ${\displaystyle k_i}$). Но сразу же заметим, что, поскольку решение ищется с точностью по постоянного слагаемого ${\displaystyle C_i}$, то коэффициент Фурье на нулевой частоте игнорируется.

Если ${\displaystyle n_z>1}$, то определим поле ${\displaystyle {\hat{u}}_i}$ от каждого среза ${\displaystyle {\hat{g}}_i}$ результирующей модели:

${\displaystyle {\hat{u}}_i={{ℱ}}^{-1}Diag\left[{ℱ}{\hat{k}}_i\right]({ℱ}{\hat{g}}_i)}$. (3)

Исходное поле должно равняться сумме полей от срезов результата:

${\displaystyle \hat{u}=\sum _{i=1..n_z}{\hat{u}}_i}$. (4)

Тогда решение обратной задачи можно сформулировать как распределение поля по "срезам" в соответствии с (4), после чего по каждому срезу задача решается однозначно по формуле (2).

Спектральный базис

Важным этапом является переход к спектральному базису. Обозначим модель поля ${\displaystyle \hat{U}={ℱ}\hat{u}}$, эквивалентное распределение физического параметра ${\displaystyle {\hat{G}}_i={ℱ}{\hat{g}}_i}$, ${\displaystyle {\hat{K}}_i={ℱ}{\hat{K}}_i}$. Найдя ${\displaystyle \hat{U}}$, легко получить ${\displaystyle \hat{u}={{ℱ}}^{-1}\hat{U}}$.

Так как ${\displaystyle {{ℱ}}^{-1}}$ - линейный оператор, то его можно вынести за сумму в (1): ${\displaystyle \hat{u}={{ℱ}}^{-1}\hat{U}=\sum _{i=1..n_z}Diag\left[{ℱ}{\hat{k}}_i\right]({ℱ}{\hat{g}}_i)}$. Тогда:

${\displaystyle \hat{U}=\sum _{i=1..n_z}Diag\left[{\hat{K}}_i\right]{\hat{G}}_i}$. (5)

Пусть в качестве начального приближения мы начали с некоторого распределения физического параметра ${\displaystyle \hat{\psi }}$, ${\displaystyle \hat{\mathrm{\Psi }}={ℱ}\hat{\psi }}$. Обозначим поле от такого распределения через ${\displaystyle \hat{V}}$:

${\displaystyle \hat{V}=\sum _{i=1..n_z}Diag\left[{\hat{K}}_i\right]{\hat{\mathrm{\Psi }}}_i}$. (6)

Центральная идея данного метода решения обратной задачи состоит в том, что целевое поле ${\displaystyle \hat{U}}$ и текущее ${\displaystyle \hat{V}}$ можно связать через вектор ${\displaystyle \hat{X}}$, представляющий собой мультипликативную невязку для начального приближения:

${\displaystyle \hat{U}=Diag[\hat{X}]\hat{V}}$. (7)

Вектор ${\displaystyle \hat{X}}$ определяется поэлементно. Для каждого номера элемента ${\displaystyle j}$:

1. Если элемент ${\displaystyle {\hat{V}}_j\ne 0}$, то ${\displaystyle {\hat{X}}_j=\frac{{\hat{U}}_j}{{\hat{V}}_j}}$.
2. Если элемент ${\displaystyle {\hat{V}}_j=0}$ и ${\displaystyle {\hat{U}}_j=0}$, то конкретное значение ${\displaystyle {\hat{X}}_j}$ не играет роли, но для определённости можно взять ${\displaystyle {\hat{X}}_j=0}$.
3. Если элемент ${\displaystyle {\hat{V}}_j=0}$, а ${\displaystyle {\hat{U}}_j\ne 0}$, то решение не определено (то есть, решение (5) может существовать, но для выбранного ${\displaystyle \hat{\mathrm{\Psi }}}$ оно - не определено).

В дальнейшем предположим, что условие (3.) никогда не выполняется. При этом, если условие (1.) выполняется всегда, то формула для вычисления ${\displaystyle \hat{X}}$ принимает следующий вид:

${\displaystyle \hat{X}=Dia{g}^{-1}\left[\sum _{i=1..n_z}Diag\left[{\hat{K}}_i\right]{\hat{\mathrm{\Psi }}}_i\right]\hat{U}}$. (8)

Обе стороны (6) умножаем на ${\displaystyle Diag[\hat{X}]}$ и используя (7) получаем:

${\displaystyle \hat{U}=Diag[\hat{X}]\sum _{i=1..n_z}Diag\left[{\hat{K}}_i\right]{\hat{\mathrm{\Psi }}}_i=\sum _{i=1..n_z}Diag\left[{\hat{K}}_i\right]Diag[\hat{X}]{\hat{\mathrm{\Psi }}}_i}$. (9)

Сравнив (9) с (5), легко заметить, что следующее выражение определяет эквивалентное распределением физического параметра:

${\displaystyle {\hat{G}}_i^{(\mathrm{\Psi })}=Diag[\hat{X}]{\hat{\mathrm{\Psi }}}_i}$. (10)

Адаптированный для ГИС ИНТЕГРО метод Приезжева состоит в следующим: начиная от некоторого распределения ${\displaystyle \psi }$, по формуле (8) расчитывается ${\displaystyle \hat{X}(\mathrm{\Psi })}$, после чего вычисляются все срезы эквивалентного распределения по формуле (10).

Корректность

Равенство нулю хотя бы одного элемента ${\displaystyle {\hat{V}}_j}$ при неравенстве нулю ${\displaystyle {\hat{U}}_j}$ приводит к тому, что ${\displaystyle \hat{X}}$ не существует. На практике это приводит к неустойчивости метода, когда значения некоторых ${\displaystyle {\hat{X}}_j}$ получаются очень большие, обычно - для высокочастотных составляющих. Если для выбранной геометрии и параметризации решение выглядит неустойчивым, то рекомендуется использовать регуляризацию.

Некоторые свойства

1. Результат не определён, если ${\displaystyle \psi }$ принадлежит ядру оператора прямой задачи A.
   1. Результат также не определён, если хотя бы один коэффициент ${\displaystyle \sum _{i=1..n_z}Diag\left[{\hat{K}}_i\right]{\hat{\mathrm{\Psi }}}_i}$ равен нулю.
2. Для того, чтобы определить результат, достаточно указать параметризацию - начальное приближение ${\displaystyle \mathrm{\Psi }}$, причём если его указывать в некотором смысле случайно, то, скорее всего, такая параметризация будет формально корректна.

Центральная спектральная проекция

Зафиксируем значения частот ${\displaystyle {\nu }_x}$ и ${\displaystyle {\nu }_y}$ (то есть, индекс ${\displaystyle j}$), отбросив все остальные значения. Тогда вектора ${\displaystyle \hat{X}}$, ${\displaystyle \hat{U}}$, ${\displaystyle \hat{V}}$ состоят из одного элемента, поэтому примем, что они - скаляры. В этом случае ${\displaystyle {\hat{G}}_i^{(\psi )}=\hat{X}{\hat{\mathrm{\Psi }}}_i}$, а ${\displaystyle {\hat{G}}^{(\psi )}=\hat{X}\hat{\mathrm{\Psi }}=\frac{\hat{U}}{\hat{V}}\hat{\mathrm{\Psi }}=\frac{\hat{U}}{{\hat{K}}^T\hat{\mathrm{\Psi }}}\hat{\mathrm{\Psi }}}$. Здесь можно узнать центральную проекцию: вектор ${\displaystyle {\hat{G}}^{(\psi )}}$ имеет то же направление, что и ${\displaystyle \hat{\mathrm{\Psi }}}$, но такую величину, что он находится на гиперплоскости эквивалентных решений, заданных уравнением ${\displaystyle {\hat{K}}^T{\hat{G}}^{(\psi )}=\hat{U}}$. Поэтому вектор ${\displaystyle \hat{\mathrm{\Psi }}}$ также иногда называется "направляющим" вектором. Но, к сожалению, этому сложно дать геолого-геофизическую интерпретацию.

Алгоритм

Алгоритм имеет простой вид:

1. Выбрать исходные данные; опционально экстраполировать исходные данные. Если данные содержат пропуски, заполнить их (например, с помощью функции экстраполяции).

2. Выбрать параметризацию ${\displaystyle \mathrm{\Psi }}$ - аналитическую (неявную) в виде одной из встроенных функций, или явную в виде свойства ТОС (первого приближения распределения физического параметра).

3. Рассчитать эквивалентное распределение физического параметра G через кнопку "Применить".

4. Осуществить проверку результата, его интерпретацию. В случае необходимости вернуться к шагу 1 или 2.

Особенности программной реализации

Программно метод основан на дискретных преобразованиях Фурье сеточной модели гравитационного потенциала, которые соответствуют непрерывным преобразованиям подлежащей непрерывной модели с известными оговорками. Одна из проблем - подразумевающаяся периодичность поля, которая выражается в виде краевых эффектов. Для устранения краевых эффектов рекомендуется выполнить экстраполяцию поля, подробнее см. [1] и Экстраполяции_поля. При этом для наибольшей вычислительной производительности рекомендуется выбирать такие величины экстраполяции, чтобы размеры сетки были равны степени двойки - например, 512 на 512, 512 на 1024 и др. В качестве альтернативы можно использовать экстраполяцию "отражение".

Другая особенность реализации состоит в том, что получаемая сеточная модель G полностью объясняет исходное поле. То есть, если на поле влияют какие-то другие факторы, - например, источники гравитационного поля вне сетки G, то будет найдено неадекватное распределение G, которое объяснит влияние и этих факторов.

В реализацию включена опциональная регуляризация, ограничивающая число обусловленности оператора эквивалентного перераспределения ${\displaystyle A_{\psi }^{-1}(u)}$.

Параметры

• Исходная ТОС 2D и исходное свойство - ТОС и имя свойства, содержащие исходное поле.
• Экстраполяция - Повтор, Специальная и Отражение: тип экстраполяции. Повтор: значение точек за краем данных равно значению точек с противоположного края. Отражение: значение точек за краем данных равно значению точек с этого же края с отступом, равным отступу за край сетки.
• Высота поля - высота, на которой измерено (на которую нормировано/пересчитано) поле, м.
• Регуляризация - выполнять ли регуляризацию оператора эквивалентного перераспределения физического параметра.
• Альфа - коэффициент регуляризации.
• Альфа, дБ - коэффициент регуляризации в логарифмических единицах.
• Аналитическая параметризация и Явная параметризация - вид направляющего вектора.
• Множитель по глубине, экспонента по Z и смещение по Z - коэффициенты аналитической параметризации.
• Распределить по глубинам в окне - дополнительное управление распределением физического параметра: Прямоугольное - распределить физический параметр строго в окне от верхней границы до нижней, задаваемых в Глубина верхней границы, м. и Глубина нижней границы, м.; Гауссиан - перераспределить физический параметр ближе к центру, задаваемым параметром Глубина центра, м., при этом разлёт контролируется параметром Ширина, м..
• Сохранить для аналитической параметризации - сохранить параметризацию как новое свойство в целевой 3D ТОС. При этом открывается новое диалоговое окно, в котором предлагается ввести имя свойство для записи в него параметризации. Сохранённое свойство можно использовать как явную параметризацию, в том числе после произведения различных трансформаций, например, в калькуляторе свойств.
• Имя свойства для явной параметризации - имя свойства в существующей целевой ТОС, содержащее явную параметризацию.
• Глубина верхней границы, глубина первого слоя, шаг по глубине, количество слоёв, глубина нижнего слоя и глубина нижней границы - параметры целевой ТОС 3D.
• Целевая ТОС 3D и целевое свойство - ТОС и имя свойства, в которые будет записано эквивалентное перераспределение. Если свойство с указанным именем уже существует, будет запрошено подтверждение за перезапись свойства. Предлагается несколько имён свойств на выбор, но можно ввести собственный вариант имени свойства.

Экстраполяция

Так как сеточная (конечноэлементная) модель - конечна, а уравнения физических законов определены на бесконечных плоскости и нижнем полупространстве, необходимо определить физические свойства в областях вне сетки. в ГИС INTEGRO это определяется через экстраполяцию.

• Экстраполяция повтором* - за границами сетки используются значения, взятые из данных сетки, отталкиваясь от её противоположной границы. Это - "естественный режим работы" для дискретных преобразований Фурье.

• Экстраполяция отражением* - за границами сетки используются значения, взятые из данных сетки, отталкиваясь от прилежащей границы.

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

Если сетка исходного поля содержит пропуски, заполнение пропусков (например, стандартной экстраполяцией) - обязательны.

Соответствие прямой и обратной задач

Важным вопросом может являться соответствие решения прямой и обратной задачи типа ${\displaystyle A{A}^{-1}u=u}$, то есть решение прямой задачи от результата обратной должно давать исходное поле. Это тождество соблюдается в следующих случаях:

1. Используется встроенная экстраполяция повтором
2. Используется встроенная экстраполяция отражением, задача - гравитационная или магнитная на полюсе (склонение = 0).
3. Используется людая другая экстраполяция - например, из утилиты Экстраполяции_поля, но область экстраполяции не удаляется из результата (то есть, результат решения обратной и прямой задач используется как есть, без удаления краёв, на которых осуществлялась экстраполяция).

В других случаях это тождество не соблюдается.

Параметризация (направляющий вектор)

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

Аналитический вид направляющего вектора

В случае аналитического вектора ${\displaystyle \mathrm{\Psi }}$ выбирается вид функции, а также параметры сетки, распределение физического параметра в которой, как предполагается, полностью объясняет исходное поле. От выбора сетки зависит окончательный результат. Сетка определяется: начальной глубиной, шагом по глубине, количеством слоёв. Другие параметры рассчитываются автоматически для наглядности.

На данный момент доступно 4 функции:

• "Функция гаусса второго порядка, латеральный вид": ${\displaystyle f(x,y,z)=(1-x^2-y^2)e^{(-x^2-y^2)/z^2}}$. Направляющий вектор - набор гауссианов с увеличивающимся диаметром по мере увеличения глубины ("юбка").

• "Функция гаусса второго порядка": то же самое, что и в предыдущем пункте, но вычисление коэффициентов направляющего вектора происходит в спектральной форме - без промежуточного этапа трансформации из латерального (пространственного) вида в спектральный вид с помощью дискретных преобразований Фурье. Основные отличия от предыдущего варианта - в скорости вычисления; отличия в эквивалентных перераспределениях - в основном, только для экстремальных значений параметров множителя по глубине и экспоненте при z.

• "Минимум ${\displaystyle L_2}$-нормы": ${\displaystyle f(x,y,z)=k(x,y,z)}$. Для значений параметров "множителя по глубине" и "экспоненте при z" по умолчанию даёт наименьшее распределение физического параметра среди класса эквивалентных. В непрерывной форме данной параметризации соответствует решения в классе единственности гармонических функций.

• "Единичный импульс": ${\displaystyle f(x,y,z)=1}$, если x=0 и y=0, иначе ${\displaystyle f(x,y,z)=0}$. В этом случае решение ${\displaystyle g(x,y,z)}$ не зависит от z: ${\displaystyle g(x,y,z)=g(x,y)}$.

Вид функции также параметризуется:

• ${\displaystyle \alpha }$ - множителем по глубине;

• ${\displaystyle \beta }$ - экспонентой при z;

• ${\displaystyle z_0}$ - смещение по z.

Если f(x,y,z) - латеральный вид выбранной функции, то в латеральном виде направляющий вектор ${\displaystyle \psi =f(x,y,\alpha (z_0+z))|z_0+z{|}^{-\beta }}$.

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

1. Обычно увеличение множителя по глубине стягивает решение наверх, но, в общем случае, эффект зависит как от характера поля, так и от выбранной геометрии сетки среды.
2. Увеличение экспоненты при z переносит массы (аномалообразующие значения физического параметра) наверх.
3. Увеличение и того, и другого коэффициентов по отдельности или вместе (переносе "основной массы" наверх) уменьшает число обусловленности (улучшает обусловленность, уменьшает некорректность). И наоборот, при переносе основной массы вниз решение может начать осциллировать (то есть, "идти волнами").
4. Увеличение значения ${\displaystyle \alpha }$ и уменьшение ${\displaystyle \beta }$ (относительно нуля; например, ${\displaystyle \alpha }$=5 и ${\displaystyle \beta }$=-3) переносит массы по глубине к середине результирующей модели (при условии, что для модели среды указан достаточно большой интервал глубин).

Явный вид направляющего вектора

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

Регуляризация

Коэффициент регуляризации ${\displaystyle \alpha }$ вводится пользователем (не путать с параметром аналитической параметризации - множителем по глубине ${\displaystyle \alpha }$). Регуляризация осуществляется отдельно от остальной параметризации. Задаётся в пределах от 0 до 1, при этом 0 - отсутствующая регуляризация, а при 1 распределение физического параметра = 0 во всех узлах сетки. Для удобства добавлено альтернативное поле для ввода коэффициента регуляризации - логарифмическое, единицы - дБ, при этом ${\displaystyle {\alpha }_{linear}=10^{{\alpha }_{logarithmic}/10}}$.

Интерпретация результата

Для гравитационной задачи результат определяет эквивалентное избыточное распределение физического параметра, единица измерения дана в Прямая задача (геопотенциальные поля) (ДПФ). Избыточность физического параметра обозначает, что его значения заданы относительно некоторого постоянного слагаемого ${\displaystyle C_i}$ в ${\displaystyle i}$-ом горизонтальном срезе. Эквивалентность распределения (в данном случае - по полю) обозначает, что полученное распределение имеет гравитационный эффект, равный исходному с машинной точностью (вычислений и представления чисел).

Слагаемое ${\displaystyle C_i}$ может происходить из градиентной модели с нулевым аномальным полем. То есть, распределение физического параметра ${\displaystyle C(z)}$, конечно, продуцирует гравитационное поле, но без аномалий.

Моделирование и расчёты для ${\displaystyle C(z)}$ представляют из себя проблему, т.к. обычно под полем понимают именно аномальное поле. Поле от ${\displaystyle C(z)}$ - также постоянное слагаемое, которое имеет константный, не зависящий от координаты на плоскости, гравитационный эффект, для которого справедливы следующие утверждения:

1. Его значение эквивалентно коэффициенту ДПФ для нулевой частоты. 2. Он является удаляемым слагаемым при процедуре удалении тренда, которая является стандартной для обработки аномальных полей.

Таким образом, участие ${\displaystyle C(z)}$ в обратной задаче имеет только излишний, шумящий характер, и поэтому при решении обратных задач приравнивается нулю.

См. также

Магнитные параметры Земли IGRF

Ссылки на литературу

[1] С.В. Мицын, Г.А. Ососков "Экстраполяция сеточных моделей геофизичеких полей методом конечных разностей". Геоинформатика, №3, 2016, сс. 29-34