Государственный комитет Российской Федерации

По высшему образованию

Красноярский Государственный Университет

В.А. Кочнев

АДАПТИВНЫЕ МЕТОДЫ

РЕШЕНИЯ ОБРАТНЫХ ЗАДАЧ ГЕОФИЗИКИ

Учебное пособие

Красноярск

РОССИЙСКАЯ ФЕДЕРАЦИЯ

ГОСУДАРСТВЕННЫЙ КОМИТЕТ РОССИЙСКОЙ ФЕДЕРАЦИИ ПО

ВЫСШЕМУ ОБРАЗОВАНИЮ

КРАСНОЯРСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

ВЫЧИСЛИТЕЛЬНЫЙ ЦЕНТР СО РАН (КРАСНОЯРСК)

В.А. Кочнев

АДАПТИВНЫЕ МЕТОДЫ

РЕШЕНИЯ ОБРАТНЫХ ЗАДАЧ ГЕОФИЗИКИ

Учебное пособие Красноярск УДК 550. ББК 26.2В К Рецензенты: С.В. Гольдин- чл.-корр. РАН, профессор, зав.кафедрой НГУ А.М. Федотов- профессор, зав.кафедрой КПИ Кочнев В.А.

К 758 Адаптивные методы решения обратных задач геофизики учеб. пособие :

Краснояр. гос. ун-т; Красноярск, 1993. 120 с.

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

Предназначена для студентов, спирантов и геофизиков.

Табл. 20 Ил. 11. Библиогр.: 57 назв.

Без объявл. ББК 26.2В К 6 Л 5(02) ©В.А. Кочнев, Красноярский Государственный Университет, ISBN 5-230-07988- Содержание Предисловие

Глава 1. Введение в метод

Глава 2. Оценки одного параметра

2.1. Адаптивный вариант оценки среднего и дисперсии по неравноточным измерениям......... 2.2 Условия сходимости оценок

2.3 Упражнения

Глава 3. Адаптивный метод решения систем линейных Уравнений

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

3. 2. Обоснование метода решения

3.3. О сходимости метода

3.4 Примеры, иллюстрирующие свойства метода

3.5 Заключение

3.6 Упражнения

Глава 4. Адаптивный метод решения обратных задач гравиметрии

4. 1. Введение

4.2. Выбор модели и метода решения прямой задачи

4.3 Особенности системы уравнений и обратной задачи

4.4 Обоснование метода решения обратной задачи

4.6 Проверка точности решения прямой задачи

4.7. Проверка метода решения обратных задач на простых модельных примерах

4.8. Решение обратной задачи с привлечением априорной информации

4. 9. Погрешность определения параметров модели при различных горизонтальных размерах блоков

4.10. Проверка метода на сложных задачах.

4.11. Результаты и свойства метода

4.12. Некоторые рекомендации по методике решения прямых и обратных задач гравиметрии с применением пакета ADG-2

4.12.1. Решение прямых задач

4.12.2. Решение обратных задач

4.12.3. Работа с избыточными и полными плотностями

4.12.4. Редуцирование кривых dg

4.13. Упражнения

Глава 5. Решение прямых и обратных двумерных задач магнитометрии

5.1. Введение

5.2. Выбор модели среды

5.3. Выбор метода решения прямой задачи магнитометрии

5.4. Разработка алгоритма решения обратной задачи магнитометрии

5.5. Пакет программ решения прямых и обратных задач магнитометрии (АDM-2 )................ 5.6. Исследование особенностей магнитных полей на моделях

5.6.1. Необходимость модельных исследований

5.6.2. Проверка правильности работы программы решения прямой задачи

5. 6. 3. Моделирование аномалий магнитных полей на простых объектах

5.6.4. Исследование на более сложных моделях

5.7. Исследование устойчивости решения обратных задач по модельным данным.................. 5.7.1. Предисловие

5.7.2. Исследование устойчивости обратных задач для однослойной модели среды............ 5.8. Решение обратных задач на реальных данных

5.9. Выводы

5.10. Рекомендации по моделированию и решению обратных задач с использованием пакета ADM-2

5.10.1. Рекомендации по моделированию магнитных полей

5.10.2. Рекомендации по решению обратных задач

5.11. Упражнения

5.11.1. Моделирование аномалий от простых объектов

5.11.2. Моделирование полей от сложных аномальных объектов

Глава 6. Обзор литературы

6.1. Адаптивные методы в теории управления

6.2. Методы решения обратных задач в геофизике

6.3. Адаптивные методы в сейсморазведке

ЗАКЛЮЧЕНИЕ

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

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

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

Одним из простейших вариантов решения обратной задачи является подбор такой модели, которая дала бы теоретическое поле, совпадающее или близкое к наблюдаемому.

Сопоставляя прямую и обратную задачи, необходимо отметить следующие их особенности.

Прямая задача, как правило, имеет единственное решение. Заданной модели при заданных условиях наблюдения соответствует единственное поле.

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

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

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

