авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 || 3 | 4 |   ...   | 6 |

«УЧЕНЫЕ ЗАПИСКИ САНКТ-ПЕТЕРБУРГСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА № 444 СЕРИЯ ФИЗИЧЕСКИХ И ГЕОЛОГИЧЕСКИХ НАУК Издается с 1958 ...»

-- [ Страница 2 ] --

+i AJ0 (ikrf )ekvs t d, f = vs k cos kzdk 0 i +i (2) BH0 (ikrs )ekvs t d, s = vs k cos kzdk 0 i +i (2) CH1 (ikrs )ekvs t d, s = vs k sin kzdk 0 i 2 v vs 1 + 2 и s = 1 + 2 vs. На основании граничных усло где s = 1+ 2, f = vf p вий (непрерывности радиальной компоненты смещения и напряжения, равенство ну лю тангенциальной компоненты напряжения на границе жидкость упругость) можно c Д. А. Александров, А. В. Бакулин, Б. М. Каштан, Распространение низкочастотных трубных волн в радиально-слоистых проницаемых средах построить систему линейных уравнений для определения коэффициентов A, B и C.

Дальнейшее решение задачи предполагает вычисление интегралов, для чего необходи мо найти полюсы и точки ветвления подынтегральных функций относительно. Для нахождения полюсов необходимо решить дисперсионное уравнение 2 (2) (2) vf 2is v f v 4 H0 (y) H0 (x1 ) s vs s J0 (x) + f J1 (x) 2 + 4s s = 0, (1) 4 2 (2) (2) vs kr0 vs H1 (y) H1 (x1 ) x = ikr0 f, y = ikr0 s, v = ivs, которое накладывает связь на фазовую скорость v и волновое число k. Вещественные корни этого уравнения относительно v (мнимые относительно ) соответствуют трубной волне, которая распространяется без затухания.

В случае низких частот (kr0 0) удается выписать явное выражение для скорости трубной волны:

vs vt =.

f vs + s vf Такое выражение для скорости низкочастотной трубной волны на основании квази статического подхода получено Уайтом [7]. Трубная волна может возбуждаться источ ником, находящимся на стенке скважины [5], полем точечного источника, находящегося вне скважины [4], а также медленной волной в жидком слое, пересекающем скважи ну [2].

Зачастую скважина оказывается заполнена не жидкостью, а двухфазной средой смесью твердых частиц и жидкости. Такую двухфазную среду можно моделировать пористой средой Био [12]. В безграничной среде Био распространяются две продоль ные и одна поперечная волны. Различие между продольными волнами заключается в колебаниях жидкой и твердой фаз вещества: в случае быстрой продольной волны колебания происходят в фазе, в случае медленной волны в противофазе. П. В. Кра уклис показал, что в модели скважины, заполненной средой Био, распространяются две трубные волны [3].

Значительный прикладной интерес представляет модель, состоящая из несколь ких цилиндрических слоев, поскольку с ее помощью можно описать распространение волн в обсаженных скважинах с цементированием или в скважинах с более сложной конструкцией. Такую модель рассмотрел в своей работе А. А. Сидоров [15], приме нив матричный метод, предложенный Л. А. Молотковым [6]. Был рассмотрен вектор W = (rr rur uz rrz )t, состоящий из нормального напряжения rr, смещения вдоль оси симметрии uz, поверхностного напряжения rrz и изменения площади поперечного сечения rur. Основываясь на уравнениях движения в упругой среде, можно показать, что такой вектор удовлетворяет следующему дифференциальному уравнению:

d W = D(r, k, )W. (2) dr Дифференцирование по времени t и координате z снято в результате повторного Фурье-преобразования и перехода к частоте и волновому числу k. Матрица D, помимо 36 Д. А. Александров, А. В. Бакулин, Б. М. Каштан радиуса, частоты и волнового числа, содержит параметры Ламе и µ. Предметом исследования в работе является фундаментальная матрица G, связывавшая значения вектора W для разных значений радиуса:

W (r) = G(r, r0 )W (r0 ).

Легко заметить, что матрица G удовлетворяет дифференциальному уравнению (2) c единичными граничными условиями:

d dr G(r, r0 ) = D(r)G(r, r0 ), G(r0, r0 ) = I.

Такой подход позволяет вместо задачи о нахождении волнового поля в системе из большого количества однородных цилиндрических упругих слоев ограничиться вычис лением матрицы Gi (r, r0 ) для каждой среды. Тогда вектор W в произвольной точке среды r может быть найден в результате перемножения разрешающих матриц Gi :

W (r) = Gn (rn, rn1 )Gn1 (rn1, rn2 )...G1 (r1, r0 )W (r0 ), (3) где ri координаты границ цилиндрических слоев. На основании системы уравнений (3), в которой вектор W (r0 ) представляет собой граничные условия, можно получить дисперсионное уравнение. В работе [8] рассмотрена модель с двумя упругими коак сиальными трубами, пространство между которыми заполнено жидкостью. Показано, что на низких частотах дисперсионное уравнение для такой модели распадается на два квадратных уравнения относительно фазовой скорости. Решением каждого урав нения являются скорости трубной и пластинчатой волн. Таким образом, в этой системе с чередующимися жидкими и упругими слоями распространяются две трубные и две пластинчатые волны.

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

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

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

Зазор между витками обмотки 0.2 мм трубу с пластиковой обмоткой. Для моделирования непроницаемого экрана внутренняя труба заменяется сплошной алюминиевой трубой без перфораций. Основной функци ей экрана и гравийного фильтра является предотвращение выноса в скважину песка из резервуара.

Кроме того, гравийный фильтр рассеивает энергию высокоскоростных потоков неф тесодержащей жидкости, поступающей через перфорации, тем самым защищая экран от повреждений. Поэтому важно знать, насколько качественно был установлен гра вийный фильтр, а также сохраняет ли он свою структуру в процессе добычи. Если из гравийного фильтра будет вымыт песок, появятся участки, где скорость потока жид кости достаточно высока, чтобы повредить экран. В свою очередь, это повлечет вынос в скважину песка из резервуара [17]. Кроме того, вещество окружающей породы может перемешаться с гравием фильтра и частично его заменить. При этом проницаемость 38 Д. А. Александров, А. В. Бакулин, Б. М. Каштан смеси уменьшится со 100–500 до 1 Д. Уменьшение проницаемости на одних участках скважины приведет к перераспределению потоков жидкости и, возможно, к увеличе нию скорости потоков через другие перфорации. Слишком высокая скорость, в свою очередь, может привести к вымыванию гравийного фильтра, разрушению экрана и остановке добычи [17].

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

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

Результаты численного моделирования Скважина без гравийного фильтра. Рассмотрена идеализированная модель скважины со свободной внешней границей, состоящая из четырех цилиндрических сло ев: двух жидких слоев, проницаемого экрана и обсадной колонны. Численные расчеты реализованы при помощи алгоритма численного решения уравнения, описывающего среду Био [12]. Для моделирования проницаемого экрана используется пористая среда Био. Центральная частота сигнала составляет 500 Гц. Конечно-разностным алгорит мом получены значения радиальной компоненты смещения внешней стенки обсадной колонны, поскольку оптоволоконные приемники в эксперименте измеряют именно эту компоненту волнового поля.

В статье А. В. Бакулина [8] показано, что в такой четырехслойной модели, где пори стый экран заменен упругим слоем, в низкочастотном диапазоне распространяются че тыре волны: две трубные и две пластинчатые. Согласно численным расчетам, в модели с проницаемостью экрана 1 мД также наблюдаются две трубные волны. По скоростно му признаку их можно разделить на быструю и медленную. Таким образом, при низких значениях проницаемости пороупругий экран эквивалентен упругому непроницаемому веществу, для которого обе трубные волны распространяются без затухания.

Эмпирическим путем установлено, что в эксперименте проницаемость экрана, со стоящего из перфорированной алюминиевой трубы с пластиковой обмоткой, составля ет 250–1000 Д. В то же время закупоривание пор может уменьшить проницаемость до нуля. Численные расчеты выполнены для моделей, в которых проницаемость экра на варьируется от 1 мД до 100 кД. С увеличением проницаемости экрана обе волны начинают затухать по мере распространения вдоль скважины, однако зависимость за тухания от проницаемости экрана для этих мод различна. Для оценки характера этой зависимости используются двумерные спектры частота медленность:

u(z, t)eit eipz dpd, U (, p) = где u(z, t) сейсмические данные;

частота;

p медленность, скоростной спектр (рис. 3). Скоростной спектр получается в результате интегрирования первого спектра по частоте и перехода от медленности к скорости.

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

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

Поведение медленной трубной волны существенно отличается от быстрой волны.

Как видно на рис. 4, г, затухание медленной волны увеличивается по мере роста прони цаемости. При значении проницаемости более 100 мД медленная трубная волна стано вится неразличима на фоне быстрой трубной волны. Скорость медленной волны увели чивается с увеличением проницаемости (рис. 4, в), в то время как в простых моделях с колонной жидкости, окруженной пористой средой, наблюдается обратный эффект [16].

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

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

На рис. 5, а отображена данная зависимость для модели, в которой проницаемость экрана внутренней трубы была равна 10 мД. Такое распределение смещения практи 40 Д. А. Александров, А. В. Бакулин, Б. М. Каштан Рис. 4. Зависимость скорости и амплитудного коэффициента от проница емости экрана для быстрой (а, б ) и медленной (в, г) трубных волн в модели без гравийного фильтра чески совпадает с моделью, в которой пороупругий экран заменен упругим слоем. Это сходство еще раз подтверждает тот факт, что при низких значениях проницаемости пороупругий экран эквивалентен упругому непроницаемому веществу.

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

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

Проницаемость экрана для моделей: а 10 мД;

б 150 мД;

в 10 Д колебаний. Как только проницаемость экрана увеличивается, разница давлений при водит к обмену жидкостью через экран. Наиболее интенсивным обмен становится при значениях проницаемости около 300 мД, что сопровождается быстрым выравниванием профиля аксиальных смещений, а также резким падением скорости (см. рис. 4, а). При высоких значениях проницаемости вертикальная компонента смещения перестает ме няться с радиусом (рис. 5, в). Профиль радиальной компоненты становится линейным, что характерно для модели с одной трубой, заполненной жидкостью.

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

Скважина с гравийным фильтром. Рассмотрим модель скважины с гравийной набивкой. Основным отличием от предыдущей модели является наличие гравийного фильтра в пространстве между трубами. Для моделирования гравийного фильтра ис пользуется пороупругая среда. Модуль сдвига скелета этого пористого вещества был равен 0.01 ГПа. А. В. Бакулин и др. [8], основываясь на изучении каналовой волны, заключили, что скорость поперечных волн в такой среде должна быть низкой, но не нулевой, порядка 30–80 м/с. Численные расчеты показывают, что в модели с такими параметрами распространяются одна трубная и одна пластинчатая волна, как буд то источник находился в жидкости, окруженной несколькими упругими радиальными слоями. Однако эксперимент демонстрирует наличие как быстрой трубной волны, так и медленной. На основании экспериментальных данных, можно сделать вывод о том, что скорость поперечных волн в смеси песка и воды пренебрежимо мала. Поэтому 42 Д. А. Александров, А. В. Бакулин, Б. М. Каштан гравийная набивка моделируется как пористый слой с почти нулевым модулем сдвига скелета.

Гравийный фильтр и непроницаемый экран. На рис. 6 отражены зависимо сти скорости трубных волн и амплитудных коэффициентов от проницаемости песка.

Несмотря на тот факт, что обе трубы и экран, и обсадная колонна непроницаемы, наблюдается значительное влияние проницаемости песка на характер распространения волн, особенно на медленную трубную волну. Как показано на рис. 6, б и г, обе трубные волны имеют максимум поглощения в модели, в которой проницаемость песка равна 300 мД. В случаях, когда проницаемость песка мала или, напротив, высока, обе вол ны распространяются практически без затухания. Этому может быть дано следующее объяснение. Когда проницаемость песка низкая, жидкая и твердая фазы гравийного фильтра двигаются синхронно и ведут себя как эффективная жидкость из-за нулевого модуля сдвига.

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

Как видно на рис. 7, радиальные профили компонент смещения для модели с высоко проницаемым песком и эффективной жидкостью вместо него совпадают.

Распространение низкочастотных трубных волн в радиально-слоистых проницаемых средах Рис. 7. Зависимость компонент смещения от расстояния до оси скважины для модели с гравийным фильтром (а) и эффективной жидкостью вместо фильтра (б ) Рис. 8. Отношение вертикального смещения поровой жидкости (Vz ) к вер тикальной компоненте смещения скелета (Uz ) пористой среды, моделировавшей гравийный фильтр При промежуточных значениях проницаемости относительное движение жидкости в порах приводит к затуханию обеих трубных волн и особенно медленной волны. Вяз кие силы (рис. 8) удерживают жидкость в пористом веществе песка до тех пор, пока проницаемость не достигнет значения около 10 Д. Дальнейшее увеличение проницае мости ведет к значительному росту амплитуды колебаний жидкой фазы по сравнению с твердой фазой. При этом мы переходим в высокочастотный режим, где центральная частота сигнала (500 Гц) превышает частоту Био для данной модели (80 Гц).

При частоте колебаний выше критической частоты Био движение жидкости в порах перестает быть ламинарным [11].

44 Д. А. Александров, А. В. Бакулин, Б. М. Каштан Гравийный фильтр и проницаемый экран. В последней рассмотренной моде ли среда Био используется для моделирования и гравийного фильтра, и экрана. На рис. 9 представлены зависимости скоростей и амплитудных коэффициентов быстрой и медленной трубных волн в модели с гравийным фильтром и экраном, проницаемость которого составляет 100 Д. Варьируется проницаемость песка. При высоких и низких значениях проницаемости затухания быстрой волны не наблюдается. Низкая проница емость вещества, моделирующего песок гравийного фильтра, не позволяет жидкости двигаться относительно скелета пористого вещества. Поэтому, несмотря на высокую проницаемость экрана, сообщение между жидкостью в порах и жидкостью внутри экра на отсутствует. При увеличении проницаемости песка начинается обмен жидкостью, что приводит к значительному затуханию быстрой волны и полному поглощению мед ленной волны. Высокая проницаемость песка позволяет жидкости двигаться по порам, практически не встречая сопротивления, и быстрая волна перестает затухать. Поэтому в пределе высокой проницаемости систему песок гравийного фильтра экран можно рассматривать как один слой, близкий по своим параметрам к жидкому слою.

Рис. 9. Зависимость скорости и амплитудного коэффициента от проницаемо сти песка для быстрой (а, б ) и медленной (в, г) трубных волн в модели с гравий ным фильтром и проницаемым экраном (100 Д) Распространение низкочастотных трубных волн в радиально-слоистых проницаемых средах Сравнение численных расчетов с результатами эксперимента На рис. 10 представлены сейсмограммы, полученные в ходе эксперимента для трех моделей: а) скважина с непроницаемым экраном (сплошной алюминиевой трубой) и без гравийного фильтра;

б) скважина с проницаемым экраном (см. рис. 2) и без гравийного фильтра;

в) скважина с проницаемым экраном и гравийной набивкой. Как видно на сейсмограммах рис. 10, а и б, при замене сплошной алюминиевой трубы экраном с про ницаемостью 250–1000 Д амплитуда быстрых трубных волн практически не меняется, как и предсказано моделированием (см. рис. 4, б ). Кроме того, моделирование предпо лагает падение скорости быстрой трубной волны на 10% (см. рис. 4, а). В эксперименте скорость трубной волны также падает, однако изменение скорости составляет порядка 20% от значения в модели с непроницаемым экраном.

Рис. 10. Результаты измерений в эксперименте: а с непроницаемым экраном без гравийного фильтра;

б с проницаемым экраном без гравийного фильтра;

в с проницаемым экраном и гравийным фильтром В эксперименте возможность менять проницаемость гравийного фильтра не преду смотрена. Однако измерения проводятся в процессе установки гравийного фильтра.

