Анализ сингулярного разложения линеаризованного оператора динамической теории упругости для случая вертикального сейсмического профилирования Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 12, № 6, 2007

АНАЛИЗ СИНГУЛЯРНОГО РАЗЛОЖЕНИЯ

ЛИНЕАРИЗОВАННОГО ОПЕРАТОРА ДИНАМИЧЕСКОЙ ТЕОРИИ УПРУГОСТИ

ДЛЯ СЛУЧАЯ ВЕРТИКАЛЬНОГО СЕЙСМИЧЕСКОГО ПРОФИЛИРОВАНИЯ*

И. Ю. Сильвестров Институт нефтегазовой геологии и геофизики СО РАН,

Новосибирск, Россия e-mail: ilya_s@uiggm.nsc.ru

The singular value decomposition analysis of linearized forward operator of dynamic elasticity is presented. This operator is used in Newton method for solving one of the inverse problems of petroleum seismic. Such analysis allows to determine main features of the problem that should be taken into account when the numerical algorithm for solving the inverse problem is developed. Using this analysis the optimal parametrization for inversion was founded out. The principal impossibility of determining the smooth components of solution was also shown. Finally, a significant influence of the source wavelet to the solution of the inverse problem has been demonstrated.

Введение

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

* Работа выполнена в сотрудничестве с Московским научным центром фирмы "Шлюмберже", а также при частичной финансовой поддержке Российского фонда фундаментальных исследований (грант № 05-05-64227) и СО РАН (Молодежный Лаврентьевский грант 2007).

( Институт вычислительных технологий Сибирского отделения Российской академии наук, 2007.

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

В качестве примера можно привести самый первый шаг решения обратной задачи — выбор параметризации среды. Как известно, различные тройки параметров, описывающие упругую изотропную среду, не эквивалентны для обращения. Еще в работе [4] с привлечением диаграмм рассеяния от точечных объектов было показано, что упругие импедансы являются более подходящими для обращения отраженных волн, чем параметры Ламе или скорости продольных и поперечных волн. Эти рассуждения проводились на физическом уровне строгости, поэтому представляется небесполезным более строгое обоснование этого факта именно для интересующей нас системы наблюдений. С этой целью в настоящей работе привлекается анализ сингулярного разложения производной оператора, отображающего параметры среды в данные наблюдений, что позволяет обосновать и другие известные из практики явления, возникающие при решении обратной задачи. Отметим, что частично аналогичный анализ для поверхностной системы наблюдения был выполнен в работе [5].

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

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

Р

Р~

д / ^ ди1 А шу и + 2^—— дх1

дх1 д

д2и1 д2и2

дг2 дх1 У чдх2 ' дх1

дх2

ди1 ди2

дх2

ди1 ди2, .

М^-т + ^гМ -

дх2

А и + 2^

дх1 ди2

дх2

д

^(г) —8(х - х3), дх1

д

^(г)—5(х - х3),

дх2

и

4=0

ди д

0,

(1)

4=0

где х3 — вектор координат источника. Здесь мы полагаем, что ось х1 находится на поверхности земли, а х2 перпендикулярна к ней и направлена в глубину. В качестве зондирующего сигнала использовался импульс Рикера:

^(г) = (1 - 2(п^о(г - 1/^о))2) ■ ехр(-(п^ - 1/^о))2)

с доминирующей частотой щ = 30 Гц.

По отношению к системе (1) рассмотримобратную задачу. По информации

и