Для преодоления неединственности и неустойчивости используются методы регуляризации [41]. Основным способом регуляризации является поиск решения вблизи априорно заданного начального приближения параметров модели.

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

На основе адаптивного подхода созданы методы решения различных обратных задач сейсмометрии, гравиметрии и магнитометрии, которые положены в основу многих пакетов программ.

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

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

Для первого знакомства достаточно прочесть главу 1- «Введение в метод».

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

В главах 4-5 показаны возможности адаптивных методов при решении задач гравиметрии и магнитометрии. В процессе изучения этих глав предполагается использование пакетов программ на персональных ЭВМ.

Обзор литературных источников и список литературы приведены в шестой главе.

Со стороны автора глубокая признательность всем, кто способствовал подготовке и улучшению пособия.

Глава 1. Введение в метод Вы приступили к изучению адаптивных методов и Вам хочется побыстрее узнать, что такое адаптивные методы и чем они отличаются от неадаптивных?

Прежде чем ответить на этот вопрос, начнем с простого примера.

Представим себе, что имеется некоторый постоянный, но неизвестный нам параметр x*, который можно определить или измерить с погрешностью. результат i-го измерения обозначим через xi. модель процесса измерения примем Будем считать, что i – является случайно нормально распределенной помехой [17, 22] с математическим ожиданием (М), равным 0 и с дисперсией (D=2).

При этих условиях оптимальной оценкой (имеющей наименьшую дисперсию) будет среднее, т.е.

Эта хорошо известная оценка среднего найдена неадаптивным методом.

Представим себе, что мы провели еще одно k+1 измерение и решили узнать среднее из k+1 измерений. Мы можем воспользоваться формулой (1.2), высчитав снова сумму и поделив ее на k+1, но можем уточнить среднее, построив следующую формулу Убедитесь, что оценка среднего, полученная по формуле (1.3), точно совпадает с оценкой Формула (1.3) имеет следующие преимущества перед формулой (1.4):

1. значение среднего вычисляется после каждого измерения, что очень важно в системах реального времени.

2. не может возникнуть переполнения при больших k и x.

Оба эти преимущества имеют решающее значение при включении их в системы реального времени, в которых необходимо знание усредненных оценок какого-либо параметра в любой момент времени. Совместно с оценкой среднего можно оценивать и дисперсию, по следующей формуле где x k +1 = x k +1 x k, т.е. разность, стоящая в скобках формулы (1.3).

Далее будем называть ее невязкой, а множитель, стоящий перед скобкой, корректирующим коэффициентом. В формуле (1.5) корректирующий коэффициент k +1 = отличается от такового в формуле (1.3), где k +1 = и это связано с тем, что если для оценки параметров x* можно иметь одно измерение, приняв x1 = x1, то для оценки дисперсии необходимо, как минимум, иметь два значения.

Квадрат разности между ними и будет первой оценкой дисперсии. При больших k корректирующие коэффициенты в формулах (1.3) и (1.5) будут мало различимы, их можно брать равными.

Алгоритмы (1.3) и (1.5) являются адаптивными, поскольку в них оценка параметров x и D ведется через уточнение предыдущих значений по невязке. Знание таких параметров, как мы увидим дальше, очень важно для построения оптимальных более сложных алгоритмов. Два простых примера оценки параметров показывают нам идею последовательного пополнения недостающих априорных сведений о модели в ходе поступления текущей информации. Накопление и немедленное использование текущей информации для устранения неопределенности из-за недостаточной априорной информации с целью оптимизации избранного показателя качества и является наиболее характерной чертой адаптации [48].

Задачу оценки одного параметра более подробно обсудим во второй главе.

Там же покажем условия сходимости оценок.

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

При однократном измерении погрешность получаемой оценки может определяться (оцениваться) точностью прибора или другими соображениями, связанными с условиями измерения или эксперимента. Продолжая тему оценки одного параметра, будем считать, что измерения являются неравноточными и что для каждого из них нам известна оценка погрешности i. Тогда оптимальной оценкой (имеющей наименьшую дисперсию) будет Выбор веса обратно пропорционального дисперсии измерения дает оценку x с минимумом дисперсии [22]. Адаптивный вариант оценки, будет иметь вид где Pk = Pi -вес среднего значения на k-ом шаге.

Перейдя от весов к 2 и используя выражения Получим Вес k+1 оценки x будет равен сумме весов P k +1 = P k + Pk +1. Записав последнее выражение через дисперсию, получим Значение k2+1 иногда может служить оценкой погрешности среднего, если k2 точны. Реально же такое бывает редко и поэтому возникает необходимость в поиске независимой оценки дисперсии, получаемой численно через невязки. Для оценки дисперсии соответственно получим где - = если по аналогии с оценкой x ввести первый шаг, в котором задается априорное значение дисперсии D0 и погрешности дисперсии DD. В частном случае можем принять D0 = DD = 3 14 [17].

