Применение модификационного преобразования Фурье для моделирования вариаций поля силы тяжести

Введение

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

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

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

Целью аппроксимации является построение математической модели, максимально приближенной к «идеальной» модели, которая воспроизводит измеренные данные с минимальными ошибками и высокой достоверностью.

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

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

Постановка задачи

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

Sg(t + T) = T Sg(t). (1)

Здесь Sg(t) определяет исходный случайный сигнал, T — период, совпадающий с Range(t), где размах Range(t) = max(t) — min(t) определяет приближенно длину исходного сигнала Sg(t).

Во многих случаях это предположение недоказуемо, и поэтому многие исследователи пытаются преодолеть этот недостаток различными методами и способами. Если мы захотим заменить чистый периодический сигнал на некую апериодическую копию, то мы столкнемся со следующей проблемой. Проблема подгонки дискретных апериодических сигналов не решается с помощью интегрального преобразования Фурье апериодических сигналов, поэтому его нельзя использовать для предсказания дискретного апериодического сигнала за пределами заданного интервала временного окна. Дискретные представления многих аналоговых сигналов играют важную роль в их обработке. Они содержат необходимую информацию, связанную со свойствами сигналов, и допускают их дальнейшую обработку [18]. В традиционной схеме сигналы могут быть представлены в виде рядов Тейлора-Макларена, Дирихле, Лорана, Лежандра, Паде, Прони и Фурье. Хотелось бы подчеркнуть, что эти разложения в ряды используются без какого-либо математического обоснования.

Разложение в классический ряд Фурье (РФ) является наиболее простым и часто используемым инструментом в области обработки сигналов. Однако он не позволяет выделить как субгармонические, так и интергармонические компоненты заданного сигнала и во многих ситуациях имеет серьезные недостатки [5][8][16][26].

Предлагаемый метод позволяет преодолеть ограничения Фурье-анализа. Фурье-анализ основан на временно-частотных методах [9], которые использовались в последние несколько десятилетий, а именно: дробный [4][24][25], кратковременный [15][21][23], оконный FT [11][13], Габора [22][29][31], вейвлет [14][17][28], Гильберта — Хуанга [10][20][30], преобразования Фурье — Бесселя [6][7][12] и даже разложение по эмпирическим модам [27]. Эти ссылки лишний раз показывают, что существует множество различных приближений для преодоления основного недостатка (1). Однако внимательный анализ позволяет сформулировать следующую задачу:

Существует ли простое и универсальное преобразование исходного сигнала, позволяющее преобразовать его в другой цифровой сигнал, имеющий строго период равный 2π? Оказывается, этот вопрос может быть решен положительно, если записать следующее соотношение

Sg(t) = a · cos(F(t)) + b, (2)

и попытаться анализировать вместо исходной функции Sg(t) аргумент F(t) функции косинуса для получения желаемого результата. Функция F(t) представляет собой объединенный ЧФМ-сигнал и обеспечивает желаемый интервал [-1,1] для cos(F(t)), а искомый аргумент F(t) попадает в периодический интервал [0, π]. Поэтому конечный результат разложения любого сигнала сохраняет структуру вида (2), если к нему добавить разложение F(t) в форме подгоночной функции Yft(t, K) в виде конечного отрезка РФ.

(3)

Разложение (3) учитывает свойство аргумента F(t)=F(t ± π), определяющее полу-периодическую функцию. Таким образом, эти два простых выражения (2) и (3) решают задачу разложения любой случайной функции Sg(t) в РФ, базируясь на разложении аргумента F(t).

Описание гравиметрических данных

Для реализации этой постановки задачи были проведены эксперименты по регистрации вариаций силы тяжести одновременно несколькими гравиметрами на базе фундаментального гравиметрического пункта «Ледово» [2].

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

Гравиметры нового поколения SCINTREX CG-6 оснащены блоками компьютерной обработки, которые позволяют вводить вычисленные поправки в показания прибора в режиме реального времени, существенно снизив влияние помех и обеспечив точность измерений порядка 0,001 мГал. Автоматическое снятие показаний и непрерывная дискретизация данных минимизируют ошибки оператора

Гравиметр ГНУ-КВ был адаптирован для регистрации вариаций гравитационного поля с использованием видеокамеры в качестве ключевого элемента процесса цифровизации показаний. Это инновационное решение позволяет с высокой точностью фиксировать и детально анализировать вариации гравитационного поля, в том числе амплитудно-частотные характеристики (АЧХ) [1].