= и°5(х2,г); 0 <г<т; х2 е {хь ,...,хМз\хи < ... < х^Л

Х1=0

(2)

о режиме колебаний, записанной в приемниках сейсмических волн, расположенных в вертикальной скважине, восстановить параметры среды А,ц,р ниже горизонтальной линии х2 = х^, на которой находится последний приемник, т. е. ниже забоя скважины. Здесь и далее мы полагаем, что источник колебаний находится на поверхности земли. Систему (1)-(2) можно рассматривать как нелинейное операторное уравнение:

где т — параметры среды, и°Ьз — данные наблюдений, В — нелинейный оператор, действующий из пространства моделей в пространство данных. Для решения этого нелинейного уравнения возможно использование метода Ньютона:

где О В — производная Фреше оператора В; гпк — модель среды на к-м шаге.

Для рассматриваемого класса задач оператор ОВ является компактным (см. [6]), поэтому задача (4) является условно-корректной, по Тихонову. При решении такого рода задач могут применяться методы регуляризации, основанные, например, на усечении сингулярного разложения компактного оператора (прилож. А.1).

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

В (т) = и°Ьз,

(3)

ОВ(тк) < тк+1 - тк >= и° 3 - В(гйк),

(4)

Источник

с

Скважина с приемниками

Плоская Р-волна

1000

Забой скважины

Целевая область

2000

0 М 400

Рис. 1. Модель среды, используемая для анализа сингулярного разложения

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

2. Анализ сингулярного разложения

Сингулярное разложение матрицы ОБ производилось при помощи пакета программ линейной алгебры ЬАРАСК 1.1. При этом фактически мы искали собственные числа и собственные векторы матрицы ОБ*ОБ.

На рис. 2 показаны сингулярные числа матрицы ОБ. Видно, что они быстро стремятся к нулю. Это и следовало ожидать, так как данная матрица — аппроксимация компактного оператора. Таким образом, задача действительно является плохо обусловленной. Отметим, что перегиб графика в районе 2500-го сингулярного числа отвечает пределу используемой машинной точности. В данном случае расчеты проводились с одинарной точностью. Всем следующим сингулярным числам и векторам доверять нельзя. Это означает, что даже в идеальном случае отсутствия помех в данных устойчиво может быть восстановлена только проекция истинного решения на 2500 старших сингулярных векторах. Естественно, что в практических задачах такой точности входных данных добиться невозможно, поэтому будем считать допустимым числом обусловленности V(А) ~ 10-4, что соответствует числу старших сингулярных векторов, привлекаемых при усечении сингулярного разложения г = 1500 (см. прилож. А.1).

2.1. Выбор параметризации среды

Первым этапом при решении обратной задачи (4) является выбор параметров среды, оптимальных для обращения. Как известно, упругая изотропная среда описывается набором из трех параметров. В дополнение к плотности (р) ими могут быть: параметры Ламе (А,^); скорости продольной и поперечной волны (ур = л/(А + 2^)/р, у3 = л/р>/р); упругие импедансы (1Р = рур, 1Б = ру3). Хорошо известно, что они не эквивалентны с точки зрения обращения [4]. Поэтому возникает вопрос: какие параметры следует пытаться восстанавливать ниже забоя скважины при интересующей нас системе наблюдений?

Рис. 2. Сингулярные числа оператора БВ в логарифмической шкале

Отметим, что так как О В — производная Фреше оператора В, то его матричное представление меняется в зависимости от параметризации:

ОВ ОВ дт

дт' дт дт''

где т' и т — соответственно новый и старый набор параметров.

Наиболее эффективной параметризацией для обращения будет та, которая обеспечит для одной и той же размерности матрицы ОВ наименьшее отклонение решения обратной задачи от его проекции на правые сингулярные векторы, соответствующие большим сингулярным числам (прилож. А.1). Для описания этого отклонения применялось разложение искомого решения в ряд Фурье. Чтобы учесть двумерность среды, мы использовали гармоники, наклоненные под углом 5° к горизонтали, так как именно такой наклон границ характерен для практических приложений. В качестве характеристики отклонения выступает угол между гармониками Фурье и линейной оболочкой старших сингулярных векторов. На рис. 3 показаны линии уровня этих углов в зависимости от числа привлекаемых сингулярных векторов. Видно, что для фиксированного числа векторов углы, соответствующие Р-импедансу, существенно меньше, чем для других параметров. Аналогичный эффект имеет место и для Б-импеданса. Таким образом, можно сделать вывод о том, что упругие импедансы являются параметрами,

Vо

>4

е

10

15

8 К

0

1

и

5Г 20 §

О

К

25

500 1000 1500 2000 2500

Число старших сингулярных векторов Р-импеданс

500 1000 1500 2000 2500

Число старших сингулярных векторов Плотность (в базисе импедансов)

И с Л 5

е § 10

в и о

I"

1-н &

120 К

25

500 1000 1500 2000 2500

Число старших сингулярных векторов

500 1000 1500 2000 2500

Число старших сингулярных векторов

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

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

Из рис. 3 видно, что при использовании Р-импеданса (впрочем, как и скорости) углы уменьшаются с увеличением пространственной частоты. Значит, только высокочастотные компоненты решения могут быть восстановлены, а это есть проявление известной проблемы определения гладкой ("трендовой") составляющей скорости. Поэтому только первый шаг итерационного процесса Ньютона (4) имеет смысл. Последующие шаги требуют уточнения скоростной модели с использованием так называемых кинематических подходов, в которых в качестве входных данных для обращения используются именно времена прихода волн.

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

1100

8щ>0

2000

8т2=0

100 М 400

у-V-'

8т 1

юо м

~у—

ётг

ёт,=0

б 1100

1500

400 1

100

8т,

400 )