Песок в пространство между экраном и обсадной колонной намывается с одного тор ца экспериментальной установки [9]. В процессе установки фильтра скважину можно условно разделить на три зоны: область без песка;

область смеси воды и песка;

область, где песок полностью заполнил межтрубное пространство (200–300 Д). Предполагая, что первые две зоны покрывают диапазон значений проницаемости от 200 Д до беско нечности, мы можем качественно сравнить результаты эксперимента с теоретическими расчетами. На рис. 11 представлены сейсмограммы, полученные в процессе установ ки гравийного фильтра. Прямоугольником на рисунке обозначен участок установки, где песок достиг максимальной высоты. Перед этим участком имеется область, лишь частично заполненная песком. На сейсмограмме отчетливо видно, что именно в этой области наблюдается уменьшение амплитуды волн. Эта амплитудная аномалия переме щается по мере заполнения скважины песком и исчезает, как только процесс установки фильтра окончен.

Таким образом, мы наблюдаем хорошее согласие между экспериментом и численны ми расчетами. Значения скоростей трубных волн, а также характер затухания для мо делей с гравийным фильтром и без него согласуются с численными расчетами (рис. 11, а и б ). Кроме того, между предельными значениями бесконечной проницаемости для во 46 Д. А. Александров, А. В. Бакулин, Б. М. Каштан Рис. 11. Экспериментальные сейсмограммы, полученные в процессе уста новки гравийного фильтра.

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

Основные расхождения между экспериментом и моделированием связаны прежде всего с поведением медленной трубной волны. Согласно численным расчетам, медлен ная мода полностью исчезает при значении проницаемости экрана более 0.1 Д (см.

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

рис. 9, г). Тем не менее медленная волна в эксперименте наблюдается в обеих моделях (см. рис. 10, б и г). Видно, что наибольшие различия между экспериментом и числен ными расчетами возникают при наличии проницаемого экрана. Таким образом, исполь зование модели пористой среды Био при моделировании экрана не позволяет корректно описать распространение медленной трубной волны в системе. Экран состоял из трех слоев: алюминиевой трубы с перфорациями диаметром 9.5 мм, направляющих, рас положенных с зазором 17 мм, и пластиковой обмотки с зазором 0.2 мм (см. рис. 2).

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

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

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

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

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

Работа выполнена при поддержке гранта АФГИР RUG-1-30005-ST- и ГК 02.740.11.0331.

Указатель литературы 1. Крауклис П. В., Крауклис Л. А. Волновое поле точечного источника в скважине // Во просы динамической теории распространения сейсмических волн. 1976. № 16. С. 41–53.

2. Крауклис П. В., Крауклис Л. А. Возбуждение трубной волны в скважине медленной вол ной, распространяющейся в жидком слое // Зап. науч. сем. ПОМИ. 1995. № 230. С. 115– 124.

48 Д. А. Александров, А. В. Бакулин, Б. М. Каштан 3. Крауклис П. В., Крауклис Л. А. Возбуждение трубной волны радиальным и вертикаль ным источниками, прикрепленными к стенке скважины // Зап. науч. сем. ПОМИ. 2007.

№ 297. С. 154–161.

4. Крауклис П. В., Крауклис Л. А. Трубная волна от точечного источника, находящегося вне скважины // Зап. науч. сем. ПОМИ. 2006. № 332. С. 99–122.

5. Крауклис П. В., Крауклис Л. А. Возбуждение трубной волны радиальным и вертикаль ным источниками, прикрепленными к стенке скважины // Зап. науч. сем. ПОМИ. 2007.

№ 342. С. 153–163.

6. Молотков Л. А. Матричный метод в теории распространения волн в слоистых упругих и жидких средах. Л.: Наука, 1984. 270 c.

7. Уайт Дж. Э. Возбуждение и распространение сейсмических волн / пер. с англ.

О. В. Павловой и С. В. Гольдина. М.: Недра, 1986. 261 c.

8. Bakulin A., Sidorov A., Kashtan B., Jaaskelainen M. Real-time completion monitoring with acoustic waves // Geophysics. 2008. Vol. 73. E15–E33.

9. Bakulin A., Sidorov A., Kashtan B., Jaaskelainen M. Downhole acoustic surveillance of deep water wells // The Leading Edge. 2008. Vol. 27. P. 518–531.

10. Bakulin A., Alexandrov D., Sidorov A., Kashtan B. Acoustic waves in sand-screened deep water completions: Comparison of experiments and modeling // Geophysics. 2009. Vol. 74.

P. 45–56.

11. Biot M. A. Propagation of elastic waves in cylindrical bore containing a uid // Journal of the Acoustical Society of America. 1952. Vol. 23(9). P. 997–1005.

12. Biot M. A. Theory of propagation of elasticwaves in a uid-saturated porous solid // Journal of the Acoustical Society of America. 1956. Vol. 28. P. 168–191.

13. Kirkendall C. K., Dandridge A. Overview of high performance bre-optic sensing // Journal of Physics D: Applied Physics. 2004. Vol. 37. R197–R216.

14. Lamb H. On the velocity of sound in a tube as aected by the elasticity at the walls // Manchester Memories. 1898. Vol. 42. P. 1–16.

15. Sidorov A., Bakulin A., Kashtan B., Ziatdinov S., Alexandrov D. Low-frequency symmetric waves in uid-lled boreholes and pipes with radial layering // Geophysical Prospecting.

2008. Vol. 57. P. 863–882.

16. Tang X. M., Cheng A. Quantitative borehole acoustic methods. Elseiver, 2004.

17. Wong G. K., Fair P. S., Bland K. F., Sherwood R. S. Balancing act: Gulf of Mexico sand control completions, peak rate versus risk of sand control failure // SPE. 2003. P. 844–850.

Вопросы геофизики. Выпуск 44. СПб., 2011 (Ученые записки СПбГУ;

№ 444) Д. В. Аникиев, Б. М. Каштан, А. С. Благовещенский, В. А. Мулдер ТОЧНЫЙ ДИНАМИЧЕСКИЙ МЕТОД РЕШЕНИЯ ОБРАТНОЙ ЗАДАЧИ СЕЙСМИКИ НА ОСНОВЕ ИНТЕГРАЛЬНЫХ УРАВНЕНИЙ ГЕЛЬФАНДА ЛЕВИТАНА Введение При определенных условиях постановка обратной задачи сейсмики в качестве дина мической обратной задачи математической физики становится математически коррект ной. Термин динамическая означает, что речь идет о задаче, в которой используются непосредственно значения волновых полей. Это отличает такую задачу от кинематиче ской обратной задачи, в которой используются лишь времена пробега волн через изуча емую среду. Динамическое рассмотрение обладает рядом преимуществ по сравнению с кинематическим, например, отсутствует необходимость выделения (пикирования) вре мен пробега, удаления многократных отражений.

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

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

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

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

c Д. В. Аникиев, Б. М. Каштан, А. С. Благовещенский, В. А. Мулдер, 50 Д. В. Аникиев, Б. М. Каштан, А. С. Благовещенский, В. А. Мулдер 1. Решение обратной динамической задачи сейсмики на основе линейных интегральных уравнений Гельфанда Левитана и преобразования Радона Математической предпосылкой развития теории динамических обратных задач по служила работа И. М. Гельфанда и Б. М. Левитана [2]. В работе решение обратной за дачи сводится к решению линейного интегрального уравнения, названного впослед ствии уравнением Гельфанда Левитана. Важные результаты были также получены М. Г. Крейном [3] и В. А. Марченко [4]. После данных работ появилось большое количе ство литературы, посвященной использованию их результатов в различных уравнениях и системах уравнений.

Обратная задача для систем уравнений теории упругости может быть сформули рована как в спектральной области (аналогично теории квантового рассеяния), так и во временнй. Иногда термин динамическая обратная задача относят именно к вре о меннй постановке (см., например, [5]). Решение обратной динамической задачи тео о рии упругости (так называемой обратной задачи Лэмба) в спектральной области было предложено А. С. Алексеевым [6] в предположении, что плотность среды и упругие па раметры Ламе зависят только от одной пространственной переменной (вертикальной).