Нетрудно убедиться, что оценки x, получаемые по формулам 2.1 и 2.3, дают одинаковые значения при равноточных измерениях, т. е. при одинаковом i придем к формулам 1..2 и 1.3.

Рассмотрим еще одно обобщение, предполагая, что ведется измерение значения функции U i, которое связано с x * следующим образом:

Тогда для x k +1 будем иметь k - оценка дисперсии среднего на k-ом шаге.

Произведение в ( 2.6 ) a x k будем называть прогнозным значением U i, а разность - невязкой между фактическим и прогнозным значениями ( U ). Т.е. невязка на k+1 шаге: U = U k +1 a x k Для оценки дисперсии получим что следует из свойства дисперсии DU=a2D. В частном случае при а=1 выражения 2.6 и 2.7 эквивалентны 2.3 и 2. 4.

На примере оценки одного параметра видны некоторые особенности адаптивных методов и главные из них:

1.Неизвестный параметр уточняется по невязке между очередным измерением и прогнозным значением измеряемого параметра 2.Переход от невязки к параметру происходит путем умножения невязки на корректирующий коэффициент k +1, зависящий в свою очередь от k2+1 и k и коэффициента а (2.6).

3.Независимо от оценки k +1 осуществляется численная оценка дисперсии Dk+1.

Рассмотрим условия сходимости к точной оценке х приведенных выше алгоритмов.

2.2 Условия сходимости оценок А. Дворецкий предложил для исследования сходимости рассматривать любую процедуру стохастической аппроксимации как обычный детерминированный метод, но с наложением на него случайной составляющей [42]. Следуя этому подходу, формулу 2. запишем в виде где k - ошибки среднего значения на k-ом шаге;

k +1 - ошибка к+1 измерения.

В этом выражении случайная составляющая сосредоточена в последних двух слагаемых. Обозначим каждое из слагаемых через ri. На этом примере проиллюстрируем условия сходимости, доказанные А. Дворецким. В [42] показано, что для сходимости процедур стохастической аппроксимации в среднеквадратическом и с вероятностью единица, требуется выполнение двух условий;

- помеховая составляющая ( ri ) должна быть несмещенной:

- сумма дисперсий случайной составляющей должна быть конечной при любой бесконечной процедуре поиска:

Для выполнения первого условия необходимо, чтобы второе и третье слагаемые выражения (2.9) стремились к 0, если k +1 - конечно и k +1 стремится к 0. Второе слагаемое стремится к 0, если выполняется условие M [ i ] = 0. При его невыполнении оценка, среднего окажется смещенной.

Для анализа второго условия, следуя Дж. Уайлду [42] будем считать, что. суммарr. Если помехи образуют, последоная ошибка, внесенная всеми помехами равна i вательность независимых случайных величин с дисперсией 2, то дисперсия суммарной ошибки равна математическому ожиданию суммы квадратов отдельных ошибок Учитывая сказанное, условие (2. 11) для (2. 9) будет Если при k и k 0, то второе слагаемое конечно и тоже стремится к0.

Для анализа первого слагаемого при 0 л +1 1 запишем.

Правая часть неравенства, будет конечной лишь в том случае, если i2 <.

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

шаге не производится. Получите результат при ураганной помехе и при 5. Напишите программу и исследуйте процедуру для оценки амплитуды и дисперсии реальной части спектра для разных w, для выборки, формирующейся по следующей модели Используя формулы 2. 1 и 2. 7, рассчитайте значения x k и Dk для w = 0.1w0 ; 0.2w0 ;....;0.9w0 ; w0 ;1.1w0 ;....;2w0. Постройте и проинтерпретируйте графики x k = f ( w) и Dk = f ( w) 6. То же, что и в 5-ом пункте, сделайте для функции 7. То же, что и в 5-ом сделайте, для выделения реальной и мнимой части проинтерпретируйте графики x( w) и D x ( w) - для реальной и y ( w) и D y ( w) те и объясните отличие графиков x( w) и D( w), полученных в 7 и 8 пунктах; Выберите правильный алгоритм.

9. То же, что и в 7-ом пункте, посчитайте для функции xi = a cos w + t i.

Постройте графики x( w) и D x ( w), y ( w) и D y ( w) и сделайте анализ.

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

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

3.1 Постановка задачи Предположим, имеется система линейных уравнений где. [aij ]- матрица коэффициентов системы.

[U i ] - вектор- столбец, являющийся результатом физических измерений или некоторых оценок.

Допустим, что известны:

1. Средняя квадратическая погрешность каждого i-того измерения, дающая вектор Ui 2. Вектор начальных приближений неизвестных - x 0 j .

