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

6

Размещено на http://www.allbest.ru/

6

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

Имомназаров Х.Х.,

Холмуродов А.Э.

Аннотация

ламе сейсмический рекуррентный пористый

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

ВВЕДЕНИЕ

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

На практике, как правило, характеристики среды являются неизвестными. Поэтому в приложениях наибольший интерес представляют обратные задачи об определении переменных коэффициентов дифференциального уравнения, связаные с физическими характеристиками среды, в которой наблюдаются сейсмические явления. Например, система уравнений теории упругости для изотропной среды содержит три параметра: плотность вещества и параметры Ламе, характеризующие упругие свойства вещества. Очень часто эти параметры являются функциями координат. Так как непосредственно их измерить зачастую невозможно, то задача об определении свойств вещества является обратной. При этом часто в качестве информации для ее решения задаются характеристики физического процесса, измеряемые на границе области. Наиболее ярко это проявляется в задачах геофизики, основной задачей которой является изучение внутреннего строения Земли на основе наблюдений, осуществляемых на ее поверхности или в приповерхностном слое. Именно потому, что непосредственные измерения физических параметров осуществимы только в очень узком приповерхностном слое, единственным путем изучения внутреннего строения Земли является метод, основанный на рассмотрении обратных задач. А.С. Алексеевым в 1962 году впервые был рассмотрен ряд постановок обратных задач теории распространения волн для модели упругих сред [1]. Более полную информацию по теории обратных задач можно найти в [2], [3], [4].

Основой данной работы послужили статьи в журналах [7], [8], [9], в которых рассматриваются алгоритмы решения одномерных прямой и обратной задач распространения сейсмических волн в упругой неоднородной среде. Но, на самом деле, Земля не является упругой средой, и поэтому было решено рассмотреть эти же задачи для случая распространения сейсмических волн в пористых средах, которые дают более точное приближение математической модели к реальности. В данной работе представлен численный метод интегрирования вдоль характеристик, который позволяет эффективно вычислять коэффициент Ламе по заданному сейсмическому сигналу, измеренному на поверхности с произвольной исходной возмущающей функцией. Метод характеристического интегрирования заключается в преобразовании дифференциального уравнения в каноническую форму. Интегрирование вдоль соответствующих характеристик и приближенное вычисление интегралов методом трапеций приводит к явной разностной схеме. Это позволяет восстановить профиль коэффициента Ламе внутри среды, одновременно восстанавливая и волновые поля смещений упругого пористого тела и жидкости.

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

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

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

1. ПОСТАНОВКА ПРЯМОЙ ЗАДАЧИ

Одномерное уравнение распространения сейсмических волн в пористой среде с учетом поглощения энергии, обусловленного коэффициентом межкомпонентного трения , имеет вид [5], [6]

(1)

. (2)

где и - компоненты вектора смещения частиц упругого пористого тела и жидкости с парциальными плотностями и соответственно, - коэффициент Ламе.

Пусть система уравнений (1), (2) справедлива при . Предположим, что пористая среда покоится при :

(3)

(4)

Пусть на границе приложена сила:

для (5)

где - скорость волны при .

Требуется по этой информации и заданным дважды непрерывно дифференцируемым функциям , , непрерывным функциям , определить волновые поля , из (1)-(4). Такую задачу будем называть прямой динамической задачей распространения сейсмических волн в пористой среде.

2. ПОСТАНОВКА ОБРАТНОЙ ЗАДАЧИ

Рассмотрим уравнения (1), (2) с начальными условиями (3), (4) и граничным условием (5). Пусть задано смещение упругого пористого тела на поверхности :

(6)

Обратная задача состоит в том, что по условию (6) требуется определить коэффициент из (1)-(5) (остальные функции , и заданы).

3. РЕШЕНИЕ ОБРАТНОЙ ЗАДАЧИ

Проведем преобразования сначала для уравнения (1.1). Будем использовать систему координат времени пробега с новой переменной

,(7)

Где

- скорость волны.

В новой системе координат волновое уравнение преобразуется к виду

(8)

Где

_ известная функция,

_ полное сопротивление среды.

Так как уравнение (8) не меняется при умножении на константу, нормируем так, чтобы . Заметим, что при введенной таким образом функции вычисление коэффициента Ламе становится эквивалентно вычислению .

Положим , и . С новыми независимыми переменными волновое уравнение второго порядка (8) становится гиперболической системой уравнений первого порядка

(9)

Где

Граничные условия преобразуются к виду

 .(10)

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

Определяя

имеем

Умножая систему (9) слева на и полагая, что _ дифференцируемая функция, получаем каноническую форму системы (9):

где . В покомпонентной записи система выглядит так

(11)

(12)

где _ производная .

Уравнения (11), (12) можно также записать в интегральной форме. Для произвольной фиксированной точки из области и любого можно проинтегрировать уравнение (11) вдоль характеристики от до . Тогда получим следующее соотношение

(13)

Аналогично проинтегрировав уравнение (12) вдоль характеристики от до получаем второе соотношение:

(14)

Полученные уравнения (2.7) и (2.8) уже подходят для численного обращения.

Предположим, что и из (5) и (6) соответственно известны для . Зная и на границе для и то, что на линии , можно восстановить , и в области , показанной на рис. 1.

Рис. 1 Область определенности интервала .

Положим

.

Пусть , но достаточно мало. Тогда можно приблизить под интегралом выражением

(15)

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

(16)

(17)

