CC BY
6(Г
dx
dt
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ И
ПРОЦЕССЫ УПРАВЛЕНИЯ N 4, 2009 Электронный журнал, рег. N П2375 от 07.03.97 ISSN 1817-2172
Г
http://www. neva. ги/journal http://www. math. spbu. ru/diffjournal/ e-mail: [email protected]
Теория обыкновенных дифференциальных уравнений
Дифференциально-разностные уравнения
УДК 517.958+517.962.2
Некоторые постановки задач интерпретации временных последовательностей на основе линейного моделирования
Драница Ю.П.*, Драница А.Ю.**
*МГТУ г. Мурманск, **ЗАО "Ланит"г. Москва
Рассмотрен один из возможных подходов к интерпретации экспериментальных данных на основе теории линейных динамических систем. Анализируемые данные рассматриваются как выход некоторой линейной динамической системы, описываемой обыкновенным дифференциальным уравнением (ОДУ) 1-ого порядка. Предложен способ идентификации параметров ОДУ по текущим данным, без привлечения какой-либо априорной информации.
Рассматриваются возможные направления применения и развития метода. Проанализированы методы решения обратных и некорректно поставленных задач на примере принципов локации, установлены причины их низкой разрешающей способности. Выдвинуты предложения, позволяющие на 1-2 порядка увеличить разрешение этих методов.
One of possible ways of the interpretation of analyzed data based on the theory of linear dynamical systems is described in this report. Analyzed data are considered as the result of a linear dynamical system
determined by an ordinary differential equation (ODE) i-th level. The way of identification of ODE parameters by present data without using any initial information is given. Possible ways of using and development of this method are described here.
On the example of principals of location methods of solution of opposite and ill-posed tasks were analyzed, reasons of their low resolution were found. Methods for increasing the resolution on one-two order are proposed.
Ключевые слова: Обыкновенные линейные дифференциальные уравнения, обратные задачи, некорректно поставленные задачи, задачи локации, сейсморазведка, высокое разрешение.
Keywords: ordinary linear differential equation, opposite tasks, ill-posed tasks, tasks of location, seism research, high resolution.
Введение
Ценность результатов интерпретации обрабатываемых данных возрастает, если она опирается на физическую сущность изучаемого явления, которая, в свою очередь, формализуется адекватным математическим аппаратом. Взгляд авторов заключается в том, что разрабатываемые методы должны быть не только формально-математическими, но и иметь хорошее методологическое обоснование. В этом контексте нами предлагается формализовать на едином методическом уровне некоторые постановки и решения задач на основе линейного моделирования. Интерпретация обрабатываемых данных выходом линейной системы позволяет на принципиально новой основе организовать все этапы обработки информации, от ее предобработки до окончательной интерпретации.
Теория линейных систем, основанная на теории линейных ОДУ, в настоящее время представляет мощный и развитый математический аппарат. Наибольшее применение этот аппарат получил в технических приложениях: электротехнике, электронике и других, и, как правило, для описания достаточно простых технических устройств. Однако попытки формального применения идей, развитых в этих приложениях, для исследования более сложных природных объектов, экономических систем и т.д. натолкнулись на серьезные проблемы. Самыми существенными недостатками прямого применения этих методов к сложным объектам является неустойчивость и низкая разрешающая способность получаемых решений.
В работе проанализирована возникшая ситуация на примере задач локации, к которым относится и сейсморазведка. Установлено, что классические постановки этих задач являются обратными и некорректно поставленными с традиционным решением в рамках так называемых амплитудных методов, разрешающая способность которых принципиально ограничена длиной волны зондирующего импульса.
Динамическое моделирование позволяет ввести в линейную постановку нелинейный элемент - фазу сигнала. Нелинейный компонент, введенный в линейную модель, дает возможность отойти от амплитудного метода обработки информации. В результате значительно повышается разрешающая способность моделирования. Работа посвящена задаче перехода от статического к динамическому моделированию.
1. Общая постановка проблемы интерпретации данных линейной моделью
Рассмотрим стационарную во времени относительно статистик первого и второго порядков линейную систему, на вход которой воздействует сигнал /(¿) и описываемую линейным неоднородным ОДУ с постоянными коэффициентами 1-ого порядка
а{)х(Г)(Ь) + а1х{1~1)(Ь) + а2х( _2)(£) + ■■■ + а _ хх(Ь) + щ = / (£), (1)
где х(п)(£) - пронзводпая п-ого порядка выходного сигнала линейной системы; £ - параметр типа времени. Решение этого дифференциального уравнения для различных значений входной величины / (£) дает соответствующую выходную величину х(£).
Описание с помощью ОДУ имеет тот недостаток, что его постоянные коэффициенты трудно измерить, к тому же выходную характеристику проще вычислить методами операционного исчисления. Поэтому в технических приложениях для описания зависимости (1) применяют так называемую импульсную переходную характеристику (ИПХ) или весовую функцию к(Ъ) доступные непосредственному измерению. Для этого на вход системы подают специальные тестовые сигналы /(£), измеряют выходной сигнал системы и по этой информации строят к(Ъ) [19].
С точки зрения теории ОДУ построенная таким образом весовая функция к(Ъ) является частным решением неоднородного уравнения (1). В том случае, если на вход линейной системы подается короткий импульс, в идеальном случае дельта-функция Дирака, то весовая функция к(Ъ) в теории
ОДУ получила название функции Грина [17]. Будем считать, что в результате специальных экспериментов для линейной системы построена функция Грина Тогда в установившемся режиме сигнал выхода линейной системы можно выразить через так называемый интеграл свертки [20]
ь ь 1
к(Л)/(г - Л)(Л = ! /(Л)к(г - Л)(\, (2)
а а 1
где Л - переменная интегрирования типа времени. Уравнение (2) устанавливает связь между входным, выходным сигналами и ИПХ линейной системы. Так как мы имеем дело с экспериментальными данными, т.е. обладающими измерительными ошибками, то более реалистичным описанием линейной системы будет следующее выражение
ь ь1
к(Л)/(г - Л)(Л + /(Л)к(г - Л)(Л + п(г), (з)
а а1
где п(г) - ошибки измерения и моделирования. В зависимости от конкретной задачи, пределы интегрирования могут быть константами, переменными, или несобственными числами На уравнении (3) удобно продемонстрировать различные варианты постановок линейных задач и особенности их решения, в частности, постановки прямых и обратных задач, фильтрации данных и другие.
1.1. Некоторые классические постановки задач
1.1.1. Задача линейной оптимальной фильтрации данных
Отметим, что уравнение (3) в теории линейных систем часто рассматривается как некоторый фильтр с системной функцией к(Ъ) [20, 27]. Особое внимание заслуживают так называемые оптимальные фильтры, например, фильтр Колмогорова-Випера. Оптимальная фильтрация в широком смысле предназначена для решения обширного круга задач [20, 22]: фильтрации шумовой компоненты сигнала; фильтрации в заданной полосе частот; сглаживание функций; конструировании так называемых инверсных фильтров и других.
При построении оптимальных фильтров основная проблема заключается в необходимости априорного задания корреляционных свойств шумовой компоненты и полезного сигнала. Классическая теория оптимальной фильтрации
в настоящее время позволяет построить фильтры только для достаточного узкого класса шумов близких к белому или некоррелируемому [22]. Это связано с тем, что в классической постановке не существует методов идентификации шумов с другими законами распределения.
1.1.2. Формулировка прямой задачи
Для иллюстрации прямых и обратных постановок рассмотрим типичную задачу сейсмического зондирования Земной коры. Зондируемую толщу Земли обычно представляют как последовательность слоев геологических пород с неизвестными наклонами границ и толщинами, плотностями вещества и скоростями звука. На поверхности Земли возбуждается сейсмический импульс (взрыв), или зондирующий импульс (ЗИ), который распространяется внутрь земли, частично отражается от ее слоев, вновь достигает поверхности земли и регистрируется системой датчиков, расставленных на поверхности.
Допустим, что акустические (и любые другие) свойства земных пород известны в любой точке трехмерного пространства г и описываются функцией /(г). В точке г0 Земной поверхности возбуждается ЗИ. При этом известна форма ЗИ и ее эволюция во времени, т.е. задана функция ^(г0,£). Закономерность распространения ЗИ в пространстве тоже известна и описывается, например, волновым уравнением.
Для указанного примера прямая задача заключается в нахождении отклика среды у(г, £) на сейсмический импульс в любой момент времени для любой точке пространства г, если известно, что распространение ЗИ в пространстве подчиняется волновому уравнению.
Таким образом, классическая постановка прямой задачи заключается в задании дифференциального или интегрального уравнений, по которому развивается процесс в пространстве и во времени, а также начальных и граничных условий. В рассматриваемом случае роль начальных или граничных условий выполняют функции /(г) ^(г0,£). Очевидно, что к прямой задаче можно отнести и задачу фильтрации данных (при известной функции фильтра
В частности, из приведенного примера видно, что для решения прямой задачи требуется задание большого объема априорной информации. Отметим, что для многих практически важных постановок эта априорная информация фактически и является искомым решением. Вполне очевидно, что прямая задача, с точки зрения сейсмической разведки, особого интереса не представляет.
1.1.3. Обратные и некорректно поставленные задачи
Обратная задача заключается в построении свойств объекта по наблюденным данным и априорной информации, имеющейся у исследователя [24]. Для приведенного выше примера требуется, по зафиксированному на поверхности Земли отклику слоистой среды на ЗИ, определить структуру и параметры зондируемой среды: толщины и глубины залегания отдельных прослоев, их плотности, акустические скорости, вещественный состав и другую информацию, т.е. восстановить функцию начальных условий /(г). Первичной или основной задачей сейсморазведки является выявление акустически неоднородных слоев Земной толщи и глубин их залегания - так называемая задача отбивки границ.
Рассмотрим особенности постановки обратной задачи на примере решения уравнения (3). Как указывалось выше, в технических приложениях весовая функция к(Ъ) может быть измерена непосредственно. Следовательно, оценке подлежит только функция /т.е. одна из подынтегральных функций. В этом случае уравнение (3) называется интегральным.
При решении обратных задач важную роль играет их устойчивость. Задача устойчива, если малым флуктуациям левой части уравнения (3), соответствуют малые флуктуации решений. Если задача изначально является неустойчивой, то решать ее нет смысла, поскольку погрешности алгоритмов, накапливающиеся в ходе решения численными методами, неизбежно приведут к тому, что будет найдено неверное решение.
Как правило, любые экспериментальные данные задачи характеризуются наличием шумов, которые коренным образом меняет идеологию решения обратных задач [24]. Если сама задача является устойчивой, то существование шума может устойчивость нарушать. В результате, двум различным сигналам /1 и /2, искаженных шумом, будут соответствовать очень похожие измерения у1 и у2 (и наоборот). Поэтому встает вопрос, можно ли извлечь из измерений полезную информацию о сигнале, если наличие шума делает задачу неустойчивой? Такие постановки называются задачами, поставленными некорректно, и для их решения развиты специальные методы, основанные на привлечении дополнительной априорной информации о решении.
Следует также отметить, что класс некорректно поставленных задач шире класса обратных, т.к. неустойчивость могут иметь и прямые задачи. Поскольку экспериментальные данные по необходимости случайны, естественная постановка обратных задач достигается в рамках теории статистического оценивания неизвестных параметров. Именно стохастичность модели форми-
рования данных обусловливает фундаментальную трудность, возникающую при обращении информации - неустойчивость решения.
Неустойчивость означает, что в пределах естественных флуктуации шума с наблюдаемыми данными примерно в равной мере согласуется множество возможных оценок исходного объекта, включая и существенно отличные от него. Несколько утрируя, можно сказать, что основная проблема связана не с нахождением подходящего решения некорректной задачи, а, напротив, с их обилием. Необходимо найти способ выбора из этого множества решения, заведомо близкого к оригиналу.
1.2. Формальная интерпретация некорректно поставленных задач
Интегральное уравнение (3) удобно представить в следующем операторном виде
Л/(¿,а ) = у(£) + п(£) (4)
где а - вектор параметров модели, А - оператор отображения, п - вектор ошибок наблюдений и моделирования.
В этом случае постановка задачи решения уравнения (4) состоит в следующем [3]. Пусть задан непрерывный оператор, который взаимно однозначно отображает элементы /(£, а) метрического пространства Е1 в элементы у(£) метрического пространства Е2. Требуется найти решение операторного уравнения (4) в классе функций /(£, а), если функция /(£) неизвестна, но имеются измерения у1, у2,..., у/ функции у(£) в I точках £1? £2,..., £/.
Отметим, что понятие "измерения"в контексте данной проблемы трактуются расширительно и могут включать не только какие-либо измерения физических величин, но и, например, результаты вычислительного эксперимента, решения прямых задач, какие-либо эмпирические зависимости, графики, законы физики и любую другую априорную информацию, которую предполагается использовать при построении зависимостей. Оператор А уравнения (4) может представлять любое правило, по которому осуществляется отображение функциональных зависимостей. Мы же будем рассматривать отображения, осуществляемые посредством дифференциальных или интегральных уравнений.
Говорят, что решение операторного уравнения (4) устойчиво, если малая вариация правой части у(£) приводит к малому изменению решения. Это означает, что для всякого достаточно малого £ найдется такое 5(е), что как
только выполняется условие
Р2(у(г),у(г + е)) < в (е),
окажется выполненным и неравенство
Р1(/(г),/(х + е)) ^ е.
Здесь Р1? Р2 означают расстояния, определяемые в метриках Е1 и Е2 соответственно. Говорят также, что задача решения операторного уравнения (4) поставлена корректно по Адамару, если его решение существует, единственно и устойчиво. Если не выполняется хотя бы одно из перечисленных требований, то говорят, что задача решения операторного уравнения (4) поставлена некорректно.
Типичной некорректно поставленной задачей является, например, интегральное уравнение (3). Уравнение (3), в зависимости от пределов интегрирования, является интегральным уравнением Фредгольма, или Вольтера первого рода, отображающим множество непрерывных на отрезке [0,т] функций / (г) на множество непрерывных на этом же от резке функций у (г). Как показано в [3], решение уравнения (3) является некорректно поставленным в
у( г)
могут соответствовать сколь угодно большие изменения зависимости /(г). Природа этого факта заключается в том, что интеграл в правой части уравнения (3) является сглаживающим оператором. Это сглаживание приводит к тому, что большие по амплитуде возмущения функции /(г)7 но сосредото-
Е1
функции у (г) (тем меньшие, чем меньше этот участок). При решении обрат-
у( г)
предполагать значительные отклонения функции /(г).
Современная практика имеет два подхода к решению этой проблемы [3]. Первый метод заключается в разложении функций по заранее фиксированному конечному базису. В этом случае решение зависит от конечного числа параметров, и если число наблюдений за функцией превышает это число, то решение находится методом наименьших квадратов. Во втором случае рассматривается класс равномерно ограниченных функций. Тогда решение может быть найдено введением регуляризирующего функционала [25].
1.3. Особенности обратных постановок для сложных объектов и процессов
В теории кибернетики существует методический подход, основанный на так называемом принципе "черного ящика", предполагающим моделирование
по наблюдению входов и выходов. С точки зрения разработчика или пользо-
"
"
являются ОДУ или ИПХ, входом - функция /(£), выходом - наблюденное
"
вать следующим образом: по замеренным входам и выходам линейной системы требуется аппроксимировать ее ИПХ. Знание ИПХ позволяет однозначно связать входные и выходные сигналы системы.
Однако получить ИПХ для естественных природных, экономических или гуманитарных систем, а также в случае сейсмических измерений, подачей тестовых сигналов на изучаемые объекты, не представляется возможным. Для сейсморазведки, например, определение формы ЗИ затруднено сложностью законов распространения волнового импульса в Земной толще. Геологические породы являются своеобразным фильтром высоких частот с переменными параметрами, в результате форма ЗИ меняется в пространстве и во времени, т.е. ^(г,£). И эти изменения могут быть значительными.
"
хотя бы потому, что даже формальное определение входов и выходов таких объектов является нетривиальной задачей. В этих условиях мы вынуждены предполагать, что любые измеренные данные являются отображением некоторых выходов абстрактной линейной системы, на входы которой поступают некоторые внешние сигналы. Эти входные сигналы преобразуются линейной системой в соответствии с ее ИПХ (или ОДУ) и поступают на ее выходы.
Таким образом, из трех объектов, определяющих поведение сложных систем, два могут быть принципиально неизмеряемы. Поэтому для сложных
"
постановка задачи: по наблюденным данным восстановить как ИПХ системы, так и ее входные сигналы. Из этих рассуждений следует, что для интерпретации сложных систем и объектов требуется априорная оценка ИПХ, основанная только на наблюденных "выходных"данных. Метод такой оценки предлагается ниже.
2. Предлагаемые принципы линейного моделирования
2.1. Задача оценки авторегрессионной модели
Пусть имеются измерения временной последовательности некоторых фи-
зических величин, выполненные в п точках, с постоянной дискретностью Д^ Эти измерения образуют временной ряд
х= {хЛ,х2,. ■ ■ ,хп). (5)
Для упрощения дальнейших выкладок будем считать, что ряд (5) нормирован следующим образом
хг = {хг — х)/ах, (6)
где х, ах- соответственно среднее значение и среднее квадратическое отклонение ряда (5). Нормированные по формуле (6) данные являются безразмерными, имеют нулевое среднее значение и единичную дисперсию. В дальнейшем будем предполагать, что данные (5) уже преобразованы по формуле (6) поэтому значок нормировки указывать не будем. Подчеркнем, что измерения в точке могут представлять как скалярную величину, так и вектор. Не ограничивая общности рассматриваемого подхода, будем считать, что измерения в точке представляют скаляр. Переход к векторным данным будем оговаривать особо.
Эмпирическая оценка функции автокорреляции (ФАК) g, определенная на нормированных данных (5), в зависимости от временного сдвига г, вычисляется по формуле
1 п-|т I
д(т) =-—г У^ хк хк+т ,т = 0 ± 1, ± I, (7)
п — |т| 1 1 к=1
где I - максимальная задержка (или лаг) сигнала, используемая при оценке ФАК. Отметим, что ФАК на нулевом сдвиге представляет дисперсию ряда (5) и для заданной нормировки д(0) = ах = 1, а коэффициенты (7) называются коэффициентами ФАК. В статистике показано, что оценки (7) являются несмещенными и состоятельными. Альтернативой несмещенной оценке ФАК вида (7) является так называемая смещенная оценка [4, 24]
1 п—1т 1
г{т) = ХкХк+т, (8)
п
к=1
которая связана с предыдущей соотношением
г(т) = П-Ш д(т) (9)
п
и, следовательно, смещение оценки (8) будет равно
т
п
При п ^ то это смещение стремится к нулю и оценка (8) становится асимптотически несмещенной. Сравнительный анализ дисперсий оценок (7) и (8) показывает [4], что для типичных приложений дисперсия оценки (7) будет больше, чем для оценки (8), поэтому в дальнейшем будем пользоваться оценкой ФАК вида (8). Ввиду свойств симметрии, ФАК (8) является четной функцией г(т) = г ( т), поэтому для ее вычисления достаточно ограничиться
нахождением г(т), при т ^0. Определим па данных (5) авторегрессионную ()
/
хь = ^ ак х_ + = 1, 2,...,п, (10)
к=1
где аак - постоянные коэффициенты авторегрессии; I - лаг модели ; £ - дискретный параметр времени. Зависимость от собственной предыстории процесса и обуславливает приставку "авто" в названии модели. При авторегрессионном моделировании принято считать [20, 23], что процесс представляет белый шум с пулевым средним значением и дисперсией о2. Функция включает как приборные помехи в данных, так и ошибки моделирования, связанные с неадекватностью используемой модели данным.
Обозначим вектором а = (аа1, аа2,..., аа/)т коэффициенты модели (10), где (*)т - оператор транспонирования; г(т) = гт - коэффициенты ФАК. В этом случае уравнение (10) можно рассматривать как свертку последовательности данных с коэффициентами регрессии а. Домножим уравнение (10) на ххь и усредним, в результате получим
го = а1 г1 + а2г2 +-----+ а/г/ + . (И)
Если умножить (10) на хь_ при г = 1, 2,... , I и провести усреднение, то получим так называемую систему уравнений Юла-Уокера [24]
п = а,1Г{) + а2Г1 + а3г2 + ... + щ г— г2 = а1Г1 + а2го + а3Г1 + ... + аг г— \ Г3 = а1Г2 + а2Г1 + азг0 + ... + агг—
г'1 = ам-1 + а2 гг-2 + азГ1_ з + ... + аг го
(12)
Введем для удобства векторы г0 = (г0,г1,... ,г—1)т, г1 = (г1,г2,... ,гг)т и матрицу размера I х I
Щго) =
г0 г1
г1 г0
г—1
Г1-2
Г'1—1 Г'1—2 ... Го
(13)
Для принятых обозначений уравнение (11) и система (12) в матричном виде записываются соответственно следующим образом:
го = атг 1 + а2,
(14)
В.(го)а = п. (15)
Матрица (13) представляет так называемую матрицу Теплица [18], а вектор г1 - это сдвинутый на единичный шаг вектор г0. Матрица Теплица является симметричной, положительно определенной с одинаковыми элементами на субдиагоналях, параллельных главной.
Уравнения Юла-Уокера (12) представляют собой систему из I линейных алгебраических уравнений с I неизвестными и позволяют оценить вектор параметров авторегрессии а, например, по методу наименьших квадратов. После этого формула (11) даст оценку дисперсии а2 ошибок моделирования и данных.
Мы привели классическую схему оценок коэффициентов авторегрессии (10). Основным недостатком этой методики является большие вычислительные затраты, связанные с обращением матрицы Теплица И при решении системы уравнений (12). Однако специфическая структура матрицы Теплица позволяет сократить затраты на ее обращение методом Левинсона [2]. Другой способ минимизации вычислительных затрат при расчете коэффициентов
авторегрессии заключается в использовании так называемых быстрых алгоритмов [21].
Не вызывает принципиальных затруднений решение данной задачи и для векторных временных рядов [5, 18]. По сути дела уравнение (10) является одношаговым прогнозом текущих значений временного ряда (скалярного или векторного) по его нескольким предыдущим значениям. Авторегрессия (10) является частным случаем линейного прогноза временного ряда на основе статистических методов. Наш подход заключается в принципиально иной, не статистической, интерпретации полученного уравнения авторегрессии.
2.2. Переход к разностному уравнению
С другой стороны, соотношение (10) можно проинтерпретировать несколько иным образом [7, 12, 16]. Для этого перепишем уравнение авторегрессии в следующем виде
/
хь ак хг-к = V*, £ = 1, 2,..., п, (16)
к=1
где V* - некоторая временная функция. Введем следующие обозначения: 60=1 - коэффициент при первом члене левой части выражения (16), ^ = —а^. Тогда выражение (16) можно представить так
/
хь—к = V*, £ = 1, 2,...,п. (17)
к=0
Уравнение (17) теперь можно интерпретировать [17] как обыкновенное неоднородное разностное линейное уравнение с постоянными коэффициентами (ОРУ) 1-ого порядка. Таким образом, мы поставили в соответствие АР(/) разностное уравнение с постоянными коэффициентами порядка I. Но с точки зрения линейного моделирования это уравнение описывает динамическую линейную систему /-ого порядка с дискретными сигналами на выходе.
Сигнал V* в теории линейных систем принято называть внешним возмущением линейной системы. Отметим, что уравнение (17) формально является просто иной записью уравнения (10). Однако интерпретацию процесса типа белого шума еь мы будем выполнять в принципиально другой постановке по сравнению с авторегрессионной моделью. Этот переход методически обоснован в рамках предполагаемого линейного моделирования, т.к. в обоих случаях связываются текущее значение ряда с несколькими его предыдущими. При
этом, разумеется, модули соответствующих коэффициентов а моделей (10) и Ь (17) совпадают.
Установленная связь между авторегрессией и ОРУ и является основной идеей настоящей работы. Но путь к пониманию этого факта оказался совсем не тривиальным и потребовал много усилий и времени. Сложность проблемы перехода от авторегрессии к ОРУ, по нашему мнению, заключается, в том, что, во-первых, существенно отличаются классы задач решаемых этими методами и, во-вторых, теоретической основой, на которой методики построены.
Авторегрессионные методы являются частным случаем регрессионного моделирования. В свою очередь, регрессионный анализ является прикладной наукой, описывающий методы обработки наблюдений на основе математической статистики. Основная идеология, используемая в регрессионном анализе, заключается в минимизации отклонений между натурными данными и модельными. В свете статистического подхода, параметры регрессионной модели должны удовлетворять ряду требований: состоятельности, несмещенности и эффективности. Требование несмещенности модели (10) приводит, в частности, к тому, что среднее значение ошибки моделирования должно равняться нулю.
ОРУ является частным случаем разностных методов обработки информации, которые возникли в результате развития численных подходов к решению задач. Основной задачей разностного моделирования является замена непрерывных по своей природе методов решения, например, дифференциальных или интегральных уравнений, дискретными вычислительными схемами и выяснение условий, при которых эта замена будет обладать требуемыми качествами: точностью аппроксимаций, устойчивостью и сходимостью к решению, минимальностью вычислительных затрат и т.д. Методическим аппаратом решения этих вопросов является классический математический анализ.
На поверхностный взгляд эти методы никак не связаны между собой. Однако проведенные исследования показывают, что между авторегрессией и ОДУ существует тесная связь. Эта связь немедленно обнаруживается при попытке интерпретации абстрактных математических конструкций на содержательном уровне. В данном случае мы попытались представить обрабатываемые данные не абстрактным набором чисел, никак не учитывающим природы данных, а как временную последовательность, которая является выходным сигналом, некоторого гипотетического генератора.
Эта смена акцентов позволяет теперь рассматривать данные (5) как дис-
кретный выход некоторой линейной динамической системы. Формирование функции выхода осуществляется, во-первых, внешним шумом V*, воздействующим на вход линейной системы, и, во-вторых, параметрами линейной системы, определяемыми левой частью разностного уравнения (17). Такой подход дает возможность воспользоваться теорией разностных уравнений [17]. Рассмотрим однородную часть уравнения (17)
/
^ Ькхь-к = 0,£ = 1, 2,...,п. (18)
к=0
Будем искать решение уравнения (18) в виде пробной функции хь = Аь. Подставив эту функцию в уравнение (18) и проведя элементарные преобразования, получим
Аь—/(&0А/ + М/—11 + ... + к- 1А + Ь/) = 0, £ =1,2,..., п. (19)
Рассмотрим выражение (19) более подробно. При конечном аргументе (£ — /) первый множитель этого выражения никогда не равен нулю, следовательно, для выполнения равенства необходимо, чтобы сумма в скобках обращалась в ноль
Ь)А' + Ь1 А" + ••• + Ь/ = 0. (20)
В теории ОРУ выражение (20) называется характеристическим уравнением, которое является полиномом /-ого порядка. Как известно, любой полином /-ого порядка имеет ровно / корней, включая кратные. Каждый &-ый
Аьк
называемую фундаментальную систему решений (ФСР) однородного ОРУ. Общее решение уравнения (18) выражается через фундаментальную систему в следующем виде [17]
хь = С1(А1 )ь + С2(А2)ь + • • • + С (А)ь, (21)
где С - некоторые константы (г = 1, 2,..., /), определяемые начальными условиями. Решение (21) предполагает, что все корни А{ различны. Если некоторый корень, скажем А1? имеет кратность т, то соответствующий член в решении (21) будет иметь вид (С1 + £С2 + • • • + £т—1СТО)А1-
При действительных коэффициентах уравнения (20) корни могут быть или реальными, или комплексно сопряженными. Рассмотрим &-ый комплекс-
но сопряженный корень Хк = ак ± гЬк этого уравнения. Так как модуль рк и аргумент Шк этого корня вычисляются соответственно по формулам
Рк = \/а\ + Ь2к, Шк = &т^(±Ьк/ак), (22)
то его показательная форма будет иметь вид Хк = рк ех.р(±шк). В последнем выражении учитываем только главное значение функции arctg. В теории линейных систем параметр Шк называется собственной круговой частотой линейной системы. При единичной дискретности данных, собственные частоты занимают диапазон О = [—ж; ж] рад./с. Отметим, что значение частоты \шк\ = ж соответствует частоте Найквиста. Если уравнение (18) имеет действительные коэффициенты, то два члена, соответствующих простым комплексно-сопряженным корням Хк = Рк ех.р(±Шк), можно заменить выражением [17]
р{ (Ак сов(шк 1) + Вк вт(ик ^), (23)
где ААк, ВВк - неизвестные константы. Отметим, что отрицательные частоты не несут дополнительной информации по сравнению с положительными, поэтому обычно за эффективный частотный диапазон принимают интервал О = [0; ж] рад./с. Для реального корпя решение будет иметь вид (т.к. Ьк =0, то Шк = 0)
АкРк. (24)
На основании принципа суперпозиции линейных систем и замен (22-24), общее решение однородного уравнения (21) будет иметь следующий вид
з г
хг = ^2 А'Р + р\(Аг соб(ш$) + Вг ът^)), (25)
1=1 г=з'+1
где 2 ~ число реальных корней. Отметим, что решения (21, 25) устойчивы, если р{ < 1, для вс ех г = 1, 2,...,1. Здесь мы не рассматриваем задачу о выборе оптимального порядка авторегрессионного уравнения (16), которая является самостоятельной проблемой. По этому вопросу существуют многочисленные рекомендации, некоторые из них можно найти в [21].
Таким образом, преобразования (22-24) показывают, что ФСР однородного ОРУ являются система затухающих показательных функций и и затухающих по показательному закону гармоник различных частот. Эти частоты
никак не связаны с периодом наблюдений, как, например, в разложении Фурье, и являются естественным следствием линейных связей в данных.
Задача аппроксимации частотного состава дискретного процесса на основе авторегрессионной модели в литературе хорошо известна [21]. Прямое решение этой задачи затрудняется тем, что оцениваемые параметры входят в функции (25) нелинейно. Для решения этой задачи могут использоваться итеративные процедуры градиентного спуска или метод Ньютона. Однако все эти алгоритмы требуют больших вычислительных затрат.
Вычислительные затруднения привели к разработке субоптимальных процедур оценивания, основанных на решении линейных уравнений и получивших название метода наименьших квадратов Прони. Метод Прони сводит нелинейные методы оценивания к процедуре факторизации полиномов. Однако как традиционные алгоритмы оценок, так и метод Прони, имеют чисто формальный характер и не учитывают природу обрабатываемых данных.
Предложенная нами процедура оценки динамических параметров ряда (5) на основе линейного моделирования имеет четкий физический смысл, позволяющий на содержательном уровне решить эту задачу. Но, что является более важным, нам удалось связать чисто статистическую задачу с теорией дифференциальных уравнений. Установленная адекватность между статистическим подходом и классическим математическим анализом позволяет на принципиально новом уровне осуществлять постановки и решения большого класса практически важных задач.
2.3. Переход к непрерывному времени
Будем рассматривать неоднородное ОРУ (17) как дискретный аналог неоднородного линейного ОДУ (1), описывающего некоторую линейную си/
этого уравнения
Ь0х(|)(£) + Ь1х( 1—11)(£) + Ь2х^—2)(£) + • • • + Ь ^х(£) + Ь = 0. (26)
х(£) =
ехр(&£). Вычислив необходимые производные пробной функции и подставив их в соотношение (26), после элементарных преобразований, по аналогии с выражением (20) получим следующий характеристический полином
Ь0^ + Ь^ + • • • + Ь/—+ Ь/ = 0. (27)
Согласно общей теории [17], если полином (27) имеет корни к17 к2^.к/, то общее решение однородного уравнения (26) будет иметь вид
хг = С1 ехр^) + С2 ехр(Ы) + • • • + Сехр(к/г). (28)
Решение (28) предполагает, что все корпи кг различны. Если некоторый к1
(28) будет иметь вид (С1 + Ю2 + • • • + гт-1Ст) ехр(к1г). Если коэффициенты уравнения (26) действительны, то комплексные корни встречаются сопряженными парами (ак ± гшк)- Соответствующая пара решений также будут комплексно сопряжена и может быть заменена действительными членами [17]
ех.р(а^)(Ак сов(шкг) + Вк вт^и^)). (29)
Для реального корня решение будет иметь вид (т.к. Шк=0)
Ак ехр(акг), (30)
где Ак, Вк - неизвестные константы. После внесения замен (29, 30) в уравнение (28), общее решение однородного ОДУ (26) будет иметь вид
з/
хг = ^2 Аг ехр(агг) + ^ ехр(агг)(Аг сов(игг) + Вг зт(ш$)), (31)
1=1 г=з^1
где 2 - число реальных корней. Неизвестные константы Ак, В к должны быть
аг
ОДУ получили название постоянных затухания. Отметим, что система затухающих экспонент и затухающих по экспоненте гармоник также образуют ФСР однородного ОДУ.
При расчете частот и постоянных затухания неявно предполагалось, что временная дискретность данных (5) равна единице Д1=1. Переход к реальным значениям параметров (арк7 ирк)7 при другой дискретизации данных, осуществляется по следующим формулам
арк = ак/Дгс-1, Шрк = Шк/Дгс-1.
Очевидно, что для существования устойчивого решения (31) требуется выполнение условия аг ^ 0, для всех г = 1, 2, ...,1. Таким образом, переход от
авторегрессии к ОРУ позволил на физически обоснованном уровне получить частотную и динамическую информацию об анализируемых данных.
2.4. Связь между дискретной и непрерывной моделями
Для установление связи между дискретной и непрерывной моделями временной последовательности рассмотрим &-ый модуль характеристического уравнения (20), описываемый выражением (22). Учитывая очевидное тождество /(х) = е/п^получаем следующую цепочку преобразований
р* = у]а* + Ь* = ехр(1п^а* + Ь* ). (32)
Введем следующее обозначение а* = 1п(у/а*+&! ), тогда &-ый корень характеристического уравнения (20) можно записать в следующем виде
р* = ехр(а*). (33)
Подставляя выражение (33) в решение (25), получаем решение однородного ОРУ в следующем виде
з /
хг = ^ Л ехр(аг£) + ^ ехр(аг£)(Аг со8(ыг£) + Д (34)
г=1 г=з'+1
где ^ - число реальных корней, £ - по-прежнему дискретный параметр времени. Выражение (34) по форме совпадает с решением однородного ОДУ (31). Таким образом, любому ОРУ /-ого порядка формально может быть сопоставлено ОДУ также /-ого порядка с решениями вида
з /
хг = ^ Л ехр(аг£) + ^ ехр(аг£)(Аг сов(ыг£) + Д (35)
г=1 г=з'+1
где параметры «¿и ¡^ введены для того, что бы подчеркнуть непрерывность полученной модели.
Итак, чисто формально установлена связь между однородным ОРУ (18) и однородным ОДУ (26). Эти модели описывают две линейных системы с одинаковыми ИПХ с единственным отличием: выходной сигнал в первом случае дискретный и непрерывный - во втором. Отметим, что дискретность выходного сигнала не является свойством линейной системы, а только отображает способ регистрации данных.
Преобразование (32) позволяет непосредственно перейти от модели с дискретным выходом к модели с непрерывным выходом. Естественно, что этот переход является только некоторым приближением, т.к. процесс дискретизации должен определенным образом сказываться на оценках параметров линейной системы. Поэтому поставим задачу найти более точные оценки {аг,Шг) однородного ОДУ по аналогу ОРУ (18).
Согласно общей теории [17], выражение (35) можно интерпретировать как результат наложения решений, вырабатываемых 2 подсистемами первого и (I - ]) второго порядков, т.е. работа непрерывной линейной системы может быть представлена следующей суммой ОДУ
з I
- акХк + х'к - 2акхк + (а' + и1)хк = 0, (36)
к=1 к=Л+1
где (•)', (•)" - соответственно первая и вторая производные. О соответствии между теперь уже непрерывным соотношениями (35) и (36) можно убедиться прямой подстановкой компонент решений в соответствующие ОДУ. Рассмотрим, например, ОДУ второго порядка
х' - 2акхк + (ак + Шк)хк = 0 (37)
и выясним влияние на него процесса дискретизации. Представим решение, соответствующие этому уравнению, в виде Хк (г) = е^ак+Шкгде г - мнимая единица. Истинные значения первой и второй производных этого решения равны
х'к (г) = (ак + гШк )е(&к+г*к
х'(г) = (ак + гШк )2е^&к+г^к(38)
Аппроксимация первой производной центральной разностью имеет вид
(&к+г^к)(г+1) - е(ак+г(Ък)^г-1))/2 =
= [(е&к - е-&к)сов(Шк) + г(е&к + е-&к)вт(Шк)]
е(&к)г/2. (39)
При выводе выражения (39) использовались известные формулы Эйлера. Отношение аппроксимированного и истинного значения первых производных будет равно
х'карг (г)
х*аРг(£) = (е^ - е-"*)со5(ы*) + ¿(еА* + е-"*)^п(ы*) = х* (£) 2(а* + ¿ы*)
Л (а* ,ы*). (40)
Аналогичные вычисления для второй производной дают следующие выражение
х*арг(£) = е 2е + е =
= [(е(А* + е-"*)сой(Ы*) - 2 + ¿(е^ - е-"*)згп(ы*)]
Отношение аппроксимированного и истинного значения для вторых производных будет равно
х*У (£) = (еА* + е-"*)с05(ы*) - 2 + ¿(еА* - е-"*)^п(ы*) = х*(£) (а* + ¿ы* )2
/2(а* ,ы*). (42)
С точки зрения линейного подхода, аппроксимация производных разностями представляет линейную фильтрацию с передаточными функциями (39) и (41). В свою очередь, отношения (40) и (42) являются комплексными коэффициентами усиления производных за счет разностной аппроксимации.
/
стем первого и второго порядков. Подсистемы второго порядка, в частности, имеют вид
х*+1 + Ь*х хг + Ь*0 х^1 = 0, (43)
ы* а *
позволяющими представить уравнение (37) в виде
х*сарг - 2а* х*арг + (а2 + ы2)х = 0 (44)
Учитывая соотношения (40) и (42), уравнение (44) преобразуется к следующему виду
х" - 2^к(ак,Шк)акх' + /к(ак,Шк)(а' + и')х = 0.
//
(45)
Так как функции /1(ак,Шк) и /2(ак,Шк) являются вычисляемыми, то выражения (37) и (45) дают решением поставленной выше задачи. В частности
Отметим, что формально коэффициенты ОДУ могут быть оценены без использования уравнения авторегрессии. Действительно, с помощью сплайн аппроксимации можно вычислить производные нужных порядков во всех точках кривой. Используя оптимизирующую процедуру, например, метод наименьших квадратов, затем можно оценить коэффициенты уравнения (26), т.е. решить поставленную выше задачу. Однако такой подход является крайне непродуктивным.
Самым существенным недостатком этой схемы является невозможность оценок точности аппроксимации производных сплайн методами и, следовательно, точности оценок коэффициентов модели (26). Напротив, аппроксимация производных разностями позволяет выполнить такие оценки в частотной области, более того, появляется возможность коррекции оценок. Таким образом, авторегрессионная модель является неявным методом оценок производных наблюденных данных.
2.4.1. Экспериментальная оценка связей между дискретной и непрерывной моделями
Как видно из вышеизложенного, дискретизация данных выхода не меняет структуры линейной модели, а вносит погрешность в оценки ее параметров. Так как вносимые погрешности вычисляются аналитически, они могут быть скомпенсированы. Как видно из выражений (46), погрешности определяются функциями /1(а',Шк) и /2(ак, Шк) которые являются комплексными. Так как оцениваемые параметры являются чисто реальными числами, выделим из этих функций реальную компоненту.
На рис.1 изображены зависимость реальной части функций /1(ак,шк) и /к(ак,Шк) от частоты при различных значениях параметра а. Как видно из этих рисунков, аппроксимация производных разностями приводит к подавлению высоких частот. При этом аппроксимация первой производной подавляет
ак = ~т?—'—г ак }1(ак ,Шк)
/к(ак ,Шк)
Ш1 = 12(ак,Шк)(а\ + Ш2к) - ак.
(46)
высокие частоты сильнее, чем аппроксимация второй. Другими словами, аппроксимация производных разностями занижает значения производных тем сильнее, чем ниже порядок производной. Параметр затухания а влияет на зависимость аппроксимаций производных значительно слабее, чем частота.
1 Д
О -,-,-,-,-,-,-,-,-,-,-,-,-,-,-,-,-,-,-,-,-,-,-,-,-,-,-,-,-,-,-
о.оо о;га 0.41 г' 0.Я1 1,01 1 1.42 1 1.В2 2рз з.гз здз г.ез з
ш
Рис. 1. Коэффициенты усиления производных разностями в зависимости
аа аа
Как видно из формул (46), важно не само значение параметров /\(ак, ¿к) и /к(ак,Шк), а их отношение. Зависимость этого отношения от частоты, при
а
из этого рисунка, отношение на первой трети частотного диапазона имеет практически линейную зависимость, затем крутизна линии начинает сильно возрастать. Это отношение является коррекцией коэффициентов затухания линейной системы, следовательно наибольшему исправлению должны подвергаться высокие частоты. Отношение функций на низких частотах слабо зависит от коэффициента затухания.
Зависимость точности аппроксимации частоты разностной моделью более сложная, чем коэффициента затухания. Учитывая, что по натурным данным коэффициенты затухания близки к нулю, вторая формула соотношения (46) может быть упрощена
- /2(ак,Шк)и2к. (47)
Из формулы (47) и рис. 1 видно, что аппроксимация производных разностями завышает линейные оценки частот. В целом влияние дискретизации
на оценки динамических параметров непрерывной линейной модели следующие. На низких частотах это влияние незначительно, но по мере роста частоты аппроксимация ухудшается. При этом оценки коэффициентов затухания (по модулю) занижаются, а частот - завышаются. Влияние дискретизации на оценки частот более значительное, чем на коэффициенты затухания. Сделанные выводы подтверждаются рис. 9, на котором изображены оценки динамических параметров кривой рис. 3 и их откорректировонных по формулам (46) значения.
о -I-;-,-;-,-,-,-,-,-,-,-;-,-,-,-,-,-,-,-,—,-,-,-,-,-,-;-,-;-,—
О.ОП 0.20 0.41 О.В1 О.В1 1.П1 1.22 1 .42 1 .Б2 1 ,В2 2.ПЗ 2,23 2.43 2.63 2.В4
О)
Рис. 2. Зависимость отношения функций {2{а ш)/ а(а^ш ) от частоты. а=0
а
Исследование на натурных данных На рис. 3 изображены отсчеты некоторого геофизического параметра в 256 точках временной сетки с единичной дискретностью. Данные пронормированы по формуле (6).
По изложенной выше методике определялись частоты (размерностью рад/с) и коэффициенты затухания процесса (размерностью 1/с). В результате обнаружено, что процесс имеет одну нулевую и тридцать ненулевых частот, заключенных в интервале О = [0,^]. Все данные были упорядочены по возрастанию частоты. На рис. 4 изображены, в зависимости от порядкового номера данных: эмпирические оценки частот (1), коэффициенты затухания, увеличенные в 20 раз (2) и относительная энергия (умноженная на 4) составляющих сигнала, связанных с отдельными частотами (3).
Анализ распределения энергии составляющих сигнала по частотам (рис.4, кривая 3) показывает, что в процессе имеется две частоты, несущих наибольшую энергию сигнала. Первая энергонесущая частота (и ~ 0.111 рад/с, а ~ -0.036 1/с) является самой низкой из обнаруженных ненулевых частот
с вкладом в общую энергию сигнала ~58%. Вторая одиночная частота имеет следующие параметры (ш ~ 0, 24 рад/с, а « —0,103 1/с) с вкладом в общую энергию сигнала ~20%.
■1 13 25 (37 /45 73 85 197 109 ,121 «3 145 157 169 181 Л93 205 2171229 2к'|д25В
Рис. 3. Исследуемая экспериментальная временная зависимость.
1\ 2 Д 4 ¡6^-6, 7 8 Э 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
V ^ V \
Рис. 4. Оценки: частот (1), коэффициентов затухания (2) и энергии составляющих сигнала (3).
Следовательно, на двух частотах сосредоточено более 75% общей энергии сигнала. Более детальные расчеты показывают, что на первой трети частотного диапазона сосредоточено более 97 % общей энергии сигнала, а оставшийся диапазон содержит практически белый шум. Анализ рис. 2 и 9 показывает, что на первой трети частотного диапазона влияние процесса дискретизации на оценки динамических параметров разностями является минимальным. Следовательно, во многих случаях за оценку коэффициентов
ОДУ можно принять коэффициенты соответствующего ОРУ. Более точные оценки коэффициенотов ОДУ могут быть получены за счет их коррекции по формулам (46-47), рис. 9.
3. Некоторые постановки задач на основе линейной модели
3.1. Постановка задач предобработки
Переход от авторегрессии к разностному уравнению, а затем ОДУ позволяет эффективно и на твердой теоретической основе решать многие задачи предобработки, анализа и интерпретации данных. Как видно из вышеизложенного, функции ФСР аппроксимируются достаточно простыми аналитическими выражениями, что позволяет продлить их на любой временной интервал. А так как по построению ФСР являются линейно независимыми, они образуют естественный базис разложения.
Предобработку данных по предлагаемой методике можно формализовать следующим образом. Пользуясь, например, рекомендациями, изложенными в [21], подбирается оптимальная по сложности регрессионная модель для заданной временной последовательности. По изложенной выше методике вычисляется набор динамических параметров записи (частоты, коэффициенты затухания). Сообразуясь с поставленными целями предобработки, формируется один или несколько базисов на основе ФСР и выполняется обработка, использующая эти базисы. Некоторые постановки и решения, основанные на указанном принципе, изложены ниже.
Задача аппроксимации ФАК на основе базиса из функций ФСР рассмотрена и решена в работах [7, 12]. Кроме того, в [7] поставлены и решены задачи сглаживания, фильтрации в полосе частот и экстраполяции ФАК временной последовательности. Там же рассмотрена задача аппроксимации ФАК коррелируемого шума. Полученные результаты могут быть использованы для надежного оценивания ФАК, конструировании оптимальных фильтров, например, Колмогорова-Винера.
Основной проблемой при расчете оптимального фильтра является необходимость априорного задания корреляционных свойств сигнальной части данных и помехи. Без этой информации выполнить расчет фильтра невозможно. Классическая теория оптимальной фильтрации [20, 22] позволяет успешно фильтровать только достаточно узкий класс сигналов типа белого шума. Это связано с тем, что в классической теории нет методов, позволяющих выполнить соответствующие корреляционные оценки для для сигналов с распределением, отличным от белого шума.
Как видно из вышеизложенного, динамическое моделирование позволяет рассматривать любой временной ряд как суперпозицию сигналов с определенными динамическими параметрами, которые могут быть оценены. Полученная информация позволяет затем оценить корреляционные свойства любой группы составляющих смеси. Последнее означает, что эта группа компонент сигнала процедурой фильтрации может быть выделена из смеси или отфильтрована. В частности, в работе [13] поставлена задача разделения (сепарация) сигнала на основе ФСР для произвольной комбинации составляющих сигнала. Решение задачи сепарация сигналов позволяет, при необходимости, увеличить разрешение данных, что особенно актуально для сейсморазведки.
В работе [12] поставлена и решена задача высокоразрешающей оценки спектра мощности временной последовательности. Высокое спектральное разрешение достигнуто за счет аналитического продолжения ФАК анализируемого процесса и оцененной на интервале [0,г], на временной интервал (—то, то) . Полученные результаты могут использоваться для изучения тонкой спектральной структуры исследуемой временной последовательности.
Следующий шаг в развитии динамического анализа заключается в переходе к многоканальным (векторным) временным рядам. Такого рода данные возникают при измерении нескольких физических величин в каждой временной точке. Векторные временные ряды содержат о процессе значительно больше информации, чем скалярные. Здесь автоматически возникают такие, например, понятия, как взаимная корреляция, взаимный спектр, частная корреляция [1] и другие, которые несут дополнительную информацию о процессе.
Однако для векторных временных рядов из-за роста вычислительных затрат особенно актуальной становится задача оптимизации вычислений. В связи с этой проблемой в работе [5] поставлена и решена задача создания быстрого алгоритма оценки функций взаимных корреляции и коэффициентов авторегрессии для векторных временных рядов. Полученное решение является шагом к созданию теории динамической обработки векторных временных рядов.
Таким образом, предлагаемый подход к моделированию позволяет: идентифицировать статистические параметры шумов, присутствующих в данных; построить оптимальные фильтры; организовать полосовую фильтрацию; сепарацию сигналов и многие другие задачи.
3.2. Постановка обратных задач
Рассмотренная выше задача сейсмической разведки по идеологии во мно-
гом близка к довольно обширному классу задач, например, гидро и радиолокации, томографии, оптической и электронной микроскопии и некоторым другим. Все эти постановки мы объеденим общим термином - задачи локации.
Принципиально задачу локации в общем случае можно сформулировать в виде следующей схемы. Будем считать, что имеется некоторая станция, которая посылает в пространство зондирующий импульс. При этом форма ЗП описывается уравнением (31). Предположим, что распределение отражающих способностей объектов вдоль луча, по которому распространяется ЗП, отображается функцией /а зависимость этого распределения от времени объясняется разными расстояниями от станции до объектов.
При сделанных предположениях сигнал, регистрируемый на приемниках зондирующей станции, будет определяться сверткой функций (31) и/(1)7 которая и является выходом линейной системы, т.е. регистрируемыми данными. По измеренным данным требуется оценить распределение отражающих объектов в пространстве, их размеры, форму и, возможно, их материальный состав. Примем, что теоретической основой решения данной задачи является интегральное уравнение (3).
Если ЗП представляет короткий импульс, в идеале дельта-функцию Дирака, то регистрируемые данные от разных отражающих объектов будут разделены между собой и никакой проблемы с их идентификацией не возникает, т.е. задача является корректно поставленной. Однако любой реальный ЗП имеет конечную временную длительность, в результате интегрирование, согласно выражения (3), приводит к сглаживанию (размытию) выходной последовательности. Это размытие является следствием интерференции сигналов, отраженных от соседних объектов. Суперпозиция сигналов приводит к тому, что задача становится некорректно поставленной. Таким образом, разрешение зависит от временной продолжительности ЗП. С другой стороны, из волновой оптики следует, что при размерах объектов меньших или сравнимых с длиной волны наблюдается явление дифракция. Таким образом, длина волны является естественной мере разрешающей способности методов локации, построенных на анализе величины сигнала.
Другой аргумент зависимости разрешающей способности методов локации от длины волны ЗП заключается в следующем. По мере развития волнового процесса во времени, он начинает охватывать все большую область пространства. В результате огибающая волновой энергии во фронтальной зоне приобретает характер затухающей функции (рис. 5). В этом случае основная энергия ЗП сосредотачивается на участке 0.5^1 периода видимой длины вол-
ны ЗИ, который можно принять за эффективную длину импульс-посылки.
Из приведенных выше рассуждений следует, что задачи локации в традиционной постановке относятся к обратным и некорректно поставленным в случае, если длина импульс-посылки (или длина волны) ЗИ больше расстояния между отражающими объектами. Именно такая ситуация, например, возникает в сейсмической разведке, томографии близко расположенных объектов и других областях. При расположении объектов на расстояниях, меньших длины импульс-посылки ЗИ, наступает перекрытие сигналов от разных отражающих объектов, т.е. наблюденные данные (5) представляют интерференцию сигналов от различных источников.
В общем случае может быть неизвестной и меняющаяся во времени форма сигнала ЗИ Следовательно, наиболее полная постановка задачи локации заключается в решении уравнения (3) при неизвестной форме сигнала ЗИ В [14] рассматривается постановка задачи анализа сложных динамических систем. В этой же работе решена задача вычисления функции Грина, которая может быть использована для оценки частных решений неоднородного ОДУ. Функция Грина является аналогом ИПХ к{Ъ) и позволяет непосредственно решить уравнение (3), использовав, например, его регуляризацию по Тихонову [25].
Решение поставленной выше задачи локации является достаточно сложной проблемой в случае ее некорректной постановки. Трудности, возникающие при интерпретации интерферирующих записей хорошо видны по работе [26], в которой рассмотрен ряд алгоритмов решения этой задачи. Один из этих алгоритмов основан на итерационном подходе с уточнением на каждом этапе формы ЗИ, амплитудно-кинематических параметров сейсмической записи с использованием метода сингулярного разложения. Однако, не смотря на разносторонний методический материал, используемый в алгоритме, тщательную проработку этапов его работы, в результате проведения численных оценок было установлено, что разрешающая способность метода не выше полупериода ЗИ.
4. Развитие методики
4.1. Ревизия классических постановок задач локации
Выполненные нами оценки функции Грина позволяют непосредственно решить уравнение (3), воспользовавших, например, методом регуляризации по Тихонову [25]. И хотя авторами не выполнено соответствующее моделирование, можно предположить, что получаемые решения будут не очень хоро-
шими.
Дело в том, что любая регуляризация, как правило, приводит к сглаживанию решений. Но при зондировании малоразмерных по сравнению с длиной волны ЗИ объектов, их распределение в пространстве носит хаотический характер с распределением, близким к белому шуму. Следовательно, и решение должно иметь шумоподобный вид, а любое сглаживание таких зависимостей немедленно приводит к их полному разрушению, или к потере разрешающей способности.
Особенно критична низкая разрешающая способность для методов сейсморазведки. Поэтому, несмотря на колоссальную избыточность сейсмических данных, высокую стоимость их получения, в геологии метод считается вспомогательным. Чтобы сейсморазведка приобрела статус самодостаточного метода, требуется повысить его разрешающую способность на 1-2 порядка. Это обстоятельство имеет объективные причины и требует основательной ревизии всей методологической основы обработки сейсмической информации для выявления причин неудач.
Несмотря на актуальность проблемы, иная постановка задач сейсморазведки до сих пор не предложена. Все попытки увеличить разрешение сейсморазведки фактически сводились к усовершенствованию технической стороны измерений. В результате средства для сейсмических измерений стали чрезвычайно сложны, а сами измерения очень дороги. Однако прогресс в повышении разрешения при этом имеет величины порядка процентов, или, в лучшем случае, десятков процентов.
Малый рост качества интерпретации сейсмических измерений, по нашему мнению, заложен в первоначальной и неудачной постановке проблемы. Анализ показывает, что задача сейсморазведки может быть сформулирована, поставлена и решена на принципиально другой методической основе. При этом разрешающая способность метода увеличивается на 1-2 порядка и более, т.е. фактически лимитируется частотой Найквиста. Аналогичная ревизия, вероятно, должна быть проведена и для других сфер применения методов локации.
Существующие способы обработки локационной информации (в том числе традиционное решение уравнения (3) и и алгоритм [26]) относятся преимущественно к так называемым амплитудным методам, в которых основным источником информации является значение (амплитуда) измеренного сигнала. Амплитудные методы практически не используют фазовую информацию, заключенную в измеренном сигнале. Эти методы пытаются восстано-
вить структуру последовательности отражений или анализом экстремальных свойств сигнала, или построением различного рода фильтров, либо аппроксимацией ИПХ за счет введения различной априорной информации о форме ЗИ и прямого решения уравнения свертки (3). Все это удовлетворительно работает для простейших случаев: редкие отражения на расстояниях, сравнимых или больших длины волны ЗИ.
Но в случае интерференции не величина сигнала, а именно его фазовые соотношения несут наибольшую информацию. Поэтому амплитудные методы имеют принципиальное ограничение в разрешающей способности, которая определяется длиной волны зондирующего импульса. Смещение акцентов моделирования с амплитуды на фазу и является потенциалом развития линейного моделирования. Увеличение на 1-2 порядка разрешения методов обработки для интерферирующих записей можно достичь только используя корректный фазовый анализ данных. При этом задача из разряда некорректно поставленных переводится в категорию вполне корректно поставленных, т.к. фаза не обладает столь патологическими свойствами, как величина сигнала. У авторов есть определенные теоретические и реализованные практически проработки в этом направлении [6, 8, 16]. Проделанная работа показывает, что эти проблемы решаемы как в теоретическом, так и практическом плане. Модельные эксперименты показали, что увеличение разрешения методов локации на 1-2 порядка реально.
4.2. Моделирование задачи сейсмической разведки
Чтобы продемонстрировать возможности фазовых методов решений обратных задач, прокомментируем выполненные нами некоторые модельные расчеты с точки зрения задач сейсморазведки. На рис. 7 (кривая 1) изображена модель акустически неоднородной Земной толщи, где по горизонтальной оси отложены условные глубины залегания геологических слоев, а по оси ординат - так называемые коэффициенты отражения границ акустических неоднородностей. На эту пачку сверху падает звуковая волна - сейсмический ЗИ. Форма некоторых ЗИ, используемых при моделировании изображены на рис. 5. Моделью ЗИ является выражение (34) в котором частоты, коэффициенты затухания и амплитуды были получены с помощью датчика случайных чисел.
Акустическая неоднородность среды приводит к тому, что падающая звуковая волна частично отражается от этих слоев пропорционально коэффициентам отражения. Поэтому рис. 7 можно интерпретировать как глубинную последовательность коэффициентов отражения среды. Кривая содержит 256
отсчетов глубин, на каждом дискрете генерировалось случайное число ( с равномерным законом распределения в интервале —0.5 ^ ( ^ 0.5, которое и принималось за коэффициент отражения. Отметим, что положительные значения сигнала означают отражения звуковой волны без изменения фазы, а отрицательные - изменяют ее фазу на п.
^—*
1 3 51 7 а/
Рис. 5. Формы некоторых ЗИ, используемые при моделировании.
Часть энергии падающей звуковой волны проникает через границу раздела сред и распространяется в нижележащие слои Земной толщи, снова отражается от нижележащего слоя и т.д. В результате часть энергии отраженных волн достигает поверхности Земли и регистрируется датчиками. Измеренная датчиками звуковые волны изображены на рис. 6. Эти измерения представляют свертку сейсмических 311 (рис. 5) и коэффициентов отражения слоев среды (рис.7).
Постановка основной задачи сейсморазведки заключается в построении глубинного распределения коэффициентов отражения слоев среды (рис. 7, кривая 1) по измеренным данным (рис. 6). Очевидно, что по своей поставке задача является обратной и поэтому актуальной становится оценка устойчивости решений. Для этих целей на результат свертки 311 и коэффициентов отражения накладывался белый шум различной мощности, который моделировал ошибки измерений.
Отметим, что восстановить довольно хаотический сигнал (рис. 7, кривая 1) с достаточной точностью амплитудными методами принципиально невозможно, поэтому сгенерированные данные, с точки зрения апробирования методики, является наиболее интересными. Как уже отмечалось, любая регуляризация, присущая амплитудным методам, сглаживает решение. Очевидно,
что если применить любые сглаживающие операции к данным (рис. 7, кривая 1), то решения немедленно будут разрушены, или полностью потеряют разрешающую способность.
Рис. 6. "Измеренные"сигналы, используемые для оценки коэффициентов отражения.
Тем не менее, сделаем оценку разрешающей способности амплитудных методов. Если исходить из 1/2 периода волны ЗИ (рис. 5), то разрешение метода будет не более ~130 м. Другую оценку разрешения можно сделать на основе анализа особых точек кривой, например, положения ее экстремумов. В результате получаем, что амплитудные методы позволяют идентифицировать всего около десятка отражающих объектов.
Отметим, что модельная кривая (рис. 7, кривая 1) отражает так называемую тонкую слоистость геологической среды и довольно близка к натурным данным, полученных в скважинах. Для справки отметим, что реальная дискретизация информации, используемая в сейсморазведке, составляет ДЪ = (0.5 ^ 2)10_3с, следовательно, расстояние между отсчетами глубин составляет АН « 4 м (при ДЪ = 1.10_3с и скорости звука V = 4.103 м/с). Длина модельной волны 314 также согласована с реальными данными, получаемых от глубинных отражений. Итак, предполагается, что имеются только "измерения" рис. 6., и никакой другой априорной информации нет. Требуется оценить глубинное распределение коэффициентов отражения.
На рис. 7 (кривая 2) изображена оценка глубинного распределения коэффициентов отражения среды с использованием фазы сигнала. Анализ рис. 7 показывет, что модельные данные и их оценки довольно хорошо согласованы с коэффициентом корреляции гсь > 0.95. На рисунке хорошо видно, что алго-
ритм четко выделяет все особенности модельной кривой, например, моменты смены фазы звуковой волны. Видно, что наибольшие погрешности оценок возникают в экстремальных точках кривой. Таким образом, все отражающие объекты выделяются практически полностью, наблюдается погрешность в оценке величины их отражающей способности. Некоторые дополнительные сведения о поведение оценок сигнала можно получить по рис. 8, на котором изображен некоторый фрагмент рис. 7.
О,Б
■О ,6
Рис. 7. Модельные данные (ряд 1) и оценки (ряд2).
Рис. 8. Фрагмент рис. 7.
Отметим, что на рис.7 изображена некоторая типическая оценка. В результате проведенных вычислительных экспериментов установлено, что для конкретных данных коэффициент корреляции гсь может меняться в интервале [0.9 ^ 0.98]. Причина флуктуации качества оценок пока не установлена. Однако выполненный анализ показал, что чисто формальными методами
разброс оценок может быть уменьшен. Однако эта задача требует более тщательной теоретической проработки и экспериментирования.
Из рисунка рис. 7 видно, что разрешение метода лимитируется в данном случае дискретностью данных (АН ~ 4 м), т.е. частотой Найквиста. Следует ожидать, что при техническом уменьшении дискретности данных, разрешение метода будет увеличиваться, что не противоречит разработанной методике. Отметим, что расчет коэффициентов отражения позволяет автоматически выполнить оценки скорости звука в слоях, их плотности и глубины залегания, т.е. решить основную задачу сейсморазведки.
--
—- 2
1\2/\ 4 7 г- Э^Ш 1112 ]Д14 15 1В 1 7 18 1 9 20 21 22 23 24 25 26 27 28 2Э 30 31
^^ ^—3 4
Рис. 9. 1 - оцененные частоты; 2 - исправленные частоты; 3 - оцененные коэффициенты затухания; 4 - исправленные коэффициенты затухания.
В результате экспериментов по изучению устойчивости получаемых решений установлено, что ошибки измерений до определенного уровня не сказываются на качестве решения. После перехода энергии ошибок измерения некоторого порога, качество решения (коэффициент гсъ) начинает постепенно деградировать. Постепенное снижение качества оценок по мере роста энергии измерительных ошибок, говорит об устойчивости получаемых решений. Другая серия экспериментов по проверке устойчивости решений заключалась в свертывании коэффициентов отражения с 314, имеющими различную форму и видимую частоту (рис. 5). В результате установлено, что постановленные кривые для разных форм 314 практически идентичны, что также говорит о устойчивости получаемых решений.
Следует отметить, что полученные результаты являются предварительными, т.к. теоретические проработки в области фазового анализа данных в настоящее время еще полностью не закончены. Тем не менее, проведенное
моделирование вселяет определенный оптимизм.
Очевидно, что возможность увеличения на несколько порядков разрешающей способности, выводят методы локации на принципиально новый уровень. Например, в настоящее время в мире накоплен гигантский архив сейсмической информации. Его переобработка по новой методике с повышением разрешения на 1-2 порядка является вполне реальной. Экономические, научные и практические выгоды такой переобработки вполне очевидны.
Мы вполне отдаем себе отчет в том, что для локации мелкомасштабных объектов требуется соответствующая энергетика отражений. Однако фазовые методы являются более помехоустойчивыми и чувствительными, чем амплитудные. Вероятно, и здесь речь может идти о величинах в 1-2 порядка. Однако этот вопрос требует более тщательной теоретической проработки и экспериментирования.
Мы рассмотрели задачу активной локации с использованием ЗИ. Аналогичным образом могут быть поставлены и решены задачи так называемой пассивной локации. Известно, что различные объекты, например, автомобили, самолеты и др. являются источниками широкополосных электромагнитных излучений. Другим источником информации может быть акустический шум, генерируемый, например, движением судов в воде. Даже неподвижные объекты, находящиеся в водной среде, изменяют естественный акустический фон окружающей среды, что также является источником информации. Всю эту информацию также можно использовать с целью опознания объектов и оценке их некоторых геометрических и материальных характеристик.
Пассивными методами могут, например, изучаться и так называемые естественные шумы горных пород. Увеличение разрешения позволит, например, проанализировать внутреннюю структуру естественных шумов, выявлять те или иные их особенности. Вполне вероятно, что полученная информация позволит на принципиально новом уровне взглянуть на процессы, происходящие в Земных недрах. К пассивным методам можно также отнести обработку данных о землетрясениях и в этом направлении, вероятно, методика также позволит получить новую информацию.
5. Заключение
Анализ теории показывает, что классический математический анализ по своей идеологии является амплитудным. Аналогичное заключение можно сделать и относительно его дискретных реализаций - численных методов вычислений. Действительно, основой математического анализа являются про-
изводные, т.е. анализ значения функции в точке и ее малой окрестности. Это автоматически накладывает на анализируемые данные довольно серьезные ограничения - функции должны быть непрерывными, достаточно гладкими и т.д. Свойствами гладкости и непрерывности по классическому подходу, должны обладать и получаемые решения.
Если экспериментальные данные не обладают этими свойствами, то обычно выдвигается предположение о наличии в измерениях некоторых шумов. При этом ссылаются на ошибки измерений, ошибки, возникающие при преобразовании данных при переходе к цифровому формату и т.д. Информация сглаживается и после этого производится собственно обработка по классическому сценарию. Действительно реальные экспериментальные данные часто имеют достаточно "рваный" (нерегулярный, негладкий) вид. Но можно ли объяснить такой характер информации только действием измерительных ошибок?
Как правило, реальные измерительные ошибки имеют энергию значительно меньшую, чем им приписывают. В большинстве случаев нерегулярные компоненты сигналов имеет совершенно другую природу. Например, хаотическую флуктуацию данных при непрерывном измерении скорости ветра, температуры воздуха и т.п. высокочувствительными приборами можно объяснить проявлением турбулентности атмосферы, т.е. тонкой внутренней структуры наблюдаемого динамического процесса, а не действием каких то ошибок наблюдений.
Более того, любые природные и другие сложные процессы в соответствующем масштабе времени не являются гладкими (в этом легко убедится, проанализировав, например, различные биржевые индексы, геофизические измерения и т.д.). Чтобы учесть такое поведение информации, современная практика иногда прибегает к нелинейному моделированию. Мы не будем рассматривать достоинства и недостатки такого подхода. Отметим только, что при разработке любой теории или модели должен действовать принцип Оккама, образно выраженный А. Эйнштейном: "Все нужно делать просто, насколько это возможно, но не проще".
В данном случае этот принцип можно сформулировать так: "зачем прибегать к нелинейному, плохо формализованному моделированию, если такие же результаты (а может быть и лучшие) можно достичь более простым и хорошо теоретически обоснованным линейным подходом". Тем более, что процедура вычисления фазы является очень нелинейной операцией. Следовательно, в линейную модель вносится элемент нелинейности, что усиливает как поста-
новку задач, так и ее решение.
Здесь, вероятно, просматривается некоторая аналогия с принципами ней-роподобного моделирования [11]. Из этой теории, в частности, следует, что с помощью линейных операций и последовательного соединения произвольных нелинейных элементов, можно получить любой требуемый результат с любой, наперед заданной точностью. Другими словами, соединение ОДУ с нелинейной операцией вычисления фазы в принципе гарантирует получение точной аппроксимации любой наблюденной зависимости. В цитируемой выше работе, однако, отмечается, что далеко не всегда эта аппроксимация будет отвечать требуемым качествам: однозначности, устойчивости получаемых результатов, адекватными экстраполяционными и обобщающими качествами. Однако сама по себе обнаруженная связь представляется очень интересной и, возможно, объясняет полученные выше результаты.
Таким образом, часто "негладкость"экспериментальных данных является их принципиальным свойством, поэтому любое сглаживание должно быть хорошо обоснованным. Нужно не сглаживать данные, если для этого нет достаточно серьезных причин, а разрабатывать адекватную имеющейся информации математическую модель. В частности, экспериментами установлено, что сглаживание приводит к ухудшению разрешения фазовых методов. Это можно объяснить тем, что сглаживание фактически приводит к той или иной деформации данных, что искусственно меняет фазовые соотношения сигнала.
Отметим, что предложение использовать фазу сигнала при интерпретации данных не представляет совершенно новую и неизвестную процедуру. Так, например, в технических приложениях имеется ряд примеров успешного применения фазовых методов анализа информации. Например, в физике - это голография, в измерительной технике - измерительные системы, основанные на принципах когерентной оптики [8, 9, 10]. Эти казалось бы далекие друг от друга приложения объединяет общая концепция - они используют фазовую информацию. Особенно впечатляют успехи, достигнутые за последние 10-15 лет когерентными измерительными системами. За это время пройден путь от опытных образцов до серийных изделий, которые на порядок и более увеличили потребительские качества измерительных приборов и систем по сравнению с амплитудными методами измерений.
Предпринятая нами попытка использования фазовых методов для решения ряда задач линейного анализа данных показала, что это направление исследований, вероятно, является достаточно перспективным. Метод представ-
ляет ряд этапов преобразования информации, при этом каждый этап имеет собственную математическую модель, в той или иной мере опирающуюся на динамический подход. Начальным этапом этой цепочки преобразований является аппроксимация данных выходом некоторой абстрактной линейной системы. Математическое и методологическое обоснование этой идеи и было изложено в настоящей работе.
Список литературы
[1] Бендат Дж., Пирсол А. Применение корреляционного и спектрального анализа. М.: Мир, 1983.
[2] Блейхуд Р. Быстрые алгоритмы цифровой обработки сигналов. М.: Мир, 1989.
[3] Вапник В.Н., Глазкова Т.Г. и др. Алгоритмы и программы восстановления зависимостей. М.: Наука, 1984.
[4] Губанов В.С. Обобщенный метод наименьших квадратов. СПб.: Наука, 1997.
[5] Драница Ю.П. Об одном методе моделирования нестационарных динамических систем и процессов. //Вестник МГТУ. Тр. Мурм. гос. технич. ун-та. Т. 3, N 1, 2000.
[6] Драница Ю.П. Об одном методе решения некорректно поставленных задач. //Вестник МГТУ. Тр. Мурм. гос. технич. ун-та. Т. 3, N 1, 2000.
[7] Драница Ю.П. Моделирование одномерных динамических процессов с целью предварительной обработки результатов. //Вестник МГТУ. Тр. Мурм. гос. технич. ун-та. Т. 4, № 1, 2001.
[8] Драница Ю.П., Жеребцов В. Д., Слипченко В. А. Частотно-временной метод обработки фазовых измерений в геофизике. //Журнал Измерительная техника, N6, 2001.
[9] Драница Ю.П., Жеребцов В. Д. Об одном методе обработки фазовых измерений в геофизике. // Доклад на VI Международной научно-технической конференции Современные методы и средства океанологических исследований, часть 2, М.: 2001,.
[10] Драница Ю.П., Жеребцов В. Д. и др. Использование волоконно-оптических технологий в геофизике. // Журнал "Геофизика № 6,2002.
[11] Драница Ю.П. Принципы нейроподобного моделирования геофизических объектов и процессов. // Вестник МГТУ. Тр. Мурм. гос. технич. ун-та. Т. 5, N 2, 2002.
[12] Драница Ю.П., Драница А.Ю. Корректный метод оценки спектральной плотности на основе линейной модели. // Матер. Междунар. научно-технич. конф. Наука и образование-2007. Мурманск, МГТУ, 2007. НТЦ Информрегистр 0320700491 от 05.03.2007.
[13] Драница Ю.П., Драница А.Ю. Фильтрация и декомпозиция сигнала на основе метода коллокации. // Матер. Междунар. научно-технич. конф. Наука и образование-2007. Мурманск, МГТУ, 2007. НТЦ Информрегистр 0320700491 от 05.03.2007.
[14] Драница А.Ю., Драница Ю.П. Применение линейной модели для количественной диагностики экономических систем. // Матер. Междунар. научно-технич. конф. Современнные проблемы экономики, управления и юриспруденции-2007. Мурманск, МГТУ, 2007. НТЦ Информрегистр 0320700489 от 05.03. 2007.
[15] Драница Ю.П., Драница А.Ю. Об одной постановке обратной задачи локации. // Матер. Междунар. научно-технич. конф. Наука и образование-2008. Мурманск, МГТУ, 2008. НТЦ Информрегистр 0320800238 от 21.01. 2008.
[16] Драница Ю.П., Драница А.Ю. Некоторые аспекты интерпретации экспериментальных данных на основе теории линейных динамических систем. // Вестник МГТУ. Тр. Мурм. гос. технич. ун-та. Т.12, № 1, 2009.
[17] Г. Корн, Т. Корн Справочник по математике (для научных работников и инженеров). М.: Наука, 1977.
[18] Ф. Клаербоут Джон Теоретические основы обработки геофизической информации. М.: Недра, 1981.
[19] Краус М., Вошни Э. Измерительные информационные системы. М.: Мир, 1975.
[20] Кулханек О. Введение в цифровую фильтрацию в геофизике. М.: Недра, 1981.
[21] С.Л. Марпл-мл. Цифровой спектральный анализ и его приложения. М.: Мир, 1990.
[22] Никитин A.A. Теоретические основы обработки геофизической информации. М.: Недра, 1986.
[23] Рожков В.А., Трапезников Ю.А. Вероятностные модели океанологических процессов. Л.: Гидрометеоиздат, 1990.
[24] Теребиж В.Ю. Введение в статистическую теорию обратных задач. М.: Физматлит, 2005.
[25] Тихонов А.Н., Гончарский A.B., Степанов В.В., Ягола А.Г. Численные метода решения некорректных задач. М.: Наука, 1990.
[26] Троян В.Н, Соколов Ю.М. Методы аппроксимации геофизических данных на ЭВМ. Л.: Из-во Ленинградского Университета, 1989.
[27] Хемминг Р.В. Цифровые фильтры. М.: Недра, 1987.