В данной работе анализ АЧХ ограничивается только для пары приборов ГНУ-КВ и SCINTREX CG-6, так как SCINTREX CG-6 является базовым прибором при регистрации вариаций гравитационного поля на фундаментальном гравиметрическом пункте «Ледово».

Алгоритм обработки данных

Так как данные от двух гравиметров представляют собой прямоугольные матрицы с числом строк N = 1000 и столбцов М = 55, то (для иллюстрации предлагаемого метода и сохранения приемлемого объема страниц данной статьи) можно выделить только три характерных сигнала с максимальным, средним и минимальными значениями соответственно. Значения этих трех столбцов исходных матриц находятся из так называемой трубы доверия (ТД), которая показывает распределения этих трех значений по всей совокупности столбцов исходных матриц. Для первичного анализа выберем матрицу данных гравиметра GNUKV. Труба доверия для этой матрицы показана на рисунке 1а, б.

На рисунке 1а, б показаны распределения экстремальных значений ТД по всем 55 столбцам. Максимальные и минимальные значения выделены соответственно красным и зеленым цветами. Распределение средних значений выделено синим цветом и ввиду малости их значений показано отдельно на правом рисунке. Максимальное и минимальные значения соответствуют столбцу матрицы под нулевым и 14-м номером и явно видны на левом рисунке.

На рисунке 2а, б показаны три самых характерных сигнала, выделенных по характерным экстремальным точкам ТД и подлежащих подгону. Напрямую такой подгон с высокой точностью практически невозможен, но расчет их аргументов F_p(t), допускающих строгое разложение в РФ, позволяет решить поставленную задачу. Для дальнейших целей удобно использовать в качестве независимой переменной нормированную на единицу переменную x_j = j/N.

На рисунке 3 показаны все три функции аргументов, подлежащих подгону. На левом рисунке (a) красным и зеленым цветом соответственно показаны функции F_0,2(t) для максимальных и минимальных значений, а на правом рисунке (б) показана функция F₁(t) для распределения средних значений, показанная на рисунке 2 справа.

Эти три функции строго периодические (нетрудно заметить, что все они лежат в интервале [0, π]), и поэтому к ним применимо разложение (3). Для более точных расчетов мы применим к ним метод Неортогонального Комбинированного Фурье Анализа Сглаженных Сигналов (НОКФАСС) [19], предложенного одним из авторов (РРН) этой статьи. Суть этого комбинированного метода состоит в том, что весь спектр частот сдвигается к центру на угол π и из всего избыточного спектра частот используется только их малая часть, расположенная справа от резонансной частоты.

На рисунке 4 объяснена суть НОКФАСС. Слева (а) показана функция F₀(t), а справа — (б) ее Ф-спектр. Резонансная частота располагается ровно посередине на N/2 = 500. Для подгоночных целей достаточно взять безразмерную полосу частот [500, 650].

Так как спектр периодической функции известен и приемлемая полоса частот для подгона определена, то дальнейшие действия осуществляются по стандартному МНК.

На рисунке 5 слева (а) показан укороченный спектр частот, применимый для подгоночных целей. Он расположен в полосе [2—400]. Справа (б) показан подгон функции F₀(t). Точность подгона весьма высокая и величина относительной ошибки меньше чем 0,2%. Она рассчитывается по известной формуле

 (4)

Более точные значения приведены в таблице 1.

На рисунке 6 показано распределение модуля амплитуд Amd_k = √Ac_k² + As_k² слева (а) и распределение фаз, вычисляемое по формуле arctg(As_k/Ac_k) для функции F₀(t) справа (б). Заметим, что распределение фаз практически равномерное, так как последовательность ранжированных амплитуд (ПРА), выделенная сплошной синей линией на правом рисунке, практически совпадает с отрезком прямой линии.

На рисунке 7а, б, в проиллюстрирован переход от функций аргумента тригонометрической функции F₀(t) к исходной функции распределения максимумов, расположенных на правом рисунке вместе с ее подгоночной функцией. Для нахождения окончательных параметров a и b используются выражения (2) и (5), показанное ниже. Эти параметры приведены в таблице 1.

Ex(t) = a · cos(Yft(t)) + b. (5)