2000

8Х> 0;

8Х

8/1=0;

100 м 400 100 м

-у-

8ц

8р=0

400 )

100 М 400

у-у-'

8р

в

1юо 81)р>0;

1500 ■ г

2000 Ь

100 М

<51Л=0

400

8Ур

81Р> 0;

1100

1500 г

-- 2000 400 100 М

) V.

-у—

81Р

400 1

Ш=0;

100 М у—

8Ш

8р=О

400 1

100 М 400

У-V-'

8р

Рис. 4. Поведение г-решения в зависимости от выбранной параметризации: а — местоположение разрыва одного из параметров; б — г-решение, соответствующее исходной неоднородности по Л и однородной среде по у и р; в — г-решение, соответствующее исходной неоднородности по ьр; г — г-решение, соответствующее исходной неоднородности по 1Р

100 М 400

П> Ш р к ц р Ьр V, Р

Рис. 5. Второй правый сингулярный вектор оператора ББ в зависимости от выбранной параметризации

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

На рис. 4, б показано г-решение для параметров Ламе, соответствующее исходной неоднородности по А (подчеркнем, что по остальным параметрам среда однородна). Такое решение содержит возмущения по всем параметрам, что делает невозможным оценить, какой же тип неоднородности имеет место на самом деле. Тот же эффект характерен и для скоростей продольных волн (рис. 4, в). Более того, можно показать, что поведение г-решения и для остальных параметров из данных двух типов параметризации будет аналогичным. В противоположность этому г-решение, соответствующее неоднородности по Р-импедансу, состоит из ненулевой компоненты только по этому параметру (рис. 4, г). То же будет происходить и для ¿"-импеданса.

Такое поведение г-решения можно объяснить, рассматривая структуру одного из старших сингулярных векторов для различных параметризаций. Этот вектор — базисный в пространстве моделей, поэтому он состоит из трех частей, отвечающих каждому параметру. На рис. 5 изображен второй правый сингулярный вектор оператора ОБ для изучаемых трех наборов параметров. Отчетливо видно, что компоненты, соответствующие ^ и р, похожи между собой. Учитывая, что та же картина наблюдается и для других старших векторов в этой параметризации, получается, что г-решение, являющееся линейной комбинацией этих векторов (см. [7]), будет состоять из похожих компонент по всем параметрам, что мы и наблюдали на рис. 4, б. Тот же эффект обнаруживается и при использовании скоростей. В противоположность этому компоненты старшего вектора, соответствующего упругим импедансам, отличаются друг от друга, поэтому линейная комбинация таких векторов позволяет приблизить искомое решение с достаточной точностью. Таким образом, действительно подтверждается, что оптимальными параметрами для обращения являются упругие импедансы.

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

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

1000 г

