CC BY
54
ИИІІІІІІІІІІІІ^^ ..<^1111111 і іе І
СЕМИНАР 3
ДОКЛАД НА : СИМПОЗИУМЕ "НЕДЕЛЯ ГОРНЯКА -2000” : ....і.;.;.;.;.;.;.;.;.;.;.;.;.;.;.; МОСКВА, МГГУ, 31 января - 4 февраля 2000 года
І !і&ь *
:::: ^ Л.С. Загорский, 2000
і! і! УДК 550.342:51:622.831
У Л.С. Загорский
ОБ ИНТЕРПОЛЯЦИИ В ОБРАТНЫХ ЗАДАЧАХ
■'■'■'■'■'■'■'■СЕЙСМИКИ'"''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Перед проведением сейсмических измерений мы располагаем редкой сеткой геологических скважин. Задача интерполяции на данном этапе исследований превращается в задачу аппроксимации, т.к. степень аппроксимирующего многочлена т «П - степени интерполяционного. Следующий этап решения обратной задачи- проведение сейсмических измерений методом сейсмической томографии. Построение аппроксимирующего многочлена в задаче разложения искомой функции медленности по базисным функциям составляет следующий этап решения. В этом случае решение задачи сейсмической томографии сводится к нахождению коэффициентов разложения по этим функциям. Задача решается методом последовательных приближений [1,2]. В этом случае удобна схема с использованием ортонор-мированных на множестве точек X систем функций , которые являются полиномами Чебышева [3]:
п А!!]
Рк,п(1) = £ (~^сС+>—Т, (1)
s = 0 п[ 1
где к = 0,1,2,.т - степень полинома, т < п; Скс,Скс+s -биномиальные коэффициенты; число точек разбиения, уменьшенное на единицу; t = (х - Хо )/h; h - шаг;
= t(t - 1 )(t - 2)...(t - s + 1) ;
п[= п(п -1 )(п - 2)...(п - s +1).
Решение задачи ищется в виде обобщенного полинома вида:
м
т(х) = £ Ь} ф }(х) (2)
}=1
где т(х) - медленность среды, величина, обратная скорости распространения сейсмической волны V(x). Основная же система уравнений выглядит так [1,2]:
di= | т(x)ds , (3)
4
— время волны в лучевом приближении, at — координаты
источника, c — координаты
приемника, ds — длина дуги.
M
Или dt = ^ Gjbj , где матрица G j=1
определена следующим выражением:
Gj = |фj(x)ds . (4)
Т.е. имеем систему уравнений d = Gb для M неизвестных bj ,которая решается итерационными
методами SIRT либо ART [1,2]. Примеры практического применения метода SIRT содержатся в работах [2,4].
При решении задачи в плоскости, перпендикулярной пласту угля, удобна следующая схема [5,6,7]:
a) в обратной задаче Штурма- Лиувилля для волн типа Лява определяется по первой моде нулевое приближение к распределению Vs(z) и интервал его поиска;
b) по высшим модам результат уточняется на основе решения возмущенного интегрального уравнения Гельфанда- Левитана -Марченко выше и ниже пласта (задача на всей оси).
После разделения переменных в волновом уравнении для волн типа SH задача сводится к обратной задаче Штурма-Лиувилля:
— д2У/dZ2 +Q(Z)Y = XY (—ю< Z <х>) (5)
В общем случае обратная задача Штурма-Лиувилля на всей оси состоит в восстановлении функции Q(Z) по матрице-функции RM(X) [8].
Решение обратных динамических задач сейсмики началось с работ А.С. Алексеева [9], который впервые использовал динамические параметры сейсмической записи и разработал методы решения двух- и трехмерных задач сейсмики. Однако, если не учитывать неидеальную упругость среды, то это существенно влияет на точность получаемых сейсмических разрезов.
Оставался открытым вопрос о нахождении таких методов решения обратных задач, которые позволяли бы использовать минимальное количество параметров реальной сейсмической записи для получения результатов с высокой степенью разрешенности.
Представляется целесообразным для этой цели использовать интерференционные волны типа Лява и головные волны.
Поэтому был разработан метод решения обратной задачи на базе интегрального уравнения Гельфанда-Левитана с использованием волн типа Лява [5,6,7]. Определена спектральная матрица-функция уравнения Штурма-Лиувилля
Наиболее просто обратная задача решается в случае доминирования первой моды в спектре волн типа Лява.
Начальные и граничные условия ставятся в соответствии с методикой работы [11].
Для шахтной сейсморазведки: начальные условия:
f = f(z,r,t) + U о (r,z)b'(t) + Ul(r,z)b(t);f = 0,t < 0; граничные условия:
N
і і -U
U n = U n ;ц(г)—
lz=0+s \z=0—s 7 -z
z=0+s
-U
= ц( z)—
dz
z=0—s
Для поверхностной сейсморазведки: начальные условия:
f = f(z,r,t) + U0(r,z)b' (t) + Ul(r,z)b(t);f = 0,t < 0;
граничное условие:
I n / OU CT zz=0 = 0;v(z)—
z=0
dU
: 0, откуда-------
dz
= 0.
(7)
(8)
(9)
z=0
После разделения переменных и подстановки по методу Трикоми уравнение оператора для волн типа Лява записывается в следующем виде:
2
д 2U
L(U) =—^- + (—d(z) +
— k2 )U = 0,
dz
VS (z)
1 дц
1 д 2ц
где d(z) = —- (—/ц) +- — / ц 4 dz 2 dz2
(10)
(11)
Необходимым условием является сопряжение и в точке z=0. Здесь z- координата в глубь массива. Плотность p(z) и модуль сдвига ц( z) связаны со скоростью поперечных волн
соотношением V, = —
'р
Решение одномерной обратной задачи ищется при условии непрерывности функции скорости поперечной волны от координаты z (или аппроксимации к ней) и далее экстраполируется результат решения одномерной обратной задачи на трехмерную с помощью концепции локальных дисперсионных кривых фазовых и групповых скоростей и амплитуд волн типа Лява.
Построение интервалов поиска решения сводится к аналитическому продолжению непрерывной функции V, (2) с использованием формулы следов:
N
q(t) = Еak + ^k — 2%k(t)
(12)
k=1
при условии p(z) = const, либо ()z / д/ц = o( 1).
Можно для этого использовать программу «D3»[13], где использовано условие §l(zi) = ai ,Pi, либо данные
геологоразведочных скважин, пробуренных на исследуемом участке пласта.
Из работ Б.М. Левитана и анализа вещественных
амплитуд волн Лява a j прямо следует, что все точки £ j (t) удовлетворяют дифференциальному уравнению:
dt
± 2a ^ — R( % j(t))
~N
П ( %к —% j(t)) k * j
j = 1....N,
(13)
R( %) = %П (%—*k)( %—fo), k=1
где
%j(z) =
V2 (z)
-—к
j
вертикальное волновое
число, которое меняется в процессе углубления в массив и определяет поворот сейсмических лучей в волноводе, ^ сдвиг.
Если N=1, то знаменатель дроби равен 1. Точка ^
всегда принадлежит лакуне, что обеспечивается сменой знака перед корнем, иначе частота пропускается.
Если N=1, то можно использовать следующее уравнение:
M M
— V/(t)X( X^i) — (2nF(i)/V(i)) — (щ(і) + ^(i))) — Е( 2nF(i))
д Vs(t) _ _______i=1__________________________________________i=1_________
dt
MVs(t)
(14)
Здесь а-і(і)и Рі(і^ - границы 1 лакуны в спектре волн типа Лява; F(i)- частота волн типа Лява; і - номер частоты, полученной после спектрально- временного анализа основного тона волн типа Лява; V(i) - фазовая скорость волны типа Лява.
Если N Ф1, то следует использовать следующее уравнение:
(^тт _ N ,а , в „, ----Г— = Eak + ^k — X k — %k(t)
V ц(t) k=1
(15)
Если Хк <<1, а также точно известно значение скорости поперечной волны V, (И) на глубине первого шага h от поверхности, то достаточно использования лишь одной моды, при этом ошибка представления функции R(Хк)
составляет о( Х\).
При использовании не более трех сдвигов на частоту и известной скорости поперечной волны Vs(И) также можно использовать одну моду, что следует из теории интерполяции.
Основная идея метода состоит в аналитическом продолжении функции скорости поперечной волны, построении интервалов поиска решения и далее решения интегрального уравнения Г ельфанда-Левитана.
Из анализа представления для функций Вейля, полученных Б.М. Левитаном,
V—R( Лk) ,
Я lS(X)l(X — Z)^ ' “ы \s'(4k)\(4k — Z)
Z 1f JRX) m1(Z) ----------------
Z 1t 4RX) . 2 N 4~ R(%)
m2 (Z) =-i^^T^dX+ 2 ------
я |S(X)(X- Z) к=1 S'(лк)\(- Z)
(16)
(17)
(m1 определяет потенциал на интервале (—ж,0) ; m2 — на интервале (0,ж);
N
R( X ) =ХП ( X — ak)( X — $k); k=1
2
СО
аки^к - границы лакуны в спектре; X — собственное значение оператора Штурма-Лиувилля,
Х = ю2 /V2 (0) — к2 ) следует, что:
Б() = 0 соответствует каналовым и поверхностным волнам (дискретный спектр), а нули радикала
±.{рйїї)8(Х)——ЩХ) = Q(Z1,X) (19)
соответствуют боковым волнам, непрерывный спектр. Собственное значение оператора имеет вид:
Х =
(20)
Элементы Х),с(Х), г[(Х) спектральной матрицы-функции КМ (Х ) выражаются через вещественные многочлены Р(?1 ,Х),Б(Х) и функцию К(Х) , введенные Б.М. Левитаном [8]:
Р(21,А) = П (А - £ (21 )),£ (0) = 4,5(А) = (А - По )П (А - Ъ (0)),
_ 1 P(z1 ,Х) 1 Б(X)
дХ 2п ±Л[Я(Х)’ дХ 2п ±Л[Я(Х)
дц_ 1 Q(Zl ,Х)
дХ 2% ±Л[Я(Х)
^ Х , п, Х , N ±УІ— К($к^1))
Q(Zl,Х) = P(Zl,Х) ■ 2--------Ї-------;---------.
к=1(Х — ^к(21)) ■Р (^К(21))
Функция Грина возмущенного уравнения Штурма-Лиувилля (А.Б. Хасанов) [10]
G± (21 ,22 ,Х) = -1 [ г± (Х_)у± (21 ,Х/)ф± (22 ,Х)р(Х)дХ + 2% Е
Е N
N 1 ~ ~
+ 2І Е—У± (21,Х к)Ч± (22,Х к) к=0 тк
где т± = f±(21 ,Хк)\ - нормировочные числа,
(21)
і і^і 12 ' '
\а(Х)\1 — \Ь(Х)\1 = 1,Ш(/^) = / • g — / • g, (25)
здесь W - Вронскиан пары фундаментальных решений возмущенного уравнения Штурма-Лиувилля, которое сводится к следующей системе интегральных уравнений [12]:
21
— А(г1 ,х2) — С(г1 ,х2) — |С(ї,22) ■ А(г1 ¿)д = 0, (26)
—21 —21
А(21,22) + в(21,22) — | 0(ї,22) ■ А(21 ,ї)дї = 0,
21
Поправка р(21 ) = 2дА(21,21) / д21.
Система интегральных уравнений решается относительно ядра А(21,21) и находятся Vs(2) и р(2) из условия:
Ь г 2
2'\\2дА(к+1 )(21,21,1)/д21 — р<к)(21,1)\ }1 /2 /Ь < Ае (27)
I=1 Iі 1
где к- номер итерации, р(к) (21) — вычисленная поправка на предыдущей итерации, 1- номер частоты квантования
ю1.
Поскольку на больших расстояниях R от источника вкладом волн непрерывного спектра
0( 1 /Я) можно пренебречь по сравнению с волнами
дискретного спектра 0( 1 / т/я ) , то функция Грина задана с точностью до 0( 1 /Я) [5]:
G( 21,22 ,г,у,і) =
N ^ 1 д
= 2І 2І ~Ч±(21,Хк)у±(22,ХкМ1(Ккг)Е-ш к=0 0 тк 2%
Способ нахождения размеров лакун в спектре волн типа Лява изложен в работах автора [5-7].
Определение размеров первой лакуны необходимо для вычисления нулевого приближения к решению задачи и границ его поиска.
Значения ОС1 и (31 определяются следующим образом:
(28)
(Х1 —^1(Ю) — Х2 +Х1 •р1
А2 ^Х1
/(—Х1 +Р1) (29)
ч?(21,Х) =Л1 Р^ЕХРІ ±І^ЩХ) I1
ди 0 Р(и,Х)
т± = ^Щ~к) ,Р(Х) = Р(Х; /4ЩХ),
Р( Х к)
(22)
г - (Х) = Ь( Х )/а(Х ),Ь( Х) __ Ш{Г—,/+ },а( Х; =
2І,[ЩГ)
Р( Х )
ІІ^ЩХ)
(23)
№{/+,/— },
г + ( Х ^
а( Х )
2
со
к=1
к=1
(І1 =
1
где А1 - амплитуда первой моды нормальной волны Лява;
(30)
р? (—§1+ Х1) + &§! + (Х1 —\1(к)) — х2 (Х1 -^1(к))
12Х1
А2Х
А2§1
А12 Х1
А12 §1
-§1 + х2§1 = 0.
где Р1 =
— Ь ±^Ь2— 2а
4 ас
(31)
причем знак перед корнем выбирается так, чтобы значение Р1 было неотрицательным.
Для обеспечения устойчивости вычислений в знаменателях формул (29) и (31) следует прибавить Є .
Из вышеизложенного следует, что если амплитудный ряд имеет вид: А1 Д0...0, то (,Р1 вычисляются точно,
если все члены Ак ряда не равны нулю, то ((1 ,Р1
вычисляются приближенно.
При этом в качестве амплитуды первой моды следует
/ А1
брать А1 —---------
N
П (Х1 —§к) к=2
А
А
для выполнения условия:
(32)
где р (Х1) =(Х1 -^1),
N
Р( Х1) = П (Х1 -1к). к=1
Погрешность определения скоростного разреза составляет Р = 20% при доверительной вероятности 0.95.
При поверхностной сейсморазведке интервальные фазовые скорости находятся по методике СВАН для соседних сейсмоприемников.
Рис. 1. Вертикальный сейсмический разрез для пласта угля ш.«Алмазная» АО «Гуковуголь» при удалении 140 м от источника возбуждения
асимптоты дисперсионной кривой. Здесь для каждой частоты необходимо построить распределение групповых скоростей волн, а затем, используя весь диапазон частот, локальные
дисперсионные кривые в каждом блоке массива. Этого можно добиться с использованием программы полосовой фильтрации FLTRTMG, разработанной в УкрНИМИ и программы сейсмической томографии SOND, разработанной в МГУ им. М.В.
Ломоносова, определение амплитуд
нормальных волн возможно с использованием программ МИТП РАН. Амплитуду первой моды определяет наша программа «Б3» [13].
Итак, для каждого блока массива определены
фазовые скорости, которые являются входными данными для нашей программы «SUPERS21» [14].
Для сейсмопросвечивания массива с целью
определения вертикального сейсмического разреза
следует использовать цифровую сейсмостанцию SSS-1 («Дружба») с относительной погрешностью единичного измерения времени At = 1% . Цифровая сейсмостанция необходима для обработки получаемой информации на ЭВМ. Метод восстановления вертикального
сейсмического разреза по спектральной матрице-функции уравнения Штурма-Лиувилля для волн типа Лява был опробован для условий АО «Эстонсланец», «Ленинградсланец» и «Гуковуголь», для границ штатов Оклахома и Техас, верхней мантии в Индии, а также для слоя между полупространствами и трехслойной модели среды по опубликованным литературным данным.
Анализ результатов расчетов по разработанной программе «SUPERS21» показал, что максимальное отклонение рассчитанного сейсмического разреза от модельного не превосходило 13,4 %, при этом
восстанавливается строение всего сложного волновода, состоящего из N вложен-ных волноводов. Указанные выше программы внедрены в АО «Гуковуголь», где используются для решения задач прогноза горногеологических условий залегания угольного пласта.
Пример практического применения метода (с использованием программы автора «Б3») приведен на рис.1-2. Показаны, для условий шахты <Алмазная» АО «Гуковуголь»: вертикальный сейсмический разрез (скорость поперечной волны)
При шахтной сейсморазведке определить интервальные групповые скорости можно лишь методом сейсмической томографии при сейсмопросвечивании массива с учетом существования высокочастотной
2
2
+
Рис. 2. Скорость поперечной волны на глубине +1, +2, +3 м от оси пласта ш.»Алмазная» при различных удалениях от источника возбуждения
при удалении 140 м от источника возбуждения колебаний (рис. 1); скорость поперечной волны на глубине +1, +2,+3 м при различных удалениях от источника возбуждения (рис. 2).
Удалению 90 м соответствует тектоническое нарушение, длина лавы составляла 180 м.
ш ш 0 п ш % * ш ш а а а в « в ф а в '
1. Иванссон С. Сейсмическая скважинная томография- теория и методы вычислений //ТИИЭР.-1986. - Т.74.- № 2. С. 99-111.
2. Ефимова Е.А., Рудерман Е.Н.
Возможности применения цифровой томографии для интерпретации
геофизических данных. - М.: ВИЭМС, 1982. - 55с. (Обзор ВИЭМС).
3. Демидович Б.П., Марон И.А.,
Шувалова Э.З. Численные методы анализа. Приближение функций, дифференциальные и интегральные уравнения.-М.: Наука,
1967.- 368 с.
4. Рубан А.Д., Черняков А.Б., Загорский Л.С., Новиков А.Н. Томографические методы восстановления волновых полей в задачах геомеханического контроля //ФТПРПИ.- Новосибирск, Изд. Наука Сибирское отд., 1993, № 3. — С.37-44.
5. Загорский Л. С. Восстановление вертикального сейсмического разреза по спектральной матрице-функции уравнения Штурма-Лиувилля //Доклады академии наук. -М.: Наука. - 1998. - Т. 358. № 5.
6. Загорский Л. С. Практическое применение метода восстановления вертикального сейсмического разреза горного массива // Горный информационноаналитический бюллетень. - М.: МГГУ, 1997.-- № 2.-С. 84-86
7. Загорский Л.С. Условия
разрешимости обратной задачи
определения строения массива по дисперсии основного тона волны 8Н. //Сб. трудов V сессии РАО. Проблемы геоакустики: Методы и средства /Под ред. проф. В.С. Ямщикова.- М.:МГГУ.- 1996.-С.42-45.
8. Левитан БМ. О разрешимости обратной задачи Штурма-Лиувилля на всей оси // Дифференциальные уравнения с частными производными.- Новосибирск: Наука, 1980. - С.234-239.
9. Алексеев А.С., Кремлев А.Н., Жерняк Г. Ф. Обратная задача дифракции акустических волн на малых скоростных неоднородностях.- Геология и геофизика, 1981. —№ 1.- С. 111-118
10. Хасанов А.Б. Обратная задача
рассеяния для возмущенного
конечнозонного оператора Штурма-Лиувилля /ДАН СССР, 1991.-Т.318.— № 5.-С.1095-1098.
11. Владимиров В.С. Уравнения математической физики.- М.: Наука, 1981. 512 с.
12. Загорский Л.С. Определение
строения горного массива по спектру интерференционных волн: Автореф.
дисс. докт. физ.-мат. наук/ ОИФЗ им. О.Ю. Шмидта, М., 1998, 35 с.
13. Загорский Л. С. Восстановление трехмерного сейсмического разреза по спектральным характеристикам первой моды волны SH // Информационный бюллетень. Алгоритмы и программы (ГосФАП)— М.: ВНТИЦ , 1997. -№ 1. - С.3-
4.
14. Загорский Л. С. Определение вертикального сейсмического разреза по дискретному спектру волн SH // Информационный бюллетень. Алгоритмы и программы (ГосФАП)- М.:ВНТИЦ,1997.-№ 1. — С.3.
Загорский Лев Сергеевич. — доктор физико-математических наук. Институт горного дела им. А.А. Скочинского.
