ОЦЕНИВАНИЕ СВОЙСТВ ПОЧВЫ С ИСПОЛЬЗОВАНИЕМ ВОЛНОВЫХ СИГНАЛОВ СЕЙСМИЧЕСКИХ ПОВЕРХНОСТНЫХ ВОЛН
ESTIMATION OF SOIL PROPERTIES USING WAVE SIGNALS OF SEISMIC SURFACE WAVES

Перекрестные ссылки на родственные заявки
По этой заявке испрашивается преимущество приоритета предварительной заявки №61/087933 на патент США, которая была подана 11 августа 2008 года.
Область техники, к которой относится изобретение
В общем это изобретение относится к области геофизической разведки, а более конкретно к обработке сейсмических данных. По существу, изобретение представляет собой способ обращения данных о сейсмической поверхностной волне для получения упругих свойств почвы, таких как скорость поперечной волны в почве, или модуль сдвига, или затухание поперечной волны. В частности, можно получать упругие свойства как функцию глубины от поверхности или глубинный разрез свойств. Способ также можно использовать для оценивания местоположений аномалий в почве, таких как пустоты или погребенные объекты.
Предпосылки создания изобретения
Сейсмические поверхностные волны, также называемые поверхностными волнами-помехами, или волнами Релея, или волнами Лява, заключены в области вблизи земной поверхности, и поэтому распространение их зависит от упругих свойств верхней части разреза, особенно от скорости поперечной волны как функции глубины от поверхности. Скорость поперечной волны прямо связана с жесткостью почвы, определяемой модулем сдвига почвы (Скорость поперечной волны равна корню квадратному из модуля сдвига, деленному на плотность). Скорость поперечной волны или профиль модуля упругости и другие упругие характеристики почвы можно непосредственно использовать для инженерных или других целей или можно использовать косвенно для повышения качества геофизической разведки под почвой или областью верхней части разреза. В дополнение к модулю сдвига или скорости волны упругое затухание или спад амплитуды с расстоянием также является полезной информацией, которую можно использовать для инженерных изысканий и геофизической разведки.
Скорость поперечной волны в почве или модуль сдвига можно определять, обращая кривые дисперсии (зависимости фазовой скорости от частоты) поверхностных волн-помех для получения профиля скорости поперечной волны в среде. Вследствие эффектов уплотнения скорости обычно ниже вблизи земной поверхности и возрастают с глубиной. Высокочастотные составляющие поверхностных волн заключены вблизи поверхности и следуют по более медленным слоям почвы. С другой стороны, низкочастотные составляющие следуют по более глубоким, более быстрым слоям. Поэтому скорость поверхностной волны изменяется с частотой, то есть является дисперсионной. В частности, скорость снижается с повышением частоты. Форму дисперсионной кривой как функцию частоты можно сравнивать с расчетными дисперсионными кривыми, вычисленными для многослойного профиля скорости, и после этого свойства профиля, то есть толщину слоя и модуль сдвига, можно обновлять для лучшего соответствия измеренным дисперсионным кривым. Многослойность приводит к резонансным эффектам и захвату различных мод поверхностной волны-помехи. Более высокую точность можно получать, обращая дисперсионные кривые для основной и более высокого порядка мод поверхностной волны-помехи.
После того как поверхностные волны преобразованы в глубинные разрезы скорости поперечной волны в почве или модуля сдвига, затухания поперечной волны или другого свойства, информацию можно использовать непосредственно в качестве важного инженерного параметра при расчете конструкций, таких как здания и мосты. Другие применения, связанные с непосредственным использованием поверхностных волн для определения характеристик верхней части разреза, включают в себя оценивание реакции участка на землетрясение, контроль уплотнения почвы, картирование пологой поверхности, оценивание прочности подземных материалов, оценивание дорожного покрытия, обнаружение погребенных искусственных сооружений или аномалий, оценивание пустот вокруг канализационной сети и обнаружение глубины залегания коренной породы. Профили скорости в верхней части разреза можно косвенно использовать, чтобы повышать качество определения физической структуры или физического свойства более глубоких подземных областей для оценивания запасов или извлечения углеводородов. Поскольку верхняя часть разреза характеризуется низкой скоростью и является неоднородной, она оказывает сильное влияние на сейсмические волны, которые проходят через нее, и может ограничивать возможность определения структуры и свойств глубоких областей. Профиль скорости в верхней части разреза можно использовать для ввода временных поправок и статических поправок за отражения сейсмических волн от более глубоких зон, или информацию о скорости можно использовать для построения изображения, выполнения миграции или для обращения сейсмических данных.
Проблема, связанная со способами из предшествующего уровня техники, в которых используют поверхностные волны, волны Релея или волны Лява, заключается в том, что трудно разрешать или различать разные дисперсионные кривые для различных мод поверхностной волны. Одна характерная трудность обусловлена неопределенностью частоты и скорости, и обращение является компромиссным, когда идентифицируют неверно выбранную моду или неправильно выбирают моду. Вторая трудность обусловлена отсутствием сведений о фазе источника и неопределенностями изменений фазы при значениях выше 2π. Кроме того, в таких способах при анализе используют усредненные свойства поверхностной волны в пределах протяженности расстановки приемников, и это усреднение ограничивает разрешение. Наконец, интерференция между модами и шумом и затухание поверхностных волн могут искажать амплитуды сейсмических волн, что делает трудной идентификацию отдельных дисперсионных кривых. Далее до некоторой степени подробно будут рассмотрены традиционные способы определения модуля сдвига почвы или скорости поперечных волн в почве.
Современные способы использования поверхностных волн для описания характеристик модуля сдвига почвы или скорости волны включают в себя регистрацию сейсмических данных, за которой следуют две стадии обработки: (1) измерение дисперсионных кривых как функции частоты и затем (2) обращение дисперсионных кривых для получения модуля сдвига как функции глубины. Аналогичные способы можно использовать для получения других свойств помимо модуля сдвига, таких как характеристики затухания. Способы измерения дисперсионных кривых сильно различаются по количеству источников и количеству приемников при регистрации. В старых способах используют один источник и одну пару приемников. В новейших способах используют один источник и многочисленные приемники (порядка 20 или больше), разнесенные с образованием регулярных интервалов. В большей части способов используют источник продольных волн, который возбуждает волны Релея, но также можно использовать источник поперечных волн для возбуждения волн Лява. Способ использования волны любого типа является одним и тем же.
Стадия 1: способы с парой приемников
Старые способы с использованием волн Релея для описания модуля сдвига почвы представляют собой «способы установившегося состояния», включающие применение сейсмического вибратора для создания вибрации грунта на одной частоте или на медленно изменяющейся частоте (патент США №3864667, Bahjat, 1975). На каждой частоте измеряют разность фаз между откликами двух геофонов. На основании этих измерений получают свойства верхней части разреза между приемниками. Однако для регистрации таких измерений затрачивают значительное время. В 1980-годах был разработан способ спектрального анализа поверхностных волн (САПВ) (Nazarian et al., “Use of spectral analysis of surface waves method for determination of moduli and thickness of pavement systems”, Transport. Res. Record 930, 38-54 (1983)). В способе определяют дисперсионную кривую, сначала вычисляя взаимный спектр мощности между сигналами, регистрируемыми двумя датчиками, и после этого развертывая фазу.
Для способа спектрального анализа поверхностных волн и способа установившегося состояния характерны одинаковые проблемы. Поскольку в данный момент времени используют только пару приемников, трудно провести различие между влияниями различных мод поверхностной волны и влияниями любого регистрируемого шума. Расстояние между приемниками и между парой приемников и источником изменяют, чтобы минимизировать, но не исключить влияния мод высокого порядка. Характерной проблемой, связанной с этими измерениями, является неоднозначность около фазовой постоянной 2π. Приемники должны находиться достаточно близко друг к другу с тем, чтобы фаза не изменялась больше чем на 2π. Важно суметь провести различие между изменением фазы на Δ и изменением на Δ+π; последнее соответствует более низкой скорости. Иногда используют несколько дополнительных приемников для содействия развертыванию фазы.
Как указывают Holschneider и соавторы в “Characterization of dispersive surface waves using continuous wavelet transforms”, Geophys. J. Int. 163, 463-478 (2005), в последнее время вейвлет-преобразование используют в качестве способа, позволяющего повысить ослабление одной моды поверхностной волны-помехи и минимизировать вклады других мод. Вейвлет-преобразование представляет собой частотно-временное преобразование, которое может обеспечить лучшую локализацию отдельных мод. Способом снижают интерференцию мод, но шум все еще остается проблемой вследствие 2π-неоднозначности фазы. Holschneider и соавторы выполнили построение математической модели эффектов распространения поверхностной волны в области вейвлет-преобразования и использовали ее для нахождения сначала дисперсионных кривых, а затем кривых затухания одной моды за один раз.
Стадия 1: многочисленные приемники
Способ многоканального анализа поверхностных волн (МАПВ) был разработан Kansas Geological Survey. В этом способе размещают 20-65 или большее количество приемников и регистрируют данные от одного импульсного или вибрационного источника (Park et al., “Multichannel analysis of surface wave”, Geophysics 64, 800-808 (1999); и Park et al., “Multichannel analysis of surface waves (MASW) active and passive methods”, The Leading Edge 26, 60-64 (2007)). Данные обрабатывают как одну сейсмограмму общего пункта возбуждения, амплитуды нормируют или масштабируют, а затем преобразуют в частотно-волночисловую (f-k) область или в область частота-медленность (f-p). Затем дисперсионные кривые для одной или нескольких мод выбирают в точках максимальных амплитуд в областях f-k или f-p. Дисперсионные кривые для одной или нескольких мод используют для последующего обращения свойств верхней части разреза (Beaty et al., “Repeatability of multimode Rayleigh-wave dispersion studies”, Geophysics 68, 782-790 (2003)). Регистрацию можно повторять при каждом новом возбуждении, направляемом в новую расстановку приемников, и каждую сейсмограмму общего пункта возбуждения отдельно анализировать, чтобы для каждой расстановки получать одномерный профиль скорости в верхней части разреза. Затем сочетают каждый одномерный профиль со схемой интерполяции для образования двумерного профиля верхней части разреза.
Многоканальный способ является усовершенствованием способа спектрального анализа поверхностных волн. Одно преимущество заключается в том, что при использовании разнесенных на небольшие расстояния приемников происходит минимизация 2π-неоднозначности изменений фазы; однако невозможно оценивать фазу источника по одной сейсмограмме (Hermann and Ammon, “Surface Waves, Receiver Functions, and Crustal Structure: Version 3.3”, in Computer Programs in Seismology, Saint Louis University, http:www.eas.slu.edu/People/RBHermann/CPS330.htl. (2004)). Его второе преимущество заключается в том, что по определению преобразование в f-k- или f-p-области включает в себя объединение или суммирование по трассам, которое повышает разрешающую способность по частоте и уменьшает проблемы, связанные с шумом. Однако взамен этого происходит потеря разрешающей способности в поперечном направлении; при этом невозможно обнаруживать изменения скорости по ширине расстановки. Как и в случае двух каналов, можно выбирать диапазон удалений для подчеркивания или ослабления различных мод (Xia et/ al., “Utilization of high-frequency Rayleigh waves in near-surface geophysics”, The Leading Edge 23, 753-759 (2004)). В дополнение к этому приемники необходимо разносить равномерно, располагая близко друг к другу, чтобы не иметь ложных частот пространственной волны-помехи. Кроме того, поскольку амплитуды пространственной волны-помехи сильно ослабляются, они изменяются от трассы к трассе, и при суммировании амплитуды искажаются. Используют нормализацию или уравновешивание амплитуд, но все же трудно выбирать и различать многочисленные интерферирующие моды поверхностной волны-помехи. Lefebvre и Benhassen (заявка №2005/0143924 А1 на патент США) сообщили об использовании вейвлет-преобразования для повышения возможности различения разных мод. Forbridger (“Inversion of shallow-seismic wavefields: I. Wavefield transformation”, Geophys. J. Int. 153, 719-734 (2003)) указали на трудности, связанные с выделением многомодовых дисперсионных кривых, и проблемы, связанные с последующим обращением, в случае когда дисперсионные кривые неправильно выделены или неправильно идентифицированы. Forbridger, а также Ryden и Park (“Fast simulated annealing inversion of surface waves on pavement using phase-velocity spectra”, Geophysics 71, R49-R58 (2006)) исключили выделение дисперсионных кривых, а вместо него выполняют непосредственное обращение результатов f-p-преобразования. Эти способы включают в себя математическую модель распространения поверхностной волны в f-p-области, для которой требуются несколько предположений и приближений.
Стадия 1: многочисленные приемники и многочисленные источники
Одновременное использование многочисленных положений источников и многочисленных приемников для получения кривых изменяющейся в поперечном направлении фазовой скорости рассмотрели Ernst и соавторы в “Tomography of dispersive media”, J. Acoust. Soc. Am. 108, 105-115 (2000) и в “Removal of scattered guided waves from seismic data”, Geophysics 67, 1240-1248 (2002). Применением не является описание характеристик верхней части разреза, а ослабление рассеиваемой поверхностной волны-помехи при сейсмической разведке на нефть. Их способ включает в себя ряд последовательных операций. Сначала они обращают изменяющуюся в поперечном направлении фазовую скорость как функцию частоты, используя томографический способ, основанный на обобщенных временах пробега. Они предполагают, что поперечные изменения являются небольшими и можно получать ослабление одной моды во временном окне. Для справедливости этого последнего предположения требуется, чтобы источник и приемники были на достаточно большом расстоянии друг от друга для хорошего разделения мод во времени, но вследствие сильного затухания мод поверхностной волны трудно или невозможно получать данные для описания характеристик верхней части разреза. Использование обобщенных времен пробега включает вычисление производной фазы данных, и проблемы фазовой неопределенности при определении фазовой скорости поверхностной волны-помехи возрастают.
Стадия 1: затухание в зависимости от кривых частоты
Хотя большинство способов из предшествующего уровня техники сосредоточено на дисперсионных кривых фазовой скорости, использование кривых затухания как функции частоты рассмотрели Xia и соавторы (работа цитировалась). Наряду с модулем сдвига качественный показатель (Q) как функция глубины также является важной технической величиной, но обращение данных о затухании дает меньшую стабильность. Обычно предполагают, что затухание не зависит от частоты (Ernst et al., “Removal of scattered guided waves from seismic data”, Geophysics 67, 1240-1248 (2002); и Kulesh et. al., “Modeling of Wave Dispersion Using Continuous Wavelet Transforms II: Wavelet Based Frequency-Velocity Analysis”, Pure & Applied Geophysics 165, 255-270 (2008)). Однако в случае использования поверхностных волн для описания характеристик верхней части разреза это предположение является ограниченным. Поскольку затухание обычно уменьшается при возрастании глубины почвы, затухание поверхностных волн должно уменьшаться таким же образом в зависимости от частоты, поскольку фазовая скорость является дисперсионной.
Стадия 2: обращение дисперсионных кривых для профилей скорости в верхней части разреза
Имеются несколько алгоритмов, пригодных для нахождения профиля скорости в верхней части разреза на основании дисперсионных кривых, но успех всех таких способов зависит от точности входных дисперсионных кривых. Обращение представляет собой задачу нелинейной оптимизации модели, в которой моделью является профиль скорости в верхней части разреза. Параметры включают в себя глубину слоя и модуль сдвига слоя. Алгоритмы включают в себя линеаризированное обращение методом наименьших квадратов, алгоритм Левенберга-Марквардта, квазиньютоновский, а в последнее время имитированного отжига (Beaty et al., “Simulated annealing inversion of multimode Rayleigh wave dispersion curves for geological structure”, Geophys. J. Int. 151, 622-631 (2002)). Доступное программное обеспечение включает в себя свободно распространяемое программное обеспечение (Hermann and Ammon, “Surface Waves, Receiver Function, and Crustal Structure: Version 3.3”, in Computer Programs in Seismology, Saint Louis University, http:www.eas.slu.edu/People/RBHermann/CPS330.htl. (2004)) и коммерческое программное обеспечение (SeisOpt ReMi™, http:www.optimsoftware.com; и Kansas Geological Survey: http:/www.kgs.ku.edu/software/surfseis./index.html.).
Остается необходимость в усовершенствованном способе преобразования поверхностных волн в глубинные профили свойств верхней части разреза путем получения высокоразрешающих, изменяющихся в поперечном направлении кривых многомодовой дисперсии и затухания. В частности, в способе должны минимизироваться неоднозначности фазы источника. Настоящее изобретение удовлетворяет этим нуждам.
Краткое изложение сущности изобретения
В одном осуществлении изобретением является способ оценивания модуля сдвига, скорости поперечной волны, затухания поперечной волны или другого физического свойства области верхней части разреза геологической среды на основании трасс сейсмических данных, соответствующих по меньшей мере одному положению источника и множеству положений приемников при сейсмическом исследовании области, содержащий этапы, на которых:
(а) разделяют область верхней части разреза на одну или несколько ячеек;
(b) для каждой трассы и соответствующих мест нахождения сейсмического источника и приемников вычисляют протяженность траектории луча через каждую имеющуюся промежуточную ячейку от места нахождения источника до места нахождения приемника;
(с) для каждой трассы одновременно находят решения для по меньшей мере двух поверхностно-согласованных составляющих, каждая из которых характеризует эффекты фильтрации распространения сейсмических поверхностных волн через ячейку или связь соответствующего сейсмического источника или приемника с грунтом, при этом указанным решением итерационно оптимизируют составляющие, сравнивая спрогнозированную поверхностную волну, вычисленную с использованием предполагаемых или итерационно обновленных составляющих, наряду с информацией о траектории луча из этапа (b), с соответствующей трассой данных из исследования;
(d) выбирают одну или несколько из поверхностно-согласованных составляющих и используют их для вычисления численным обращением модуля сдвига или другого свойства области верхней части разреза.
Краткое описание чертежей
Настоящее изобретение и его преимущества можно лучше понять при обращении к нижеследующему подробному описанию и сопровождающим чертежам, на которых:
Фигуры 1А-В - иллюстрации регистрации сейсмических данных с использованием многочисленных приемников в случае линейного (1А) и площадного (1В) применений;
Фигуры 2А-В - иллюстрации регистрации сейсмических данных с использованием многочисленных источников в случае линейного (2А) и площадного (2В) применений;
Фигуры 3А-В - иллюстрации регистрации сейсмических данных с использованием многочисленных пар источник-приемник при различных удалениях, охватывающей область поверхности в случае линейного (3А) и площадного (3В) применений;
Фиг. 4 - иллюстрация физических процессов, которые влияют на волновой сигнал, когда поверхностные волны распространяются от источника к приемнику вдоль земной поверхности;
Фиг. 5 - схема небольшого участка сейсмического исследования на картографическом виде, показывающая местоположения источников и приемников и разделение на области или ячейки и показывающая примеры траекторий лучей от источников к приемникам через различные области;
Фиг. 6 - блок-схема последовательности основных этапов в одном осуществлении способа настоящего изобретения;
Фиг. 7 - блок-схема последовательности этапов, предназначенных для выполнения части оптимизации модели из фиг. 6 в одном осуществлении изобретения, при этом образуют поверхностно-согласованные составляющие (передаточные функции), определяющие распространение поверхностной волны;
Фиг. 8 - вид четырех трасс для четырех расстояний (удалений) источник-приемник при компьютерной имитации полного волнового уравнения для поверхностных волн в случае единственного тонкого подземного слоя в пределах полупространства;
Фигуры 9A-D - результаты оптимизации модели для примера из фиг. 8, где изменения волновых сигналов разложены в комплексные (амплитудные и фазовые) составляющие источника и составляющие распространения как функции частоты;
Фиг. 10 - спрогнозированные волновые сигналы, вычисленные по составляющим из фигур 9A-D; волновые сигналы можно сравнивать с входными сейсмическими данными из фиг. 8;
Фигуры 11А-В - данные (11А) трасс и f-k-спектр (11В) при компьютерной имитации полного волнового уравнения для поверхностных волн в случае 20 тонких подземных слоев в пределах полупространства, при этом можно видеть интерференцию 6 мод поверхностных волн-помех;
Фиг. 12 - отобранные волновые сигналы для небольшого количества удалений (протяженностей трасс до приемников) из данных трасс на фиг. 11А; при этом волновые сигналы не имеют простой формы, а образованы суперпозицией многочисленных мод;
Фиг. 13 - спрогнозированные волновые сигналы для таких же удалений, как и для фиг. 12; в этом прогнозе использованы параметры для 6 мод;
Фиг. 14A-F - виды одной трассы (14А) данных из фиг. 12 и прогнозов, в которых спрогнозированные волновые сигналы ограничены одной (14В), двумя (14С), тремя (14D), четырьмя (14Е) и шестью (14F) модами; включение дополнительной моды повышает соответствие трассе (14А) данных; и
Фиг. 15 - кривые дисперсии скоростей (скорости как функции частоты) для 6 различных мод поверхностных волн-помех из фигур 11А-В, вычисленные способом настоящего изобретения.
Изобретение будет описано применительно к примерам осуществлений. Однако в том смысле, что нижеследующее описание является специфическим для конкретного осуществления или конкретного использования изобретения, оно предназначено только для иллюстрации и не подразумевается ограничивающим объем изобретения. Наоборот, оно предполагается охватывающим все варианты, модификации и эквиваленты, которые могут быть включены в объем изобретения, определяемый прилагаемой формулой изобретения.
Подробное описание примеров осуществлений
Настоящее изобретение представляет собой способ получения по сейсмическим данным кривых дисперсии скоростей и затухания как функции частоты для поверхностных волн. Предпочтительно регистрировать данные, используя многочисленные источники и многочисленные приемники. После этого кривые или их модификации можно использовать для обращения свойств верхней части разреза, таких как скорость поперечной волны, или модуль сдвига, или затухание поперечной волны как функции глубины, известными численными методами. Оптимизацию модели и избыточность данных используют, чтобы найти решение для параметров поверхностно-согласованной модели, которые наилучшим образом представляют волновые сигналы в зарегистрированных сейсмических данных и изменения волновых сигналов, когда поверхностные волны распространяются вдоль земной поверхности от источников к приемникам. Параметры связаны с отдельными составляющими или передаточными функциями фильтров для каждого местоположения источника, каждого местоположения приемника и распространения через каждую область поверхности. Параметры модели могут включать в себя многочисленные моды поверхностной волны и поверхностно-согласованные поперечные вариации.
Поэтому в изобретении предпочтительно использовать сейсмические данные, регистрируемые при наличии множества источников и приемников, применяя способ обращения из двух стадий. На стадии 1 вариации волновых сигналов поверхностных волн разлагают в поверхностно-согласованные передаточные функции, предпочтительно, для каждого источника, каждого приемника и каждой небольшой области поверхности. В дальнейшем на стадии 2 передаточные функции для каждой области обращают, чтобы определить свойства почвы или свойства верхней части разреза (такие как скорость поперечной волны) как функцию глубины. Способом можно находить решения в случае сложного многомодового характера поверхностных волн для сред с изменениями свойств по вертикали и в поперечном направлении. При этом исключаются погрешности и ограничения разрешающей способности традиционных способов, возникающие вследствие ошибочной идентификации мод поверхностной волны-помехи или вследствие предположения равномерности свойств почвы в поперечном направлении.
Регистрация сейсмических данных
Сначала сейсмические данные должны быть зарегистрированы или получены. Обычно можно использовать сейсмические данные, зарегистрированные для других целей, например для построения изображения геологической среды, см. ниже обсуждение предпочтительных параметров регистрации. Если такие данные отсутствуют, то должно быть проведено специальное исследование для получения свойств поверхностных волн. Если информация о верхней части разреза необходима только вдоль линии поверхности, то данные можно регистрировать при использовании источников и приемников вдоль линии, как при стандартной регистрации двумерных сейсмических данных. Если необходима информация с площади, то источники и приемники должны быть распределены по всей площади, при этом источники и приемники должны охватывать исследуемую площадь. Как и при любой регистрации сейсмических данных, местоположения источников и приемников должны измеряться и должны иметься наряду с зарегистрированной трассой данных. Можно использовать сейсмический источник любого вида, такой как взрывчатое вещество, вибраторы, воздушные пушки, падающий груз, ударные воздействия, сейсмические пушки и т.д. Предпочтительно возбуждать источники среди многочисленных приемников (например десяти или большего количества), как показано на фигурах 1А-В для линейной (фиг. 1А) и площадной (фиг. 1В) регистрации. На фиг. 1А показаны один источник (10) и несколько приемников (11). Понятно, что между источником и каждым приемником имеется траектория (12) прямой волны. На фиг. 1 возможные траектории лучей проведены только в небольшом количестве. Расстояние между источником и приемником называют удалением. С другой стороны, как показано на фигурах 2А и 2В, в динамике во времени, когда источник перемещают на различные места, на каждом приемнике (20) следует регистрировать данные от многочисленных возбуждений (21) (предпочтительно, десяти или большего количества), используя различные траектории (22) лучей. Эти два требования легче всего выполнить, образовав расстановку приемников из большого количества приемников (60-1000) и затем возбуждая источники, один за другим, по всей расстановке. Нет необходимости в том, чтобы количество возбуждений было таким же, как количество приемников. Для снижения затрат количество возбуждений может быть меньше, чем количество приемников, или наоборот. Если имеющимися приборами можно регистрировать только ограниченное количество каналов, то регистрацию можно ограничить приемниками с максимальным расстоянием (максимальным удалением) от каждого источника.
Кроме того, предпочтительно располагать источники и приемники так, чтобы многочисленные траектории лучей приемник-источник (предпочтительно, 10-20 или больше) проходили через каждую исследуемую область в соответствии с различными удалениями, находящимися в диапазоне от заданного минимального удаления до заданного максимального удаления. На фиг. 3А область показана прямоугольником 30. На чертеже показаны некоторые из траекторий (33) лучей, которые проходят через прямоугольник 30 от источника, такого как 31, к приемнику, такому как 32. Имеется некоторое количество таких траекторий лучей (не все показаны), соответствующих различным удалениям. На фиг. 3В показаны те же самые характеристики при площадном исследовании. Минимальное удаление должно быть достаточно большим, чтобы обеспечивалось распространение плоской волны, которое устанавливается на расстоянии около 1/2 длины поверхностной волны. Это расстояние примерно равно глубине проникновения. Предпочтительно, чтобы максимальное удаление было достаточно небольшим для исключения излишнего затухания поверхностных волн. Обычно оно в 2-3 раза больше минимального удаления. Можно выполнять регистрацию на больших и меньших удалениях. Предпочтительно, чтобы разнесение приемников было не больше величины заданного разрешения в поперечном направлении, но, если возможно, целесообразно иметь половину этой величины дискретизации. Источники и приемники можно укладывать согласно регулярной сетке или нерегулярной сетке или можно использовать псевдослучайную выборку. Может возникнуть необходимость планировать исследование с учетом использования дорог для размещения источников или с учетом обхода зданий и сооружений. Кроме того, может быть полезным размещение источников и приемников за пределами заданной площади, подлежащей охвату, чтобы гарантировать достаточную избыточность данных на краях.
Для оптимизации параметров исследования, таких как минимальное удаление, максимальное удаление и дискретизация приемников и источников, полезно иметь информацию о площади. Может оказаться целесообразным выполнение некоторого предварительного испытания для определения частотного диапазона источника и диапазона скоростей поверхностных волн-помех. Это позволит вычислять длины поверхностных волн. Кроме того, может быть полезно выполнять регистрацию в центре площади, возбуждение в двумерный профиль при частой дискретизации приемников, чтобы f-k- или f-p-анализ можно было использовать для получения исходных параметров. Для этого особого профиля дискретизация (то есть разнесение) приемников должно быть меньше, чем максимальная частота, умноженная на наименьшую скорость, чтобы данные не были искажены вследствие недостаточной частоты выборки.
Стадия 1: оптимизация модели
Одна характеристика, помимо различий траекторий лучей, по которой поверхностные волны отличают от других сейсмических волн (таких как нисходящие волны, которые отражаются от границ раздела), представляет собой волновой сигнал. Волновые сигналы поверхностных волн являются высокоамплитудными, низкочастотными и осциллирующими. С другой стороны, сейсмический импульс источника обычно представляет собой высокочастотный импульс или сейсмический импульс небольшой длительности. Физические процессы при таком изменении волнового сигнала пояснены на фиг. 4, на которой показана часть многослойной геологической среды 110. Когда мода поверхностной волны распространяется от источника, расположенного в точке 101 на поверхности, к приемнику, расположенному в точке 102 на поверхности, вдоль земной поверхности (по траектории 103 луча), она подвергается значительной фильтрации средой и при этом изменяется от импульса (104) источника небольшой длительности до более продолжительного, осциллирующего низкочастотного волнового сигнала (105) поверхностной волны. Поверхностная волна задерживается и затухает. Чем больше расстояние (103) распространения, тем больше изменение. При наличии оценки волнового сигнала (104) источника и регистрируемого выходного сигнала (105) это изменение можно выразить количественно, вычислив передаточную функцию для каждой трассы. Ее можно назвать передаточной функцией полной трассы, и она характеризует суммарную фильтрацию средой, которая приводит к изменению волнового сигнала от волнового сигнала, генерируемого на источнике, до волнового сигнала, регистрируемого на геофоне.
Фильтры обычно характеризуют их передаточными функциями или их импульсными характеристиками. Передаточную функцию определяют как результат деления в частотной области выходного сигнала фильтра на входной сигнал фильтра. Для схемы регистрации сейсмических данных из фиг. 4 передаточная функция T(f) имеет вид
<mrow><mi>T</mi><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo><mtext> </mtext><mo>=</mo><mtext> </mtext><mfrac><mrow><mtext>выходной сигнал(</mtext><mi>f</mi><mo stretchy="false">)</mo></mrow><mrow><mtext>входной сигнал(</mtext><mi>f</mi><mo stretchy="false">)</mo></mrow></mfrac><mtext> </mtext><mo>=</mo><mtext> </mtext><mfrac><mrow><mi>D</mi><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo></mrow><mrow><mi>S</mi><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo></mrow></mfrac></mrow>
<mrow><mtable><mtr><mtd><mrow /></mtd><mtd><mrow /></mtd><mtd><mrow /></mtd><mtd><mrow /></mtd></mtr></mtable></mrow>
,
(1)
где D(f) является преобразованием Фурье трассы 105 данных и S(f) является преобразованием Фурье сейсмического импульса 104 источника. Передаточную функцию можно также получить по взаимной корреляции сейсмического импульса и трассы данных, разделив ее на автокорреляцию сейсмического импульса источника.
<mrow><mi>T</mi><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo><mtext> </mtext><mo>=</mo><mtext> </mtext><mfrac><mrow><msup><mi>S</mi><mo>∗</mo></msup><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo><mi>D</mi><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo></mrow><mrow><msup><mi>S</mi><mo>∗</mo></msup><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo><mi>S</mi><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo></mrow></mfrac></mrow>
<mtext> </mtext>
,
(2)
где звездочка обозначает комплексно сопряженный элемент. Передаточная функция является комплексной; чтобы полностью характеризовать эффекты фильтрации, она должна включать в себя амплитуду и фазу или действительную и мнимую части. Вычисление обратного преобразования Фурье дает импульсную характеристику фильтра. Передаточная функция и импульсная характеристика являются эквивалентными понятиями в частотной и временной области, соответственно.
Хорошо известно, что можно произвести свертку друг с другом отдельных линейных фильтров (или выполнить умножение в частотной области), чтобы получить комбинированный фильтр. На фиг.4 можно видеть, что общее воздействие на волновой сигнал поверхностной волны складывается из различных физических процессов, каждый из которых зависит от местоположения на поверхности или в области. Например, из связи энергии источника, преобразуемой в моду поверхностной волны-помехи в точке 101, распространения через область 106, распространения через область 107 и связи приемника с грунтом в точке 102. Кроме того, фильтрация вследствие распространения через область 106 должна быть функцией расстояния (108), пробегаемого волной по области 106, и фильтрация вследствие распространения через область 107 должна быть функцией расстояния (109), пробегаемого волной по области 107. Фильтрация благодаря действию распространения включает в себя влияние задержки или скорости как функции частоты (дисперсии) и ослабление амплитуд как функции частоты. Поэтому для одной моды поверхностной волны полную передаточную функцию T(f) трассы можно разложить на отдельные передаточные функции или фильтры для каждого из упомянутых выше физических процессов, или
<mrow><mi>T</mi><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo><mtext> </mtext><mo>=</mo><mtext> </mtext><msub><mi>T</mi><mrow><mn>101</mn></mrow></msub><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo><msub><mi>T</mi><mrow><mn>106</mn></mrow></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><msub><mi>d</mi><mrow><mn>108</mn></mrow></msub><mo stretchy="false">)</mo><msub><mi>T</mi><mrow><mn>107</mn></mrow></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><msub><mi>d</mi><mrow><mn>109</mn></mrow></msub><mo stretchy="false">)</mo><msub><mi>T</mi><mrow><mn>102</mn></mrow></msub><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo></mrow>
<mtext> </mtext>
.
(3)
Используя уравнения (1) и (3), можно получить модельное выражение M(f), которое описывает регистрируемый волновой сигнал поверхностной волны-помехи.
101
106
108
107
109
102
<mtext> </mtext>
M(f)=S(f)T(f)T(f,d)T(f,d)T(f).
(4)
Уравнение (4) является моделью волнового сигнала для единственной моды поверхностной волны-помехи. Дополнительно предположим, что данные состоят из линейной суперпозиции многочисленных мод. Кроме того, предположим, что на начальной стадии отдельные моды не взаимодействуют или не связаны. Поэтому для N мод модель комплексного волнового сигнала поверхностной волны-помехи можно представить в виде
<mrow><mi>M</mi><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo><mtext> </mtext><mo>=</mo><mtext> </mtext><mi>S</mi><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mrow><msub><mi>T</mi><mrow><mn>101,</mn><mi>j</mi></mrow></msub><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo><msub><mi>T</mi><mrow><mn>106,</mn><mi>j</mi></mrow></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><msub><mi>d</mi><mrow><mn>108</mn></mrow></msub><mo stretchy="false">)</mo><msub><mi>T</mi><mrow><mn>107,</mn><mi>j</mi></mrow></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><msub><mi>d</mi><mrow><mn>109</mn></mrow></msub><mo stretchy="false">)</mo><mo>,</mo><msub><mi>T</mi><mrow><mn>102,</mn><mi>j</mi></mrow></msub><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo></mrow></mstyle></mrow>
<mtext> </mtext>
,
(5)
где суммирование производится по всем отдельным модам от j=1 до N. Количество мод обычно небольшое, порядка 2-6.
В уравнении (5) отдельные физические процессы представлены в виде передаточных функций или фильтров в частотной области. Предпочтительно включать все такие фильтры, даже если для нахождения свойств верхней части разреза необходимы только эффекты распространения. Передаточные функции являются комплексными и включают в себя амплитуду и фазу или действительную и мнимую части как функцию частоты. В случае распространения они включают в себя эффекты скорости и затухания. Несколько различных процессов можно включать в модельное выражение. Некоторые из них являются известными или могут быть измерены. Например, один фильтр является фильтром системы регистрации. Импульсная характеристика системы регистрации может быть измерена специалистом в данной области техники. Еще одна может быть инструментальной характеристикой датчиков, которую можно получить от изготовителей. Третья является функцией разброса амплитуд, относительно которой известно, что она обратно пропорциональна корню квадратному из полного расстояния от источника до приемника для каждой трассы. Может быть полезно сосредотачивать несколько связанных процессов в одном члене. Например, если имеется группа приемников, эффекты группирования можно сосредотачивать в члене связи приемника с грунтом. В обоих случаях они связаны с конкретным положением приемника на земной поверхности. Кроме того, может быть полезно включать один член, который можно полагать известным, и затем иметь второй переменный поправочный член. Например, можно оценивать усредненный сейсмический импульс источника, а затем в член связи источника с грунтом можно включать связь и вариации сейсмического импульса от возбуждения к возбуждению.
При линейном или площадном исследовании в настоящем изобретении используют избыточность трасс данных, которые имеют различные траектории лучей и расстояния источник-приемник, для разложения полных изменений волнового сигнала поверхностной сейсмической волны на отдельные поверхностно-согласованные составляющие или фильтры, которые представляют отдельные физические процессы. Это использование избыточности данных показано на фиг. 5 для площадных данных. На фиг. 5 показано картографическое изображение небольшого разреза из возможного сейсмического исследования. На показанном разрезе имеется источник, местоположение которого обозначено 201. На самом деле источник может быть расположен ниже поверхности в сейсмической скважине. На схеме имеются 24 таких источника, при этом каждый обозначен символом солнечной вспышки. Имеются 49 приемников, обозначенных значком трапеции, один из которых конкретизирован цифрой 202. В предположении, что регистрацию выполняют для каждого приемника при каждом возбуждении источника, можно иметь 24×49 или 1176 трасс данных. В этом изобретении площадь исследования разделяют на области или ячейки. В этом случае прямоугольные ячейки показаны пунктирными линиями, а одна центральная обозначена на чертеже позицией 203. Каждую пару источник-приемник соединяет траектория луча, которую в первом приближении можно считать прямой линией. Показаны только некоторые примеры траекторий лучей, такие как 204, 205 и 206, и с различными путями пробега в выбранной области 203. В предположении только одной моды поверхностной волны каждую частоту трасс данных в частотной области можно использовать для нахождения неизвестных компонентных поверхностно-согласованных передаточных функций. В этом примере одной моды неизвестные величины включают в себя 24 члена связи источников с грунтом, 49 членов связи приемников с грунтом и 9 членов распространения (поскольку имеются 9 дискретных ячеек), всего 82 неизвестные величины. (Этим иллюстрируется оптимизация на этапе 303 из фиг. 6, который рассмотрен ниже). Поскольку имеются 1176 блоков данных и 82 неизвестные величины, задача является переопределенной. В случае, когда имеются две моды, могут быть 164 неизвестные величины. Ключевой частью этого изобретения является то, что различные неизвестные величины, связанные с различными физическими процессами, определяют одновременно. Имеются 49 трасс, которые несут информацию об источнике в 201, имеются 24 трассы, которые несут информацию о приемнике 202, и имеются по меньшей мере 460 траекторий лучей через область 203, каждый из которых должен быть снабжен весом в соответствии с отличающейся протяженностью пути через область 203. Специалисты в данной области техники должны распознать эту задачу как двумерную задачу томографии на прямолинейных лучах для поверхности.
В способе настоящего изобретения предпочтительно использовать модельное выражение с параметрами, которые характеризуют влияния важных физических процессов, и функцию стоимости, которая показывает согласие между данными и модельным выражением. Целесообразно включать как можно больше известных эффектов, то есть инструментальных характеристик системы регистрации. При конкретной регистрации некоторые составляющие могут быть неважными и могут не учитываться. Например, если используют одиночные приемники и если они в основном хорошо связаны с грунтом, связь с грунтом можно не учитывать или заменить простой скалярной величиной. Как минимум должно быть по меньшей мере два параметра, например, член распространения через одну область и член источника для распределения энергии в моду поверхностной волны-помехи. Функция стоимости для оптимизации методом наименьших квадратов (например) представляет собой сумму по количеству трасс квадрата разности между фактическими данными и модельным выражением для каждой трассы. Сумму по всем трассам данных используют при нахождении параметров. Функцию стоимости можно выражать в частотной области, во временной области или в некоторой области преобразования, в том числе, например, в f-k-области, Радона, вейвлет-области, Габора, комплексной трассы, Гилберта, или других областях, которые должны быть известны специалистам в данной области техники. Кроме того, к данным и модели можно применять фильтр, такой как взаимно-корреляционный фильтр. В дополнение к этому сейсмические записи можно обрабатывать интерферометрическими методами по активным и пассивным источникам. Предпочтительно использовать частотную область и связывать функцию стоимости с разностями между действительными и мнимыми представлениями, то есть представлениями амплитуд и фаз из данных и модельного выражения. Пример функции стоимости (для этапа 303) с использованием действительной и мнимой частей имеет вид
<mtable columnalign="left"><mtr><mtd><mi>min</mi><mtext> </mtext><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mrow><mtext>Трасса N</mtext></mrow></munderover><mrow><mrow><mo>[</mo><mrow><msup><mrow><mrow><mo>|</mo><mrow><mi>r</mi><mi>e</mi><mi>a</mi><mi>l</mi><mo stretchy="false">(</mo><msub><mi>D</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>s</mi><mi>x</mi><mo>,</mo><mi>s</mi><mi>y</mi><mo>,</mo><mi>r</mi><mi>x</mi><mo>,</mo><mi>r</mi><mi>y</mi><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mtext> </mtext><mo>−</mo><mtext> </mtext><mi>r</mi><mi>e</mi><mi>a</mi><mi>l</mi><mo stretchy="false">(</mo><msub><mi>M</mi><mn>1</mn></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>s</mi><mi>x</mi><mo>,</mo><mi>s</mi><mi>y</mi><mo>,</mo><mi>r</mi><mi>x</mi><mo>,</mo><mi>r</mi><mi>y</mi><mo>,</mo><mover accent="true"><mi>T</mi><mo>→</mo></mover><mo stretchy="false">)</mo><mo stretchy="false">)</mo></mrow><mo>|</mo></mrow></mrow><mn>2</mn></msup></mrow><mo>]</mo></mrow></mrow></mstyle><mo>+</mo></mtd></mtr><mtr><mtd><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mrow><mtext>Трасса N</mtext></mrow></munderover><mrow><mrow><mo>[</mo><mrow><msup><mrow><mrow><mo>|</mo><mrow><mi>i</mi><mi>m</mi><mi>a</mi><mi>g</mi><mo stretchy="false">(</mo><msub><mi>D</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>s</mi><mi>x</mi><mo>,</mo><mi>x</mi><mi>y</mi><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mtext> </mtext><mo>−</mo><mtext> </mtext><mi>i</mi><mi>m</mi><mi>a</mi><mi>g</mi><mo stretchy="false">(</mo><msub><mi>M</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>s</mi><mi>x</mi><mo>,</mo><mi>s</mi><mi>y</mi><mo>,</mo><mi>r</mi><mi>x</mi><mo>,</mo><mi>r</mi><mi>y</mi><mo>,</mo><mover accent="true"><mi>T</mi><mo>→</mo></mover><mo stretchy="false">)</mo><mo stretchy="false">)</mo></mrow><mo>|</mo></mrow></mrow><mn>2</mn></msup></mrow><mo>]</mo></mrow></mrow></mstyle><mtext> </mtext><mo>,</mo><mtext> </mtext><mo stretchy="false">(</mo><mn>6</mn><mo stretchy="false">)</mo></mtd></mtr></mtable>
<mover accent="true"><mi>T</mi><mo>→</mo></mover>
<mover accent="true"><mi>T</mi><mo>→</mo></mover>
<mrow><msub><mi>M</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mover accent="true"><mi>T</mi><mo>→</mo></mover><mo stretchy="false">)</mo></mrow>
1
2
3
i
трасс
где символ
представляет собой вектор, составленный из всех отдельных компонентных передаточных функций (T(f),T(f),T(f) и т.д.). Значения вектора
определяют путем минимизации квадрата разностей между модельной трассой
и трассой D(f) данных, просуммированного по количеству N трасс. В уравнении (5) для разностей используют показатель степени 2, который соответствует задаче оптимизации методом наименьших квадратов, но можно использовать другие показатели степени, называемые нормами, такие как показатель 1, норма L1. Имеются многочисленные различные способы, называемые способами оптимизации модели, для нахождения параметров модели, которые наилучшим образом согласуются с данными. Все они находятся в объеме способа настоящего изобретения.
Стадия 2: обращение свойств верхней части разреза
После уже описанного этапа оптимизации модели полученные параметры модели для различных поверхностно-согласованных составляющих как функции частоты можно преобразовать в кривые скорости или затухания. Затем эти кривые, по отдельности или совместно, можно обратить (на этапе 308, рассмотренном дополнительно ниже) для параметров верхней части разреза, используя любой алгоритм оптимизации модели, такой как наименьших квадратов, Левенберга-Марквардта, квазиньютоновский и имитированного отжига. Как вариант можно использовать само модельное выражение, передаточные функции или фильтрованную версию параметров модели. В любом случае профиль верхней части разреза принимают в качестве исходной модели, дисперсионные кривые или другие функции вычисляют, используя принятую модель, и затем сравнивают их с найденными по сейсмическим данным. После этого модель профиля обновляют для лучшего согласования с полученными параметрами, при этом скорость поперечной волны в слое является одним из параметров модели, который корректируют в процессе итерационного обновления или оптимизации. В тот момент, когда несоответствие между моделируемыми и полученными результаты сейсмических исследований уменьшается ниже приемлемого допустимого отклонения или достигается другое условие остановки, итерационный процесс заканчивают, и параметры модели, в том числе скорость поперечной волны, которые включают в последнюю обновленную версию модели профиля скорости в верхней части разреза, представляют собой выходные результаты способа настоящего изобретения. Термин «верхняя часть разреза», используемый в этой заявке, не может быть точно определен, но иногда может считаться означающим нахождение в пределах одной длины поверхностной волны-помехи от поверхности, которое зависит от частоты.
Этапы способа
Основные этапы в одном осуществлении способа настоящего изобретения показаны на блок-схеме последовательности действий из фиг. 6. В качестве входных данных для способа требуются трассы сейсмических данных и местоположения источника и приемников для каждой трассы. Трассы могут быть однокомпонентными или многокомпонентными трассами, где в данном случае слово «компонента» относится к детекторным блокам. На этапе 301 определяют местоположения на поверхности для членов распространения. Это можно сделать, разделив площадь исследования на одну или несколько областей, таких как ячейки, подобные показанной позицией 203 на фиг. 5. Может быть выгодно начинать с больших областей, а затем подразделять их при последующих итерациях (см. этап 306). Области могут быть регулярными или нерегулярными. Как вариант сетку можно выбирать на основании другой информации, такой как геологические объекты, типы почвы, поднятие суши или результаты анализов скорости и частоты мод поверхностных волн-помех. На этапе 302 вычисляют протяженности траекторий лучей через каждую область для каждой трассы и сохраняют для использования в будущем (В случае практических применений многие или все этапы способа выполняют на компьютере, имеющем память для хранения данных в запоминающем устройстве компьютера или вспомогательные устройства хранения данных). Достаточно использовать прямолинейный луч между местами нахождения источника и приемника, но если информация о скорости для различных областей доступна из предшествующих итераций, то можно использовать криволинейный луч с траекторией, определенной по закону Снелла. Для вычисления траекторий лучей можно использовать любой прибор построения лучей. Протяженность траектории луча представляет собой расстояние вдоль земной поверхности и может включать в себя изменения возвышения поверхности.
На этапе 303 все сейсмические данные трасс или блоки данных используют, чтобы одновременно найти решение для параметров двух или большего количества поверхностно-согласованных составляющих или передаточных функций. Решение представляет собой разложение данных в различные составляющие. Каждая составляющая должна представлять физические процессы, которые обладают эффектом фильтрации или ограничения полосы относительно сейсмической волны, таким как эффекты излучения источником, или обнаружение приемником, или распространения через область поверхности. Эффекты источника или приемника в первую очередь обусловлены степенью переменной связи этих приборов с геологической средой, но включают в себя любые другие эффекты, связанные с излучением сейсмической энергии и ее распределением в моду поверхностной волны-помехи (поверхностной волны) или с обнаружением сигнала из грунта приемником, если такие эффекты изменяют сейсмический импульс. Высокая степень неизменной связи источников и приемников с грунтом во всех случаях является желательной, но никогда не получается идеальной. Параметры поверхностно-согласованных составляющих не являются параметрами сред, такими, как скорость и плотность для отдельных глубинных слоев, но для физического процесса представляют изменения волнового сигнала на различных частотах. Как минимум следует включать одну составляющую источника и одну составляющую распространения или две составляющие распространения. Например, используют фильтр источника, оцененный для сейсмических трасс при всех местоположениях источников, и используют фильтр распространения для определенной моды распространения поверхностной волны, оцененный во всех ячейках, пересекаемых сейсмическими трассами. В еще одном примере минимума можно включать два фильтра распространения, каждый из которых соответствует отличающейся моде распространения поверхностной волны. Составляющие определяют при наилучшем соответствии между волновыми сигналами поверхностной волны в данных и волновыми сигналами поверхностной волны, прогнозируемыми согласно модельному выражению (то есть, согласно математической модели волнового сигнала для поверхностных волн после распространения от заданного места нахождения источника до заданного места нахождения приемника через ячеистую сетку, определенную на этапе 301), которое включает в себя описанные выше параметры. Соответствие между измеряемыми и моделируемыми волновыми сигналами предпочтительно определять с помощью функции стоимости, основанной на любой норме, описанной выше. Любой способ оптимизации модели, такой как метод сопряженных градиентов, скорейшего спуска или Левенберга-Марквардта, можно использовать для нахождения параметров модели. Целесообразно корректировать данные для известных или аппроксимированных физических процессов или включать их в модельное выражение. Например, можно включать характеристику прибора или входной сейсмический импульс вибратора. Для решения может потребоваться (в зависимости от модельного выражения) выбор исходного набора параметров модели, и любую априорную информацию можно использовать для определения исходных значений. Кроме того, может быть целесообразно выполнять оконную обработку, уплощение, масштабирование или нормирование данных. Кроме того, для данных и модели можно использовать преобразования, а соответствие между данными и моделью сравнивать в области преобразования. Кроме того, может быть целесообразно придавать веса для различных удалений. Например, более высокие частоты теряются или затухают при более значительных удалениях, и поэтому этим удалениям можно придавать веса при вкладе их в функцию стоимости. В зависимости от конкретного способа оптимизации этап 303 может быть итерационным процессом оптимизации несмотря на отсутствие петли, такой как 306, показанной на фиг. 6.
Можно ограничивать оптимизацию на этапе 303 выбором параметризации на этапе 301. На этапе 305 оценивают согласие между вычисленными и реальными волновыми сигналами, предпочтительно во временной области, и в зависимости от согласия оценивают возможность разделения областей для получения дополнительного разрешения в поперечном направлении. Если согласование с данными не считают достаточным, то процесс можно повторить (ответвление 306) при другой параметризации, начиная с этапа 301. Согласие может указывать на поперечную изменчивость, и области можно разделить. В дополнение к этому можно включать дополнительные составляющие, например, большее количество мод поверхностных волн. Кроме того, можно вводить поправки за эффекты высшего порядка, такие как криволинейные лучи или фокусировка, дефокусировка амплитуд, потери на прохождение и рассеяние. При повторении этапа 303 параметры модели из предшествующих решений можно использовать для получения исходных значений для следующего этапа оптимизации.
Когда на этапе 305 оптимизацию считают хорошей и при желании, если заданное разрешение в поперечном направлении получено, параметры модели оптимизируют способом настоящего изобретения, и стандартные способы можно использовать для вычисления скорости поперечной волны и других параметров верхней части разреза, указанных на этапах 307 и 308. На этапе 307 полученные параметры вводят в форму для обращения свойств в верхней части разреза. Например, кривые дисперсии скоростей для каждой моды можно извлекать и использовать для обращения. Как вариант также можно получать другие функции параметров. Примеры включают в себя составную передаточную функцию, решение по плоским волнам или прогнозируемые, свободные от шума волновые сигналы. Другие модификации включают в себя выполнение взаимных корреляций или преобразований, таких как вейвлет-преобразование, f-k-преобразование или f-p-преобразования.
Наконец, на этапе 308 полученные в заданной форме параметры используют для описания характеристик верхней части разреза. Для исходных свойств верхней части разреза, таких как скорость поперечной волны, делают предположение. Оно может включать в себя задание количества слоев и назначение каждому слою толщины и свойства, такого как скорость поперечной волны. Затем форму полученных параметров, то есть дисперсионные кривые для каждой моды, вычисляют для исходного набора свойств верхней части разреза и сравнивают с полученными на этапе 307. Это сравнение может включать в себя вычисление разности в некоторой степени, например в степени 2 (наименьшие квадраты), между результатом из этапа 307 и значениями исходных параметров верхней части разреза. Затем параметры верхней части разреза обновляют для получения лучшего согласования. Можно использовать стандартные алгоритмы оптимизации модели.
Этап 303
j
<mrow><msub><mi>M</mi><mi>v</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>s</mi><mi>x</mi><mo>,</mo><mi>s</mi><mi>y</mi><mo>,</mo><mi>r</mi><mi>x</mi><mo>,</mo><mi>r</mi><mi>y</mi><mo>,</mo><mover accent="true"><mi>T</mi><mo>→</mo></mover><mo stretchy="false">)</mo></mrow>
Далее в качестве примера осуществления изобретения некоторые признаки настоящего изобретения описываются более подробно. Прежде всего описывается предпочтительная параметризация отдельных поверхностно-согласованных составляющих или передаточных функций. Начинают с трассы с заданными координатами (sx, sy) источника на поверхности и координатами (rx, ry) приемника на поверхности. Траекторию луча между (sx, sy) и (rx, ry) определяют для пробега через каждую область j, и при этом расстояние пробега в области j составляет x. В таком случае модельное выражение для вертикальной составляющей
имеет вид
<mtable columnalign="left"><mtr><mtd><msub><mi>M</mi><mi>v</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>s</mi><mi>x</mi><mo>,</mo><mi>s</mi><mi>y</mi><mo>,</mo><mi>r</mi><mi>x</mi><mo>,</mo><mi>r</mi><mi>y</mi><mo>,</mo><mover accent="true"><mi>T</mi><mo>→</mo></mover><mo stretchy="false">)</mo><mo>=</mo></mtd></mtr><mtr><mtd><mo>=</mo><msup><mo stretchy="false">(</mo><mi>a</mi><mi>o</mi><mi>f</mi><mi>f</mi><mi>s</mi><mi>e</mi><mi>t</mi><mo stretchy="false">)</mo><mrow><mo>−</mo><mn>1</mn><mo>/</mo><mn>2</mn></mrow></msup><mi>S</mi><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo><mi>R</mi><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo><mi>I</mi><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo><msub><mi>C</mi><mi>v</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>r</mi><mi>x</mi><mo>,</mo><mi>r</mi><mi>y</mi><mo stretchy="false">)</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>n</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mrow><msub><mi>P</mi><mrow><mi>n</mi><mo>,</mo><mi>v</mi></mrow></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>s</mi><mi>x</mi><mo>,</mo><mi>s</mi><mi>y</mi><mo stretchy="false">)</mo><msup><mi>e</mi><mrow><mo>−</mo><mi>i</mi><mn>2</mn><mi>π</mi><mi>f</mi><mstyle displaystyle="true"><munderover><mo>∑</mo><mi>j</mi><mrow><mtext>траектория</mtext></mrow></munderover><mrow><msub><mi>s</mi><mrow><mi>n</mi><mo>,</mo><mi>f</mi></mrow></msub><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mstyle></mrow></msup></mrow></mstyle></mtd></mtr></mtable>
<mtext> </mtext>
,
(7)
<mrow><msub><mi>M</mi><mi>h</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>s</mi><mi>x</mi><mo>,</mo><mi>s</mi><mi>y</mi><mo>,</mo><mi>r</mi><mi>x</mi><mo>,</mo><mi>r</mi><mi>y</mi><mo>,</mo><mover accent="true"><mi>T</mi><mo>→</mo></mover><mo stretchy="false">)</mo></mrow>
а соответствующая горизонтальная составляющая
имеет вид
<mtable columnalign="left"><mtr><mtd><msub><mi>M</mi><mi>v</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>s</mi><mi>x</mi><mo>,</mo><mi>s</mi><mi>y</mi><mo>,</mo><mi>r</mi><mi>x</mi><mo>,</mo><mi>r</mi><mi>y</mi><mo>,</mo><mover accent="true"><mi>T</mi><mo>→</mo></mover><mo stretchy="false">)</mo><mo>=</mo></mtd></mtr><mtr><mtd><mo>=</mo><msup><mo stretchy="false">(</mo><mi>a</mi><mi>o</mi><mi>f</mi><mi>f</mi><mi>s</mi><mi>e</mi><mi>t</mi><mo stretchy="false">)</mo><mrow><mo>−</mo><mn>1</mn><mo>/</mo><mn>2</mn></mrow></msup><msup><mi>e</mi><mrow><mi>i</mi><mi>π</mi><mo>/</mo><mn>2</mn></mrow></msup><mi>S</mi><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo><mi>R</mi><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo><mi>I</mi><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo><msub><mi>C</mi><mi>h</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>r</mi><mi>x</mi><mo>,</mo><mi>r</mi><mi>y</mi><mo stretchy="false">)</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>n</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mrow><msub><mi>P</mi><mrow><mi>n</mi><mo>,</mo><mi>h</mi></mrow></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>s</mi><mi>x</mi><mo>,</mo><mi>s</mi><mi>y</mi><mo stretchy="false">)</mo><msup><mi>e</mi><mrow><mo>−</mo><mi>i</mi><mn>2</mn><mi>π</mi><mi>f</mi><mstyle displaystyle="true"><munderover><mo>∑</mo><mi>j</mi><mrow><mtext>траектория</mtext></mrow></munderover><mrow><msub><mi>s</mi><mrow><mi>n</mi><mo>,</mo><mi>j</mi></mrow></msub><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mstyle></mrow></msup></mrow></mstyle></mtd></mtr></mtable>
(8)
Предшествующие выражения представлены для геофонов; можно получить аналогичные выражения для других датчиков, таких как акселерометры или гидрофоны. В уравнениях (7) и (8) некоторые составляющие являются детерминированными, то есть их можно измерить или аппроксимировать. Другие параметры являются неизвестными и должны находиться при оптимизации модели на этапе 303.
Детерминированные составляющие:
-1/2
<mrow><msqrt><mrow><msup><mrow><mo stretchy="false">(</mo><mi>s</mi><mi>x</mi><mtext> </mtext><mo>−</mo><mtext> </mtext><mi>r</mi><mi>x</mi><mo stretchy="false">)</mo></mrow><mn>2</mn></msup><mtext> </mtext><mo>+</mo><mtext> </mtext><msup><mrow><mo stretchy="false">(</mo><mi>s</mi><mi>y</mi><mtext> </mtext><mo>−</mo><mtext> </mtext><mi>r</mi><mi>y</mi><mo stretchy="false">)</mo></mrow><mn>2</mn></msup></mrow></msqrt></mrow>
(aoffset) - коэффициент разброса амплитуд или поправка за расхождение для поверхностной волны, распространяющейся вдоль двумерной поверхности (Berkholt, Applied Seismic Wave Theory, Elsevier Science Publisher, p.142 (1987)). Величина aoffset представляет собой абсолютное значение полного расстояния от источника до приемника и имеет вид
.
S(f) - сейсмический импульс источника в частотной области: он представляет собой усредненный или характеристический сейсмический импульс для всего исследования. Способы измерения его описаны ниже.
R(f) - импульсная характеристика приемника в частотной области: ее можно получить от изготовителя. Она представляет собой усредненную характеристику датчиков, но если отдельные кривые чувствительности измерены, то их можно использовать.
I(f) - импульсная характеристика системы регистрации в частотной области. Ее можно получить от изготовителя или измерить, подавая импульсное напряжение на вход системы регистрации.
jπ/2
e - поворот фазы на 90° для горизонтальной составляющей.
Неизвестные составляющие:
v
h
C(f, rx, ry) или C(f, rx, ry) - связь приемника с грунтом для вертикального или горизонтального датчика на месте (rx, ry) нахождения приемника на поверхности. Этот член представляет вариацию чувствительности датчика и связь приемника с грунтом на месте (rx, ry). Связь приемника или геофона с грунтом можно описать комплексной передаточной функцией (Krohn, “Geophone ground coupling”, Geophysics 49, 722-731 (1984)). Она также может включать в себя эффект группирования, для которого может потребоваться азимутальная зависимость источник-приемник. В этой реализации предполагается независимость от мод, но можно использовать член связи приемника с грунтом, который изменяется в зависимости от моды.
n,v
n,h
P(f, sx, sy) или P((f, sx, sy) - связь источника с грунтом и распределение в моду n для вертикального и горизонтального датчиков, соответственно. Это важный сосредоточенный параметр. Он включает в себя вариации от источника к источнику и связь источника с грунтом на месте (sx, sy). Он также включает в себя распределение амплитуд энергии источника в каждую моду. Каждую моду можно считать интерференцией плоских волн, движущихся в многослойной структуре, образующей волновод. В многослойной среде различные частоты распространяются избирательно и поэтому каждая мода будет иметь отличающийся амплитудный спектр. Кроме того, основная мода будет начинаться с нулевой частоты, но для последующих мод нижние критические частоты будут повышаться (Aki and Richards, Quantitative Seismology: Theory and Methods, volume I, W.H. Freeman and Co., 259-318 (1980)). Различия между горизонтальными и вертикальными составляющими приводят к различным поляризациям разных мод.
n,j
j
s(f) - комплексная медленность при распространении через область j моды n. Комплексную медленность умножают на протяженность x траектории луча через ячейку j. В этой реализации медленность принимают не зависящей от азимута или направления распространения. Если на месте исследования имеется достаточная анизотропия, азимутальную зависимость можно использовать.
ikx
<mrow><mi>k</mi><mtext> </mtext><mo>=</mo><mtext> </mtext><mfrac><mrow><mn>2</mn><mi>π</mi><mi>f</mi></mrow><mrow><mi>v</mi><mo>'</mo></mrow></mfrac></mrow>
<mrow><mfrac><mn>1</mn><mrow><mi>v</mi><mo>'</mo></mrow></mfrac><mtext> </mtext><mo>=</mo><mtext> </mtext><mfrac><mn>1</mn><mi>v</mi></mfrac><mtext> </mtext><mo>+</mo><mtext> </mtext><mfrac><mi>i</mi><mrow><mn>2</mn><mi>v</mi><mi>Q</mi></mrow></mfrac></mrow>
n,j
n
n
Экспоненциальный член в уравнениях (7) и (8) представляет собой форму распространения e плоской волны, где k является комплексным волновым числом (Aki and Richards, упомянутая выше работа, с.173). Волновое число k можно записать как
, где
, а Q является качественным показателем эффективного затухания. Для каждой области j определяют член s(f) медленности, который включает в себя действительную часть, обратную величину фазовой скорости v(f) для поверхностной моды n и мнимую часть, связанную с частотно-зависимым качественным показателем Q
<mrow><mi>r</mi><mi>e</mi><mi>a</mi><mi>l</mi><mo stretchy="false">(</mo><msub><mi>s</mi><mrow><mi>n</mi><mo>,</mo><mi>j</mi></mrow></msub><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mtext> </mtext><mo>=</mo><mtext> </mtext><mfrac><mn>1</mn><mrow><msub><mi>v</mi><mi>n</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo></mrow></mfrac></mrow>
<mtext> </mtext>
<mrow><mi>i</mi><mi>m</mi><mi>a</mi><mi>g</mi><mo stretchy="false">(</mo><msub><mi>s</mi><mrow><mi>n</mi><mo>,</mo><mi>j</mi></mrow></msub><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mtext> </mtext><mo>=</mo><mtext> </mtext><mfrac><mn>1</mn><mrow><mn>2</mn><msub><mi>Q</mi><mi>n</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo><msub><mi>v</mi><mi>n</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo></mrow></mfrac></mrow>
<mtext> </mtext>
,
.
(9)
Поэтому комплексная медленность является параметром каждой из составляющих распространения для каждой области j в уравнениях (7) и (8).
В настоящем изобретении находят решение для фильтра геологической среды, зависящего от области. Для полного учета влияний фильтрации геологической среды на волновой сигнал поверхностной волны предпочтительно включать параметры скорости и затухания как функцию частоты. Кроме того, в случае более чем одной моды влияния фаз и амплитуд являются связанными и не могут быть найдены независимо. На этой стадии не требуется оценивать вертикальный профиль скорости или затухания в геологической среде. Вместо этого находят только частотную зависимость скорости и затухания. После того как они найдены, их можно использовать на этапе 308 для определения профиля скорости и затухания.
Возможность аппроксимации сейсмического импульса S(f) источника в уравнениях (1), (4), (5), (7) и (8) позволяет сделать оптимизацию модели более робастной. В литературе авторы обсуждают проблемы и трудности при распознавании фазы сейсмического импульса источника, особенно в случае вибрационных источников (Ziolkowski, “Why don't we measure seismic signatures”, Geophysics 56, 190-201 (1991); и Gibson and Lamer, 1984, “Predictive deconvolution and the zero-phase source”, Geophysics 49, 379-397 (1884)). Вследствие этих осознаваемых трудностей член источника не включают в задачи оценивания скорости поверхностной волны-помехи и рассмотрение передаточной функции (уравнение (1)) не используют для ослабления поверхностной волны-помехи. Например, в обеих ранее упомянутых статьях Ernst et al. член источника не включен в первое решение для поля фазовых скоростей.
Эти трудности могут быть исключены в настоящем изобретении способами аппроксимации сейсмического импульса источника и с помощью члена связи источника с грунтом и распределения излучения, описанного выше. Членом связи источника с грунтом корректируются неадекватности аппроксимации сейсмического импульса источника. Trantham (“Controlled-phase acquisition and processing”, Soc. Explor. Geophy. Expanded Abstracts 13, 890 (1994)) показал, что регистрация сейсмического импульса может быть управляемой. Кроме того, Krohn (Международная заявка WO2004/095073) изложил способ обработки данных вибратора путем выполнения деконволюции с использованием сигнатуры вибратора, определяемой на основе измерений на вибраторах, и формирования данных с получением желаемой импульсной характеристики или сейсмического импульса. Полученный обработкой сейсмический импульс вибратора, образованный методом Krohn, является сейсмическим импульсом, предпочтительным для использования в этом изобретении при регистрации вибросейсмических данных. С другой стороны, можно использовать автокорреляцию свип-сигнала. В случае данных, образуемых с использованием группы воздушных пушек, можно использовать сигнатуру воздушной пушки в дальнем поле. Сигнатуру в дальнем поле можно измерять или можно моделировать на основании характеристик отдельных воздушных пушек. Наконец, в случае данных от взрыва или других данных без сигнатуры источника, сейсмический импульс источника можно образовать, используя ближние трассы, путем вычисления взаимного спектра мощности или автоспектра мощности при временном окне вокруг первых вступлений. Корень квадратный из спектра мощности используют в качестве амплитуды сейсмического импульса. Фазовый спектр эквивалентен минимальной фазе, при которой амплитудный спектр может быть вычислен стандартными способами. Один способ получения фазы вычислением представляет собой метод преобразования Гилберта логарифма амплитудного спектра сейсмического импульса.
Упомянутые выше неизвестные составляющие представляют в виде функций частоты. Значения этих составляющих оценивают для каждой частоты в пределах всех частот поверхностной волны-помехи. Поскольку полоса частот поверхностной волны-помехи ограничена со стороны низких частот, частотный диапазон, который необходимо охватывать, является небольшим, например от 3 до 25 кГц. Один путь нахождения неизвестных составляющих заключается в независимом нахождении отдельных значений частот, то есть одной частоты за один раз. В этом осуществлении изобретения индивидуальными параметрами на каждой частоте являются выделенные амплитуда и фаза источника, амплитуда и фаза связи приемника с грунтом и действительная и мнимая части фазовой медленности. Проблема, связанная с этим подходом, заключается в том, что решение для одной частоты может в значительной степени отличаться от решения для следующей частоты, особенно на краях полосы, где амплитуды являются небольшими. Предпочтительный способ заключается в наложении ограничения на этапы оптимизации с тем, чтобы составляющие были гладкой функцией частоты. Один путь осуществления этого заключается не в нахождении отдельных значений частот, а коэффициентов полинома или кривой сплайна, которые представляют составляющие (см. ранее упомянутые статьи Ernst et al.). В этом случае инверсию выполняют одновременно в диапазоне частот. Еще один путь уменьшения количества неизвестных заключается в связывании параметров скорости и затухания при использовании причинно обусловленного ограничения.
Теперь предпочтительный способ выполнения этапа 303 будет описан и детализирован на фиг. 7. В этой реализации входными данными для этапа 303 нахождения решения являются сейсмические трассы с местоположениями источников и приемников и протяженности траекторий через каждую область для каждой трассы. В дополнение к этому как часть этапа 303 получают результаты измерений или приближенные значения для системы регистрации, датчика и импульсные характеристики (сейсмические импульсы) источников. Способ можно использовать только для вертикальных трасс, только для горизонтальных трасс или для обеих. Равным образом также можно получить выражения для датчиков давления или акселерометров.
На этапе 401 из фиг. 7 поверхностные волны ослабляют с помощью временного окна. Цель заключается в ограничении помехи от других волн, таких как отраженные волны или первые вступления. Кроме того, можно ограничить диапазон удалений или используемые трассы. При более значительных удалениях могут теряться высокие частоты, и измерениями можно выявить, что они являются менее полезными. Кроме того, небольшие удаления могут быть бесполезными, поскольку они соответствуют ближнему полю источника. Можно использовать другие способы ослабления или усиления поверхностных волн, такие как фильтрация, взаимная корреляция с опорной трассой, интерферометрические операции взаимной корреляции и суммирования и т.д.
На этапе 402 выполняют преобразование Фурье данных и присоединенных импульсных характеристик в частотную область. Как вариант данные можно оставить во временной области или преобразовать в некоторую другую область. Например, в f-k-область, вейвлет-область, область Габора или область Радона.
-1/2
-π/2
На этапе 403 выполняют деление трасс D(f) данных в частотной области на детерминированные составляющие ((aoffset), S(f), I(f) и R(f), e) и сохранение результата в виде D'(f). Это не является обязательным; цель заключается в исключении умножения на эти члены каждый раз, когда выполняют вычисления по уравнениям (7) и (8) трасс модели.
Теперь уравнения (7) и (8) можно переписать, избавившись от известных параметров
<mrow><msub><msup><mi>M</mi><mo>′</mo></msup><mi>v</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>s</mi><mi>x</mi><mo>,</mo><mi>s</mi><mi>y</mi><mo>,</mo><mi>r</mi><mi>x</mi><mo>,</mo><mi>r</mi><mi>y</mi><mo>,</mo><mover accent="true"><mi>T</mi><mo>→</mo></mover><mo stretchy="false">)</mo><mtext> </mtext><mo>=</mo><mtext> </mtext><msub><mi>C</mi><mi>v</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>r</mi><mi>x</mi><mo>,</mo><mi>r</mi><mi>y</mi><mo stretchy="false">)</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>n</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mrow><msub><mi>P</mi><mrow><mi>n</mi><mo>,</mo><mi>v</mi></mrow></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>s</mi><mi>x</mi><mo>,</mo><mi>s</mi><mi>y</mi><mo stretchy="false">)</mo><msup><mi>e</mi><mrow><mo>−</mo><mi>i</mi><mn>2</mn><mi>π</mi><mi>f</mi><mstyle displaystyle="true"><munderover><mo>∑</mo><mi>j</mi><mrow><mtext>траектория луча</mtext></mrow></munderover><mrow><msub><mi>s</mi><mrow><mi>n</mi><mo>,</mo><mi>f</mi></mrow></msub><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mstyle></mrow></msup></mrow></mstyle></mrow>
(10)
и
<mrow><msub><msup><mi>M</mi><mo>′</mo></msup><mi>h</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>s</mi><mi>x</mi><mo>,</mo><mi>s</mi><mi>y</mi><mo>,</mo><mi>r</mi><mi>x</mi><mo>,</mo><mi>r</mi><mi>y</mi><mo>,</mo><mover accent="true"><mi>T</mi><mo>→</mo></mover><mo stretchy="false">)</mo><mtext> </mtext><mo>=</mo><mtext> </mtext><msub><mi>C</mi><mi>h</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>r</mi><mi>x</mi><mo>,</mo><mi>r</mi><mi>y</mi><mo stretchy="false">)</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>n</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mrow><msub><mi>P</mi><mrow><mi>n</mi><mo>,</mo><mi>h</mi></mrow></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>s</mi><mi>x</mi><mo>,</mo><mi>s</mi><mi>y</mi><mo stretchy="false">)</mo><msup><mi>e</mi><mrow><mo>−</mo><mi>i</mi><mn>2</mn><mi>π</mi><mi>f</mi><mstyle displaystyle="true"><munderover><mo>∑</mo><mi>j</mi><mrow><mtext>траектория луча</mtext></mrow></munderover><mrow><msub><mi>s</mi><mrow><mi>n</mi><mo>,</mo><mi>f</mi></mrow></msub><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mstyle></mrow></msup></mrow></mstyle></mrow>
(11)
0
<mrow><msup><mi>e</mi><mrow><mi>i</mi><mn>2</mn><mi>π</mi><mo stretchy="false">(</mo><mi>o</mi><mi>f</mi><mi>f</mi><mi>s</mi><mi>e</mi><mi>t</mi><mo stretchy="false">)</mo><msub><mi>s</mi><mn>0</mn></msub></mrow></msup></mrow>
На этапе 404 целесообразно дополнительно улучшить данные и параметры. Полезно выполнить уплощение данных, используя единственное значение усредненной медленности s, или обратной величины скорости поверхностной волны-помехи. Данные умножают на фазовый член
, чтобы выполнить уплощение данных, и соответственно корректируют все параметры медленности. Кроме того, целесообразно использовать 10 км в качестве единицы расстояния и для масштабирования амплитуд данных с тем, чтобы все неизвестные параметры имели подобные численные значения, то есть в диапазоне 1-10, так что на этапе оптимизации модели можно лучше осуществлять обновление. Могут быть полезными другие способы регуляризации данных.
Оптимизация модели является линейной, когда имеется одна мода; когда имеются несколько мод она является нелинейной. Как известно специалистам, хорошо знакомым с оптимизацией модели, для нелинейной оптимизации могут требоваться ограничивающие и демпфирующие члены (этап 405) и исходный набор параметров (этап 406). Целесообразно ограничивать амплитуды связи источника и приемника с грунтом положительными числами. Кроме того, чтобы функция стоимости исходно не была чересчур малой, целесообразно нормировать функцию стоимости на каждой частоте делением ее на усредненную среднеквадратическую амплитуду трасс данных. Кроме того, как рассматривалось выше, может быть предпочтительно применять демпфирование или перепараметризацию составляющих как сплайнов, чтобы составляющие плавно изменялись в зависимости от частоты. Веса для нормирования различных удалений или амплитуд трасс также могут быть полезными. Предпочтительно делать более значительные удаления более важными, например, производя деление функций стоимости на корень квадратный из удаления. Моды более высокого порядка являются преобладающими при более значительных удалениях, так что зависимый от моды вес удаления может быть полезным. Если амплитуды трасс при исследовании являются сильно изменчивыми, функции стоимости можно нормировать, используя амплитуды трасс, чтобы сделать более равным вклад всех трасс в решение. Начальные значения (этап 406) можно задавать, используя результаты измерений данных, при этом применяют f-k-анализ, или результаты других измерений, или другую информацию. Если этап оптимизации повторяют, то можно использовать предшествующие значения.
Этап 407 оптимизации выполняют для одной частоты за один раз или для одной группы частот за один раз. Конкретные способы реализации этапа 407 описаны ниже. В конце каждой группы оценивают величину полного вычета (значение функции стоимости), чтобы определить, имеется или не имеется (этап 408) сходимость к решению. Если сходимость отсутствует, новые значения можно выбрать для диапазона удалений, масштабирования, нормирования и т.д. Если получается хорошее решение, оптимизацию выполняют для дополнительных групп частот (этап 409).
Этап 407: решения для одной моды
Когда специалист-практик в данной области техники на основании сейсмических данных способен определить, что имеется только одна мода поверхностной волны-помехи, то задача оптимизации является линейной относительно параметров, а амплитудные и фазовые составляющие можно разделять и находить независимо. Даже в случае, когда имеются несколько мод, может быть целесообразно найти решение для одной моды, чтобы облегчить определение набора исходных параметров, особенно на низких частотах до нижней критической частоты мод более высокого порядка. Когда имеется только одна мода уравнение (8) превращается в
<mrow><msub><msup><mi>M</mi><mo>′</mo></msup><mi>v</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>s</mi><mi>x</mi><mo>,</mo><mi>s</mi><mi>y</mi><mo>,</mo><mi>r</mi><mi>x</mi><mo>,</mo><mi>r</mi><mi>y</mi><mo>,</mo><mover accent="true"><mi>T</mi><mo>→</mo></mover><mo stretchy="false">)</mo><mtext> </mtext><mo>=</mo><mtext> </mtext><msub><mi>C</mi><mi>v</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>r</mi><mi>x</mi><mo>,</mo><mi>r</mi><mi>y</mi><mo stretchy="false">)</mo><msub><mi>P</mi><mi>v</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>s</mi><mi>x</mi><mo>,</mo><mi>s</mi><mi>y</mi><mo stretchy="false">)</mo><msup><mi>e</mi><mrow><mo>−</mo><mi>i</mi><mn>2</mn><mi>π</mi><mi>f</mi><mstyle displaystyle="true"><munderover><mo>∑</mo><mi>j</mi><mrow><mtext>траектория луча</mtext></mrow></munderover><mrow><msub><mtext>s</mtext><mtext>j</mtext></msub><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mstyle></mrow></msup></mrow>
<mrow><mtable><mtr><mtd><mrow /></mtd><mtd><mrow /></mtd></mtr></mtable></mrow>
(12)
Аналогичное выражение можно получить для других составляющих, то есть передаточных функций. Поскольку отсутствует суммирование по модам, амплитудные и фазовые составляющие можно разделять и находить независимо. Используя А для амплитуды и беря логарифм амплитуд, и Φ для фазы, и развертывая фазу, получаем следующее:
<mtable columnalign="left"><mtr><mtd><mi>log</mi><mi>A</mi><mo stretchy="false">(</mo><msub><msup><mi>M</mi><mo>′</mo></msup><mi>v</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>s</mi><mi>x</mi><mo>,</mo><mi>s</mi><mi>y</mi><mo>,</mo><mi>r</mi><mi>x</mi><mo>,</mo><mi>r</mi><mi>y</mi><mo>,</mo><mover accent="true"><mi>T</mi><mo>→</mo></mover><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mo>=</mo><mi>log</mi><mi>A</mi><mo stretchy="false">(</mo><msub><mi>C</mi><mi>v</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>r</mi><mi>x</mi><mo>,</mo><mi>r</mi><mi>y</mi><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mo>+</mo></mtd></mtr><mtr><mtd><mo>+</mo><mi>log</mi><mi>A</mi><mo stretchy="false">(</mo><msub><mi>P</mi><mi>v</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>s</mi><mi>x</mi><mo>,</mo><mi>s</mi><mi>y</mi><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mo>+</mo><mrow><mo>(</mo><mrow><mn>2</mn><mi>π</mi><mi>f</mi><mstyle displaystyle="true"><munderover><mo>∑</mo><mi>j</mi><mrow><mtext>траектория луча</mtext></mrow></munderover><mrow><mi>i</mi><mi>m</mi><mi>a</mi><mi>g</mi><mo stretchy="false">(</mo><msub><mi>s</mi><mi>j</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo><mo stretchy="false">)</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mstyle></mrow><mo>)</mo></mrow><mo>,</mo></mtd></mtr></mtable>
(13)
<mtable columnalign="left"><mtr><mtd><mi>Φ</mi><mo stretchy="false">(</mo><msub><msup><mi>M</mi><mo>′</mo></msup><mi>v</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>s</mi><mi>x</mi><mo>,</mo><mi>s</mi><mi>y</mi><mo>,</mo><mi>r</mi><mi>x</mi><mo>,</mo><mi>r</mi><mi>y</mi><mo>,</mo><mover accent="true"><mi>T</mi><mo>→</mo></mover><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mo>=</mo><mi>Φ</mi><mo stretchy="false">(</mo><msub><mi>C</mi><mi>v</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>r</mi><mi>x</mi><mo>,</mo><mi>r</mi><mi>y</mi><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mo>+</mo></mtd></mtr><mtr><mtd><mo>+</mo><mi>Φ</mi><mo stretchy="false">(</mo><msub><mi>P</mi><mi>v</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>s</mi><mi>x</mi><mo>,</mo><mi>s</mi><mi>y</mi><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mo>−</mo><mrow><mo>(</mo><mrow><mn>2</mn><mi>π</mi><mi>f</mi><mstyle displaystyle="true"><munderover><mo>∑</mo><mi>j</mi><mrow><mtext>траектория луча</mtext></mrow></munderover><mrow><mi>r</mi><mi>e</mi><mi>a</mi><mi>l</mi><mo stretchy="false">(</mo><msub><mi>s</mi><mi>j</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo><mo stretchy="false">)</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mstyle></mrow><mo>)</mo></mrow><mo>.</mo></mtd></mtr></mtable>
<mtext> </mtext>
(14)
Член стоимости для амплитудного и фазового членов становится
<mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mrow><mtext>n трасс</mtext></mrow></munderover><mrow><msup><mrow><mrow><mo>|</mo><mrow><mi>log</mi><mo stretchy="false">(</mo><mi>A</mi><mo stretchy="false">(</mo><msub><mi>D</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>s</mi><mi>x</mi><mo>,</mo><mi>s</mi><mi>y</mi><mo>,</mo><mi>r</mi><mi>x</mi><mo>,</mo><mi>r</mi><mi>y</mi><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mtext> </mtext><mo>−</mo><mtext> </mtext><mi>log</mi><mo stretchy="false">(</mo><mi>A</mi><mo stretchy="false">(</mo><msup><mi>M</mi><mo>′</mo></msup><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>s</mi><mi>x</mi><mo>,</mo><mi>s</mi><mi>y</mi><mo>,</mo><mi>r</mi><mi>x</mi><mo>,</mo><mi>r</mi><mi>y</mi><mo>,</mo><mover accent="true"><mi>T</mi><mo>→</mo></mover><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mo stretchy="false">)</mo></mrow><mo>|</mo></mrow></mrow><mn>2</mn></msup><mtext> </mtext><mo>=</mo><mtext> </mtext><mi>min</mi><mo>,</mo></mrow></mstyle></mrow>
<mrow />
(15)
<mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mrow><mtext>n трасс</mtext></mrow></munderover><mrow><msup><mrow><mrow><mo>|</mo><mrow><mi>Φ</mi><mo stretchy="false">(</mo><mi>A</mi><mo stretchy="false">(</mo><msub><mi>D</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>s</mi><mi>x</mi><mo>,</mo><mi>s</mi><mi>y</mi><mo>,</mo><mi>r</mi><mi>x</mi><mo>,</mo><mi>r</mi><mi>y</mi><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mtext> </mtext><mo>−</mo><mtext> </mtext><mi>Φ</mi><mo stretchy="false">(</mo><mi>A</mi><mo stretchy="false">(</mo><msup><mi>M</mi><mo>′</mo></msup><mo stretchy="false">(</mo><mi>f</mi><mo>,</mo><mi>s</mi><mi>x</mi><mo>,</mo><mi>s</mi><mi>y</mi><mo>,</mo><mi>r</mi><mi>x</mi><mo>,</mo><mi>r</mi><mi>y</mi><mo>,</mo><mover accent="true"><mi>T</mi><mo>→</mo></mover><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mo stretchy="false">)</mo></mrow><mo>|</mo></mrow></mrow><mn>2</mn></msup><mtext> </mtext><mo>=</mo><mtext> </mtext><mi>min</mi><mo>.</mo></mrow></mstyle></mrow>
<mrow />
(16)
Решение можно получать, используя стандартные линейные способы, такие как обращение матрицы, методы сопряженных градиентов или разложение по особым значениям. При наличии только одной области и в пренебрежении членами связи приемников с грунтом линейную регрессию можно использовать для нахождения амплитудных и фазовых параметров источника и одной области.
Этап 407: решения для многочисленных мод
<mrow><mover accent="true"><mi>T</mi><mo>→</mo></mover><mo stretchy="false">(</mo><mi>f</mi><mo stretchy="false">)</mo></mrow>
Когда имеются несколько мод, оптимизация является нелинейной. Амплитудные и фазовые члены связаны и не могут быть разделены. В этом случае можно предположить существование конкретного количества мод (двух или большего количества) и использовать полную функцию стоимости, основанную на действительной и мнимой частях данных, задаваемых уравнением (6). Задача заключается в нахождении решения для полного многомодового вектора
. Поскольку задача является нелинейной, требуется код нелинейной оптимизации, такой как метод Левенберга-Марквардта, скорейшего спуска или Ньютона. Метод Левенберга-Марквардта может быть предпочтительным; он является быстрым и устойчивым. Он является особенно полезным для получения аналитических форм производных функции стоимости или якобиана. Эти методы являются итерационными. Для них требуется исходная модель, а решением итерационно обновляют модель для улучшения согласования с данными. Все из этих итераций находятся в рамках этапа 303 на блок-схеме последовательности действий из фиг. 6. Может быть полезно начинать с небольшого количества мод, находить решение для небольшого количества мод, затем переходить (306) обратно к началу, использовать параметры, ранее определенные, и добавлять новые параметры для мод более высокого порядка при следующем решении на этапе 303. При необходимости весь процесс можно повторять много раз для достижения наилучшего согласования со всеми модами поверхностной волны.
ПРИМЕРЫ
В этой заявке показана пара примеров, относящихся к оцениванию кривых дисперсии скоростей и затухания. Они являются примерами полученных моделированием на компьютере данных с использованием полного решения трехмерного вязкоупругого волнового уравнения для горизонтально залегающих слоев. Моделирование было выполнено с использованием сейсмического импульса Рикера на источнике частотой 20 Гц. Поскольку модель представляет горизонтально залегающие слои, изменение типов почвы в поперечном направлении отсутствует, поэтому имеется только одна область, для которой требуется только один источник. Кроме того, отсутствуют изменения связи приемника с грунтом. В этом случае при параметризации сокращаются члены распространения и источника для каждой моды. Это означает, что две составляющие (также называемые передаточными функциями или фильтрами) будут использоваться в этих примерах, при этом обе являются комплексными величинами. Однако в случае поперечной неоднородности в предпочтительном способе будут использоваться многочисленные источники.
Первым примером является одиночный слой толщиной 5 м в пределах полупространства. Имеется только одна мода поверхностной волны, которая является дисперсионной, показанная на фиг. 8 для четырех расстояний (удалений) источник-приемник, находящихся в диапазоне от 500 м (81) до 2000 м (82) с шагом 500 м. Представлено несколько трасс из входных сейсмических данных для этапа 301. Выходные данные этапа 303 оптимизации модели показаны на фигурах 9A-9D. Они включают в себя распределение энергии источника (то есть часть энергии источника, которая переходит в поверхностные волны, а не в объемные волны, которые создают сейсмические отражения) в амплитуду (9А) и фазу (9В), и действительную (9С) и мнимую (9D) части медленности. Каждый из этих параметров является функцией частоты. Кроме того, можно заметить, что параметром распределения энергии источника является не сейсмический импульс Рикера частотой 20 Гц, а поправка за распределение энергии в моду поверхностной волны-помехи. Прогнозируемые волновые сигналы поверхностной волны-помехи вычислены с использованием оптимизированных параметров из этапа 303 и показаны на фиг. 10. Их можно сравнить с входными данными из фиг. 8. Они полностью согласованы с частью поверхностной волны из данных.
Вторым примером является результат моделирования сложной верхней части разреза, составленного из 20 слоев, с большим градиентом скорости. Волновые сигналы (фиг. 11А) свидетельствуют о сложной интерференции многочисленных мод, а f-k-спектры (фиг. 11В), полученные в результате двумерного преобразования Фурье данных, указывают на присутствие до 6 мод наряду с дисперсией скоростей и поведением амплитуд. Самая медленная мода или основная мода 801 начинается при 0 Гц, но последующие моды имеют возрастающие нижние критические частоты. Нижняя критическая частота для третьей моды обозначена 802. Скорости и амплитуды были выбраны по f-k-спектрам и использованы в качестве исходных моделей. Для первой моды использовался линейный способ, а затем последующие решения были получены нелинейным способом с добавлением двух мод одновременно.
На фиг. 12 показаны подлинные трассы данных для пяти удалений в диапазоне от 100 м (верхняя) до 500 м (нижняя) равными шагами по 100 м, а на фиг. 13 показаны полученные прогнозированием волновые сигналы для пяти трасс. Они находятся в хорошем соответствии. Волновыми сигналами иллюстрируется сложность волновых сигналов и интерференция многочисленных мод. Модельное выражение с шестью модами может воспроизвести всю сложность волновых сигналов. Для точного представления волновых сигналов требуется нахождение многочисленных мод. На фигурах 14A-F показано для сравнения влияние включения одной, двух, трех, четырех и шести мод при вычислении прогнозируемых волновых сигналов. На фиг. 14А показана одна трасса данных, а на фигурах 14B-F показаны результаты прогнозирования этой трассы данных с использованием одной (14В), двух (14С), трех (14D), четырех (14Е) и шести (14F) мод. При включении только одной моды согласуется задняя, медленная часть волнового сигнала. При четырех модах имеется хорошее соответствии, а при шести модах очень хорошее соответствие. Наконец, на фиг. 15 показаны выходные данные этапа 307, кривые дисперсии скоростей для шести мод поверхностной волны-помехи. Эти кривые могут быть входными данными этапа 308 для определения характеристик верхней части разреза.
В изложенной выше заявке для иллюстрации приведены конкретные осуществления настоящего изобретения. Однако специалисту в данной области техники должно быть понятно, что возможны многочисленные модификации и варианты осуществлений, описанных в этой заявке. Все такие модификации и варианты предполагаются находящимися в объеме настоящего изобретения, определяемом прилагаемой формулой изобретения.