Разрыв _ Р-импеданса

2000

М 400

Рис. 6. Влияние формы импульса зондирующего сигнала на решение обратной задачи: а — местоположение разрыва Р-импеданса; б — истинный разрыв (штриховая линия ) и восстановленный разрыв (сплошная линия)

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

Заключение

Приведенный анализ сингулярного разложения линеаризованного оператора, отображающего параметры среды в данные НВСП, — важный предварительный шаг перед непосредственным численным решением обратной задачи. Хотя использовалась достаточно простая модель среды, удалось строго обосновать явления, которые наблюдаются при реализации полного обращения волновых полей [2]. Получено подтверждение оптимальности выбора импедансов для параметризации упругой среды ниже забоя скважины. Плотность при этом определить невозможно. Также было показано проявление проблемы определения гладкой составляющей скорости, что приводит к необходимости дополнительного использования кинематических подходов. И, наконец, выявлено существенное влияние импульса исходного сигнала на решение обратной задачи. Поэтому необходимо либо проводить дополнительную деконволюцию получающегося решения, либо рассматривать постановку задачи, в которой импульс определяется наравне с параметрами среды [8].

Автор выражает глубокую признательность и благодарность В.А. Чеверде и В.И. Костину за постоянное внимание и консультации во время выполнения работы.

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

[1] Tarantola A. Inversion of seismic reflection data in the acoustic approximation // Geophys. 1984. Vol. 49, N 8. P. 1259-1266.

[2] Mora P. Nonlinear two-dimensional elastic inversion of multioffset seismic data // Geophys. 1987. Vol. 52, N 9. P. 1211-228.

[3] Pratt R.G., Shin C., Hicks G.J. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion // Geophys. J. Int. 1998. N 133. P. 341-362.

[4] Tarantola A. A strategy for nonlinear elastic inversion of seismic reflection data // Geophys. 1986. Vol. 51, N 10. P. 1893-1903.

[5] Lebrun D., Richard V., Mace D., Cuer M. SVD for multioffset linearized inversion: Resolution analysis in multicomponent acquisition // Geophys. 2001. Vol. 66, N 3. P. 871882.

[6] Cheverda V.A., Clement F., Khaidukov V.G., Kostin V.I. Linearized inversion of data of multi-offset data for vertically inhomogeneous background //J. Inv. Ill-Posed Problems. 1998. Vol. 6, N 5. P. 453-484.

[7] Kostin V.I., Tcheverda V.A. r-Pseudoinverse for compact operators in Hilbert space: existence and stability //J. Inv. Ill-Posed Problems. 1995. Vol. 3, N 2. P 131-148.

[8] Воронина Т.А., Чеверда В.А. Оптимизационный подход к обработке данных метода вертикального сейсмического профилирования // Геология и геофизика. 1994. № 5. С. 127-139.

[9] Годунов С.К., Антонов А.Г., Кирилюк О.П., Костин В.И. Гарантированная точность решения систем линейных уравнений в евклидовых пространствах. Новосибирск: Наука, 1988.

[10] Трехмерные задачи математической теории упругости и термоупругости / Под ред. В.Д. Купрадзе. М: Наука, 1976.

A. Приложения

A.1. Сингулярное разложение прямоугольной матрицы

Напомним основные факты, касающиеся сингулярного разложения прямоугольных матриц [9].

Для любой матрицы A размером m х n (в дальнейшем будем считать, что m > n) существуют ортонормированные базисы {Xi,... xn} и {yi,... ym} и числа s1 > s2 > ... > sn > 0 такие, что

AXi = Sitji. (5)

Выражение (5) называется сингулярным разложением матрицы A, векторы Xi и yi — соответственно правые и левые сингулярные векторы, а si — сингулярные числа. Число обусловленности v(A) выражается через сингулярные числа:

V (A) = *.

sn

Рассмотрим теперь систему линейных уравнений:

Ах = у. (6)

Ошибку 8х в ее решении можно оценить через ошибку 8у во входных данных, как