3, Вектор погрешностей неизвестных x j По поводу этих допущений необходимо заметить, что при решении конкретных задач, связанных с оценкой параметров модели, указанные сведения обычно имеются. Мерой погрешности правых частей в первом приближении может быть точность прибора. Например, при определении времени прихода отраженной волны погрешность сильно варьирует и определяется как функция от соотношения сигнал-помеха на конкретном участке трассы [13]. При таком подходе каждое уравнение системы (1) становится неравноточным [14], что более адекватно физической сути задачи, чем подход, который вытекает из решения равноточных алгебраических уравнений. Следовательно, на поиск путей оценки значений вектору U при предлагаемой постановке задачи должно обращаться серьезное внимание. Приступая к оценке параметров физической модели, исследователь обычно знает о приближенных значениях параметров и может задавать начальное приближение вектора неизвестных. На минимизации некоторого функционала с учетом этого вектора и построен принцип регуляризации по Н. А. Тихонову [41]. Из физической сути задачи следует, что степень достоверности начальных значений параметров известна не одинаково точно. И это очень важно учесть при решении, для характеристики достоверности и введен вектор погрешностей неизвестных x j. В качестве первого приближения значений x j могут быть взяты возможные пределы изменения неизвестного.

В самой неопределенной ситуации значения вектора неизвестных могут быть приняты равными 0, а значения x j одинаковые, равные какому-либо большему числу, которое будет характеризовать возможную наибольшую разницу между истинными x *j, и начальными значениями.

Определив условия задачи, перейдем к конструированию алгоритмической части процедуры решения систем линейных уравнений.

3. 2. Обоснование метода решения В соответствии с адаптивным подходом (14) уточнение неизвестных ведется по каждому алгебраическому уравнению раздельно. Очередность поступления уравнений может быть нормальной или определяться величиной невязки, как в методе релаксаций [24]. Запишем i--тое уравнение системы (1) Возьмем в качестве прогнозных значений неизвестных x j начальные значения x0 j, т.е. x j = x0 j. Подставим прогнозные значения в левую часть уравнения, получим некоторое значение U i, которое назовем прогнозным.

Найдем разность между фактическим значением (правая часть) и прогнозным и будем называть её невязкой Величина невязки в общем случае обусловлена следующими причинами:

1. Отклонениями прогнозных значений x j от истинных x *j 2. Ошибками измерения 3. Неточностями коэффициентов [aij ] или, другими словами, неадекватностью модели физической, из которой получены измерения, модели математической, из которой получены прогнозные значения. Неадекватность обычно выявляется и устраняется на этапе отладки алгоритма. Дальнейший разговор будем вести предполагая, что система уравнений адекватна и что доминирующими в создании невязок являются первые два фактора.

Следуя [27], представим невязку в виде что в линейной системе дает Обозначив слагаемые через z j, получим Задача уточнения n неизвестных по одной невязке, в общем случае, имеет множество решений. Все слагаемые являются случайными величинами, распределенными, по нормальному закону, что следует из введенных ранее предположений и свойств случайных величин [17]. Найдем совместную плотность вероятности распределения параметров (от которых зависит U i.) в n+1 - мерном пространстве.

Значения z j выберем таким образом, чтобы плотность вероятности была максимальной (метод максимума апостериорной плотности вероятности) [17].

откуда, исключив первую сумму констант и поменяв знаки, получим новый вид целевой функции Если принять то локальная целевая функция будет иметь вид что по форме соответствует регуляризующему функционалу А.Н. Тихонова [41 с. 117], т.е. в анализируемом методе имеет место локальная регуляризация, и параметру регуляризации для каждого неизвестного и в каждом уравнении определяются величинами z2 и U. Продифференцировав по каждой из неизвестных, получим систему из n уравнений с n неизвестными. Она имеет следующий вид:

Решая систему, получим Таким образом, невязка разбрасывается пропорционально дисперсии неизвестного.

В знаменатель входит супца дисперсий неизвестных и дисперсии правой части. При U >> z2 неизвестные практически уточняться не будут. При обратном неравенi j стве невязка будет полностью распределена на уточнение параметров.

Получив оценки zj, перейдем непосредственно к оценкам неизвестных х, получим:

А это известный математикам метод Качмажа [2 1]. Таким образом создаваемый адаптивный метод является обобщением известного итерационного метода.

Нам удалось построить процедуру уточнения неизвестных по одному уравнению и она схожа с процедурой оценки одного параметра ( формулы 2.2 и 2.6 ). Следовательно, мы можем использовать условия сходимости, полученные ранее. Главное из ник k должно стремиться к 0 при k, стремящемся к бесконечности. Дисперсии уточняемых параметров будем находить по формуле Анализируя видим, что дисперсия может увеличиваться (при положительной разности т. к. произведение а на всегда положительно). Следовательно, она не может использоваться как вычисления корректируемых коэффициентов. Для этих целей будем вычислять оценку 2, через сумму весов, т.е.

Обобщим приведенный в главе 2 подход для изменения параметров, измеренных не непосредственно, а оцениваемых через функциональные связи. В этом случай оценка x ( k +1) зависит от точности измерения, от слагаемых и коэффициентов, входящих в уравнение. Допустим, измеряется U :

