Вариации поля силы тяжести, регистрируемые на поверхности земной коры, являются продуктом процессов естественного и техногенного происхождения, которые имеют различный характер (резкие мгновенные, периодические или квазипериодические, вековые). В зависимости от пространственной протяженности они могут быть локальными, региональными или глобальными. Изучение и анализ таких временных вариаций поля силы тяжести обеспечивает наблюдателя данными о кинематике и динамике Земли для решения ряда фундаментальных и прикладных задач геофизики. Исследование и разработка таких моделей в этой области еще только начинаются [
В данной работе предлагается новый подход к построению математических моделей вариаций гравитационного поля на основе использования модифицированного Фурье-преобразования инструментальных данных группы гравиметров при регистрации естественных колебаний земной поверхности.
Для понимания механизма влияния внешних процессов на вариации гравитационного поля необходимо найти математическую функцию, значения которой максимально приближены к экспериментальным гравиметрическим данным. Это достигается с помощью аппроксимации данных эмпирическими функциями. Аппроксимация позволяет изучать числовые характеристики и свойства объекта, сводя задачу к анализу более простых математических моделей. Процесс аппроксимации заключается в построении функции, которая наилучшим образом соответствует исходным данным с минимальной погрешностью.
Целью аппроксимации является построение математической модели, максимально приближенной к «идеальной» модели, которая воспроизводит измеренные данные с минимальными ошибками и высокой достоверностью.
При выборе подходящей аппроксимирующей функции необходимо учитывать параметры, наиболее точно описывающие зависимость между изучаемыми величинами. Из нескольких подходящих вариантов выбирается модель на основе баланса между точностью аппроксимации и простотой вычислительной конструкции модели.
Теория предлагает модели, гипотезы, основанные на некоторых проверяемых предположениях и постулатах. Эксперимент, в свою очередь, проверяет эти гипотезы, стараясь максимально исключить влияние неконтролируемых факторов и искажений (помех), вносимых измерительной аппаратурой (обычно эти факторы определяются как аппаратная функция или функция прибора) и учета случайных «гостевых» процессов — совокупности внешних и локальных факторов, не выясненной и до конца не изученной природы, но влияющих на изучаемые вариации гравитационного поля. При этом необходим проверенный или верифицируемый принцип, которому подчиняются практически все измерения, из математической формулировки которого можно будет вывести некую «универсальную» подгоночную функцию, которая позволила бы описать все измерения. Таким образом, на основе этого принципа можно будет получить общую математическую модель, которой должны будут подчиняться все измерения, удовлетворяющие этому принципу.
Постановка задачи
Фактически каждый исследователь, применяющий Ф-преобразование, осведомлен об основном его недостатке. Этот недостаток кроется в предположении, что исходный случайный сигнал является строго периодическим, т.e:
Sg(t + T) = T Sg(t). (1)
Здесь Sg(t) определяет исходный случайный сигнал, T — период, совпадающий с Range(t), где размах Range(t) = max(t) — min(t) определяет приближенно длину исходного сигнала Sg(t).
Во многих случаях это предположение недоказуемо, и поэтому многие исследователи пытаются преодолеть этот недостаток различными методами и способами. Если мы захотим заменить чистый периодический сигнал на некую апериодическую копию, то мы столкнемся со следующей проблемой. Проблема подгонки дискретных апериодических сигналов не решается с помощью интегрального преобразования Фурье апериодических сигналов, поэтому его нельзя использовать для предсказания дискретного апериодического сигнала за пределами заданного интервала временного окна. Дискретные представления многих аналоговых сигналов играют важную роль в их обработке. Они содержат необходимую информацию, связанную со свойствами сигналов, и допускают их дальнейшую обработку [
Разложение в классический ряд Фурье (РФ) является наиболее простым и часто используемым инструментом в области обработки сигналов. Однако он не позволяет выделить как субгармонические, так и интергармонические компоненты заданного сигнала и во многих ситуациях имеет серьезные недостатки [
Предлагаемый метод позволяет преодолеть ограничения Фурье-анализа. Фурье-анализ основан на временно-частотных методах [
Существует ли простое и универсальное преобразование исходного сигнала, позволяющее преобразовать его в другой цифровой сигнал, имеющий строго период равный 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).
Описание гравиметрических данных
Для реализации этой постановки задачи были проведены эксперименты по регистрации вариаций силы тяжести одновременно несколькими гравиметрами на базе фундаментального гравиметрического пункта «Ледово» [
В комплексе геофизической аппаратуры используются гравиметрические приборы различных типов, действующие на одних и тех же физических принципах и обладающие общим частотным диапазоном. Расположение всех гравиметров на одном геофизическом постаменте позволяет считать зарегистрированные инерционные помехи идентичными по происхождению. Давление, температура, дрейф нуля при данных измерениях не оказывают значительного влияния на результаты измерений.
Гравиметры нового поколения SCINTREX CG-6 оснащены блоками компьютерной обработки, которые позволяют вводить вычисленные поправки в показания прибора в режиме реального времени, существенно снизив влияние помех и обеспечив точность измерений порядка 0,001 мГал. Автоматическое снятие показаний и непрерывная дискретизация данных минимизируют ошибки оператора
Гравиметр ГНУ-КВ был адаптирован для регистрации вариаций гравитационного поля с использованием видеокамеры в качестве ключевого элемента процесса цифровизации показаний. Это инновационное решение позволяет с высокой точностью фиксировать и детально анализировать вариации гравитационного поля, в том числе амплитудно-частотные характеристики (АЧХ) [
В данной работе анализ АЧХ ограничивается только для пары приборов ГНУ-КВ и SCINTREX CG-6, так как SCINTREX CG-6 является базовым прибором при регистрации вариаций гравитационного поля на фундаментальном гравиметрическом пункте «Ледово».
Алгоритм обработки данных
Так как данные от двух гравиметров представляют собой прямоугольные матрицы с числом строк N = 1000 и столбцов М = 55, то (для иллюстрации предлагаемого метода и сохранения приемлемого объема страниц данной статьи) можно выделить только три характерных сигнала с максимальным, средним и минимальными значениями соответственно. Значения этих трех столбцов исходных матриц находятся из так называемой трубы доверия (ТД), которая показывает распределения этих трех значений по всей совокупности столбцов исходных матриц. Для первичного анализа выберем матрицу данных гравиметра GNUKV. Труба доверия для этой матрицы показана на рисунке 1а, б.
На рисунке 1а, б показаны распределения экстремальных значений ТД по всем 55 столбцам. Максимальные и минимальные значения выделены соответственно красным и зеленым цветами. Распределение средних значений выделено синим цветом и ввиду малости их значений показано отдельно на правом рисунке. Максимальное и минимальные значения соответствуют столбцу матрицы под нулевым и 14-м номером и явно видны на левом рисунке.
На рисунке 2а, б показаны три самых характерных сигнала, выделенных по характерным экстремальным точкам ТД и подлежащих подгону. Напрямую такой подгон с высокой точностью практически невозможен, но расчет их аргументов Fp(t), допускающих строгое разложение в РФ, позволяет решить поставленную задачу. Для дальнейших целей удобно использовать в качестве независимой переменной нормированную на единицу переменную xj = j/N.
На рисунке 3 показаны все три функции аргументов, подлежащих подгону. На левом рисунке (a) красным и зеленым цветом соответственно показаны функции F0,2(t) для максимальных и минимальных значений, а на правом рисунке (б) показана функция F1(t) для распределения средних значений, показанная на рисунке 2 справа.
Эти три функции строго периодические (нетрудно заметить, что все они лежат в интервале [0, π]), и поэтому к ним применимо разложение (3). Для более точных расчетов мы применим к ним метод Неортогонального Комбинированного Фурье Анализа Сглаженных Сигналов (НОКФАСС) [
На рисунке 4 объяснена суть НОКФАСС. Слева (а) показана функция F0(t), а справа — (б) ее Ф-спектр. Резонансная частота располагается ровно посередине на N/2 = 500. Для подгоночных целей достаточно взять безразмерную полосу частот [500, 650].
Так как спектр периодической функции известен и приемлемая полоса частот для подгона определена, то дальнейшие действия осуществляются по стандартному МНК.
На рисунке 5 слева (а) показан укороченный спектр частот, применимый для подгоночных целей. Он расположен в полосе [2—400]. Справа (б) показан подгон функции F0(t). Точность подгона весьма высокая и величина относительной ошибки меньше чем 0,2%. Она рассчитывается по известной формуле
(4)
Более точные значения приведены в таблице 1.
На рисунке 6 показано распределение модуля амплитуд Amdk = √Ack2 + Ask2 слева (а) и распределение фаз, вычисляемое по формуле arctg(Ask/Ack) для функции F0(t) справа (б). Заметим, что распределение фаз практически равномерное, так как последовательность ранжированных амплитуд (ПРА), выделенная сплошной синей линией на правом рисунке, практически совпадает с отрезком прямой линии.
На рисунке 7а, б, в проиллюстрирован переход от функций аргумента тригонометрической функции F0(t) к исходной функции распределения максимумов, расположенных на правом рисунке вместе с ее подгоночной функцией. Для нахождения окончательных параметров a и b используются выражения (2) и (5), показанное ниже. Эти параметры приведены в таблице 1.
Ex(t) = a · cos(Yft(t)) + b. (5)
Так как новый подход был подробно объяснен на примере подгона функции F0(t), то авторы избавлены от необходимости дублировать метод, используемый для других функций F1,2(t). Мы приведем только окончательные результаты фитинга этих функций, необходимых для сравнительного анализа, сохраняя по возможности прежние настройки вычислительной программы.
На рисунке 8 слева (а) показан укороченный спектр частот, применимый для подгоночных целей. Он расположен в той же полосе [ 2—400]. Справа (б) показан подгон функции F1(t). Точность подгона весьма высокая и величина относительной ошибки меньше чем 0,2%.
На рисунке 9 показано распределение модуля амплитуд Amdk = √Ack2 + Ask2 слева (а) и распределение фаз, вычисляемое по формуле Phasek = arctg(Ask/Ack) для функции F1(t), справа (б). Заметим, что распределение фаз, как и в первом случае, практически равномерное, так как последовательность ранжированных амплитуд (ПРА), выделенная сплошной синей линией на правом рисунке, как и ранее, практически совпадает с отрезком прямой линии.
Для экономии места мы опускаем рисунки, объясняющие переход к исходному усредненному сигналу. Соответствующие параметры a и b, осуществляющие этот переход по формуле (2), приведены в таблице 1.
Аналогичным образом для полноты картины приводим 4 рисунка для аргумента F2(t).
На рисунке 10 слева (а) показан укороченный спектр частот, применимый для подгоночных целей [2—400]. Как и для прежних случаев, он расположен в той же полосе [ 2—400]. Справа (б) показан подгон функции F2(t). Точность подгона весьма высокая и величина относительной ошибки, как и ранее, меньше чем 0,2%. Точные значения приведены в таблице 1.
На рисунке 11 показано распределение модуля амплитуд Amdk = √Ack2 + Ask2 слева (а) и распределение фаз, вычисляемое по формуле Phasek = arctg(Ask/Ack) для функции F2(t), справа (б). Заметим, что распределение фаз, как и в первом случае, практически равномерное, так как последовательность ранжированных амплитуд (ПРА), выделенная сплошной синей линией на правом рисунке, как и ранее, практически совпадает с отрезком прямой линии.
Общий вывод, который следует из предварительного анализа рисунков 6—11, таков: усеченные спектры Ωk для всех трех функций Fp(t) одинаковы и расположены в интервале [ 2—400]. Основные отличия для всех трех Fp(t) лежат в распределениях модулей амплитуд (рис. 6а, 9а и 11а соответственно).
Совершенно аналогичным образом обрабатываются данные по гравиметру CG-6. Приведем (в целях экономии места) только основные рисунки, выделяющие различия между ними. По ТД определяются следующие колонки. Вторая для распределения максимумов, а минимальная совпадает с нулевой колонкой. Поэтому экстремальные файлы представлены следующими рисунками.
На рисунке 12а, б показаны три самых характерных сигнала для гравиметра CG-6, выделенных по критерию ТД и подлежащих подгону. Мы сохраняем те же обозначения для аргументов Fp(t), допускающих строгое разложение в РФ, что позволяет решить поставленную задачу. Как это было проделано ранее, удобно использовать в качестве независимой переменной нормированную на единицу переменную xj = j/N.
На рисунке 13 показаны все три функции аргументов для гравиметра CG6, подлежащих подгону. На левом рисунке (a) красным и зеленым цветом соответственно показаны функции F0,2(t) для максимальных и минимальных значений, а на правом рисунке (b) показана функция F1(t) для средней функции, изображенная на прежнем рисунке справа. Заметим, что все они автоматически приводятся к интервалу [0, π].
Так как усеченный спектр для второго гравиметра также сохраняется, то мы приведем только подгон функции F0(t) и распределение модуля амплитуд. Фазовые распределения также носят равномерный характер и поэтому в целях экономии места могут быть опущены.
Рисунок 14 показывает распределение модуля амплитуд Amdk = √Ack2 + Ask2 слева (а) и подгон функции F0(t) справа (б). Распределение фаз, как и для первого гравиметра, практически равномерное, и поэтому этот рисунок мы не приводим как малоинформативный.
На рисунке 15 показано распределение модуля амплитуд Amdk = √Ack2 + Ask2 слева (а) и подгон функции F1(t) справа (б). Распределение фаз, как и для первого гравиметра, практически равномерное, и поэтому этот рисунок мы также не приводим.
На рисунке 16 показано распределение модуля амплитуд Amdk = √Ack2 + Ask2 слева (а) и подгон функции F1(t) справа (б). Распределение фаз, как и для первого гравиметра, практически равномерное, и поэтому этот рисунок мы опускаем. Точные значения параметров подгонки для гравиметра CG-6 приведены в таблице 2.
Рис. 1. Распределение максимальных значений ТД. Слева на рисунке (а) показаны распределения максимальных, средних и минимальных значений по всей совокупности столбцов исходной матрицы измерений. Ввиду того что среднее значение принимает малые значения, то это распределение выделено на отдельном рисунке справа
Fig. 1. Distribution of maximum TD values. On the left, Figure (a) shows the distribution of maximum, average, and minimum values across the entire set of columns of the original measurement matrix. Due to the fact that the average value takes on small values, this distribution is highlighted in a separate figure on the right
Рис. 2. Наиболее характерные сигналы, выделенные по экстремумам ТД и подлежащие дальнейшему подгону
Fig. 2. The most characteristic signals identified by the extremes of TD and subject to further adjustment
Рис. 3. Функции аргументов Fp(t), выделенных теми же цветами, что и на прежнем рисунке, и подлежащие подгону
Fig. 3. Functions of arguments Fp(t), highlighted in the same colors as in the previous drawing and subject to adjustment
Рис. 4. Слева (а) показан аргумент F0(t), соответствующий максимальному распределению. Справа (б) — А-преобразование, сдвинутое на угол π. Для подгоночных целей достаточно взять набор частот от центральной резонансной частоты до 650
Fig. 4. On the left (a), the argument F0(t) is shown, corresponding to the maximum distribution. On the right (b) is A transformation shifted by an angle π. For fitting purposes, it is sufficient to take a set of frequencies from the central resonant frequency to 650
Рис. 5. Слева (а) показан укороченный спектр, расположенный в полосе [ 2—400]. Справа (б) показан подгон распределения F0(t). Точность подгона весьма высокая: относительная ошибка не превышает 0,2%
Fig. 5. The shortened spectrum located in the band [ 2—400] is shown on the left (a). On the right (b), the fitting of the distribution F0(t) is shown. The accuracy of the fitting is very high: the relative error does not exceed 0.2%
Таблица 1. Основные подгоночные параметры, которые были использованы для подгона аргументов функций Fp(t) — гравиметр GNUKV
Table 1. The main fitting parameters that were used to adjust the arguments of the Fp(t) functions. — GNUKV gravimeter
| Function | Fres | Ph0 | RelErr(%) | K | a | b |
| F0(t) | 1078,131 | 0,1041 | 0,183 | 400 | 0,38594 | 0,61406 |
| F1(t) | 950,547 | 0,1038 | 0,176 | 400 | 0,0021 | 0,55684 |
| F2(t) | 944,062 | 0,5810 | 0,114 | 400 | 0,4955 | 0,4955 |
Рис. 6. Слева (а) показано распределение модуля амплитуд, справа (б) показано распределение фаз. Голубая линия выделяет последовательность ранжированных амплитуд (ПРА). Ее близость к отрезку прямой показывает, что распределение фаз — равномерное
Fig. 6. On the left (a) shows the distribution of the amplitude modulus, on the right (b) shows the phase distribution. The blue line highlights the sequence of ranked amplitudes (PRA). Its proximity to a straight line segment shows that the phase distribution is uniform
Рис. 7. Переход от функции cos(F0(t)) (рис. a) к исходной функции распределения максимумов (в) через линейное преобразование (2). Оно показано на центральном рисунке (б)
Fig. 7. Transition from the cos(F0(t)) function (Fig. (a)) to the initial maximum distribution function (b) via linear transformation (2). It is shown in the central figure (b)
Рис. 8. Слева (а) показан укороченный спектр, расположенный в полосе [ 2—400]. Справа (б) показан подгон распределения F1(t). Точность подгона весьма высокая: относительная ошибка не превышает 0,2%
Fig. 8. The shortened spectrum located in the band [ 2—400] is shown on the left (a). On the right (b), the fit of the F1(t) distribution is shown. The accuracy of the fitting is very high: the relative error does not exceed 0.2%
Рис. 9. Слева (а) показано распределение модуля амплитуд, справа на рисунке (б) показано распределение фаз. Это АЧХ относится к функции F1(t). Голубая линия выделяет последовательность ранжированных амплитуд (ПРА). Ее близость к отрезку прямой показывает, что распределение фаз — равномерное
Fig. 9. On the left (a) shows the distribution of the amplitude modulus, on the right in Figure (b) shows the phase distribution. This frequency response refers to the F1(t) function. The blue line highlights the sequence of ranked amplitudes (PRA). Its proximity to a straight line segment shows that the phase distribution is uniform
Рис. 10. Слева (а) показан укороченный спектр, расположенный в той же полосе [ 2—400]. Справа (б) показан подгон распределения F2(t). Точность подгона весьма высокая: относительная ошибка не превышает 0,2%
Fig. 10. The shortened spectrum located in the same band [ 2—400] is shown on the left (a). On the right (b), the fitting of the F2(t) distribution is shown. The accuracy of the fitting is very high: the relative error does not exceed 0.2%
Рис. 11. На рисунке слева (а) показано распределение модуля амплитуд, справа на рисунке (б) показано распределение фаз. Это АЧХ относится к функции F2(t). Голубая линия выделяет последовательность ранжированных амплитуд (ПРА). Ее близость к отрезку прямой показывает, что распределение фаз — равномерное
Fig. 11. The figure on the left (a) shows the distribution of the amplitude modulus, and on the right (b) shows the phase distribution. This frequency response refers to the function F2(t). The blue line highlights the sequence of ranked amplitudes (PRA). Its proximity to the straight line segment shows that the phase distribution is uniform
Рис. 12. Распределение максимальных значений ТД. Слева на рисунке (а) показаны распределения максимальных, средних и минимальных значений по всей совокупности столбцов исходной матрицы измерений. Ввиду того что среднее значение принимает малые значения, это распределение выделено на отдельном рисунке справа
Fig. 12. Distribution of maximum TD values. Figure (a) on the left shows the distribution of maximum, average, and minimum values across the entire set of columns of the original measurement matrix. Due to the fact that the average value takes on small values, this distribution is highlighted in a separate figure on the right
Рис. 13. Функции аргументов Fp(t), выделенные теми же цветами, что и на прежнем рисунке, (3) и подлежащие подгону
Fig. 13. The functions of the arguments Fp(t), highlighted in the same colors as in the previous figure (3) and subject to adjustment
Рис. 14. Распределение модуля амплитуд слева (а) и подгон функции F0(t) справа (б)
Fig. 14. The distribution of the amplitude modulus on the left (a) and the fitting of the function F0(t) on the right (б)
Рис. 15. Распределение модуля амплитуд слева (а) и подгон функции F1(t) справа (б)
Fig. 15. The distribution of the amplitude modulus on the left (a) and the fit of the F1(t) function on the right (б)
Рис. 16. Распределение модуля амплитуд слева (а) и подгон функции F1(t) справа (б)
Fig. 16. The distribution of the amplitude modulus on the left (a) and the fit of the F1(t) function on the right (б)
Таблица 2. Основные подгоночные параметры, которые были использованы для подгона аргументов функций Fp(t) — гравиметр CG-6
Table 2. The main fitting parameters that were used to adjust the arguments of the Fp(t) functions. — gravimeter CG-6
| Function | Fres | Ph0 | RelErr(%) | K | a | b |
| F0(t) | 1075,869 | 0,035 | 0,315 | 400 | 0,1266 | 0,70176 |
| F1(t) | 1006,841 | 0,95054 | 0,364 | 400 | 0,0112 | 0,6938 |
| F2(t) | 762,871 | 0,9103 | 0,241 | 400 | 0,4840 | 0,5161 |