Аникиев Д.В. Методы обращения сейсмических волновых полей. Удк методы обращения сейсмических волновых полей
![]()
|
ФУНКЦИОНАЛА НЕВЯЗКИ Наиболее простым способом было бы вычисление градиента с помощью производной Фреше ∂ ( ) ∂ u t m i m x x s g , , , моделированных волновых полей в приемниках: ∂ ∂ = ∂ ∂ ( ) ∂ ( ) ∂ ∑ J m J u u m t t i i t m m , , , , , , , , x x x x s g s g x x s g Данный способ предполагает решение N прямых задач, где N – число упругих параметров среды в каждой точке сетки. Это делает затруднительным использование такого подхода для реалистичных моделей, поскольку число N может достигать десятков миллиардов в трехмерном случае. Применение теории оптимального управления и использование метода множителей Лагранжа позволило определять не только значение, но и градиент функционала невязки с вычислительными затратами порядка времени одного прямого моделирования. В задачах сейсморазведки уравнения движения сред изучены достаточно хорошо, но решение обратной задачи вызывает значительные затруднения. Метод множителей Лагранжа в данном случае применяется не для составления уравнения, которому удовлетворяет распределение физических параметров, доставляющее абсолютный минимум функционалу невязки, а для определения приближенного направления на этот минимум. Данный подход к вычислению градиента функционала невязки в англоязычной литературе получил название adjoint-state method (метод сопряженных уравнений) [Plessix, 2006]. Алгоритм действий для вычисления градиента функционала невязки выглядит следующим образом. Решить прямую задачу в опорной модели для каждого источника ]( ) u x x = ∈ 0, ∀ Ω 2. В каждый приемник в качестве источника поместить производную от функционала невязки полей по моделированному полю в данном приемнике. В случае квадратичного функционала в приемнике эта производная совпадает с разностью моделированного и наблюденного волновых полей. Построить продолженное в обратном времени волновое поле из приемников, решив сопряженное уравнение движения в опорной среде (в большинстве случаев сопряженное уравнение совпадает с обычным уравнением ∂ ∂ = J t u L * , ( , T) ( , T) λ λ λ 0 0 x x 4. Вычислить градиент по формуле: ∂ ∂ = ∂ ∂ ∫ J m m dt i i T λ L u 0 Алгоритм вычисления градиента является важнейшей составляющей метода полного обращения волновых полей. Благодаря ему существует возможность эффективной реализации решения обратной динамической задачи, так как существенно сокращается количество прямых задач, которые необходимо решать для построения решения. ОПРЕДЕЛЕНИЕ НАПРАВЛЕНИЯ ПОИСКА Простейшим методом локальной минимизации является метод наискорейшего спуска. Данный метод предполагает спуск в направлении, противоположном градиенту функционала невязки для начальной модели до достижения минимума функционала вдоль этого направления. К сожалению, скорость сходимости данного метода (количество спусков для достижения приемлемой модели) даже в случае двух параметров невысока [Wolfe, 1969]. Метод сопряженных градиентов существенно сокращает число необходимых направленных спусков. В приложениях чаще всего применяется либо нелинейная версия метода сопряженных градиентов [Polak, Ri- biere, 1969], либо алгоритм Бройдена–Флетчера–Голь- дфарба–Шанно с ограниченной памятью (в англоязычной литературе – метод L-BFGS) [Nocedal, 1980]. В работе [Brossier et al., 2009] авторы сравнили эти два метода применительно к задачам сейсмики и пришли к выводу, что L-BFGS позволяет значительно сократить количество линейных поисков. Тем не менее каждая итерация метода требует дополнительных вычислений, а также сохранения направлений нескольких предыдущих итераций. Метод сопряженных градиентов до сих пор применяется во многих современных реализациях техники полного обращения волновых полей [Butzer et al., 2013; Mao et al., ПРОБЛЕМА ОПРЕДЕЛЕНИЯ ТРЕНДОВОЙ СОСТАВЛЯЮЩЕЙ Традиционно в сейсморазведке трендовая составляющая модели среды восстанавливается решением обратной кинематической задачи, после чего для полученной скоростной модели производится сейсмическая миграция с целью восстановления деталей разреза. Метод полного обращения волновых полей виде- альных условиях освещенности исследуемой среды предполагал объединение результатов этих двух методов. Однако из [Snieder et al., 1989] известно, что определение трендовой составляющей весьма затруднительно. Проблема трендовой составляющей заключается в том, что даже незначительные ошибки в плавной составляющей скоростной модели на начальном этапе приводят, как правило, к сходимости итеративного метода к ошибочной модели. Данный эффект вызван локальностью методов оптимизации, применяемых для поиска минимума целевого функционала. Малые изменения трендовой составляющей порождают существенный фазовый сдвиг регистрируемых колебаний, который ведет к наличию множества локальных минимумов целевого функционала [Jannane et al., 1989], соответствующих ошибочным моделям среды. Проблема локальных минимумов предопределила трудности при первых попытках применения метода (см [Bunks et al., 1995; Kolb et al., 1986; Mora, 1987; Pica et al., 1990; Snieder et al., 1989; Tarantola, Количество локальных минимумов функционала невязки полных полей значительно больше, чему функционала невязки времен пробега, применяемого в лучевой томографии. Как эвристически показано 47 [Bunks et al., 1995], а затем и численно подтверждено [Mulder, Plessix, 2008], наличие дополнительных минимумов обусловлено эффектом пропускания цикла. Эффект наблюдается при несоответствии времен пробега волн в зарегистрированных и смоделированных данных более чем вполовину основного периода регистрируемого сигнала (периода, соответствующего несущей частоте. Сдвиг данных на полный период доставляет локальный минимум функционалу невязки (рис. 2). Если зондирующий сигнал содержит несколько колебательных периодов, то трассы практически идентичны в большей части временных отсчетов, а моделированное поле близко к зарегистрированному, сдвинутому на целое число периодов несущей частоты. МОДИФИКАЦИИ МЕТОДА ПОЛНОГО ОБРАЩЕНИЯ ВОЛНОВЫХ ПОЛЕЙ Для решения проблемы трендовой составляющей и улучшения сходимости процедуры обращения предлагаются различные подходы. Наиболее распространенным в настоящее время является так называемый многомасштабный подход (multiscale approach), который по сути является техникой последовательной обработки данных, содержащих сначала только низкие временные частоты, а затем более высокие [Bunks et al., 1995; Silvestrov et al., 2013]. На низких частотах целевой функционал имеет более плавный характера потому вероятность спуска в локальный минимум вместо глобального значительно ниже. [Bunks et al., 1995] демонстрируют это на примере запаздывания сигнала в опорной модели относительно реально зарегистрированного сигнала. Расстояние между локальными минимумами функционала невязки соответствует кинематической ошибке, равной периоду регистрируемого сигнала. Чем ниже максимальная частота используемых данных, тем б ó льшие кинематические ошибки допустимы. Один из вариантов модификации метода полного обращения волновых полей связан с выбором функционала невязки, отличного от среднеквадратичного. Например, функционал в виде некоторой аналитической функции наблюденных и моделированных сейсмограмм приводит к изменению функции источника в уравнениях для множителей Лагранжа, что позволяет сконцентрировать обращение на определенном типе волн [Van Leeuwen, Mulder, 2008]. Различные виды функционалов обсуждаются в статье [Brossier et al., Несмотря на многочисленные модификации, для практического применения метода полного обращения волновых полей используется следующий подход [Warner et al., 2013]: сначала трендовая составляющая находится с помощью скоростного анализа, использующего только кинематические атрибуты волновых полей (времена пробега, после чего распределение локальных неоднородностей определяется с помощью локальных методов минимизации целевого функционала. При этом основными требованиями к зарегистрированным данным для успешного применения метода полного обращения волновых полей в его современной форме остаются наличие низких временных частот в спектре наблюденных данных во избежание эффекта пропускания цикла наличие сейсмических данных для значительных выносов, гарантирующих полноценное просвечивание сейсмическими волнами локальных неоднородностей на больших глубинах. Если отсутствие низких пространственных частот в зарегистрированных трассах во многих случаях может быть компенсировано подходящим выбором альтернативного функционала невязки, то данные для больших выносов, по-видимому, являются неотъемлемым условием восстановления истинного строения сред этим методом, даже в акустическом случае [Biondi, Almomin, 2012]. В качестве иллюстрации, подтверждающей возможность применения метода в отсутствие низких частотно при наличии данных для больших выносов, мы приводим фрагмент из работы [Biondi, Almomin, 2013] (рис. 3). На рис. 3, в представлен результат восстановления акустической модели среды при помощи метода полного обращения волновых полей. Стартовая горизонтально-однородная модель среды изображена на рис. 3, б, тогда как истинная модель среды – модель Мармуси (Marmousi) – на риса. Зондирующий импульс содержит лишь частоты от 5 до 10 Гц, а максимальный вынос источник приемник составляет 9 км. В результате наличие данных для больших выносов позволяет построить приемлемую модель среды даже при отсутствии низких временных частот в зарегистрированных сейсмо- граммах. При использовании алгоритмов обращения в акустическом приближении существуют очевидные физические трудности, связанные с областью применимости этих алгоритмов. Так, если восстанавливается лишь распределение скорости продольных волн в среде, то это накладывает определенные ограничения на ту информацию, которую может использовать алгоритм, поскольку коэффициенты отражения в упругой среде определяются не только скоростью продольных волн. Если же провести фильтрацию определенных типов отраженных волн, возникающих в реальной упругой анизотропной среде, то и разрешающая способность метода значительно снизится. Подобная несогласованность данных и алгоритмов приводит к трудностям в большинстве практических при- ложений. Рис. 2. Эффект пропускания цикла. При сдвиге на целое число периодов сигналы практически совпадают на протяжении большей части зарегистрированной трассы, что доставляет локальный минимум функционалу невязки ОБРАЩЕНИЕ ВОЛНОВЫХ ПОЛЕЙ ДЛЯ МОДЕЛИ ВЯЗКОУПРУГОЙ АНИЗОТРОПНОЙ СРЕДЫ Методика полного обращения волновых полей в акустическом приближении стала неотъемлемой частью обработки данных сейсморазведки в индустрии [Plessix et al., 2013], однако задача восстановления анизотропных свойств среды до сих пор находится на переходной стадии от теоретических исследований с синтетическими данными [Operto et al., 2013] к практическим реализациям метода для простейших моделей анизотропных сред [Warner et al., В модели, описываемой уравнениями идеальной упругости, сумма потенциальной и кинетической энергии среды сохраняется. Такая ситуация редко встречается в реальности, поскольку зачастую работой сил жидкого и сухого трения нельзя пренебречь. Одним из обобщений теории идеальной упругости является наследственная теория вязкоупругости Больцмана, позволяющая феноменологически описать различные механизмы поглощения, задавая определенный вид интегрального оператора релаксации. Теоретические основы обращения полных волновых полей для произвольно анизотропной вязкоупру- гой среды с причинными функциями релаксации изложены. Принципиальное отличие от случая идеальной упругости заключается в неса- мосопряженных уравнениях движения. А. Тарантола показал, что для решения сопряженной задачи необходимо использовать антипричинные функции ре- лаксации. Численное решение уравнений вязкоупругости во временной области производить неудобно из-за наличия интегрального оператора повремени, поэтому обычно вводится так называемая переменная памяти, позволяющая преобразовать интегральное уравнение в систему дифференциальных уравнений [Car- cione, 1990]. Численное решение уравнений вязкоуп- ругости в частотной области почти не отличается от соответствующего идеально-упругого случая нужно заменить вещественные параметры Ламе на комплексные. Поэтому обращение волновых полей производится, как правило, в частотной области. В общем случае анизотропной вязкоупругой среды после преобразования Фурье повремени уравнение движения записывается в виде+ ∂ ∂ ∂ ∂ = где c c c klmn klmn klmn = + = + Re Im Re Im i c ic klmn klmn . В обозначениях Фойгта тензор упругих параметров с го порядка переписывается в виде квадратной 6 × 6 матрицы. По определению, матрица факторов записывается как Q ij = Re Im C C ij ij Для производных функционала невязки по вещественной и мнимой частям элемента тензора можно получить следующие выражения − ∂ ∂ ∂ ∫ J c U x x d klmn k m l n Re , min max λ ω ω ω ∂ ∂ = − ∂ ∂ ∂ ∫ J c i U x x d klmn k m l n Im , min а в производной функционала невязки по плотности подразумевается суммирование по повторяющемуся индексу Рис. 3. Истинная скоростная модель Marmousi, использованная для получения синтетических данных (а. В зондирующем сигнале представлены частоты от 5 до Гц. Использована система наблюдений с максимальным выносом км горизонтально-однородная скоростная модель, использованная в качестве стартовой (б результат решения обратной задачи методом полного обращения волновых полей (в 49 ∂ ∂ = − ∫ J U d k k ρ ω λ ω ω ω 2 min В этих формулах λ – решение сопряженного уравнения движения с функцией источника в виде разности экспериментального и модельного полей. Вещественная и мнимая части тензора упругих параметров c klmn в соответствии с принципом причинности связаны соотношениями Крамерса–Кронига в частотной комплексной плоскости [Dahlen, Tromp, Впервые полное обращение волновых полей в вязкоакустической среде реализовано в работе [Song et al., 1995], где показано, что определение коэффициента поглощения производится с меньшей точностью по сравнению с обращением скорости. Применение вязкоакустического обращения с использованием локального метода спуска для оценки распределения поглощения и скоростей над газовой залежью описано. В работе [Watanabe et al., 2004] вязкоакустический подход применен к лабораторным данным. [Kamei, Pratt, 2008] рекомендуют осуществить сначала восстановление плотности, а затем совместное восстановление скоростей и поглощения. В статье [Long et al., 2009] предложен предобуслов- ленный градиент функционала невязки для вязко- акустического обращения волновых полей. В работе [Malinowski et al., 2011] проведен анализ чувствительности полного обращения полей в вязкоакустической среде. Среди работ по моделированию и обращению во временной области можно отметить следующие. В статье обращение полных волновых полей применено к синтетическим данным для слоистой вязкоупругой среды в предположении, что бесконечно, те. без учета поглощения продольных волн. В работе [Askan et al., 2007] вязкоакустическое обращение применено к сложной сейсмической модели среды в двумерной задаче, а в [Bai et al., 2012] – к модели Мармуси. Обращение в случае вяз- коупругой анизотропной среды поданным непро- дольного ВСП в Северном море осуществили [Barnes, Charara, 2010]. Различные виды предобусловливателей для вязкоакустического обращения использованы [Causse et al., В статье [Keers et al., 2001] подход на основе динамического лучевого трассирования применен для межскважинных данных. [Ribodetti et al., 2000] показали, что при использовании лучевого метода и приближения Борна скорость волны и фактор могут быть восстановлены одновременно при достаточно большом объеме наблюдений. В работе [Blanch, Symes, 1995] предложен алгоритм итеративного линеаризованного обращения для вязкоакустической модели среды. [Hak, Mulder, 2010, 2011] исследовали возможность построения изображений отражателей, связанных со скачками поглощения. Авторы показали, что использование причинного члена в представлении для фактора, а также геометрия наблюдений существенно влияют на восстановление этих скачков. Таким образом, переход к рассмотрению моделей среды с несколькими распределенными параметрами в задаче обращения волновых полей является сегодня одной из наиболее актуальных проблем сейсмораз- ведки. ОБРАЩЕНИЕ ИНТЕРФЕРЕНЦИОННЫХ ВОЛНОВЫХ ПОЛЕЙ В сейсморазведке при выделении положений отражающих границ поле интерференционных поверхностных волн считается вредным шумом, который маскирует полезные сигналы отраженных волн. Существуют разнообразные методы подавления этого шума. Однако при наземных сейсмических наблюдениях поверхностные волны составляют значительную часть волнового поля и могут быть использованы в качестве источника информации для детального определения параметров и свойств верхней части разреза. Интерференционные поверхностные волны характеризуются дисперсией (зависимостью фазовой и групповой скоростей от частоты. Это свойство лежит в основе метода обращения поверхностных волн для определения параметров исследуемой среды. Изучение поверхностных волн для определения внутреннего строения Земли началось в х гг., однако широко использоваться в практических целях они стали лишь в е гг., чему способствовало развитие ЭВМ, численных методов, а также совершенствование сейсмического оборудования для регистрации землетрясений. В 1950 г В. Томсон предложила в 1953 г. Н. Хаскелл развил теоретический матричный метод расчета волновых полей в горизонтально-одно- родных слоистых средах [Haskell, 1953; Thomson, 1950]. Позже [Озеров, Волин, 1959] показали, что видимый период поверхностной волны Лява характеризует толщину верхнего слоя грунта. Полученные результаты свидетельствовали о возможности исследования геологического строения приповерхностной среды поверхностными волнами. В.И. Кейлис-Борок [1960] провел фундаментальное исследование по возбуждению и распространению интерференционных волн в многослойных средах. Он впервые сформулировал в наиболее общей форме постановку обратной задачи сейсмологии и был инициатором создания методики решения обратной задачи по построению скоростного разреза Земли на основе совокупности сейсмических наблюдений, в том числе по поверхностным волнам [Левшин, Ритцволлер, 2001; Keilis-Borok, Yanovskaya, 1967]. На основе дисперсионных свойств поверхностных волн [Dormann, Ewing, 1962] создали автоматический алгоритм обращения поверхностных волн для определения структуры мантии Земли. Позже метод поверхностных волн нашел широкое применение в сейсмологических исследованиях. Применение поверхностных волн в сейсморазведке и для решения инженерных задач началось лишь с развитием методов спектрального анализа поверхностных волн, позволяющих с большей точностью оценивать их спектральные и поляризационные характеристики [Socco et al., 2010]. Вначале х гг. был предложен метод дисперсионного анализа поверхностных волн на основе оценки разницы между временами прихода сейсмических сигналов к двум разным станциям от одного точечного источника [Nazarian, Stokoe, 1984]. В настоящее время поверхностно-вол- новые сейсмические исследования проводятся с использованием большого количества приемников, что позволяет делать съемку намного быстрее и получать более точные результаты. Такие исследования часто называют многоканальным анализом поверхностных волн (MASW); данный термин введен в работах [Park et al., 1999; Xia et al., 1999]. Определение дисперсионных свойств поверхностных волн из сейсмических данных проводится на основе спектральных преобразований исходной сейсмограммы (преобразование, преобразование и их аналоги. Волновое поле при отсутствии горизонтальных вариаций может быть представлено в виде двойного интеграла Фурье по частоте и волновому числу [McMechan, Yedlin, 1981]: U x t dk d e N k D k i kx t ( , ) ( , ) ( , ) ( ) = ∫ ∫ − ω ω ω ω , (где x – расстояние между источником и приемником, t – время распространения волны, k – волновое число, N(k, ω) – функция возбуждения источника, а D(k, ω) – дисперсионное соотношение. Двойное обратное преобразование Фурье от волнового поля есть преобразование, оно позволяет построить преобразованную сейсмограмму – подынтегральное выражение поля. Амплитуда подынтегрального выражения стремится к бесконечности при значении параметров (k, ω), которые удовлетворяют дисперсионному уравнению. Пикируя положения максимумов на преобразованной сейсмограмме, можно получить значения фазовых скоростей (c = ω/k) различных мод поверхностных волн в зависимости от частоты (такие зависимости параметров от частоты называются дисперсионными кривыми. Пикирование положений максимумов зависит как от качества исходных сейсмограмм, таки от алгоритма пикирования, и может вносить существенную неоднозначность в процедуру обращения. В качестве базовой модели среды для исследования на поверхностных волнах применяется набор изотропных упругих слоев, ограниченных сверху свободной границей, а снизу упругим изотропным полупространством. В такой горизонтально-однородной среде при проведении наблюдений на сейсмограммах выделяются волны Рэлея, или нормальные моды (основная, фундаментальная и высшие, атак- же каналовые продольные волны, приходящие на более ранних временах. При наблюдениях основными составляющими волнового поля являются волны Лява. Дисперсионные свойства поверхностных волн представляются в виде дисперсионных кривых, каждая из которых соответствует определенной моде поверхностной волны. Теоретические дисперсионные кривые строятся как корни дисперсионного уравнения выбранной модели для каждого значения частоты Молотков, 1984; Петрашень и др, 1985; Aki, Richards, 2002]. Для среды, состоящей из набора упругих слоев, дисперсионное уравнение может быть записано в общем виде D(c, f, v p1 , …, v pM , v s1 , …, v sM , ρ 1 , …, ρ M , h 1 , …, h M–1 ) = = 0, где c – фазовая скорость, f – частота, M – количество слоев модели (й слой – нижнее полупространство, v pi , v si , ρ i , h i – соответственно скорость продольной волны, скорость поперечной волны, плотность и толщина го слоя. Дисперсионное уравнение имеет много решений, так как для каждого значения частоты может существовать несколько мод. Дисперсионные кривые волн Рэлея описываются вещественными корнями P-SV- дисперсионного уравнения. Они наиболее чувствительны к значениям поперечных скоростей в скоростном профиле среды и слабо зависят от значений продольных скоростей [Roth et al., 1998; Socco et al., 2010]. Волны Лява чувствительны лишь к значениям скоростей поперечных волн, при этом совсем не зависят от значений скоростей продольных волн [Aki, Richards, 2002]. Каналовые продольные волны состоят в основном из многократных отражений продольных волн в слоистой модели и чувствительны к значениям продольных скоростей в скоростном профиле среды. Дисперсионные кривые каналовых продольных волн описываются комплексными корнями дисперсионного уравнения [Boiero, 2013a; Roth et al., Метод обращения поверхностных волн заключается в восстановлении параметров выбранной скоростной модели минимизацией невязки между дисперсионными кривыми, полученными из реальных сейсмических данных и построенными методом численного поиска корней дисперсионного уравнения. На каждом шаге такой итерационной процедуры происходит обновление параметров модели до тех пор, пока различие между модельными дисперсионными кривыми и полученными из реальных данных не станет минимальным. Необходимо отметить, что для корректности решения следует правильно задавать начальную модель, выбирая подходящее число слоев в модели, а также адекватные интервалы изменения значений упругих параметров. Для этого учитывают имеющуюся априорную информацию с целью уменьшить количество искомых параметров, а также интервалы их изменений, что в итоге и приводит к уменьшению неоднозначности решения обратной задачи. Обычно метод обращения поверхностных волн применяют для восстановления значений скоростей поперечных волн и толщин слоев на основе обращения дисперсионной кривой фундаментальной моды волны Рэлея. Значения остальных параметров при этом считаются известными, а существование высших мод не рассматривается. Однако включение высших мод в процесс обращения играет важную роль в его стабилизации и положительно влияет на качество восстановления параметров среды [Socco et al., 2010; Xia et al., 2013]. В настоящее время увеличивается количество работ, в которых при обращении учитываются и высшие моды. Дисперсионные кривые волн Лява не зависят от значений скоростей продольных волн. Их обращение характеризуется меньшей неоднозначностью, чем обращение фундаментальной моды волны рэлеевского типа. Однако на практике волны Лява применяются намного реже из-за того, что для их возбуждения и приема необходимы горизонтально-поляризованные источники и приемники [Socco et al., В последнее время вследствие увеличения возможностей численного поиска комплексных корней дисперсионного уравнения широкое применение получило использование каналовых продольных волн для восстановления значений скоростей продольных волн [Boiero, Socco, 2009; Ernst, 2007, 2013; Socco et al., 2009], хотя еще совсем недавно их использование в практических целях казалось затруднительным [Roth et al., Стандартная формула функционала невязки между экспериментальными и теоретическими дисперсионными кривыми может быть записана в следующем виде 51 E N m j j j N = = ∑ 1 2 1 δ σ , (где δ j – невязка значения фазовой скорости для каждой точки дисперсионной кривой (для каждой частоты, а σ j – его экспериментальное стандартное отклонение. В зависимости от задачи, от наличия априорной информации вид функционала невязки может меняться, но это не меняет сути метода. При классическом обращении дисперсионных кривых δ j j e j t c c = − – это разница между средним экспериментальными теоретическим значениями фазовой скорости. Под средним экспериментальным значением следует понимать среднее значение из всех пикированных значений (фазовой скорости) для одного значения частоты (при проведении многократных съемок. При однократных съемках σ j принимается равной 1 или некоторому нормировочному значению. Невязки δ j определяются для каждой из выделенных дисперсионных кривых. Существует большая вероятность ошибки в сопоставлении максимумов на преобразованной сейсмограмме с конкретной теоретической дисперсионной кривой (ошибки распознавания мод. При наличии резких скачков между скоростными параметрами, а также при наличии в модели инверсии скорости (более высокоскоростной слой находится над более низкоскоростным слоем) дисперсионные кривые разных мод могут касаться друг друга. В таком случае высока вероятность неправильного распознавания мод. Кроме того, для определения теоретических дисперсионных кривых необходимо решать дисперсионное уравнение, что усложняет обращение и увеличивает вероятность итоговой ошибки – получения неверной восстановленной модели среды. Для решения проблемы распознавания мод предлагались разные методы, среди них можно выделить обращение частотно-координатного спектра приповерхностного волнового поляна основании градиентных алгоритмов [Forbridger, 2003], а также стохастический метод минимизации невязки между наблюдаемым полным спектром фазовой скорости и вычисленными в рамках слоистой модели [Ryden, Park, Новый подход в рамках анализа дисперсионного уравнения, позволяющий обойти проблему неправильного сопоставления экспериментальных и теоретических дисперсионных кривых, предложили [Ernst, 2007; Maraschini et al., 2008]. Поскольку левая часть дисперсионного уравнения (дисперсионное соотношение) стремится кв окрестности истинного решения значений дисперсионной функции в N экспериментальных точках могут быть минимизированы без определения мод и поиска корней дисперсионного уравнения. При этом в выражении формулы невязки есть значение дисперсионного выражения в точке с номером j, а σ j принимается равной 1 или некоторой константе (весовой коэффициент. Такой прямой способ минимизации невязки значительно упрощает алгоритм обращения и уменьшает вероятность итоговой ошибки в восстановлении параметров среды. Минимизация функционала невязки может быть проведена различными численными схемами и методами в зависимости от сложности задачи. Целью любого итерационного процесса минимизации является нахождение экстремума функционала невязки, в данном случае – глобального минимума выражения (8). Чаще всего функционал невязки имеет много локальных минимумов, и главная задача алгоритма минимизации найти среди них глобальный минимум. Большую роль в нахождении глобального минимума играет задание начальной модели. Если начальная модель близка к искомой, то для поиска глобального минимума можно применять градиентные методы. В более общем случае предпочтительнее использование стохастических методов, таких как методы Монте-Карло, метод отжига и др. Ограничение на использование таких методов накладывает лишь зависимость времени расчета от количества искомых параметров. В случае обращения дисперсионных кривых поверхностных волн в слоистой модели количество параметров обычно невелико, и применение стохастических методов является наиболее эффективным методом поиска глобального минимума функционала невязки. В последнее время стохастические алгоритмы минимизации применяются для обращения поверхностных волн во многих работах [Boiero, Socco, 2009; Boiero et al., 2013a; Cercato, Метод обращения поверхностных волн для восстановления параметров среды в рамках горизонталь- но-однородных слоистых сред все чаще применяется для горизонтально-неоднородных сред [Boiero et al., 2013b]. Теория распространения поверхностных волн в средах со слабой горизонтальной неоднородностью рассматривалась еще в х гг. [Бабич и др, 1976]. В дальнейшем некоторые авторы предлагали методы для определения присутствия горизонтальных вариаций параметров исследуемой среды, однако оценка этих вариаций проводилась лишь в рамках одномерной модели среды путем совмещения нескольких восстановленных слоистых моделей для разных участков профиля [Bohlen et al., В настоящее время для решения инженерных задачи задач приповерхностной сейсморазведки используется метод, основанный на использовании многократных сейсмических съемок. Вдоль сейсмического профиля получают несколько дисперсионных кривых для восстановления локальных одномерных скоростных моделей, которые затем совмещаются для восстановления псевдодвумерной структуры. Такая псевдодвумерная скоростная модель затем используется как начальная модель для совместного обращения всех дисперсионных кривых, полученных на профиле, одновременно с учетом дополнительных ограничений, определяемых из априорной информации. Данный метод получил название обращение сгори- зонтальными связями (laterally constrained inversion, LCI). Впервые он применен [Auken, Christiansen, 2004] для интерпретации электромагнитных данных, а для обращения поверхностных волн он успешно применен, В последнее время широкое распространение получил метод полного обращения волновых полей. Однако его применение для восстановления параметров приповерхностной среды часто затруднительно ввиду достаточно сложного геологического строения. В таких случаях для получения надежного результата восстановления параметров необходима начальная модель, которую может предоставить метод обращения поверхностных волн [Boiero et al., 2013a]. При работе с наземными данными, особенно в областях со сложным строением приповерхностной среды, применение полного обращения волнового поля для упругой среды является довольно сложной в вычислительном плане задачей. На основании представления волнового поля в приповерхностной среде в виде суперпозиции небольшого числа горизонтально распространяющихся мод поверхностных волн F. Ernst [2013] предложил метод, позволяющий существенно упростить данную задачу. В предложенном приближении для расчета прямой задачи используется теоретическое вычисление волнового поляна основе решения дисперсионного уравнения вида для каждой точки на поверхности модели среды. При этом учитываются наиболее значительные составляющие волнового поля и существенно экономится время решения прямой задачи по сравнению с применением конечно-разност- ного алгоритма. |