Тогда Полагая, что a и b известны точно, для оценки ошибки x запишем Дисперсия, оценки х будет равна математическому ожиданию квадрата выражения правой части уравнения (3.12) Выражение окажется близким к 0, если U b 2 y и b > 0, и будет большим, если U b 2 y и b< 0. Во избежание этих крайних ситуаций третье слагаемое в скобках отбросим, т.е. будем постулировать независимость ошибки измерения и ошибки случайной составляющей второго слагаемого. С учетом сказанного получим Преобразовав формулу, получим Аля многомерного случая формула для j-го неизвестного в i-том уравнении имеет вид Или, учитывая (2.11), найдем Здесь вывод формулы (3.17) сделан методом индукции из общих подходов метода наименьших квадратов.

В модифицированном варианте формула имеет вид при о < < 1.

Для нелинейного случая следует иметь в виду, что aij = при x *j (k ) для всех j от 1 до n.

Таким образом, в соответствии с поставленной задачей мы обосновали алгоритм итерационного уточнения параметров в одном уравнении. Он включает формулы: 3, 4, 11 и 21.

При переходе к следующему уравнению ( от k+1 к k+2 ) в качестве априорных будут использоваться уточненные значения параметров и их дисперсии. Так будет продолжаться до тех пор, пока не будут исчерпаны все уравнения системы.

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

где - малое наперед заданное число, a D ( L ) -дисперсия невязок, вычисляемая, по формуле или по более удобным рекуррентным формулам (2.2 или 2.4 ). В предельном случае, когда X X *, из (3. 24) следует т. е, D(l ) будет характеризовать средневесовую дисперсию помехи, которая от номера итерации не зависит. Следовательно, будет выполнено условие (3.23).

нова может быть принято неравенство:

Аля наглядности представим алгоритм адаптивного метода решения систем линейных алгебраических уравнений в виде блок-схемы (рис. 1).

В центре схемы видим блок памяти (БП) для хранения матрицы коэффициентов, значений правых частей, неизвестных и их дисперсий и других вспомогательных переменных. Он соединен с вычислительными блоками двойными линиями, указывающими передачу информации. В блохах А0, А1, A2 устанавливаются начальные значения констант и организуются циклы по итерациям и уравнениям. Блоки A3 – А6 реализуют вычислительные функции, указанные в них, а блоки А7-А8 переключают алгоритм на очередные уравнение или итерацию и на окончание.

3.3. О сходимости метода Созданный итерационный метод по классификации Ф.А. Фаддеева и В. Н.

Фаддеевой (Вычислительные методу линейной алгебры-М. : Физматгиз 1963) является не стационарным (циклическим), т. к. матрица Н(L) в уравнении существенно зависит от номера итерации (L). Поэтому доказательство общей сходимости построить не удается, да это и не имеет особого смысла, т. к. в реальных задачах при огромных матрицам исследовать сходимость метода для конкретной системы - задача более трудоемкая, чем ее решение и анализ на сходимость по уменьшению невязки. Поэтому в дальнейшем дли анализе сходимости будем использовать подходы, разработанные в методах стохастической аппроксимации [42]. В методах стохастической аппроксимации число независимых наблюдений предполагается неограниченным и за счет этого, при определенных ограничениях на свойства помех, возможно получение оценок со сколь угодно высокой точностью. В наших условиях мы вынуждены компенсировать недостаток независимых наблюдений многократным использованием одной и той же выборки с тем, чтобы максимально ослабить влияние априорной информации.

Условия сходимости с двумя неизвестными исследовались в работе [14]. Приведем заключительную формулу критерия сходимости на очередном k-том шаге На k-том шаге алгоритм приблизится к решению, если будет выполнено приведенное неравенство.

В частном случае при = 1;W =, т, е. имеет место сходимость независимо от соотношения других параметров. Этот частный случай соответствует методу Качмажа [21].

3.4 Примеры, иллюстрирующие свойства метода Приведем простые задачи, допускающие проверку и наглядное представление результата.

Пример 1. От начального приближения x0=0,5; y0=0.3, используя различные итерационные методы, найдем решение следующей системы Точное решение системы x * = 0; y * = 1.

Отклонение x = x0 x * = 0.5; y = y 0 y * = 2, т.е. из (7.1) o=4. Рассмотрим сначала наиболее благоприятные условия, положив априорные средние квадратические ошибки (АСКО) равными Используя адаптивный метод за 4 итерации приближаемся к точному решению (табл. 3.1 и рис. 3.2) 3.2 такое же приближение, удается получить за пять итераций.

И, наконец, предположим, что априорная информация достоверности неверна, приняв x = 2; y = 0.5 = 0.25; = 16. На рис. 3.2 видим, что для получения того же результата необходимо на две итерации больше, чем в первом варианте и на одну по сравнению со вторым.

Для сопоставления на рис. 3.2 приведены траектории результатов, полученных с помощью методов покоординатного спуска и Качмажа. В данном примере большая скорость сходимости оказалась у адаптивного метода (при первом и втором варианте АСКО), а самая медленная - у метода Качмажа.