Щ < V(А) СЩ,

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

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

X[r] = ^2(x,Xi)Xi, (7)

i= 1

где круглыми скобками обозначается скалярное произведение. Вектор Х[г] называется г-решением задачи (6). Число г привлекаемых сингулярных векторов контролирует обусловленность задачи и позволяет решить ее с приемлемой точностью, если знаешь относительную ошибку во входных данных. Можно показать (см. [7]), что построенное таким образом решение устойчиво как относительно погрешности во входных данных, так и относительно погрешности, возникающей при дискретизации оператора.

Л.2. Построение оператора ДВ

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

Для этого нам удобнее уравнения (1) записать в частотной области. В качестве параметров среды будем использовать параметры Ламе А, р и плотность р:

2 / д .. ^ ^ диЛ д ( дщ ди2 рш щ + —— А шу и + 2р—— + —— р—--+ р

дх1 \ dx1 J дх2 V dx 2 dx

д

= F (ш)~-S(x 1 - Xis, X 2 - X 2s) , (8)

dx1

2 / д / ди2\ д ( ди1 ди2

рш щ + A divu + 2p—— + ^--+

\дх2 \ ox2) дх1 \ дх2 dx1

д

= F(ш) — 8(x 1 - X1s,X2 - X2s) дХ2

+ R.C. (Radial conditions).

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

A(x1, х2) = A0 + 8A(x1, х2), р(х1,х2) = р0 + 8р(х1 ,х2), р(х 1, X2) = ро + 8р(х1, X2),

где А0,^0,р0 — константы, а ¿А^,^), ¿^(ж^^), ¿р(ж1;ж2) настолько малы, что позволяют линеаризовать уравнения (8), и будем искать решение (8) в виде 7 = 7(г) + ¿7, где

,7(0 —

7

решение задачи (8) при А = А0, ^ = р = р0. Тогда из (8) получаем:

д

и р0071 + (А0 + ^0) —-&У ¿7 + ^0Д071

дх1

д / д,(') = —^ ¿Аа1у+

джм дх1

д

и р0¿72 + (А0 + -&У ¿7 + Д ¿72

д / д,1г) д7(г)

--+ о^-

дх

дх2

дх1

2 « ри 71 ,

дх2

(9)

д / , д71г) д72г)

(¿Г

дх1

¿^^ + ) - ^ ( ¿А а1у 7« +

дх2 д^ / дх2 V дх2

2 (¿) ри 7з •

Решение системы (9), записанное с помощью матрицы Грина, имеет вид:

¿и(г, и) = / С(г,г"1 ,и)/(г! )йг!, (10)

где Пм — область возмущений параметров среды (целевая область); О — матрица Грина (Купрадзе); / — вектор правой части. Заметим, что матрица Грина системы (9) совпадает с матрицей Грина для однородной среды, в которой она выписывается явно [10]. Используя ее явный вид и предполагая, что падающая волна есть плоская продольная волна, распространяющаяся под некоторым углом к целевой области, и пользуясь приближением дальней зоны, т. е. считая, что целевая область находится достаточно далеко от источника, можно показать, что (10) переписывается в таком виде:

¿7 = ОВ < ¿т >,

где

¿7

¿71 ¿72

¿т

Ьц ¿12 ¿13 ¿21 ¿22 ¿23

а ¿у — интегральные операторы первого рода следующего вида:

< У(6, 6) >)(х2,и)

Ку (х2, и; ^2,

с ядрами, зависящими от коэффициентов матрицы Грина. Так, например,

т^ / ^ Ч/^Ч I Ш ... \ СОв ^(0 N

КИ(х2,и; ^1,^2) = Л(£)ехр —г(£) ^, F(и),

где

¿(£) =

— 2- ехр 4

«и6

и

г(е + (Х2 — б)2,

СО8 (е)

—^ г(е)

^2

ри2

А + 2^

2

п

м

р

V

р

Поступила в редакцию 7 мая 2007 г., в переработанном виде — 5 сентября 2007 г.