А. С. Благовещенским [7] было впервые рассмотрено решение динамической обратной задачи во временнй области, а также был предложен метод нелинейных вольтерров о ских систем интегральных уравнений, являющийся альтернативным по отношению к методам линейных интегральных уравнений. B. Гопинат и М. Сондхи [8, 9] независи мо предложили альтернативное интегральное уравнение (также во временнй области) о в задаче восстановления формы речевого тракта человека по акустическим измере ниям. Нелинейные интегральные уравнения во временной области также использовал В. Саймс [10]. Р. Барридж [11] исследовал применение уравнений Гельфанда Левита на и Марченко для теории упругости во временнй области и установил связь между о ними и уравнением Гопината Сондхи [9]. Подробный обзор литературы по различным точным методам решения динамической обратной задачи сейсмики для горизонтально однородной среды представлен Р. Ньютоном [12].

В большинстве работ по точным методам решения динамической обратной зада чи сейсмики строение среды предполагается горизонтально-однородным. Тем не менее существует возможность теоретически усовершенствовать модель восстанавливаемой среды. Например, в работе [13] предлагается численная схема решения обратной дина мической задачи для среды со слабой горизонтальной неоднородностью. Методы реше ния динамической обратной задачи для произвольной неоднородной среды (например, с помощью BC-метода [14]) сопряжены со значительными дополнительными требова ниями к входным данным, которые затрудняют практическое применение таких ме тодов.

Р. Кэррион [15] показал, что трехмерная обратная задача сейсмики для горизон тально-однородной среды сводится к одномерной за счет преобразования Радона [16] и решается для определенной пары параметров (, p). Кэррион [15] также представил схему численного решения обратной задачи сейсмики на основе метода характеристик во временнй области. Он показал, что использование метода характеристик позво о ляет точно восстанавливать структуру среды с кусочно-непрерывным распределением параметров.

Точный динамический метод решения обратной задачи сейсмики... Подход, рассматриваемый в данной статье, близок к работе Р. Кэрриона [15]. Оба подхода формулируются во временнй области, кроме того, для сведения многомерной о обратной задачи к одномерной в обоих случаях используется преобразование Радона.

Однако для нахождения решения одномерной задачи в настоящей статье, в отличие от работы Кэрриона, вместо метода характеристик используется метод линейных ин тегральных уравнений типа Гельфанда Левитана. Таким образом, предлагаемая в данной работе схема численной реализации является альтернативной по отношению к схеме, рассмотренной в работе [15].

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

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

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

V = V (z), = (z).

Теоретическое описание, аналогичное приведенному выше, возможно и для уравне ния упругости (см. [1]), однако в статье упругая среда не рассматривается. Для исследо В отечественной литературе (см., например, [1]) встречается термин уравнение Гельфанда Ле витана Крейна Марченко по именам отечественных математиков, исследовавших уравнения подоб ного вида в период становления теории обратных динамических задач.

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

1.2. Постановка задачи с точки зрения математической физики Математическая формулировка обратной задачи излагается здесь для трехмерного акустического случая, когда волновое поле зависит от трех пространственных пере менных. Таким образом, рассматривается уравнение акустики для трехмерного случая при x;

y +;

0 z +;

t +:

1 1 Utt = Uz + (Uxx + Uyy ). (1) V 2 z Здесь U (x, y, z, t) поле давления в зависимости от времени и пространственных координат;

(z) плотность и V (z) акустическая скорость, зависящие только от вер тикальной координаты z, так как среда предполагается горизонтально-однородной.

Математически необходимо сформулировать начальное условие, которое означает отсутствие каких-либо возмущений в среде до момента возбуждения источника:

U |t0 = 0.

Граничные условия:

U |z=0 = H(t)(x)(y), (2) Uz |z=0 = F (x, y, t). (3) Первое граничное условие (2) выполняет роль источника, (x), (y) дельта-функ ции Дирака. В силу своих свойств произведение дельта-функций в условии (2) опреде ляет местоположение источника в точке (0, 0). Функция H(t) в условии (2) определяет временню зависимость поля в точке источника. Необходимо отметить, что источник в у данной конкретной постановке задается не в виде внешней силы в правой части уравне ния (1), a в виде граничного условия (поле U при z = 0). Второе граничное условие (3) определяет характер регистрируемых данных (вертикальная пространственная произ водная поля давления). Функция F (x, y, t) в условии (3) поле, которое регистрируется приемниками, исходные данные обратной задачи, на основе которых строится ее реше ние.

В связи с предполагаемой горизонтально-однородной структурой среды обратная задача обладает цилиндрической симметрией [1]. Это означает, что зависимость функ ции F (x, y, t) от переменных x и y можно заменить зависимостью от одной радиальной координаты r. Вводя цилиндрические координаты r,, z, второе граничное условие запишем в виде Uz |z=0 = F (r, t).

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

Точный динамический метод решения обратной задачи сейсмики... Для удобства записи введем также полярный вектор r на горизонтальной плоскости с координатами (x;

y):

(x;

y) R2, r = (4) x2 + y 2 = r.

|r| = Таким образом, для построения решения обратной задачи сейсмики в случае гори зонтально-однородной среды достаточно зарегистрировать вертикальную производную давления F (r, t) вдоль некоторого горизонтального направления (заданного вектором r). В силу цилиндрической симметрии задачи такое направление произвольно.

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

В статье используется преобразование Радона, также возможно применение простран ственного преобразования Фурье. Эти преобразования и характерные особенности их применения в теории обратных динамических задач подробно обсуждаются в моно графии [1]. В предыдущих работах по теории динамических обратных задач [17–19] использовалось пространственное преобразование Фурье, тогда как в работе [20] пре образование Радона. Метод на основе преобразования Радона имеет ряд преимуществ с вычислительной точки зрения и приводит к более точным численным результатам решения обратной задачи.

Преобразование Радона широко используется в разнообразных научных исследова ния, в том числе в сейсморазведке (см. [21, 22]). При численном моделировании или обработке сейсмических данных необходим дискретный аналог преобразования Радо на [23]. В прикладной сейсмике вместо термина дискретное преобразование Радона часто применяется термин наклонное суммирование (slant stacking [24, 25]). Не ме нее распространен термин p-преобразование ( p transform [26, 27]). Дискретное преобразование Радона в отечественной литературе известно как метод регулируемо го направленного приема (РНП), который был введен независимо Ф. Рибером [28] и Л. А. Рябинкиным [29, 30].

Классическое преобразование Радона [16, 31] в трехмерном случае сопоставляет функции U (r, z, t) (r полярный вектор, введенный в (4)) функцию U(z,, p) ( и p параметры преобразования Радона) при фиксированном z:

U (z,, p) = U (r, z, + (p, r)) dxdy.

R Выражение (p, r) представляет собой скалярное произведение векторов p и r на плоскости x, y:

(p, r) = px x + py y.

54 Д. В. Аникиев, Б. М. Каштан, А. С. Благовещенский, В. А. Мулдер Тогда U(z,, p) = U (r, z, + px x + py y) dxdy. (5) R Исходя из вида подынтегрального выражения (5), величины px и py определяют обратные кажущиеся скорости распространения волн (кажущиеся медленности) вдоль осей x и y соответственно. Таким образом, вектор p c координатами (px ;

py ) имеет смысл горизонтальной компоненты вектора медленности (вектора рефракции). Максимальное значение параметра p определяется выражением pmax =, Vmax где Vmax максимальная скорость распространения волн в среде. Условие для выбора лучевого параметра p:

0 p pmax =. (6) Vmax В силу цилиндрической симметрии задачи можно полагать, что система координат выбрана по отношению к лучевому параметру удобным способом:

px = p, py = 0. (7) Это позволяет записать выражение (5) в виде U (z,, p) = rdr d U (r, z, + pr cos ).

0 Применение преобразования Радона к уравнению (1) позволяет получить следу ющее уравнение для функции U (преобразования Радона функции U ):

1 p2 V 2 U = Uz. (8) V 2 z Начальное условие примет вид U | 0 = 0.

Граничное условие U |z=0 = H( ).

Исходные данные обратной задачи:

Uz |z=0 = F (, p), где F преобразование Радона функции F.