Так как новый подход был подробно объяснен на примере подгона функции F₀(t), то авторы избавлены от необходимости дублировать метод, используемый для других функций F_1,2(t). Мы приведем только окончательные результаты фитинга этих функций, необходимых для сравнительного анализа, сохраняя по возможности прежние настройки вычислительной программы.

На рисунке 8 слева (а) показан укороченный спектр частот, применимый для подгоночных целей. Он расположен в той же полосе [ 2—400]. Справа (б) показан подгон функции F₁(t). Точность подгона весьма высокая и величина относительной ошибки меньше чем 0,2%.

На рисунке 9 показано распределение модуля амплитуд Amd_k = √Ac_k² + As_k² слева (а) и распределение фаз, вычисляемое по формуле Phase_k = arctg(As_k/Ac_k) для функции F₁(t), справа (б). Заметим, что распределение фаз, как и в первом случае, практически равномерное, так как последовательность ранжированных амплитуд (ПРА), выделенная сплошной синей линией на правом рисунке, как и ранее, практически совпадает с отрезком прямой линии.

Для экономии места мы опускаем рисунки, объясняющие переход к исходному усредненному сигналу. Соответствующие параметры a и b, осуществляющие этот переход по формуле (2), приведены в таблице 1.

Аналогичным образом для полноты картины приводим 4 рисунка для аргумента F₂(t).

На рисунке 10 слева (а) показан укороченный спектр частот, применимый для подгоночных целей [2—400]. Как и для прежних случаев, он расположен в той же полосе [ 2—400]. Справа (б) показан подгон функции F₂(t). Точность подгона весьма высокая и величина относительной ошибки, как и ранее, меньше чем 0,2%. Точные значения приведены в таблице 1.

На рисунке 11 показано распределение модуля амплитуд Amd_k = √Ac_k² + As_k² слева (а) и распределение фаз, вычисляемое по формуле Phase_k = arctg(As_k/Ac_k) для функции F₂(t), справа (б). Заметим, что распределение фаз, как и в первом случае, практически равномерное, так как последовательность ранжированных амплитуд (ПРА), выделенная сплошной синей линией на правом рисунке, как и ранее, практически совпадает с отрезком прямой линии.

Общий вывод, который следует из предварительного анализа рисунков 6—11, таков: усеченные спектры Ω_k для всех трех функций F_p(t) одинаковы и расположены в интервале [ 2—400]. Основные отличия для всех трех F_p(t) лежат в распределениях модулей амплитуд (рис. 6а, 9а и 11а соответственно).

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

На рисунке 12а, б показаны три самых характерных сигнала для гравиметра CG-6, выделенных по критерию ТД и подлежащих подгону. Мы сохраняем те же обозначения для аргументов F_p(t), допускающих строгое разложение в РФ, что позволяет решить поставленную задачу. Как это было проделано ранее, удобно использовать в качестве независимой переменной нормированную на единицу переменную x_j = j/N.

На рисунке 13 показаны все три функции аргументов для гравиметра CG6, подлежащих подгону. На левом рисунке (a) красным и зеленым цветом соответственно показаны функции F_0,2(t) для максимальных и минимальных значений, а на правом рисунке (b) показана функция F₁(t) для средней функции, изображенная на прежнем рисунке справа. Заметим, что все они автоматически приводятся к интервалу [0, π].

Так как усеченный спектр для второго гравиметра также сохраняется, то мы приведем только подгон функции F₀(t) и распределение модуля амплитуд. Фазовые распределения также носят равномерный характер и поэтому в целях экономии места могут быть опущены.

Рисунок 14 показывает распределение модуля амплитуд Amd_k = √Ac_k² + As_k² слева (а) и подгон функции F₀(t) справа (б). Распределение фаз, как и для первого гравиметра, практически равномерное, и поэтому этот рисунок мы не приводим как малоинформативный.

На рисунке 15 показано распределение модуля амплитуд Amd_k = √Ac_k² + As_k² слева (а) и подгон функции F₁(t) справа (б). Распределение фаз, как и для первого гравиметра, практически равномерное, и поэтому этот рисунок мы также не приводим.

На рисунке 16 показано распределение модуля амплитуд Amd_k = √Ac_k² + As_k² слева (а) и подгон функции F₁(t) справа (б). Распределение фаз, как и для первого гравиметра, практически равномерное, и поэтому этот рисунок мы опускаем. Точные значения параметров подгонки для гравиметра CG-6 приведены в таблице 2.

Результаты и их обсуждение

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

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