соответственно.

В случае если имеет разрыв в точке , где , но поле непрерывно на этой поверхности, уравнения (16) и (17) дают хорошую аппроксимацию при достаточно малом .

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

Рис. 2 Численная сетка, используемая в методе интегрирования вдоль характеристик. Жирные линии - характеристики гиперболической системы уравнений (9). Линия - фронт волны.

Сеточные функции для , , , и определяются следующими соотношениями:

, , ,

где нормирована так, что ,

,

,

для и . Из начальных условий получаем, что

, для

и, из граничных условий, имеем

,

.

Уравнения (13) и (14) вместе с аппроксимациями (16) и (17) подходят для обращения. Предположим, что известны для , а и вычислены для всех и . Тогда мы можем выразить и из следующих соотношений, полученных с помощью проецирования уравнений (13) и (14) на сетку с учетом соотношений (16) и (17):

(18)

(19)

для .

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

(20)

(21)

В частности вычислены и .

Чтобы вычислить , подставим в уравнение (19) вместо и . Учитывая начальные условия , получаем:

(22)

Вернемся теперь к уравнению (2). Поделив его правую и левую части на , заметим, что его можно проинтегрировать по времени . Учитывая начальные условия (3) и (4), получим:

.(23)

Проецируя это уравнение на ту же сетку, показанную на рис.2 и приближая первую производную по времени первой разностью

,

получим:

(24)

Вычитая первое уравнение из второго, получаем:

.

Теперь, аналогично приближая с помощью первой разности

,

имеем:

(25)

Заметим, что из начальных условий

, для

следует, что

.

Таким образом, соотношение (22) примет вид:

(2.20)

Выведенные соотношения (20), (21), (25) и (26) позволяют построить схему численного решения поставленной обратной задачи.

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

4. РЕШЕНИЕ ПРЯМОЙ ЗАДАЧИ

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

(27)

(28)

Где

.

Подставляя в (20) и (21) и исключая оттуда , и , а также учитывая, что , получаем соотношение

,(29)

с помощью которого вместе с соотношением (25) получаем схему решения прямой задачи.

5. ПРОГРАММНАЯ РЕАЛИЗАЦИЯ

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

Для построения алгоритма решения обратной задачи были использованы рекуррентные формулы (20), (21), (25), (26) вместе с необходимыми начальными и граничными условиями.

Схема вычислений была построена так, как показано на рис. 3

Рис. 3 Схема вычислений для решения обратной задачи

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

Для построения алгоритма решения прямой задачи были использованы рекуррентные формулы (25), (27)-(29) вместе с необходимыми начальными и граничными условиями. Схема вычислений была построена так, как показано на рис. 4.

Рис. 4 Схема вычислений для решения прямой задачи

В программе организован цикл по от 2 до N., внутри которого - цикл по от 2 до (N-k+1). На рис.4 показан шаблон для вычислений и через значения этих же сеточных функций в узлах сетки и . Полученная схема имеет первый порядок точности по , где - шаг сетки по направлению . Для выполнения всего процесса вычислений, как и для обратной задачи, требуется операций в силу того, что вычисления проводятся последовательно с помощью рекуррентных формул.

6. РЕЗУЛЬТАТЫ ТЕСТИРОВАНИЯ НА МОДЕЛЬНЫХ ДАННЫХ

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

, , .

В первом тесте использовались функции:

и .

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

Полученные результаты показаны на рис.5 и рис.6:

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

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

2. Для иллюстрации работы программы на негладких и разрывных функциях рассматривался следующий пример.

Используемые функции:

и

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

Рис. 7. График исходной функции для прямой задачи.

Рис. 8. График функции , полученной при решении обратной задачи.

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

ЗАКЛЮЧЕНИЕ

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

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

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

ЛИТЕРАТУРА

1. Алексеев А.С. Некоторые обратные задачи теории распространения волн // Изв. АН СССР. Сер. геофиз. - 1962. №11-12. с.1514-1531.

2. Лаврентьев М.М., Романов В.Г., Шишатский С.П. Некорректные задачи математической физики и анализа. М.: Наука, 1980.

3. Романов В.Г. Обратные задачи математической физики. М.: Наука, 1984.

4. Белишев М.И., Благовещинский А.С. Динамические обратные задачи теории волн. СПб.: Изд-во СПб. Ун-та, 1999.

5. Доровский В.Н., Перепечко Ю.В., Роменский Е.И. Волновые процессы в насыщенных пористых упругодеформируемых средах // ФГВ. 1993, №1. с. 100-111.

6. Imomnazarov Kh.Kh. Estimates of conditional stability of some combined inverse problems for Maxwell's equations and equations of porous media // Comp. Appl. Math., v. 20, 2001, pp. 20-34.

7. Nei-Mao Chen, Yu-Hua Chu, and John T. Kuo Analysis of the characteristic-integration method for one-dimensional wave inversion // Geophysics, 1988, v. 53 No. 8, p. 1034-1044.

8. Roberto A. Tenenbaum and Moyses Zindeluk An exact solution for the one-dimentional elastic wave equation in layered media // Journal Acoustical Society of America, 1992, v. 92, No. 6, p.3364-3370.

9. Roberto A. Tenenbaum and Moyses Zindeluk A fast algorithm to solve the inverse scattering problem in layered media with arbitrary input // Journal Acoustical Society of America, 1992, v. 92, No 6, p.3371-3378,

Размещено на Allbest.ru