Пример 2. Начиная от точки x0=y0=2, найдем решение численными итерационными методами системы Следует заметить, что решить эту систему методами покоординатного спуска без поворота системы координат нельзя. Они Не дают решения, т.к. не сходятся. Метод Качмажа (рис. 3) сходится на каждом шаге, но по мере приближения к решению скорость сходимости постепенно уменьшается. Что касается адаптивного метода (рис. 3.3 и табл. 3.2), то его траектория становится более «целеустремленной» по мере приближения к решению. Это объясняется тем, что из-за больших коэффициентов перед y, 2y сильно уменьшается и невязка в большей степени идет на уточнение x.

рации уравнения Пример 3. Начиная от точки x0 = 0.5; y 0 = 3 с использованием адаптивного метода найдем решение несовместной системы трех уравнений 0.333 x + y = 2. Имеющих различные погрешности правых частей.

U = 0. На рис.3.4 и табл.3.3 видно, что траектория решения приближается к точке пересечения 1 и2 линий, а уравнения 3 оказывает влияние только при первом (с его участием) уточнении (отрезок от 2 к 3 шагу), пока x и y являются большими.

Как видно из рис.3.5 и табл. 3.4 решение постепенно приближается к центру треугольника.

Пример 5. Используя адаптивный метод, решим нелинейную систему при следующих начальных данных:

t 0 = 0.95;

Точное решение t = 1, v = Результаты уточнения для = 0 приведены в таблице 3. Номер ите- Номер уравvk Пример 6. (Взят из [42, c.263]). От начального приближения x=10 найти корень для функции y = x 3 + x + k, где k - случайное число, которое предлагается выбирать ± c, бросая монету. Не меняя сути задания, в качестве помехи используем гармоническую функцию k = c sin(0.5k ).

На этом примере сопоставим три метода: Ньютона, стохастической аппроксимаций и адаптивной. При отсутствии помех, т.е. при с=0, метод Ньютона и адаптивный при решении этой задачи практически совпадают, отличие результатов имеет место при наличии помех, что видим из таблицы 3. 6, в которой сведены результаты поиска корня при с=2. Адаптивный метод испытывался при двух вариантах задания АСКО начального приближения:

дисперсии x20 = 10 и x20 = 100 ) первоначально одинаково быстро приближаются к решению, но, начиная с 7-го шага, видны отличия, которые отчетливо проявляются на 11 и 15 шагах. В результатах адаптивного метода флюктуации оценок корня, вызванные помехой, постепенно уменьшаются из-за постепенного уменьшения x2, а в результатах метода Ньютона они остаются на одинаковом уровне ( ± 2 ).

Сходимость метода стохастической аппроксимации [42] сильно зависит от выбора шага. Если в этом примере взять k =, то процесс уточнения будет быстро расходиться. Уже на первом шаге будет получено x1 500, а на Результаты счета сведем в таблицу 3. Из таблицы видно, что результаты метода, стохастической аппроксимации практически сильно зависят от первого или второго шага, а дальше они Приведенные выше примеры показывают нам, что адаптивный метод позволяет:

а) ускорить сходимость решения задачи не только за счет учета достоверности начальных приближений, но и за счет собственных адаптивных свойств алгоритма (примеры l, 2) б) решать несовместимые неравноточные системы (примеры 3 и 4), т. е, системы, не имеющие решения;

в) решить нелинейные системы (примеры 5 и 6).

3.5 Заключение В начале главы, мы предположили, что кроме системы уравнений нам известны:

1.Средняя квадратическая погрешность каждого i-того измерения, дающая вектор [ Ui ] 2.Вектор начальных приближений неизвестных - [x0 j ] 3.Вектор погрешностей неизвестных [ xj ] В такой постановке мы обосновали новый метод итерационного уточнения параметров в одном уравнении. Он включает формулы 3, 4, 11, 21.

При переходе к следующему уравнению ( от k+1 к k+2 :) В качестве априорных будут использоваться уточненные значение параметров и их дисперсии. Так будет продолжаться до тех пор, пока не будут исчерпаны все уравнения системы.

В связи с тем, что при выводе формул постулировалась независимость уточняемых компонент, возникает необходимость повторить уточнение по всей системе до тех пор, пока оценки не будут близки к решению, что будет определено выполнением условия (3.23).

Анализируя, мы увидели следующие свойства метода:

1. Метод относится к классу итерационных.

2. Он обобщает метод Качмажа.

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

Для дальнейшего изучения возможностей метода полезны приведенные ниже упражнения.

3.6 Упражнения 1. Выберите систему двух уравнений с двумя неизвестными. Задайте начальные приближения неизвестных и x и y. Решите систему адаптивным методом, построив таблицу решения (см. пример 1).

На координатную плоскость XY нанесите уравнения и траекторию решения методов: адаптивного, Качмажа и покоординатного спуска (см. рис 3.2). Сопоставьте скорость сходимости методов. Смените значения x и y и повторите решение, дайте интерпретацию полученных результатов.

