CC BY
63
УДК 517.958+517.962.2
Некоторые постановки задач на основе динамического моделирования
1 2 Ю.П. Драница , А.Ю. Драница
1 Политехнический факультет МГТУ, кафедра высшей математики и программного обеспечения электронно-вычислительных машин ЗАО "Ланит ", г. Москва
Аннотация. В работе рассматриваются постановки некоторых задач в рамках динамического моделирования, том числе обратных и некорректно поставленных. Установлено, что классическое решение этих задач осуществляется в рамках так называемых амплитудных методов, разрешающая способность которых принципиально ограничена длиной волны зондирующего импульса. Выполнен анализ, установивший причины малой разрешающей способности амплитудных методов анализа. Выдвинуты предложения, позволяющие на 1-2 порядка увеличить разрешающую способность методик решения некоторых некорректно поставленных задач.
Abstract. The paper considers some problems (tasks) in the framework of dynamic simulation including inverse and ill-posed tasks. It has been established that the classical solution of these tasks can be found in the context of so-called amplitude methods which resolution capability is in principle restricted by a length of a direct impulse wave. The reasons of weak resolution capability of amplitude methods have been analysed. Some proposals allowing to increase 1-2 orders of resolution capability of methods of some ill-posed tasks' solution have been made.
Ключевые слова: обыкновенные линейные дифференциальные уравнения, обратные задачи, некорректно поставленные задачи, задачи локации, сейсморазведка, высокое разрешение
Key words: ordinary linear differential equations, opposite tasks, ill-posed tasks, tasks of location, prospecting seismology, high resolution
1. Введение
Теория линейных систем, основанная на их описании обыкновенными линейными дифференциальными уравнениями (ОДУ), представляет мощный и развитый математический аппарат. Наибольшее применение этот аппарат получил в технических приложениях: электротехнике, электронике и других, и, как правило, для описания достаточно простых технических устройств. Однако попытки формального применения идей, развитых в этих приложениях, для исследования более сложных природных объектов, экономических систем и т.д. натолкнулись на серьезные проблемы. Самыми существенными недостатками прямого применения этих методов к сложным объектам является неустойчивость и низкая разрешающая способность получаемых решений.
В работе проанализирована возникшая ситуация на примере задач локации, к которым относится и сейсморазведка. Установлено, что классические постановки этих задач являются обратными и некорректно поставленными. Традиционно решение этих задач осуществляется в рамках так называемых амплитудных методов, разрешающая способность которых принципиально ограничена длиной волны зондирующего импульса.
2. Общая постановка проблемы интерпретации данных линейной моделью
Рассмотрим линейную стационарную во времени систему, на вход которой воздействует сигнал ft) и описываемую линейным неоднородным ОДУ с постоянными коэффициентами l-го порядка
+ fllx(l-1)(t) + ... + a(i-1)X(t) + a, = ft), (1)
где x(n)(t) - производная n-го порядка выходного сигнала линейной системы; t - параметр типа времени. Однородным ОДУ, сопряженным с (1), является следующее выражение
a0x(I)(t) + a1x(l-1)(t) + ... + a(l-1)x(t) + ai = 0. (2)
Описание с помощью ОДУ имеет тот недостаток, что его постоянные коэффициенты трудно измерить, к тому же выходную характеристику проще вычислить методами операционного исчисления. Поэтому в технических приложениях для описания зависимости (1) применяют так называемую импульсную переходную характеристику (ИПХ) или весовую функцию h(t), доступные непосредственному измерению. Для этого на вход системы подают специальные тестовые сигналы f(t), измеряют выходной сигнал системы и по этой информации строят h(t) (Краус, Вошни, 1975).
С точки зрения теории ОДУ, построенная таким образом весовая функция Н(Г) является частным решением неоднородного уравнения (1). В том случае, если на вход линейной системы подается короткий импульс, в идеальном случае дельта-функция Дирака, то весовая функция Н(Г) в теории ОДУ получила название функции Грина (Корн, Корн, 1977). Будем считать, что в результате специальных экспериментов для линейной системы построена весовая функция Н(Г). Тогда в установившемся режиме сигнал выхода линейной системы можно выразить через так называемый интеграл свертки (Кулханек, 1981)
Ь Ь1
*(0 = |И(Я)/Ц - Я)йЯ + п(Х) = |Щ - Я)/ДОА + п(Х), (3)
а а1
где Л - переменная интегрирования типа времени, п($) - ошибки измерений и моделирования. Уравнение (3) устанавливает связь между входным, выходным сигналами и ИПХ линейной системы. В зависимости от конкретной задачи, пределы интегрирования могут быть константами, переменными, или несобственными числами ±®. На уравнении (3) удобно продемонстрировать различные варианты постановок линейных задач и особенности их решения, в частности, постановки прямых и обратных задач, фильтрации данных и другие.
2.1. Некоторые классические постановки задач
2.1.1. Задача линейной оптимальной фильтрации данных
Отметим, что уравнение (3) в теории линейных систем часто рассматривается как некоторый фильтр с системной функцией Н(Г) (Кулханек, 1981; Хемминг, 1987). Особое внимание заслуживают так называемые оптимальные фильтры, например, фильтр Колмогорова-Винера. Оптимальная фильтрация в широком смысле предназначена для решения обширного круга задач (Никитин, 1986): фильтрации шумовой компоненты сигнала; фильтрации в заданной полосе частот; сглаживания функций; конструирования так называемых инверсных фильтров и других.
При построении оптимальных фильтров основная проблема заключается в необходимости априорного задания корреляционных свойств шумовой компоненты и полезного сигнала. Классическая теория оптимальной фильтрации в настоящее время позволяет построить фильтры только для достаточно узкого класса шумов, близких к белому, или некоррелируемому (Никитин, 1986; Теребиж, 2005). Это связано с тем, что в классической постановке не существует методов идентификации шумов с другими законами распределения.
2.1.2. Формулировка прямой задачи
Для иллюстрации прямых и обратных постановок рассмотрим типичную задачу сейсмического зондирования земной коры. Зондируемую толщу Земли обычно представляют как последовательность слоев геологических пород с неизвестными наклонами границ и толщинами, плотностями вещества и скоростями звука. На поверхности возбуждается сейсмический зондирующий импульс (ЗИ), который распространяется внутрь Земли, частично отражается от ее слоев, вновь достигает поверхности и регистрируется системой датчиков, расставленных на поверхности.
Допустим, что акустические (и любые другие) свойства земных пород известны в любой точке трехмерного пространства г и описываются функцией /(г). В точке г0 земной поверхности возбуждается ЗИ. При этом известна форма ЗИ и ее эволюция во времени, т.е. задана функция И(г0, /). Закономерность распространения ЗИ в пространстве тоже известна и описывается, например, волновым уравнением.
Для указанного примера прямая задача заключается в нахождении отклика среды у(г^) на сейсмический импульс в любой момент времени /, для любой точке пространства г, если известно, что распространение ЗИ в пространстве подчиняется, например, волновому уравнению. Таким образом, классическая постановка прямой задачи заключается в задании дифференциального или интегрального уравнения, по которому развивается процесс в пространстве и во времени, а также начальных и граничных условий. В рассматриваемом случае роль начальных и граничных условий выполняют функции /(г), И(г0,(). Очевидно, что к прямой задаче можно отнести и задачу фильтрации данных (при известной функции фильтра ^(/)).
В частности, из приведенного примера видно, что для решения прямой задачи требуется задание большого объема априорной информации. Отметим, что для многих практически важных постановок эта априорная информация фактически и является искомым решением. Вполне очевидно, что прямая задача, с точки зрения сейсмической разведки, интереса не представляет.
2.1.3. Обратные и некорректно поставленные задачи
Обратная задача заключается в определении свойств объекта по наблюденным данным и априорной информации, которая имеется у исследователя (Теребиж, 2005). Для приведенного выше
примера требуется по зафиксированному на поверхности Земли отклику слоистой среды на сейсмический импульс определить структуру и параметры зондируемой среды: толщины и глубины залегания отдельных слоев, их плотности, акустические скорости, вещественный состав и другую информацию. Другими словами - восстановить функцию начальных условий fr). Первичной или основной задачей сейсморазведки является выявление акустически неоднородных слоев земной толщи и глубин их залегания - так называемая задача отбивки границ.
Рассмотрим особенности постановки обратной задачи на примере решения уравнения (3). Как указывалось выше, в технических приложениях весовая функция h(t) может быть измерена непосредственно. Следовательно, оценке подлежит только функция f(t), т.е. одна из подынтегральных функций. В этом случае уравнение (3) называется интегральным.
При решении обратных задач важную роль играет их устойчивость. Задача устойчива, если малым флуктуациям левой части уравнения (3) соответствуют малые флуктуации решений. Если задача изначально является неустойчивой, то решать ее нет смысла, поскольку погрешности алгоритмов, накапливающиеся в ходе решения численными методами, неизбежно приведут к тому, что будет найдено неверное решение.
Как правило, любые экспериментальные данные характеризуются наличием шумов, которые коренным образом меняют идеологию решения обратных задач (Теребиж, 2005). Если сама задача является устойчивой, то существование шума может устойчивость нарушать. В результате двум различным сигналам f и f2, искаженным шумом, будут соответствовать очень похожие измерения x1 и x2 (и наоборот). Поэтому встает вопрос, можно ли извлечь из измерений полезную информацию о сигнале, если наличие шума делает задачу неустойчивой? Такие постановки называются задачами, поставленными некорректно, и для их решения развиты специальные методы, основанные на привлечении дополнительной априорной информации о решении.
Следует также отметить, что класс некорректно поставленных задач шире класса обратных, т. к. неустойчивость могут иметь и прямые задачи. Поскольку экспериментальные данные по необходимости случайны, естественная постановка обратных задач достигается в рамках теории статистического оценивания неизвестных параметров. Именно стохастичность модели формирования данных обусловливает фундаментальную трудность, возникающую при обращении информации - неустойчивость решения.
Неустойчивость означает, что в пределах естественных флуктуаций шума с наблюдаемыми данными примерно в равной мере согласуется множество возможных оценок исходного объекта, включая и решения, существенно отличающиеся от него.
2.2. Формальная интерпретация некорректно поставленных задач
Интегральное уравнение (3) удобно представить в следующем операторном виде
A ft) = x(t) + n(t), (4)
где А - оператор отображения, n(t) - ошибка наблюдений и моделирования. В этом случае постановка задачи решения уравнения (4) состоит в следующем (Вапник и др., 1984). Пусть задан непрерывный оператор, который взаимно однозначно отображает элементы ft) метрического пространства E1 в элементы x(t) метрического пространства E2. Требуется найти решение операторного уравнения (4) в классе функций f(t), если функция ft) неизвестна, но имеются измерения x1, x2, ... ,х/ функции x(t) в I точках t1, t2,...,4
Отметим, что понятия "измерения" в контексте данной проблемы трактуются расширительно и могут включать не только какие-либо измерения физических величин, но и, например, результаты вычислительного эксперимента, решения прямых задач, какие-либо эмпирические зависимости, графики, законы физики и любую другую априорную информацию, которую предполагается использовать при построении зависимостей. Оператор А уравнения (4) может представлять любое правило, по которому осуществляется отображение функциональных зависимостей. Мы же будем рассматривать отображения, осуществляемые посредством дифференциальных или интегральных уравнений.
Говорят, что решение операторного уравнения (4) устойчиво, если малая вариация правой части x(t) приводит к малому изменению решения. Это означает, что для всякого достаточно малого е найдется такое ¿(е), что как только выполняется условие
P2[x(t), x(t+e)] < ¿(е)
окажется выполненным и неравенство
Pift), ft+e)] < е.
Здесь P1, P2 означают расстояния, определяемые в метриках E1 и E2, соответственно. Говорят также, что задача решения операторного уравнения (4) поставлена корректно по Адамару, если его
решение существует, единственно и устойчиво. Если не выполняется хотя бы одно из перечисленных требований, то говорят, что задача решения операторного уравнения (4) поставлена некорректно.
Типичной некорректно поставленной задачей является, например, интегральное уравнение (3). Уравнение (3), в зависимости от пределов интегрирования, является интегральным уравнением Фредгольма или Вольтера первого рода, отображающим множество непрерывных на отрезке [0, г] функций ft) на множество непрерывных на этом же отрезке функций x(t). Как показано в (Вапник и др., 1984), решение уравнения (3) является некорректно поставленным в том смысле, что сколь угодно малым изменениям наблюденной функции x(t) могут соответствовать сколь угодно большие изменения зависимостиf(t).
Природа этого факта заключается в том, что интеграл в правой части уравнения (3) является сглаживающим оператором. Это сглаживание приводит к тому, что большие по амплитуде возмущения функции f(t), но сосредоточенные на малом участке пространства E1, переводятся в малые возмущения функции x(t) (тем меньшие, чем меньше этот участок). При решении обратной задачи, чтобы объяснить малые возмущения функции x(t), приходится предполагать значительные отклонения функции ft).
Современная практика имеет два подхода к решению этой проблемы (Вапник и др., 1984). Первый метод заключается в разложении функций по заранее фиксированному конечному базису. В этом случае решение зависит от конечного числа параметров, и если число наблюдений за функцией превышает это число, то решение находится методом наименьших квадратов. Во втором случае рассматривается класс равномерно ограниченных функций. Тогда решение может быть найдено введением регуляризирующего функционала (Тихонов и др., 1990).
2.3. Особенности обратных постановок для сложных объектов и процессов
В теории кибернетики существует методический прием, основанный на так называемом принципе "черного ящика", предполагающим моделирование по наблюдению входов и выходов. С точки зрения разработчика или пользователя модели, структура объекта спрятана в "черном ящике", который имитирует его поведенческие особенности.
В данном случае "черным ящиком" являются ОДУ или ИПХ, входным сигналом - функция f(t), выходным - наблюденное значение сигнала (4). Тогда принцип "черного ящика" можно сформулировать следующим образом: по замеренным входным и выходным сигналам линейной системы требуется аппроксимировать ее ИПХ. Знание ИПХ позволяет однозначно связать входные и выходные сигналы системы.
Однако получить ИПХ для естественных природных, экономических или гуманитарных систем, а также в случае сейсмических измерений, подачей тестовых сигналов на входы изучаемых объектов не представляется возможным. Для сейсморазведки, например, определение формы ЗИ затруднено сложностью законов распространения волнового импульса в земной толще. Геологические породы являются своеобразным фильтром высоких частот с переменными параметрами, в результате форма ЗИ меняется в пространстве и во времени, т.е. h(r, t). И эти изменения могут быть значительными.
В этих условиях мы вынуждены предполагать, что любые измеренные данные являются отображением некоторых выходов абстрактной линейной системы, на входы которой поступают некоторые внешние сигналы. Эти входные сигналы преобразуются линейной системой в соответствии с ее ИПХ (или ОДУ) и поступают на ее выходы.
Таким образом, из трех объектов, определяющих поведение сложных систем, два могут быть принципиально неизмеряемы. Поэтому для сложных систем, по сравнению с принципом "черного ящика", возникает общая постановка задачи: по наблюденным данным восстановить как ИПХ системы, так и ее входные сигналы. Из этих рассуждений следует, что для интерпретации сложных систем и объектов требуется априорная оценка ИПХ, основанная только на наблюденных "выходных" данных. Метод такой оценки предлагается ниже.
3. Предлагаемые принципы линейного моделирования
Пусть имеются измерения временной последовательности некоторых физических величин, выполненные в n точках, с постоянной дискретностью At. Эти измерения образуют временной ряд
x = (Х1, Х2, ... , Xn). (5)
В работах (Драница, 2001; Драница, Драница, 2007a; 2009) предложено представлять данные (скалярные или векторные) временного ряда (5) выходом некоторой абстрактной линейной динамической системы, работа которой описывается ОДУ (1). Там же поставлена и решена задача аппроксимации коэффициентов дифференциальных уравнений (1, 2). Все это позволяет выполнить оценки динамических параметров линейной системы (собственных частот и коэффициентов затухания линейной системы) и сформировать естественный базис. Этот базис представляет фундаментальную
систему решений (ФСР) однородного ОДУ (2) и состоит их затухающих экспонент и затухающих по экспоненте гармоник различных частот.
Построенная таким образом динамическая модель позволяет рассматривать измеряемые сигналы выходов как результат воздействия некоторых развивающихся во времени неизвестных сигналов, приложенных к входам линейной системы. Так как оценка параметров линейной системы выполнена, вся неизвестная (следовательно, искомая) компонента сосредоточена на ее входах. В результате мы приходим к постановке типичной обратной задачи, а именно: по измеренным выходным значениям и оцененным ОДУ выполнить оценки входных сигналов. Эта постановка является единообразной для любых типов данных, лишь бы они являлись временной последовательностью.
Такая постановка дает возможность полностью описать линейную систему, однако интерпретация полученных результатов может встретить определенные трудности. Действительно, полученная в результате моделирования и принятая к рассмотрению линейная система является довольно абстрактным объектом, для которого вполне определенный физический смысл могут иметь только его выходные сигналы. Вся другая информация, в том числе и сигналы входов, является результатом моделирования, тем не менее, эта информация является вполне объективной в том смысле, что опирается только на фактические данные и минимальные априорные предположения. Однако конкретная интерпретация модельной информации полностью ложится на разработчика или пользователя модели.
4. Некоторые постановки задач на основе линейной модели
4.1. Постановка задач предобработки
В работе (Драница, Драница, 2009) излагаются постановки и решения некоторых задач предобработки данных. Основной методикой решения поставленных задач является математический аппарат, позволяющий осуществить достаточно разнообразные операции манипулирования ФАК. Отметим, что успешность их решения во многом определяется тем, что исчерпывающей информацией о параметрах линейной динамической системы является ее ФАК (Губанов, 1997). ФАК по сути дела является физическим понятием, так как с ее помощью можно восстановить спектр мощности процесса и синтезировать динамическую модель случайной последовательности. ФАК позволяет установить распределение энергии процесса по частотам, его устойчивость, наличие периодичностей, характер их затухания и проч.
Физически и методически обоснованный подход к постановке и решению некоторых задач предобработки, в частности, позволяет: аппроксимировать параметры ОДУ (1); идентифицировать статистические параметры шумов, присутствующих в данных; организовать построение оптимальных фильтров; осуществить полосовую фильтрацию данных, а также их сепарацию и многие другие задачи.
4.2. Постановка обратных задач
Рассмотренная выше задача сейсмической разведки по идеологии во многом близка к довольно обширному классу задач, например, гидро- и радиолокации, томографии, оптической и электронной микроскопии и некоторым другим. Все эти постановки мы объединим общим термином - задачи локации.
Принципиально задачу локации в общем случае можно сформулировать в виде следующей схемы. Будем считать, что имеется некоторая станция, которая посылает в пространство зондирующий импульс. При этом форма ЗИ описывается функцией Н((). Предположим, что распределение отражающих способностей объектов вдоль луча, по которому распространяется ЗИ, отображается функцией п(0, а зависимость этого распределения от времени объясняется разными расстояниями от станции до объектов.
При сделанных предположениях сигнал, регистрируемый на приемниках зондирующей станции, будет определяться сверткой функций Н(() и ц($), которая и является выходом линейной системы, т.е. регистрируемыми данными. По измеренным данным требуется оценить распределение отражающих объектов в пространстве, их размеры, форму и, возможно, их материальный состав. Примем, что теоретической основой решения данной задачи является интегральное уравнение (3).
Если ЗИ представляет короткий импульс, в идеале дельта-функцию Дирака, то регистрируемые данные от разных отражающих объектов будут разделены между собой во времени, и никакой проблемы с их идентификацией не возникает, т.е. задача является корректно поставленной. Однако любой реальный ЗИ имеет конечную временную длительность и в результате интегрирование выражения (3) приводит к сглаживанию (размытию) выходной последовательности.
Это размытие является следствием интерференции сигналов, отраженных от соседних объектов. Суперпозиция сигналов приводит к тому, что задача становится некорректно поставленной. Таким образом, разрешение зависит от временной продолжительности ЗИ. С другой стороны, из волновой оптики следует, что при размерах объектов, меньших или сравнимых с длиной волны, наблюдается явление дифракции, которое также ограничивает разрешение длиной волны ЗИ. Таким образом, длина
волны является естественной мерой разрешающей способности методов локации, построенных на анализе величины сигнала.
Другая причина зависимости разрешающей способности методов локации от длины волны ЗИ заключается в следующем. По мере развития волнового процесса во времени он начинает охватывать все большую область пространства. В результате огибающая волновой энергии во фронтальной зоне приобретает характер затухающей функции (рис. 1 кривая 1). В этом случае основная энергия ЗИ сосредотачивается на участке порядка длины волны ЗИ, который можно принять за эффективную длину импульс-посылки.
Из приведенных выше рассуждений следует, что задачи локации в традиционной постановке относятся к обратным и некорректно поставленным в том случае, если длина импульс-посылки (или длина волны) ЗИ больше расстояния между отражающими объектами. Именно такая ситуация, например, возникает в сейсмической разведке, томографии близко расположенных объектов и других областях.
В общем случае может быть неизвестной и меняющаяся во времени форма сигнала ЗИ Н(Г). Следовательно, наиболее полная постановка задачи локации заключается в решении уравнения (3) при неизвестной форме сигнала ЗИ Н((). В (Драница, Драница, 2007б) рассматривается постановка задачи анализа сложных динамических систем. В этой же работе решена задача вычисления функции Грина, которая может быть использована для оценки частных решений неоднородного ОДУ. Функция Грина является аналогом ИПХ Н(Г) и позволяет непосредственно решить уравнение (3), применив, например, регуляризацию решения по Тихонову (Тихонов и др., 1990).
Отметим, что решение поставленной выше задачи локации является достаточно сложной проблемой в случае ее некорректной постановки. Трудности, возникающие при интерпретации интерферирующих записей, хорошо видны по работе (Троян, Соколов, 1989), в которой рассмотрен ряд алгоритмов решения этой задачи. Один из этих алгоритмов основан на итерационном подходе с уточнением на каждом этапе формы ЗИ, амплитудно-кинематических параметров сейсмической записи с использованием метода сингулярного разложения. Однако, несмотря на разносторонний методический материал, используемый в алгоритме, тщательную проработку этапов его работы, в результате проведения численных оценок было установлено, что разрешающая способность метода не выше полупериода ЗИ.
Отметим, что сложность анализа интерферирующих сигналов даже побудила авторов (Вапник и др., 1984) отнести подобного типа задачи к нелинейным постановкам. Анализ проблем, возникающих при решении некорректно поставленной задачи (3), излагается ниже.
5. Развитие методики
5.1. Ревизия классических постановок задач локации
Выполненные нами оценки функции Грина позволяют непосредственно решить уравнение (3), воспользовавшись, например, методом регуляризации по Тихонову (Тихонов и др., 1990). И хотя авторами не выполнено соответствующее моделирование, можно предположить, что получаемые решения будут не очень хорошими.
Дело в том, что любая регуляризация, как правило, приводит к сглаживанию решений. Но при зондировании малоразмерных по сравнению с длиной волны ЗИ объектов их распределение в пространстве носит хаотический характер с распределением, близким к белому шуму. Следовательно, и решение должно иметь шумоподобный вид, а любое сглаживание таких зависимостей немедленно приводит к их полному разрушению или к потере разрешающей способности.
Особенно критична низкая разрешающая способность для методов сейсморазведки. Поэтому, несмотря на колоссальную избыточность сейсмических данных, высокую стоимость их получения, в геологии метод считается вспомогательным. Чтобы сейсморазведка приобрела статус самодостаточного метода, требуется повысить его разрешающую способность на 1-2 порядка. Это обстоятельство имеет объективные причины и требует основательной ревизии всей методологической основы обработки сейсмической информации для выявления причин неудач.
Несмотря на актуальность проблемы, иная постановка задач сейсморазведки до сих пор не предложена. Все попытки увеличить разрешение сейсморазведки фактически сводились к усовершенствованию технической стороны измерений. В результате средства для сейсмических измерений стали чрезвычайно сложны, а сами измерения очень дороги. Однако прогресс в повышении разрешения при этом имеет величины порядка процентов, или, в лучшем случае, десятков процентов.
Малый рост качества интерпретации сейсмических измерений, по нашему мнению, заложен в первоначальной неудачной постановке проблемы. Анализ показывает, что задача сейсморазведки может быть сформулирована, поставлена и решена на принципиально другой методической основе. При этом разрешающая способность метода увеличивается на 1-2 порядка и более, т.е. фактически лимитируется частотой Найквиста. Аналогичная ревизия, вероятно, должна быть проведена и для других сфер применения методов локации.
Существующие способы обработки локационной информации (в том числе традиционное решение уравнения (3) и алгоритм (Троян, Соколов, 1989)) относятся преимущественно к так называемым амплитудным методам, в которых основным источником информации является значение (амплитуда) измеренного сигнала. Амплитудные методы практически не используют фазовую информацию, заключенную в измеренном сигнале. Эти методы пытаются восстановить структуру последовательности отражений или анализом экстремальных свойств сигнала, или построением различного рода фильтров, либо аппроксимацией ИПХ за счет введения различной априорной информации о форме ЗИ и прямого решения уравнения свертки (3). Все это удовлетворительно работает для простейших случаев: редкие отражения на расстояниях, сравнимых или больших длины волны ЗИ.
Но в случае интерференции не величина сигнала, а именно его фазовые соотношения несут наибольшую информацию. Поэтому амплитудные методы имеют принципиальное ограничение в разрешающей способности, которая определяется длиной волны зондирующего импульса. Смещение акцентов моделирования с амплитуды на фазу и является потенциалом развития линейного моделирования. Увеличение на 1-2 порядка разрешения методов обработки для интерферирующих записей можно достигнуть за счет использования корректного фазового анализа данных. При этом задача из разряда некорректно поставленных задач переводится в категорию вполне корректно поставленных, т.к. фаза не обладает столь патологическими свойствами, как величина сигнала. У авторов есть определенные теоретические и практически реализованные проработки в этом направлении (Драница, 2000; Драница и др., 2001; Драница, Жеребцов, 2001; Драница, Драница, 2008). Проделанная работа показывает, что эти проблемы решаемы как в теоретическом, так и практическом плане.
Временные отсчеты Временные отсчеты
Рис. 1. Форма зондирующего импульса (1) и Рис. 2. Модельный шум (1) и
его свертка с коэффициентами отражений (2) его оценка (2)
5.2. Моделирование задачи сейсмической разведки
Чтобы продемонстрировать возможности фазовых методов решения обратных задач, прокомментируем выполненные нами некоторые модельные расчеты (рис. 1-2) с точки зрения задач сейсморазведки. На рис. 2 (кривая 1) изображена часть модели акустически неоднородной земной толщи, где по горизонтальной оси отложены условные глубины залегания геологических слоев, а по оси ординат - так называемые коэффициенты отражения границ акустических неоднородностей. На эту пачку сверху падает звуковая волна, форма которой изображена на рис. 1 (кривая 1).
Акустическая неоднородность среды приводит к тому, что падающая звуковая волна частично отражается от этих слоев пропорционально коэффициентам отражения. Поэтому модель (рис. 2) можно интерпретировать как глубинную последовательность коэффициентов отражения среды. Модель содержит 256 отсчетов глубин, на каждом дискрете генерировалось случайное число £ с равномерным законом распределения в интервале £=[-0.5^0.5], которое и принималось за коэффициент отражения.
Часть энергии падающей звуковой волны проникает через границу раздела сред и распространяется в нижележащие слои земной толщи, снова отражается от нижележащего слоя и т.д. В результате часть энергии отраженных волн достигает поверхности Земли и регистрируется датчиками. Измеренная датчиками звуковая волна изображена на рис. 1 (кривая 2). Эти измерения представляют свертку сейсмического ЗИ и коэффициентов отражения границ среды.
Постановка основной задачи сейсморазведки заключается в построении глубинного распределения коэффициентов отражения слоев среды (рис. 2, кривая 1) по измеренным данным (рис. 1, кривая 2). Очевидно, что по своей поставке задача является обратной, и поэтому актуальной становится оценка устойчивости решений. Для этих целей на результат свертки ЗИ и коэффициентов отражения накладывался белый шум различной мощности, который моделировал ошибки измерений.
Отметим, что восстановить довольно хаотический сигнал (рис. 2, кривая 1) с достаточной точностью амплитудными методами принципиально невозможно, поэтому сгенерированные данные, с точки зрения апробирования методики, являются наиболее интересными. Как уже отмечалось, любая регуляризация, присущая амплитудным методам, сглаживает решение. Очевидно, что если применить
любые сглаживающие операции к данным (рис. 2, кривая 1), то решения немедленно будут разрушены или полностью потеряют разрешающую способность.
Тем не менее, сделаем оценку разрешающей способности амплитудных методов. Если исходить из 1/2 периода волны ЗИ (рис. 1, кривая 1), то разрешение метода будет не более ~130 м. Другую оценку разрешения можно сделать на основе анализа особых точек кривой, например, положения экстремумов. В результате получаем, что амплитудные методы позволяют идентифицировать всего около десятка отражающих объектов.
Отметим, что модельная кривая (рис. 2, кривая 1) отражает так называемую тонкую слоистость геологической среды. По крайней мере, замеры некоторых геофизических параметров, получаемых в скважинах, внешне довольно близки к этой кривой. Для справки отметим, что реальная дискретизация информации, используемая в сейсморазведке, составляет А/ = (0.5+2)-10-3 с, следовательно, расстояние между отсчетами глубин составляет АН ~ 4 м (при А/ = 10-3 с и скорости звука V = 4-103 м/с). Итак, предполагается, что имеются только "измерения" рис. 1. (кривая 2), и никакой другой априорной информации нет. Требуется оценить глубинное распределение коэффициентов отражения.
На рис. 2 (кривая 2) изображена оценка глубинного распределения коэффициентов отражения среды с использованием фазы сигнала. Анализ рис. 2 показывает, что модельные данные и их оценки довольно хорошо согласованы с коэффициентом корреляции гсЪ > 0.95. На рисунке хорошо видно, что алгоритм четко выделяет все особенности модельной кривой, например, моменты смены фазы звуковой волны. Видно, что наибольшие погрешности оценок возникают в экстремальных точках кривой.
Таким образом, все отражающие объекты выделяются практически полностью, наблюдается погрешность в оценке величины их отражающей способности.
Отметим, что на рис. 2 изображена некоторая типическая оценка. В результате проведенных вычислительных экспериментов установлено, что для конкретных данных коэффициент корреляции гсЪ может меняться в интервале [0.9 + 0.98]. Причина флуктуации качества оценок пока не установлена. Выполненный нами анализ показал, что чисто формальными методами разброс оценок может быть уменьшен. Однако эта задача требует более тщательной теоретической проработки и экспериментирования.
Из рис. 2 видно, что разрешение метода лимитируется в данном случае дискретностью данных (АН ~ 4 м). Следует ожидать, что при техническом уменьшении дискретности данных разрешение метода будет увеличиваться, что не противоречит разработанной методике.
В результате экспериментов по изучению устойчивости получаемых решений установлено, что ошибки измерений до определенного уровня не сказываются на качестве решения. После перехода энергии ошибок измерения через некоторый порог качество решения (коэффициент гсЪ) начинает постепенно деградировать. Постепенное снижение качества оценок по мере роста энергии измерительных ошибок говорит об устойчивости получаемых решений. Другая серия экспериментов по проверке устойчивости решений заключалась в свертывании коэффициентов отражения с ЗИ, имеющими различную форму и видимую частоту. В результате установлено, что восстановленные кривые для разных форм ЗИ практически идентичны, что также говорит об устойчивости получаемых решений.
Следует отметить, что полученные результаты являются предварительными, т.к. теоретические проработки в области фазового анализа данных в настоящее время еще полностью не закончены. Тем не менее, проведенное моделирование вселяет определенный оптимизм.
Очевидно, что возможность увеличения разрешающей способности на порядки выводит методы локации на принципиально новый уровень. Например, в настоящее время в мире накоплен гигантский архив сейсмической информации. Переобработка этого материала по новой методике с повышением разрешения на 1-2 порядка является вполне реальной. Практические выгоды такой переобработки информации с точки зрения экономики, науки и практики совершенно очевидны.
Мы вполне отдаем себе отчет в том, что для локации мелкомасштабных объектов требуется соответствующая энергетика отражений. Однако фазовые методы являются более помехоустойчивыми и чувствительными, чем амплитудные. Вероятно, и здесь речь может идти о величинах в 1-2 порядка. Однако этот вопрос требует более тщательной теоретической проработки и экспериментирования.
Мы рассмотрели задачу активной локации с использованием ЗИ. Аналогичным образом могут быть поставлены и решены задачи так называемой пассивной локации. Известно, что различные объекты, например, автомобили, самолеты и др., являются источниками широкополосных электромагнитных излучений. Другим источником информации может быть акустический шум, генерируемый, например, движением судов в воде. Даже неподвижные объекты, находящиеся в водной среде, изменяют естественный акустический фон окружающей среды, что также является источником информации. Всю эту информацию также можно использовать с целью опознания объекта и оценке их некоторых геометрических и материальных характеристик.
Пассивными методами могут, например, изучаться и так называемые естественные шумы горных пород. Увеличение разрешения позволит проанализировать внутреннюю структуру естественных шумов, выявлять те или иные их особенности. Вполне вероятно, что полученная информация позволит на принципиально новом уровне взглянуть на процессы, происходящие в земных недрах. К пассивным методам можно также отнести обработку данных о землетрясениях и в этом направлении, вероятно, методика также позволит получить новую информацию.
6. Заключение
Анализ теории показывает, что классический математический анализ по своей идеологии является амплитудным. Аналогичное заключение можно сделать и относительно его дискретных реализаций - численных методов вычислений. Действительно, основой математического анализа являются производные, т.е. анализ значения функции в точке и ее малой окрестности. Это автоматически накладывает на анализируемые данные довольно серьезные ограничения: функции должны быть непрерывными, достаточно гладкими и т.д. Свойствами гладкости и непрерывности, по классическому подходу, должны обладать и получаемые решения.
Если экспериментальные данные не обладают этими свойствами, то обычно выдвигается предположение о наличии в измерениях некоторых шумов. При этом ссылаются на ошибки измерений, ошибки, возникающие при преобразовании данных при переходе к цифровому формату и т.д. В результате информация сглаживается, и после этого производится собственно обработка по классическому сценарию. Действительно реальные экспериментальные данные часто имеют достаточно "рваный" (нерегулярный, негладкий) вид. Но можно ли объяснить такой характер информации только действием измерительных ошибок?
Как правило, реальные измерительные ошибки имеют энергию, значительно меньшую, чем им приписывают. В большинстве случаев нерегулярные компоненты сигналов имеет совершенно другую природу. Например, хаотическую флуктуацию данных при непрерывном измерении скорости ветра, температуры воздуха и т. п. высокочувствительными приборами можно объяснить проявлением турбулентности атмосферы, т.е. тонкой внутренней структуры наблюдаемого динамического процесса, а не действием ошибок наблюдений.
Более того, любые природные и другие сложные процессы в соответствующем масштабе времени не являются гладкими (в этом легко убедиться, проанализировав, например, различные биржевые индексы, геофизические измерения в скважинах и т.д.). Чтобы учесть такое поведение информации, современная практика иногда прибегает к нелинейному моделированию. Мы не будем рассматривать достоинства и недостатки такого подхода. Отметим только, что при разработке любой теории или модели должен действовать принцип Оккама, образно выраженный А. Эйнштейном: "Все нужно делать просто, насколько это возможно, но не проще".
В данном случае этот принцип можно сформулировать так: "зачем прибегать к нелинейному, плохо формализованному моделированию, если такие же результаты можно достичь более простым линейным подходом". Тем более, что процедура вычисления фазы является очень нелинейной операцией. Следовательно, в линейную модель вносится элемент нелинейности, что усиливает как постановку задач, так и их решение.
Здесь можно увидеть некоторую аналогию с принципами нейроподобного моделирования (Драница, 2002). Из этой теории, в частности, следует, что с помощью линейных операций и последовательного соединения произвольных нелинейных элементов можно получить любой требуемый результат с любой, наперед заданной точностью. Другими словами, соединение ОДУ с нелинейной операцией вычисления фазы в принципе гарантирует получение точной аппроксимации любой наблюденной зависимости. В цитируемой выше работе, однако, отмечается, что далеко не всегда эта аппроксимация будет обладать требуемыми качествами: однозначностью, устойчивостью получаемых результатов, адекватными экстраполяционными и обобщающими качествами. Однако сама по себе обнаруженная связь представляется очень интересной и, возможно, объясняет полученные выше результаты.
Таким образом, часто "негладкость" экспериментальных данных является их принципиальным свойством, поэтому любое сглаживание должно быть хорошо обоснованным. Нужно не сглаживать данные, если для этого нет достаточно серьезных причин, а разрабатывать адекватную имеющейся информации математическую модель. В частности, экспериментами установлено, что сглаживание приводит к ухудшению разрешения фазовых методов. Это можно объяснить тем, что сглаживание фактически приводит к той или иной деформации данных, что искусственно меняет фазовые соотношения сигнала.
Отметим, что использование фазы сигнала при интерпретации данных не является совершенно новой и неизвестной процедурой. В технических приложениях имеется ряд примеров, где был успешно
применен фазовый метод анализа информации. Например, в физике это голография, в измерительной технике - измерительные системы, основанные на принципах когерентной оптики (Драница и др., 2001; Драница, Жеребцов, 2002). Эти, казалось бы, далекие друг от друга приложения объединяет общая концепция - они используют фазовую информацию. Особенно впечатляют успехи, которые достигли за последние 10-15 лет когерентные измерительные системы. За это время пройден путь от опытных образцов до серийных изделий, которые на порядок и более увеличили потребительские качества измерительных приборов и систем по сравнению с амплитудными методами измерений.
Предпринятая нами попытка использования фазовых методов для решения ряда задач линейного анализа данных показала, что это направление исследований, вероятно, является достаточно перспективным. Метод представляет ряд этапов преобразования информации, при этом каждый этап имеет собственную математическую модель, в той или иной мере опирающуюся на динамический подход. Начальным этапом этой цепочки преобразований является аппроксимация данных выходом некоторой абстрактной линейной системы.
Литература
Вапник В.Н., Глазкова Т.Г., Кощеев В.А., Михальский А.И., Червоненкис А.Я. Алгоритмы и
программы восстановления зависимостей. Под ред. В.Н. Вапника. М., Наука, 816 с., 1984. Губанов В.С. Обобщенный метод наименьших квадратов. СПб., Наука, 318 е., 1997. Драница А.Ю., Драница Ю.П. Применение линейной модели для количественной диагностики экономических систем. Матер. Междунар. науч.-техн. конф. "Современнные проблемы экономики, управления и юриспруденции-2007". Мурманск, МГТУ, НТЦ "Информрегистр" 0320700489, с.73-75, 2007б. Драница Ю.П. Моделирование одномерных динамических процессов с целью предварительной
обработки результатов. Вестник МГТУ, т.4, № 1, с.97-115, 2001. Драница Ю.П. Об одном методе решения некорректно поставленных задач. Вестник МГТУ, т.3, № 1, с.67-78, 2000.
Драница Ю.П. Принципы нейроподобного моделирования геофизических объектов и процессов.
Вестник МГТУ, т.5, № 2, с.241-252, 2002. Драница Ю.П., Драница А.Ю. Корректный метод оценки спектральной плотности на основе линейной модели. Матер. Междунар. науч.-техн. конф. "Наука и образование-2007". Мурманск, МГТУ. НТЦ "Информрегистр" 0320700491, с.116-120, 2007а. Драница Ю.П., Драница А.Ю. Некоторые аспекты интерпретации экспериментальных данных на
основе теории линейных динамических систем. Вестник МГТУ, т.12, № 1, е.125-131, 2009. Драница Ю.П., Драница А.Ю. Об одной постановке обратной задачи локации. Матер. Междунар. науч.-техн. конф. "Наука и образование-2008". Мурманск, МГТУ, НТЦ "Информрегистр" 0320800238, с.88-91, 2008.
Драница Ю.П., Жеребцов В.Д. Использование волоконно-оптических технологий в геофизике.
Геофизика, № 6, с.30-39, 2002. Драница Ю.П., Жеребцов В.Д. Об одном методе обработки фазовых измерений в геофизике. Доклад на VIМеждун. науч.-техн. конф. "Современные методы и средства океанологических исследований", ч. 2, с.205-208, 2001.
Драница Ю.П., Жеребцов В.Д., Слипченко В.А. Частотно-временной метод обработки фазовых
измерений в геофизике. Измерительная техника, № 6, с.21-29, 2001. Корн Г., Корн Т. Справочник по математике (для научных работников и инженеров). М., Наука, 832 с., 1977. Краус М., Вошни Э. Измерительные информационные системы. М., Мир, 310 с., 1975. Кулханек О. Введение в цифровую фильтрацию в геофизике. М., Недра, 200 с., 1981. Никитин А.А. Теоретические основы обработки геофизической информации. М., Недра, 342 е., 1986. Теребиж В.Ю. Введение в статистическую теорию обратных задач. М., Физматлит, 376 с., 2005. Тихонов А.Н., Гончарский А.В., Степанов В.В., Ягола А.Г. Численные метода решения некорректных
задач. М., Наука, 214 с., 1990. Троян В.Н., Соколов Ю.М. Методы аппроксимации геофизических данных на ЭВМ. Л., Изд-во ЛГУ, 302 с., 1989.
Хемминг Р.В. Цифровые фильтры. М., Недра, 222 с., 1987.