Точный динамический метод решения обратной задачи сейсмики... Итак, уравнение (8) равнозначно одномерному уравнению акустики 1 U= Uz, (9) 2 tt z V где в роли скорости распространения волн выступает V V=, (10) 1 p2 V в роли пространственной переменной z, в роли временнй t (вместо ). Параметр о p достаточно мал для того, чтобы выполнялось условие (6), т. е. чтобы скорость V в выражении (10) имела смысл.

Рассмотрим решение получившейся одномерной задачи для уравнения (9), основан ное на интегральных уравнениях.

1.4. Вывод линейных интегральных уравнений типа Гельфанда Левитана Сведение трехмерной задачи (1) к одномерной задаче (9) позволяет использовать для нахождения параметров среды линейные интегральные уравнения типа Гель фанда Левитана. Одномерная обратная задача для уравнения (9) должна быть сфор мулирована для определенного значения лучевого параметра p (хотя бы для одного).

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

Для дальнейших выкладок удобно произвести замену пространственной переменной z на переменную с размерностью времени следующим образом:

z cos = (z, p) = dz, 0 +.

V (z ) Подынтегральное выражение здесь представляет собой вертикальную кажущуюся медленность, тем самым имеет смысл наклонного времени пробега волны [32] в среде со скоростью V на глубину z в направлении луча, определяемого лучевым пара метром p. Выражение для проще записать в виде z 1 p2 V 2 (z ) = (z, p) = dz.

V (z ) Введем дополнительную переменную, зависящую от параметра p,, V и от z (за счет скорости V (z) и плотности (z)):

1 p2 V (z) (p,, V, z) =. (11) (z)V (z) При p = 0 величина :

(0,, V, z) =, (z)V (z) т. е. равна величине, обратной акустическому импедансу V.

56 Д. В. Аникиев, Б. М. Каштан, А. С. Благовещенский, В. А. Мулдер Выражение для через, :

z (z )(z )dz.

= (12) Заметим, что z (z )(z )dz = (z)(z).

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

= =. (13) z z y Отметим, что значение z = 0 соответствует значению = 0. Для обозначения, V и на границе z = 0 ( = 0) введем нижний индекс 0:

0 = (0), V0 = V (0), 1 p2 V 0 =. (14) 0 V С учетом введенных обозначений и видоизмененного оператора (13) уравнение (9) имеет вид Utt U = 0. (15) Граничное условие:

U (, t)|=0 = H(t).

Исходные данные обратной задачи:

1 p U |=0 = F (t) = F (t).

V02 0 Введем 1 p = =, V02 0 а также F = F (t). (16) Точный динамический метод решения обратной задачи сейсмики... Тогда данные обратной задачи запишем в виде U |=0 = F (t).

Решение одномерной задачи для уравнения (15) далее описывается для фиксирован ного значения лучевого параметра p. Это позволяет опустить зависимость функции F от p для краткости и простоты изложения, однако стоит помнить, что различные значе ния p дают различные функции (p,, V, ). Разные значения p соответствуют разным зависимостям (z, p). Связь переменных 1 и 2 для значений p1 и p2 определяется выражением 1 d d =.

1 ( ) 2 ( ) 0 Введем вспомогательную функцию (, t), являющуюся решением уравнения (15):

tt ( ) = 0. (17) Пусть также удовлетворяет следующим граничным условиям ( = 0):

|=0 = (t), (18) |=0 = (t). (19) Здесь (t) функция Хевисайда, 1, t 0, (t) = (20) 0, t 0.

Свойства (, t) (см., например [1]):

|t = 0, |t = 1, |t = 1 ( t) + B( t)+ + C( + t)+. (21) Введем следующие обозначения:

± t, ± t 0, ( ± t)+ = 0, ± t 0.

1 0 ( ) d B =, 2 2 = C =.

2 = 58 Д. В. Аникиев, Б. М. Каштан, А. С. Благовещенский, В. А. Мулдер Если функция () постоянна в окрестности = 0 (на границе), то 1 0 ( ) d, B = 2 C = 0.

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

Наряду с решением (17) вводим дополнительно функцию :

(, t) = 1 (, t). (22) Пара функций V = ( ), (23) W = ( + ) (24) 2 t образуют фундаментальную систему решений задачи Коши для уравнения (17) со сле дующими условиями:

V |=0 = 0, V |=0 = (t), W |=0 = (t), W |=0 = 0.

С помощью этих функций любое решение уравнения, в том числе и решение, соот ветствующее физической реализации, имеет вид U (, t) = W H(t) + V F (t). (25) Символ обозначает свертку функций по времени. В последнем выражении использу ется функция F (t), которая была введена в формуле (16). Функция H(t) представляет собой произвольную гладкую временную функцию в источнике.

В силу того что U ||t| = 0, справедливо соотношение W H(t) + V F (t) = 0, |t|, или W (, )H(t )d + V (, )F (t )d = 0.

Точный динамический метод решения обратной задачи сейсмики... Подставляя выражения для V (23) и W (24), перепишем уравнение в виде d (, s) + (, s) H(t s)ds + dt + (, s) (, s) F (t s)ds = 0.

Заметим, что и (, s), и (, s) = 0 при s. Кроме того, H(t s) = F (t s) = 0, s t.

Поскольку |t|, то последнее равенство выполняется также при s. Теперь изменим пределы интегрирования с учетом носителей входящих в уравнение функций:


d (, s) + (, s) H(t s)ds + (, s) (, s) F (t s)ds = 0.

dt Заменяем обратно на по формуле (22):

d (H(t s) (, s)H(t s) + (, s)H(t s)) ds + dt + F (t s) (, s)F (t s) (, s)F (t s) ds = 0.

Делаем замену (s) (s) в интегралах, содержащих (, s):

d d H(t s)ds + (, s) (H(t s) H(t + s)) ds + dt dt + F (t s)ds (, s) F (t s) + F (t + s) ds = 0.

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

t+ t+ d F (s )ds + H(s )ds + (, s) H(t s) H(t + s) ds dt t t (, s) F (t s) + F (t + s) ds = 0.

60 Д. В. Аникиев, Б. М. Каштан, А. С. Благовещенский, В. А. Мулдер Произведем дифференцирование в первом слагаемом, третье и четвертое слагаемые перенесем в правую часть:

t+ F (s )ds = H(t + ) H(t ) + (, s) F (t s) + F (t + s) ds t (, s) H(t s) H(t + s) ds.

Поскольку |t|, то H(t ) = 0, кроме того F (s)|s0 = 0, тем самым выражение упрощается:

t+ H(t + ) + F (s)ds = = (, s)[H(t + s) H(t s) + F (t s) + F (t + s)]ds. (26) Последнее выражение справедливо при 0 +, t и является линей ным интегральным уравнением (для функции ) типа Гельфанда Левитана, т. е. явля ется близким по виду к уравнению Гельфанда Левитана, выведенному в монографии [1]. Уравнение (26) относится к интегральным уравнениям Фредгольма I рода. Числен ное решение такого уравнения относительно функции (, t) не является корректной математической задачей, и ее решение требует регуляризации (см. монографию [33]).

Чтобы добиться корректности, необходимо наличие во временнй функции источника о сингулярной составляющей (см. [1]).

Так, можно рассмотреть частный случай, когда времення функция H(t) являет а ся суммой обобщенной функции Хевисайда (t) [см. уравнение (20)] и произвольной гладкой функции h(t):

H(t) = (t) + h(t).

Отклик среды F (t), зарегистрированный на поверхности в силу условий (18), (19) и формулы (25), будет алгебраической суммой дельта-функции (t) и гладкой функ ции f (t):

F (t) = (t) + f (t). (27) Аналогично предыдущему случаю с помощью функций V и W, введенных ранее, можно записать любое решение уравнения, в том числе и физическое:

U (, t) = W H(t) + V F (t).

Однако теперь, с учетом вида H(t) и F (t):

U (, t) = W (, t)(t) + W (, t)h(t) V (, t)(t) + V (, t)f (t) = (t) V (, t).

= W (, t)(t) + W (, t)h(t) + V (, t)f Точный динамический метод решения обратной задачи сейсмики... Подставляя выражения для V (23) и W (24), получим:

d d (, s) + (, s) (t s)ds + (, s) + (, s) h(t s)ds + dt dt + (, s) (, s) f (t s)ds (, t) + (, t) = 0.

Итак, заменим обратно на, как в формуле (22), и изменим пределы интегриро вания с учетом носителей входящих в выражение функций:

t d 1 (, s) + (, s) + (1 (, s) + (, s)) h(t s)ds 1+ dt t + (, s) + (, s) + (1 (, s) (, s)) f (t s)ds = 0.

Проведем замену (s) (s) в интегралах, содержащих (, s):

t+ t 2(, s) + h(s)ds (, s)h(t + s)ds + (, s)h(t s)ds + t t+ t + f (s)ds (, s)f (t + s)ds (, s)f (t s)ds = 0.

t Перепишем интегральное уравнение:

t+ h(t + ) + f (s)ds = = (, s)[h(t + s) h(t s) + f (t + s) + f (t s)]ds 2(, t). (28) Таким образом, за счет дополнительного предположения о виде временнй зависи о мости в источнике получено уравнение (28), справедливое при 0 +, t.

Уравнение (28) относится к интегральным уравнениям Фредгольма II рода. Численное решение таких уравнений является корректной задачей (cм. [33]).

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

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

а 62 Д. В. Аникиев, Б. М. Каштан, А. С. Благовещенский, В. А. Мулдер Таким образом, представлен способ нахождения решения одномерной задачи (9) в виде функции (, t), которая связана с параметрами с помощью формулы (21).

1.5. Восстановление акустической скорости и плотности из решения интегрального уравнения Итак, функция (, t) получается решением интегрального уравнения (26) или (28).

Такое решение в общем виде нельзя выписать аналитически, однако можно решить численно. Функция (, t) уже известна. Покажем, как из известной (, t) находятся параметры среды.

С учетом формулы (21) для справедливо следующее предельное соотношение:

lim (, t) = (, ) = 1.

() t Отсюда () = 2. (29) (1 (, )) Заметим, что 0 (значение на границе, см. формулу (14)) единственная априор ная информация о реальной среде, необходимая для построения решения.

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

V () =.

2 () + p Связь z с следует из выражения (12) при 1:

d z() =. (30) () В более сложном общем случае, когда плотность зависит от глубины z, можно найти лишь зависимость функции от [см. формулу (29)], т. е. плотность и скорость в совокупности. Для независимого восстановления плотности и скорости V от глуби ны z необходимо решить задачу как минимум для двух различных значений лучевого параметра, p1 и p2. Действительно, рассмотрим два независимых решения, 1 (1 ) и 2 (2 ):

1 (1 ) = 2, (1 1 (1, 1 )) 2 (2 ) = 2.

(1 2 (2, 2 )) Кроме того, 1 p2 V 2 (1 ) 1 (1 ) =, (1 )V (1 ) Точный динамический метод решения обратной задачи сейсмики... 1 p2 V 2 (2 ) 2 (2 ) =.

(2 )V (2 ) Рассмотрим теперь выражение 1 p2 V 2 (1 ) 1 p2 V 2 (1 ) p2 p 2 22 = 22 2 1 (1 ) 2 (1 ) =.

2 (1 )V 2 (1 ) (1 )V (1 ) (1 ) После несложных преобразований p2 p 2 (1 ) =, 2 ( ) 1 1 2 (1 ) 2 1 (1 ) 2 (1 ) V (1 ) = 2 2 ( ).

2 2 ( ) p2 1 1 p1 2 Зависимость 2 (1 ) можно найти, подставив 2 (1 ) в известную 2 (2 ). Связь z c в таком случае имеет вид 1 d 1 2 ( ) 2 2 ( ) z(1 ) = = d.

