Задача деконволюции в инерциальной гравиметрии Текст научной статьи по специальности «Физика»

ФИЗИКА ЗЕМЛИ, АТМОСФЕРЫ И ГИДРОСФЕРЫ

Задача деконволюции в инерциальной гравиметрии

В. Л. Пантелеев, Т. С. Чесноковаа

Государственный астрономический институт имени П. К. Штернберга Московского государственного университета имени М.В. Ломоносова. Россия, 119991, Москва, Университетский просп., д. 13.

E-mail: а bulg33@mail.ru

Статья поступила 13.07.2010, подписана в печать 04.10.2010

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

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

УДК: 550.831(26); 528.563(26); 550.831.23. PACS: 93.85.Вс.

Введение

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

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

1. Характеристики упругой системы гравиметра

Поскольку в принятой системе координат вертикальная компонента ускорения ш2 и вектор силы тяжести § по направлению совпадают, то эти две силы действуют на чувствительный элемент совместно (аддитивно).

Запишем упрощенное линеаризованное дифференциальное уравнение упругой системы гравиметра. Оно

имеет вид уравнения колебаний горизонтального рычага (маятника) гравиметра под действием силы тяжести и вертикальных ускорений:

■п28 = п2(хюг + д-д0), (1)

• 2ES ■

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

хЦ) + 2ехЦ) + п2хЦ) = п2уЦ) представляет собой свертку функций /г(£)и г/(£):

x(t) =

h(r)y(t - т)с1т.

(2)

Здесь /г(т) — весовая функция или импульсная характеристика динамической системы [2, 3].

Чтобы получить частотные характеристики упругой системы, подчиняющейся уравнению (1), необходимо в качестве входной тестовой функции взять единичную гармонику y(t) = exp(iojt). Тогда выходной функцией (установившимся решением уравнения) будет комплексная функция x(t) = W(ujj) exp(tatf), где W(iw) носит название передаточной функции, или комплексной частотной характеристики системы. Передаточная функция, так же как и h(t), содержит полную информацию о свойствах данной линейной динамической системы.

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

Для того чтобы упругая система, подчиняющаяся уравнению (1), могла служить гравиметром, необходимо, чтобы система давала отклик на ничтожно малые изменения «входного сигнала» — аддитивной смеси приращения силы тяжести и вертикального ускорения.

Поскольку «амплитуда» ускорения вертикальной качки нередко достигает 1 м/с^2, а измеряемое приращение силы тяжести, доступное для измерителя, должно быть 1 мкм-с-2 (отношение сигнал/шум равно всего 1СГ6), то этот прибор должен иметь динамический диапазон не менее 106, чего достигнуть технически очень сложно.

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

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

2es{t) + n2s{t) = rr{Ag(t) + w?{t)).

(3)

Если разделить обе части уравнения на л2, то выражение (3) для упругой системы гравиметра примет вид

Ts(l)+s(l) = Ag(l) + w?(l), (4)

где Т = 2г/п2. Эта величина, имеющая размерность времени, получила название постоянной времени.

2. Передаточная функция гравиметра

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

Если дифференциальное уравнение упругой системы с достаточной точностью может быть изображено линейным уравнением первого порядка (4), то установившимся решением его будет