2.Выберите систему нелинейных уравнений с двумя неизвестными. Решите адаптивным методом и методом Ньютона. Постройте таблицы и графики, дайте интерпретацию.

3.То же что и в третьем примере, но в одно из уравнений введите помеху cos k, где k = 0,1,2,...,20. Соответственно примите U = 1.

4.Используя программу решения систем линейных уравнений (ADASLU) на персональных ЭВМ, решите выбранную вами систему линейных уравнений при различных параметрах x и y, распечатайте невязки и результаты. Проинтерпретируйте результаты.

5.То же, что и 4, но введите в уравнение помехи: anois = 0.1;1 при U = 0,1, U = 1.

Выпишите СКН и дайте объяснение уровню невязок (см. формулу 3. 25) 6. Выберите и решите недоопределенную систему с различными начальными данными, опишите результат.

7. Выберите и решите переопределенную систему при U = 0.01. Сделайте описание результата.

Глава 4. Адаптивный метод решения обратных задач гравиметрии 4. 1. Введение Обратные задачи гравиметрии (особенно трехмерные) относятся к классу наиболее трудных [9] в силу нескольких причин.

1. Отсутствие единственности решений. Одному, и тому же гравитационному полю, в общем случае, может быть подобран не один вариант модели среды.

2. Неустойчивость решения, даже тогда, когда имеет место единственное решение обратные задачи обладают свойством неустойчивости, т. е. небольшим изменениям гравиметрического поля могут соответствовать большие изменения в распределении масс.

3. Трудности решения алгебраических систем с большим числом уравнений и неизвестных (103 и более).

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

Эта задача решалась многими исследователями [2,4,19,34,36-40], но наиболее близкими к нашему направлению являются методы сеток, предложенные Г. Г.

Булахом и А. А.Юньклвым [4]. В дальнейшем для решения систем и учета ограничений на плотность С. В. Шалаевым было предложено использовать методы линейного программирования [49]. Позднее этот подход был продолжен в работах В. В. Ломтадзе [19]. Использование методов линейного программирования не всегда оказывается конструктивным. Его неконструктивность выявляется при необходимости ввода жестких ограничений. Например, при отрицательных аномалиях необходимо решать "зеркальную" задачу, т. е. переводить и аномалии, и избыточные плотности из отрицательных в положительные. А как быть, когда есть и отрицательные, и положительные? Кроме того, трудности в методе возникают и тогда, когда число неизвестных велико.

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

Используем для решения обратных задач гравиметрии изучаемый нами адаптивный метод. Впервые он был применен для решения кинематических задач сейсморазведки [11-16]. Накопленный опыт позволил создать методику разработки и исследования алгоритма конкретной обратной задачи [13]. В соответствии с методикой необходимо пройти последовательно следующие этапы:

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

Намеченный порядок разработки и используем для последующего описания.

4.2. Выбор модели и метода решения прямой задачи За основу примем слоистую модель среды, которая является достаточно общей при изучении осадочных толщ. Она же может использоваться и при моделировании большого числа интрузивных тел. Каждый слой по оси X разбивается на прямоугольные параллелепипеды, вертикальные размеры которых определяются мощностью слоя, а горизонтальные - заданными величинами dx и dy, соответствующими шагу прямоугольной сетки наблюдения или сканирования гравитационного поля. Точки наблюдения будем располагать на поверхности над центрами параллелепипедов.

Проследим логику вывода формулы расчета гравитационного поля для прямоугольного параллелепипеда. В соответствии с законом Ньютона имеем где m1 и m2 -точечные массы тел, r 2 - расстояние между ними, G - гравитационная постоянная, равная 6.67 10 Если F поделить на m1, то получим вместо силы напряженность гравитационного поля, которую обозначим через g. Введя вместо m2 – переменную M получим Напряженность - векторная, величина, имеющая то же направление что и сила F. Следуя [23], введем радиус вектор R, соединяющий точку наблюдения и источник. Направление по радиусу к источнику считаем положительным, тогда вектор g Для расчета напряженности поля объемного тела, с плотностью вещества, рассмотрим элементарный объем dv, а затей интегрируя по всему телу, получим Рассмотрим вертикальную составляющую гравитационного поля в декартовой системе. Формулу для нее получим из (4. 4) где z и z0- координаты точки источника и наблюдения.

Преобразуем формулу (4.5), учитывая, что в гравиметрии интерпретируется не все поле, а его аномальная составляющая g, вызванная аномальной или избыточной плотностью.