1 ( )(p2 p2 ) )( ) ( 2 0 Таким образом, представлены формулы, позволяющие восстановить независимо друг от друга акустическую скорость V и плотность в зависимости от глубины z на основе решения (, t) уравнения (26) или его частного случая (28).

2. Численная реализация 2.1. Модель среды На основе предложенного метода решения обратной динамической задачи были со зданы алгоритмы численной реализации для одномерного, двумерного и трехмерного случаев и написан программный код. Приведем краткое описание этих алгоритмов, а также численного моделирования и его результатов.

В качестве модели горизонтально-однородной среды при численном моделирова нии использовался вертикальный профиль известной синтетической модели Мармузи (Marmousi model, см. [34, 35]). Распределение акустической скорости в модели Мармузи показано на рис. 2. Модель использовалась во многих работах, связанных с численны ми экспериментами по решению обратных задач (см., например, [35–38]) и является эталонной моделью для проверки методов решения обратной задачи. Распределение акустической скорости в левой части модели Мармузи близко к горизонтально-одно родному, поэтому в качестве модели для численного эксперимента взят вертикальный профиль для горизонтальной координаты x = 1000 м (рис. 3, а).

Распределение плотности в модели Мармузи также близко к горизонтально-одно родному в окрестности x = 1000 м (рис. 3, б ). Профиль до глубины 40 м представляет собой водный слой, акустическая скорость и плотность в котором однородно распреде лены и равны соответственно 1500 м/с и 1000 кг/м3.

64 Д. В. Аникиев, Б. М. Каштан, А. С. Благовещенский, В. А. Мулдер Рис. 2. Распределение акустической скорости в модели Мармузи.

Вертикальный профиль изображен штриховой линией Рис. 3. Вертикальный профиль акустической скорости (а) и плотности (б ) в модели Мармузи для горизонтальной координаты x = 1000 м При численном моделировании была использована исходная пространственная дис кретизация модели Мармузи, составляющая 4 м как по горизонтали, так и вертикали:

dx = 4 м, dz = 4 м. (31) 2.2. Алгоритм численного решения обратной динамической задачи в одномерном случае Проверку алгоритма численной реализации предложенного метода решения обрат ной задачи удобно провести для одномерного пространственного случая, поскольку ис ключаются погрешности, связанные с неточностью численного преобразования Радона.

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

В случае одного измерения, когда z является единственной пространственной коор динатой, исходное уравнение акустики выглядит следующим образом:

1 Utt = Uz. (32) V 2 z Оно эквивалентно уравнению (9) при p = 0. Таким образом, решение обратной задачи для уравнения (32) аналогично решению для уравнения (9).

Ранее было показано, что для конкретного фиксированного значения параметра p можно восстановить только параметр (11). В одномерном случае, когда p = 0, величи на обратно пропорциональна акустическому импедансу (произведению акустической скорости и плотности):

(, V, z) =.

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


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

Высказанные соображения обосновывают предлагаемую численную реализацию ме тода решения одномерного варианта задачи для постоянной плотности, которую луч ше принять равной 1 г/см3 (или 1000 кг/м).

В этом случае уравнение акустики (32) преобразовывается в обычное волновое урав нение Utt = V 2 (z)Uzz с начальным условием U |t0 = 0.

Для проверки работоспособности метода в качестве источника удобно использовать функцию Хевисайда:

U (z, t)|z=0 = (t). (33) Исходные данные обратной задачи:

Uz (z, t)|z=0 = F (t). (34) Моделирование волнового поля давления (моделирование прямой задачи) планиро валось осуществить с помощью одномерного численного псевдоспектрального метода PSTD method Pseudo Spectral Time Domain method.

66 Д. В. Аникиев, Б. М. Каштан, А. С. Благовещенский, В. А. Мулдер Рис. 4. Сигнал Гаусса с центральной частотой 100 Гц во временнй области (далее метод PSTD) [39–41]. Однако моделирование прямой за о дачи для источника в виде функции Хэвисайда дает недостаточно точный результат.

В связи с этим было рассмотрено решение другой задачи (более точное с вычислитель ной точки зрения):

Utt = V 2 (z)Uzz + (t)(z) при z, 0 t, c начальными условиями U |t=0 = 0, Ut |t=0 = 0.

Решение такой задачи эквивалентно исходной задаче с точностью до коэффициен та 1/2 [42, 43]. Подобная эквивалентность задач позволяет при численной реализации в качестве источника использовать сигнал в виде дельта-функции Дирака. В каче стве модели дельта-функции был использован сигнал Гаусса с центральной частотой 100 Гц (рис. 4). Центральная частота определяет длительность импульса. Исходя из ми нимальной скорости в профиле на рис. 3, равной 1500 м/с, минимальная длина волны для импульса с центральной частотой 100 Гц 1500 м/с = = 15 м.

100 Гц Таким образом, на минимальной длине волны укладывается три узла (расстояние между узлами 4 м (31)), что является достаточным условием для метода PSTD [40].

Шаг по времени dt = 0.46 мс при моделировании был выбран из соображений стабиль ности одномерного метода PSTD.

Точный динамический метод решения обратной задачи сейсмики... В качестве скоростной модели среды V (z) использовался профиль акустической ско рости, изображенный на рис. 3, в интервале глубин от 0 до 1100 м с дискретизацией 4 м.

Значение акустической скорости до глубины 40 м (водный слой) при решении обратной задачи считается известным: V0 = 1500 м/с. В данном случае это единственная апри орная информация о среде, необходимая для решения обратной задачи предлагаемым методом.

Входные данные для решения обратной задачи это функция F (t) (3), вертикаль ная производная (производная по координате z) давления на границе. Поскольку ре зультат моделирования прямой одномерной задачи методом PSTD представляет собой поле давления в каждой пространственной точке, то для вычисления вертикальной пространственной производной давления в граничной точке достаточно сосчитать раз ность давлений для двух близких к точке границе узлов (и поделить на расстояние между узлами). Смоделированная таким образом вертикальная пространственная про изводная давления в граничной точке среды в зависимости от времени (функция F (t)) изображена на рис. 5.

Рис. 5. Вертикальная пространственная производная поля давле ния на границе в зависимости от времени При сформулированном ранее условии (33) справедливо интегральное уравнение (28), которое принимает вид t+ f (s)ds = (, s)[f (t + s) + f (t s)]ds 2(, t). (35) Уравнение (35) является уравнением Фредгольма II рода и не требует регуляриза ции. Функция f при p = 0 и 0 = 1 связана с входными данными [cм. формулы (16), (27)] следующим образом:

1 F (t) = V0 · F (t) = F (t) = (t) + f (t). (36) 68 Д. В. Аникиев, Б. М. Каштан, А. С. Благовещенский, В. А. Мулдер Таким образом, f (t) это входные данные F (t) за вычетом импульса, соответству ющего прямой волне от источника до приемника (т. е. входные данные за вычетом пер вого импульса), дополнительно домноженные на величину скорости в первом слое V0.

Заметим, что в случае источника с произвольной гладкой временнй функцией (не им о пульсного) вид уравнения (35) изменяется, и для нахождения его решения необходима регуляризация.

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

Дискретные аналоги временных переменных t и с учетом |t| :

j = j dt, j = 1, n, ti = i dt, i = j, j.

Дискретная функция f (t) представляет собой набор значений fk = f (tk ), а функция (, t):

ji = (j, ti ).

Уравнение (26), как и его частные случаи уравнения (28), (35), сформулировано для определенных, т. е. эти выражения описывают системы линейных интегральных уравнений относительно функции (, t) для различных значений.

Итак, уравнение (35) в дискретном виде (интегралы заменяются интегральными суммами) запишем следующим образом:

i+j j dt · fk = dt · jk [fk+j + fkj ] 2ji, j = 1, n, i = j, j. (37) k=1 k=j Введем единичную матрицу Iik :

Iik = 1, k = i, Iik = 0, k = i.

Тогда уравнение (37) примет вид i+j j fk = [fk+j + fkj Iik ]jk, j = 1, n, i = j, j.

dt k=1 k=j Индекс j можно считать фиксированным, так как он определяет значение, для каждого из которых записывается уравнение (35). Введем обозначения без индекса j:

k = jk, i+j Ri = fk, k= Точный динамический метод решения обратной задачи сейсмики... Kik = fk+j + fkj Iik. (38) dt Тогда дискретное уравнение (37) запишем в виде j Ri = Kik k, j = 1, n, i = j, j k=j или в матричном виде:

R = K.

Тогда = K 1 R.

Вид матрицы K [см. выражение (38)] предполагает существование обратной мат рицы K 1, для нахождения которой можно использовать, например, метод Гаусса Жордана (метод полного исключения неизвестных). При численном решении обрат ной задачи для обращения матрицы K использовалась численная реализация метода Гаусса Жордана [44].

Обращение позволяет получить дискретные значения ji функции (, t) для каж дого значения j, при котором 1 j n. Далее для нахождения дискретных значений j функции () необходимо использовать jj (или просто j ) дискретные значения (, ):

j = 2.

(1 j ) В конечном итоге для дискретных значений Vj скорости V ():

Vj = V0 (1 j ).

Для получения зависимости скорости V от глубины z находится зависимость z() [см. уравнение (30)] в дискретном виде zj :

j j (1 k )2.

zj = z(j ) = dt · Vk = V0 dt · k=1 k= Таким образом, численным решением рассматриваемой обратной задачи являются значения Vj и zj для j = 1, n, где каждому Vj соответствует свой zj.

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

Рассмотренный одномерный случай удобен для описания алгоритма, наглядной де монстрации и проверки работоспособности, однако не соответствует реальному сейсми ческому исследованию. Далее рассмотрена численная реализация алгоритма в случае 70 Д. В. Аникиев, Б. М. Каштан, А. С. Благовещенский, В. А. Мулдер Рис. 6. Сравнение результата численного решения одномерной об ратной задачи с исходной моделью среды двух измерений. Двумерный вариант обратной задачи затрагивает все аспекты реаль ного сейсмического исследования и удобен с точки зрения моделирования. Двумерный и трехмерный алгоритмы различаются лишь методом преобразования Радона.

2.3. Численная реализация в случае двумерного уравнения акустики Рассмотрим уравнение акустики для двух пространственных координат:

1 1 Utt = Uz + Uxx. (39) V 2 z Первое граничное условие U |z=0 = (t)(x), где (t) функция Хевисайда (20).

Второе граничное условие Uz |z=0 = F (x, t), где F (x, t) зарегистрированные на поверхности среды значения вертикальной про странственной производной давления в зависимости от положения на поверхности и времени. Эти данные являются исходными для решения обратной задачи.

Принципиальная схема расстановки приемника и источников изображена на рис. 7.

Прямая задача для такой схемы решалась численно с помощью двумерного варианта численного метода PSTD. При моделировании прямой задачи использовались следу ющие параметры: пространственная дискретизация dx = dz = 4 м, шаг по времени (исходя из стабильности метода PSTD) dt = 0.5 мс, минимальный вынос (расстояние от источника до ближайшего приемника) rmin =8 м, максимальный вынос rmax =3000 м, расстояние между приемниками dr = 8 м, количество приемников n = 375.

В качестве скоростной модели среды V (z) использовался тот же профиль акустиче ской скорости, что и для одномерного случая (см. рис. 3, а) в интервале глубин от 0 до Точный динамический метод решения обратной задачи сейсмики... Рис. 7. Принципиальная схема расстановки приемника и источников для численного моделиро вания прямой двумерной задачи: rmin минималь ный вынос;

dr расстояние между приемниками 1100 м. Как и прежде, значение акустической скорости в водном слое (до глубины 40 м) при решении обратной задачи считается известным: V0 = 1500 м/с. В отличие от одно мерного случая, для двумерной обратной задачи независимое восстановление скорости и плотности принципиально возможно, и далее, в качестве подтверждения, приводится пример такого численного результата решения обратной задачи. Остальные численные результаты восстановления приводятся для одного параметра акустической скорости.

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

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

Двумерная обратная задача для уравнения (39) сводится с помощью двумерного преобразования Радона к одномерной задаче для уравнения (8). В свою очередь, для численного решения получившейся одномерной задачи используется аналогичный рас смотренному ранее алгоритм.

Для преобразования Радона (как двумерного, так и трехмерного) существует мно жество численных методов. Наиболее рациональным из них является метод быстрого преобразования Радона. Дискретная формулировка этого метода и характерные особен ности подробно изложены в работе [45]. В рамках данной работы метод использовался при вычислении двумерного преобразования Радона.

Результат применения метода быстрого дискретного преобразования Радона [45] к входным данным для ряда значений лучевого параметра p показан на рис. 9. Диапазон значений лучевого параметра выбран от 0 до 0.667 мс/м, т. е. до максимального зна чения p, представляющего интерес при морских сейсмических наблюдениях (см. [25]).

Точное максимальное значение p ограничено условием (7), т. е. оценить его можно по максимальной скорости распространения волн.

Заметим, что решение обратной задачи по предложенному алгоритму может осу ществляться для каждого лучевого параметра p. Для различных p результаты восста новления параметров должны быть одинаковы. Тем не менее в силу различных ошибок численного характера и погрешностей, связанных с дискретным преобразованием Ра 72 Д. В. Аникиев, Б. М. Каштан, А. С. Благовещенский, В. А. Мулдер Рис. 8. Поле вертикальной пространственной производной давле ния на поверхности среды в зависимости от времени Рис. 9. Двумерное преобразование Радона входных данных Точный динамический метод решения обратной задачи сейсмики... Рис. 10. Сравнение результата численного решения двумерной об ратной задачи для p = 0.001 мс/м (а) и p = 0.3 мс/м (б ) с исходной моделью среды дона, результаты для разных p могут отличаться. Неточности восстановления могут быть уменьшены за счет усреднения результатов для многих p.

Полученные численные результаты восстановления акустической скорости для раз личных значений лучевого параметра показали, что для малых значений парамет ра p результат восстановления точнее описывает модель. Так, сравнение результа та численного решения обратной задачи для достаточно малого значения парамет ра p = 0.001 мс/м с исходной моделью среды представлено на рис. 10, а, а результат восстановления для p = 0.3 мс/м (значение, близкое к половине максимального) на рис. 10, б. В первом случае алгоритм позволяет получить решение с достаточной степенью точности (максимальная погрешность восстановления акустической скоро сти менее 2%), но меньшей по сравнению с одномерным случаем (см. рис. 6). Отличие связано с дополнительной процедурой дискретного преобразования Радона, которая вносит в данные численную погрешность. Так как моделирование производилось для 74 Д. В. Аникиев, Б. М. Каштан, А. С. Благовещенский, В. А. Мулдер Рис. 11. Сравнение результата численного восстановления акусти ческой скорости (а) и плотности (б ) для пары значений p1 = 0.15 мс/м и p2 = 0.3 мс/м с исходной моделью среды модели импульсного источника, данные моделирования содержат короткие импульсы, а дискретное преобразование Радона дает менее точный результат для таких данных.

В случае p = 0.3 мс/м результат решения обратной задачи (для акустической ско рости) (рис. 10, б ) является менее точным, чем для p = 0.001 мс/м (максимальная по грешность 5%). Это связано с искажениями, которые вносит процедура дискретного преобразования Радона.

Глубина восстановления в двумерном и трехмерном случаях определяется длиной апертуры приемной системы (rmax ) и максимальным временем записи.

На рис. 11 приведены результаты, подтверждающие принципиальную возможность независимого восстановления акустической скорости и плотности для пары значений лучевого параметра p1 = 0.15 мс/м и p2 = 0.3 мс/м.

Точный динамический метод решения обратной задачи сейсмики... 2.4. Влияние сейсмического шума на решение обратной динамической задачи Одной из неотъемлемых проблем при проведении реальных исследований является сейсмический шум. Влияние сейсмического шума на решение обратной динамической задачи подробно рассмотрено в работах Кэрриона и Паттона [46], Кэрриона [47], Кэр риона и Бьюба [48]. Рассмотрим влияние шума на точность восстановления параметров среды на основе двумерной модели.

При моделировании шум (в спектральной области) в виде случайной величины с нормальным распределением (имеется в виду вероятностное распределение случайной величины) в диапазоне частот от 50 до 1000 Гц (1000 Гц частота Найквиста для входных данных) был искусственно добавлен (равномерно распределен по частотному диапазону) в исходные данные двумерной обратной задачи (рис. 12). Частотный диа пазон шума выбран так, чтобы не вносить изменений в низкочастотную (от 0 до 50 Гц) область данных. Проблемы, связанные с искажением спектра в области низких частот, будут рассмотрены далее.

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

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

76 Д. В. Аникиев, Б. М. Каштан, А. С. Благовещенский, В. А. Мулдер Рис. 13. Сейсмическая трасса без импульса прямой волны, зарегистриро ванная приемником с выносом 8 м Рис. 14. Спектральная плотность сейсмической трассы приемника с вы носом 8 м (без импульса прямой волны) Таким образом, значение отношения сигнал/шум для отраженных волн значитель но ниже. На рис. 13 представлена сейсмическая трасса, зарегистрированная первым приемником (с выносом 8 м), из которой вырезан первый импульс, соответствующий прямой волне от источника. Видно, что средняя амплитуда отраженных импульсов при мерно в 2 раза выше средней амплитуды шума. Такое соотношение является типичным для сейсмических данных. Распределение спектральной плотности показано на рис. 14.

Результат восстановления акустической скорости на основе зашумленных данных (рис. 15, а) свидетельствует о том, что метод устойчив к шуму в высокочастотной обла сти (максимальная погрешность менее 3%). Кроме того, существует принципиальная возможность улучшить результат восстановления при наличии шума за счет усредне ния результатов для нескольких значений лучевого параметра.

2.5. Низкочастотная проблема при решении обратной динамической задачи Важнейшей проблемой при сейсмических исследованиях в целом и при решении об ратной задачи в частности является отсутствие во входных данных достоверной инфор мации о низких временных частотах. Во многом это связано с техническими возмож ностями приемных систем. Обычно при обработке сейсмической информации частоты ниже определенного низкочастотного предела отфильтровываются. Стандартные на земные сейсмические приемники позволяют достоверно зарегистрировать данные на Точный динамический метод решения обратной задачи сейсмики... Рис. 15. Сравнение результата численного решения двумерной обратной задачи с исходной моделью среды: а при наличии высоко частотного сейсмического шума;

б на основе данных, для которых применен фильтр низких частот частотах примерно от 5 Гц, морские от 10 Гц (см., например, [49]). Технически мож но достигнуть низкочастотного предела 2 Гц, и в случае морских сейсмических иссле дований для этого требуется источник в виде пневмопушки, обладающей значительно большей мощностью [49]. Так или иначе, данное технологическое направление стреми тельно развивается, смещая низкочастотный предел в меньшую сторону.

При решении динамической обратной задачи потеря низкочастотной информации оказывает сильное влияние на конечный результат. Для численного исследования низ кочастотной проблемы к двумерным данным, изображенным на рис. 9, был применен фильтр низких частот до 0.2 Гц. На рис. 15, б показан результат восстановления среды на основе таких данных.

Низкочастотная проблема при решении обратной динамической задачи детально рассмотрена в работе Кэрриона [47], где показано, что для получения приемлемого ре зультата восстановления среды при отсутствии низких частот необходимо с большой точностью знать форму сигнала. Результат на рис. 15, б получен в предположении, что 78 Д. В. Аникиев, Б. М. Каштан, А. С. Благовещенский, В. А. Мулдер источник является импульсным (дельта-функция), однако искажение, введенное филь трацией низких частот, изменяет характер регистрируемого сигнала. В соответствии с работой [47], даже такое незначительное искажение (всего лишь до 0.2 Гц) сильно вли яет на низкие пространственные частоты в конечном результате восстановления. Из рис. 15, б видно, что наклон (или тренд) восстановленной модели отличается от истин ного. В работе [15] также исследуется подобное влияние, а для решения низкочастотной проблемы используется процедура регуляризации [50], которая применима в том чис ле для алгоритма, рассмотренного в настоящей работе и может помочь в корректном восстановлении тренда исходной модели. Применение подобной регуляризации при под готовке работы не производилось. Регуляризация задачи при отсутствии низких частот является одним из возможных направлений развития предлагаемого метода.



Pages:     | 1 || 3 | 4 |   ...   | 6 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.