e^T/r(Ag(t - т) + wz(t - т)) (1т,

(5)

ветственно равны

А(Ш) =

1

Vi +W

Ф(о>) = - arctg Тш.

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

5(0 = Гр){А +

Рис. /. Структурная схема гравиметра

Согласно этой структурной схеме, передаточная функция гравиметра будет иметь вид

Г(£>) =

О + цЩё(О)Щ(О)Ж<'с0П(О)'

(6)

где ц — коэффициент гармонической (или стохастической) линеаризации, — передаточная функция упругой системы гравиметра, ЖЦО) — передаточная функция дополнительного фильтра нижних частот, Гсогг(1?) — передаточная функция корректирующего звена.

Структурная схема гравиметра может содержать и существенно нелинейные элементы, ограничивающие сигнал, как, например, это выполнено в морском гравиметре ГМН института ВНИИГеофизики, но и в этом случае в рабочем режиме дифференциальное уравнение гравиметра можно считать приближенно линейным, вернее линеаризованным. Приемы технической линеаризации известны и подробно описаны в литературе [5].

3. Деконволюдия и фильтрация полезного сигнала

Запишем дифференциальное уравнение гравиметра (линейное приближение) в общем виде:

ЛГ(Я){5(0} = М(Я){Д£(0 + ££'Д0}>

где

ЩБ) = аоО" + о,!)"-1 + ... + а„,

М(П) = Ь0От + + ... + Ьт.

Здесь О — оператор дифференцирования, а От = = (1т/сИт. Передаточной функцией нашего воображаемого гравиметра будет

М(О)

T{D) =

(7)

а передаточной функцией можно считать комплексную функцию

При этом амплитудно-частотная характеристика (АЧХ) и фазово-частотная характеристика (ФЧХ) будут соот-

ЩО)-

Решение будет устойчивым, если т<п, и неустойчивым, если т>п. Установившимся решением будет свертка или интеграл Дюамеля (2). Другое название этой операции — конволюция. С помощью операции конволюции определяется сигнал динамической системы, если известен входной сигнал. Обратная задача — деконволюция (определение входного сигнала

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

s(t) = T(D){Ag(t) + wz(t)}. (8)

Формально «решив» полученное уравнение относительно Ag, получим

Таким образом, восстановив входной сигнал, мы восстанавливаем одновременно и возмущающее ускорение.

Теперь подействуем на «восстановленный» сигнал оператором фильтра нижних частот Фо с действительной частотной характеристикой:

Agit) = Фо{г-1{«(0}-О'г(0} = Ф0Г^{^)ЬФ0{ш2(Щ.

В диапазоне частот измеряемого сигнала амплитудно-частотная характеристика фильтра близка к единице, а в диапазоне частот помех — близка к нулю. Конечно, такого идеального фильтра нижних частот (ФНЧ) не существует. Измеренное приращение силы тяжести будет отличаться от истинного. Поэтому вместо приращения Ag мы получим лишь оценку этой величины Ag, но всегда можно выбрать такой фильтр, который практически не будет пропускать шумы, обусловленные качкой корабля. То есть можно считать, что Ф0{г}иО. В общем случае вертикальная компонента возмущающего ускорения wz(t) не обязательно связана с вертикальной качкой. Она может быть обусловлена, например, кориолисовым ускорением при движении корабля или непостоянством высоты полета самолета и т. п. То есть реально существуют такие возмущающие ускорения wz, которые могут лежать и в диапазоне частот полезного (измеряемого) сигнала. Их необходимо учитывать отдельно, не прибегая к методам фильтрации. Итак,

Ag( 0 = #0r^{s(0}. (9)

Остается выбрать Ф0 — передаточную функцию фильтра нижних частот.

Мы выполнили численные оценки погрешностей фильтрации для фильтров различных типов, не обязательно физически реализуемых [5]. Оказалось, что точность фильтрации, когда диапазоны частот измеряемого сигнала и возмущающего ускорения не совпадают, прежде всего зависит от выбора полосы пропускания, в частности для фильтра нижних частот — частоты среза. Необходимо лишь, чтобы крутизна АЧХ в окрестности частоты среза была достаточно большой. Величину, обратную частоте среза, часто называют постоянной времени фильтра.

Обязательное условие для фильтра нижних частот — частотная характеристика алгоритма фильтрации должна быть действительной функцией частоты,

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

Ф0 (D) = W(D)W(-D).

(10)

Выходной сигнал этой системы также определяется интегралом Дюамеля, т. е. конволюцией — сверткой двух функций, одна из которых — входной сигнал, а вторая — четная весовая функция, а интегрирование выполняется от ^оо до +оо.

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

Приведем несколько характеристик ФНЧ четвертого порядка, пригодных для инерциальной гравиметрии:

1) Ф0М =

UJ" + UJf,

h0(t) = ^j- [exp(-u>oH)(acos/?u>oO + ß sin ßojo\t\

^j- [exp(-u>o/?K|)(asin aojo\t\ + ß cos aujot)],

• 71 а 71

« = sm-, ß = cos —;

2) Ф0М =

Wn

OJ,

MO = -¡r exp(-woKI) x

3 ÜJnt ( 3 . Л . ÜJQ\t\

3) ФоМ =

(wg+^xwä+w2)

h(t) = Y ехР(^шо|Ф (l + t

з7! exp V--2 J sm I'

4) Ф0(о;) =

OJ,

0

U + üjzy

= М ехР(^шоН)(15 + 15шо|*| + +и#3|) •

В приведенных формулах шо — частота среза ФНЧ, т. е. такой фильтр «срезает» все частоты, величина которых больше шо.

Остановимся на втором варианте частотной характеристики ФНЧ. Заменяя величину ш2 квадратом опера-

тора дифференцирования с обратным знаком, получим

Ф0 (D) =

Шп

U + or-'

(П)

Связь входного сигнала у{1) с выходным сигналом х{1) определится сверткой входного сигнала с весовой функцией /го(?):

.г(О = Фо(Д)МО} = Ло(О*0(О. В данном случае

Ао(0 = ^ехр

V2

3 wqí —= cos ——

\/2 у/2

OJ0\t\

3 \ . Lüoltl — sin —Lr1

sfi. J s/2

(12)

Графики этой весовой функции для разных величин частоты среза (1/10, 1/30, 1/50) показаны на рис. 2.

МО о.оз

0.02

0.01

1/10

/ / /, *""" 1/50 / У / \..........1/30 \\

\J 1/ -

-0.01

О 200 400 600 800 1000 г

Рис. 2. Весовые функции фильтра нижних частот

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

В общем случае для гравиметра с линейной системой оценка измеряемой величины будет

Для корректности операции необходимо, чтобы степень числителя в последнем выражении была меньше степени знаменателя.

В частном случае, когда Г '(£)) = 1 + TD, тогда

<^(1+77))

Фо(£>)Г (D) =

2 '

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

Весовая функция (импульсная характеристика), корректирующая динамические погрешности гравиметра, обладающего передаточной функцией первого порядка, будет иметь вид h(t) = Th0(t) + h0(t), а в случае применения ФНЧ с /го(0. определяемой формулой (12), весовая функция будет иметь вид

h{t) =

8у/2

(з-

ехр

шоИ

' у/2

9ч OJQÍ

ТЩ) cos -J=

(Г+/2г) cosing

у/2

(3 - Tlw20) sin

2\ „.•„ шоИ

y/2

На рис. 3 показан вид весовой функции, корректирующий динамические погрешности и не «пропускающий» инерциальную помеху. Видно, что она не обладает свойством симметрии. В данном случае эффективная постоянная времени гравиметра равна 30, 40, 50 и 100 с. Частота среза равна 1/30 с^1.

Для иллюстрации эффективности и точности работы предложенных алгоритмов на реальном практическом материале мы использовали в качестве «полез-

й(0

0.02

0.01

-0.01

/ч,Г=100

1 Y-T= 50

Т= 40 С/ \ *,

I ^ 1 ю0 = 1/30

/ щ

/ \«\ /'1 щ

!' щ Г=30 /,! Щ

\ п\ \

у /, 1

0 200 400 600 800 1000 t

Рис. 3. Весовые функции корректирующего фильтра

млГал

Рис. 4. Иллюстрация восстановления сигнала после гравиметра: / — «полезный сигнал», 2 — показания гравиметра, 3 — восстановленный сигнал

ного сигнала» гравиметрические профили, полученные в морских экспедициях сотрудниками кафедры геофизики геологического факультета МГУ (руководитель A.A. Булычев). Это большой массив значений Ag (800000 значений), полученный с непрерывного профиля через 0.1 с. Далее этот сигнал был отягощен вертикальным шумовым ускорением около 0.67 м/с^2. Согласно формуле (9), получали оценку Ag. Точность восстановления полезного сигнала оказалась равной

Problem deconvolution in inertial gravimetry

1.2-1.6 мкм/с2. При увеличении шага съема информации до одной секунды, погрешность восстановления сигнала возросла до 4 мкм/с2.

Рис. 4 иллюстрирует результат восстановления (изображен небольшой отрезок массива — 80000 значений).

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

Заключение

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

Список литературы

1. Пантелеев В.Л. Основы морской гравиметрии. М., 1983.

2. Ройтенберг Я.Н. Автоматическое управление. М., 1971.

3. Левин Б.Р. Теория случайных процессов и ее применение в радиотехнике. М., 1960.

4. Голован A.A., Горицкий А.10., Парусников H.A., Тихомиров В.В. Алгоритмы корректируемых инерциальных навигационных систем: Препринт. М., 1994.

5. Цыпкин Я.З. Основы теории автоматических систем. М., 1977.

V. L. Panteleev, Т. S. Chesnokova"

Р. К. Sternberg State Institute of Astronomy, M. V. Lomonosov Moscow State University, Universitetskii pr. 13, Moscow И9991, Russia. E-mail: " bulg33@mail.ru.

The method of allocation of slowly changing measured signal — gravity variations — from an additive mix of a useful signal and a hindrance much more surpassing measured variations of a specific gravity is discussed. The basic problem — an exception of the errors connected with dynamic characteristics gravirneter. The variant as which it is possible to consider alternative to known filter Kalmana is offered.

Keywords: specific gravity, vertical component of acceleration, parcel, weight function, transfer function, frequency characteristic, filter of bottom frequencies, convolution, deconvolution, correcting filter, signal restoration. PACS: 93.85.Bc. Received 13 July 2010.

English version: Moscow University Physics Bulletin 1(2011). Сведения об авторах

1. Пантелеев Валерий Леонтьенич — докт. физ.-мат наук, профессор, профессор; тел.: (495) 939-37-64, e-mail: bulg33(u>mail.rii.

2. Чеснокона Татьяна Сергеенна — нед. инженер; тел.: (495) 939-37-64, e-mail: bulg33(u>mail.rii.