Предположим, что точка наблюдения расположена в начале координат, источник напряженности представлен в виде параллелепипеда, ограниченного сверху горизонтальной плоскостью z1, а снизу z 2, по оси X плоскостями x1 и x 2, а по оси Y соответственно y1 и y 2. Переходя от интегрирования по объему к интегрированию по осям декартовой системы координат, получим вклад в поле g от элементарного параллелепипеда Проигнорировав, получим [2, 9, 36, 39] Для двухмерного случая при y будем иметь Вводя множество параллелепипедов (которые будем называть, для краткости, блоками или призмами) и, обозначая избыточную плотность в каждом блоке через ijk, а гравиметрический эффект от блока с = 1 через cijk, получим следующее уравнение для расчета g от блока в точке наблюдения с координатами х0,у0,z Для каждой точки наблюдения получим своё уравнение.

4.3 Особенности системы уравнений и обратной задачи Проиллюстрируем особенности систему уравнений и обратной задачи на простом двухмерном варианте с одним слоем. В этом случае, g в i--той точке наблюдения от J элементарных блоков будет Взяв совокупности из I наблюдений, получим систему из I уравнений с J неизвестными При I= J система будет иметь единственное, но не всегда устойчивое решение, что будет проиллюстрировано далее. В частном случае, когда поверхность наблюдения и границы слоя будут горизонтальны, то матрица коэффициентов обладает следующими свойствами:

1. диагональные элементу матрицы равны и симметричны относительно главной диагонали.

2. Элементы главной диагонали будут наибольшими.

3. Значения остальных элементов будут убывать по мере удаления от главной диагонали.

Все эти свойства четко объясняются из физических предпосылок. На главной диагонали находятся коэффициенты cii равные g ii при избыточной плотности, равной единице. Они отражают гравиметрический эффект от блоков, расположенных под точкой наблюдения. Вспомним, что точки наблюдения находятся на поверхности над центрами блоков (4.2). По мере удаления блоков от точки наблюдения их влияние на гравиметрическое поле будет уменьшаться.

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

Приведем для иллюстрации коэффициенты для трех слоев, кровли которых располагаются на глубинах 0, 500 и 1000 м, при горизонтальной величине блоков 500 м.

N слоя интервал В случае, когда поверхность наблюдения или границы слоев не горизонтальны, упомянутые свойства коэффициентов не сохраняются. Поэтому необходимо вычислять их значения для всех точек наблюдения, что приводит к увеличению числа коэффициентов в I раз (где I число точек). Для двухмерной задачи, включающей сто точек по профилю и пять слоев, потребуется хранить 100 100 5 = 50000 коэффициентов. Для трехмерной задачи их число увеличивается как минимум на порядок.

4.4 Обоснование метода решения обратной задачи Предположим, что априорно известны начальные приближения избыточных плотностей j в каждом блоке. В частном случае они могут, быть заданы константой, равной нулю.

Будем считать, что где *j - точное, но неизвестное нам значение избыточной плотности, а (j0 ) - погрешность начального приближения.

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

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

Предположим также, что правые части, т. е. значения g i, осложнены помехой и что нам известен стандарт ошибок правой части i, в частном случае он может быть задан константой 0, равной погрешности съемки.

В предлагаемом методе [11-16] система уравнений не решается в целом, а осуществляется многократное, по всей системе, последовательное уточнение неизвестных по каждому уравнению отдельно, как это осуществляется в некоторых пошаговых итерационных алгоритмах [24, 42].

Предположим, что неизвестные уточняются на k+l-ом шаге l-ой итерации, по i-ому уравнению, которые связаны следующим соотношением: k = i + I (l 1).

Подставим в него сформировавшиеся на k-ом шаге значения неизвестных (k ) :

Значение g i будем называть прогнозным значением гравитационного поля в точке i.

Найдем разность между фактическим и прогнозным значениями гравитационного поля Наличие ненулевой невязки в общем случае обусловлено следующими причинами:

• отклонением прогнозных значений избыточных плотностей (k ) от истинных, т. е. j являются не нулевыми;

• неадекватностью математической и физической модели.

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

Итак, предполагая, что основная причина невязки обусловлена первой и второй причинами, представим невязку в виде суммы Детально повторив рассуждения, приведенные в главе 3, придем к следующим оценкам значений j где k- номер шага, а i- номер уравнения.

Из приведенной формулы видим, что в большей степени уточняются избыточные плотности тех блоков, у которых большие коэффициенты и большие дисперсии погрешностей, отсюда уточненные избыточные плотности в i-том уравнении на k-том шаге будут равны:

В адаптивном методе на каждом шаге изменяется (уменьшается) значение параметра 2 в соответствии с обоснованием, приведенным в главе примем следующую формулу, пересчета дисперсий Из выражений (17-20) видим, что значению параметра, уточненному в большей степени, будет соответствовать и большее уменьшение 2.

После уточнения оценок неизвестных и изменения дисперсий переходим на очередное уравнение и процесс уточнения неизвестных повторяем, используя формулы (15, 16, 18-20). Продолжая до последнего уравнения, пройдем первую итерацию уточнения. Из-за того, что процесс уточнения идет приближенно, необходимо повторное уточнение на нескольких итерациях до тем пор, пока не будет выполнено одно из двух условий: k = k max или