Новые подходы в теории геостатистического моделирования

УПРАВЛЕНИЕ, ВЫЧИСЛИТЕЛЬНАЯ ТЕХНИКА И ИНФОРМАТИКА • МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ...

УДК 551

В. А. БАНКОВ,Н.К. БАКИРОВ,А. А. ЯКОВЛЕВ

НОВЫЕ ПОДХОДЫ В ТЕОРИИ ГЕОСТАТИСТИЧЕСКОГО МОДЕЛИРОВАНИЯ

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

Настоящая работа посвящена методам построения трехмерных геологических геостохастических моделей объемного распределения тех или иных геофизических параметров, то есть их интерполяции по скважинным данным. Неоднозначность, неполнота и, как часто бывает на практике, противоречивость входных и рассчитанных по ним производных данных, а также сложность построения адекватной математической модели пласта не позволяют нам однозначно строить геологическую модель. Таким образом, мы невольно становимся заложниками случайности. Отдельная ветвь математики - теория вероятностей, создана с целью ее (случайность) «приручить». Работы Жоржа Матерона из Центра геостатистики Парижской горной школы и его учеников создали прикладные методы геостохастического моделирования, которые легли в основу современных пакетов геологического конструирования [1-3]. Основные гипотезы геостохастического моделирования - это стационарность и изотропность (а также их ослабления: тренд плюс стационарный остаток и геометрическая анизотропия), эргодичность и гауссовость случайного поля. Однако ни стационарность, ни изотропность в своих строгих формах почти не встречаются при работе с реальными данными. Более того, зачастую не помогают и предложенные ослабления этих ограничений. Отметим также, что на практике широко распространен субъективный подход к определению необходимых для построения геостохастиче-ской реализации пласта параметров (радиусов вариограмм, порога, эффекта самородков...). Субъективизм выбора параметров, часто не

Контактная информация: (347) 273-77-35

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

нии (ОСМ) данного геофизического параметра в пределах трехмерного «куба», содержащего залежь полезного ископаемого, посредством распространения кривых каротажа. В итоге мы должны построить «равновероятные» реализации случайного геофизического поля:

а) совпадающие с эмпирическими значениями в пределах скважин;

б) с вероятностными характеристиками

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

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

Методы разрешения проблемы интерполяции можно разделить на детерменированные и вероятностные. Целью детерменированной модели является построение поверхности интерполируемой переменной максимально близкой к неизвестной величине. Одним из таких методов является кригинг (в честь Баше К^е, [4]), который заключается в нахождении наи-

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

о геологической системе. Целью вероятностной модели ОСМ является предсказание истинного характера неоднородностей между скважинами, при этом требуется воспроизведение статистических параметров поля (среднее, дисперсия, вариограмма) и совмещение значений построенного поля с его известными значениями в известных точках. Основными способами обусловленного стохастического моделирования являются: 1) последовательное гауссовское моделирование [5]; 2) метод моделирования с использованием цепей Маркова; 3) моделирование с использованием спектрального представления случайной переменной, которое является также ключевым методом цифровой обработки сигналов, статистической физики, теории турбулентности, оптики и т. д. (см., например, [6]). В данной работе при моделировании физических полей также, по существу, используется спектральное представление стационарных случайных полей.

1. СОСТОЯНИЕ ВОПРОСА

Пусть 0(1), ^ е Я - случайное поле геофизического параметра в пределах некоторого «куба» Т. Изучение большого количества каротажных диаграмм позволяет заключить, что зачастую они не обладают свойствами стационарности и гауссовости, что препятствует применению стандартной идеологии геостатистики и говорит о необходимости поиска новых, более универсальных моделей. Невозможность повторения построения модели и ограниченность предложенных моделей ковариационных функций говорит о необходимости использования непараметрического анализа.

2. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ

Пусть /к(к) - нормированная кривая каротажа по к-й скважине в отбитом маркерами пласте, переведенная линейным преобразованием на интервал [-1, 1], здесь к - вертикальная координата переменной ^ = (иь и2, к). Выявление характерной формы кривой каротажа осуществим посредством разложения функции /к(к) по базису многочленов Лежандра Р(к), ортонормированных на [-1, 1] с единичным весом:

/к

где Ак(к) - некоторый остаток,

5к,п (к) = £ акР* (к),

1 ^ (2) ак = | /к (кр (к)М.

-1

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

весом Р(к), то есть величину ак* в широких и

достаточно слабых предположениях можно приближенно считать гауссовским (в силу центральной предельной теоремы теории вероятностей). Без ограничения общности «куб» Т в вертикальных сечениях есть отрезок [-1, 1]. Рассмотрим векторный случайный процесс

Х( х) = (хо( х) Х1( х),..4-1( х))Т,

1

X * (х) = | Х( X, к)Р} (И)с1к,

-1

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

Х(хк ) =(аок, а1к ,..-аП-1,Т.

Ввиду центральной предельной теоремы векторное поле £,(х) считаем гауссовским и, ввиду геологических представлений о плавном характере осадконакопления в латеральном направлении, считаем поле £,(х) квазистацио-нарным, то есть

Х( х) = а( х) + с( х)^( х),

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

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

информацию о сигнале /к

п-1

0(Г) = X (а* (х) + а * (х)л * (х)) х

1=0

х Р1 (к) + А( х, к),

где п(х) = (По(х), П1(х),--, Пп-1(х))Т - векторный, стационарный, гауссовский случайный процесс с нулевым средним, а(х) = (а0(х), а^х),..., ап-1(х))Т - вектор средних, а(х) = diag(o0(x), о^х),..., оп-1(х)) - масштабирующая функция, А(х, к) - высокочастотный шум малой мощности, ортогональный функциям Р(к), * =

= 0,1,..., п-1.

Ниже мы рассматриваем возможные процедуры оценивания всех составляющих элементов модели (3).

2.1. Оценивание среднего а

Функцию а(х) можно оценить сглаживанием соседних значений поля £(х). Обозначим а( х) - сглаженную оценку среднего по методу ближайших т соседей:

XХ(х,)| шах|х-х,|- х-х

Л /-1 ^ 1' 1 г=1, ,т г| I

а( х) =

' ' г

X

*=1

шах х -

г=1,..,т

х - х,-

здесь х1, х2,...,хт - наиболее близкие к х (в евклидовой метрике |.|) экспериментальные точки - точки скважин, т, V > 0 - варьируемые параметры. Ясно, что за счет конструкции среднего, функция а(х) будет непрерывна по х. В тех областях, где мало или нет экспериментальных точек хк, возможно пользоваться теми или иными математическими моделями, выражающими представления геолога о характере убывания (возрастания) функции среднего на конкретном месторождении. Оценивание масштабирующих функций оДх) производится усреднением локальной оценки вариации поля П/х) по формуле, аналогичной (4).

2.2. Обусловленный кригинг гауссовского поля на основе спектрального представления для стационарных случайных функций

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

ции геофизического поля Z(u) (см., например, [7]). В частности, интерполяционная формула кригинга поля Z(u) с нулевым средним представляет собой линейную форму

N

г *(и) = X1 ^ (ик),

к=1

где г(ик) - значения случайного поля Х(и) в экспериментальных точках, а константы («веса») ^к = Ми) подбираются по методу наименьших квадратов, исходя из минимума среднекавадра-тичного уклонения линейной формы г (и) от случайного поля 2(и):

Е(г(и) - г*(и ))2

Обозначим через С(и) ковариационную функцию случайного поля г(и) и

1 = (1,...., 1 „ ) ,

с(и) = (С(и1 - и),С(и2 - и),...,С(uN - и)) ,

I = (г (идг (и 2),..., z(uN ))Т,

С = | С (и) - С (и* )|| ,,1 .

II 3 Иг, 1=1,..,N

Тогда оптимальное значение 1 задается формулой [7]:

1 = С-1с(и).

Для гауссовских случайных полей таким образом построенная формула кригинга имеет простое важное свойство: случайное поле невязки ВД = г(и) - г (и) и данные г(ик), к = = 1,...^ независимы, что открывает возможность моделирования поля ¥(и) независимо от исходных данных. Отметим также известный факт, состоящий в том, что вместо точного представления

г

в результате моделирования мы получаем реализацию

г1(и) = у1(и) + г * (и),

где поле %(и) есть независимая копия поля ¥(и). Получающаяся неопределенность в оценивании поля ¥(и) отражает имеющуюся неопределенность в определении рассматриваемого геофизического параметра. Все эти рассуждения сохраняют силу и в случае многомерного поля г(и).

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

V

V

х

моделирование, воспроизводящее исходные статистические данные (математическое ожидание, ковариацию) и значения оцениваемого поля в известных точках.

Используем для моделирования невязки 5(х) = п(х) - П (х) и построения функции кри-гинга п (х) спектральные представления для стационарных случайных процессов и полей с нулевым средним:

I ег(хЛ}ЛФ(1),

Я2

С

|^Ф(1)|2 = ЛБ (1),

где Ф(1), 1е Я2 - некоторое векторозначное комплекснозначное случайное поле с ортогональными приращениями и нулевым средним, см., например, [6]. В случае гауссовости поля п(х) случайное поле Ф(Х) также будет гауссовским, [6], С(и) - матрица ковариаций поля п(х), Б(Х) - матричная спектральная функция распределения.

Итак, в соответствии с (5),

8( х) = л( х) -л*( х) =

.К хкЯ)

Л

dф(1).

Ясно, что в точках х = хк имеем, что 5(хк) = = 0, поскольку в (6) подынтегральное выражение обращается в нуль (все X* в этих точках за исключением Хк равны нулю, а Хк = 1).

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

8(х) = | ( еі(хЛ) - £ 1 кеі(ХкЛ) ]^Ф1(1),

к=1

Л( х) = 8( х) + л*( х).

Ясно, что поле ^(х) является обусловленным и имеет такие же вероятностные характеристики, что и поле п(х), но, вообще говоря, с ним не совпадает. Выбирая различные варианты Ф^Х) поля Ф(Х), соответственно, получим

различные варианты поля ^(х).

Наконец, компьютерное моделирование поля 5( х) производим, заменяя интеграл в (7) соответствующими интегральными суммами:

гК хк,т ])

jеJ с г 2 V к=1

где АФ* = иг + г'У, * е J - гауссовские комплекснозначные вектора с нулевыми средними и

АФ-* = АФ *, Еи* = ЕУ* = 0,

с^ои у) = 0, "г, *

С^(Ц,и*) = 0,"г Ф *,-*',

С^у у) = 0,"г Ф *,-],

с^(и* ,и*) = с^у. ,у)=1 Г ЛБ (1),

2 В

где В* - элементы разбиения, т* е В*.

Легко видеть, что вычисления по рассматриваемым формулам принципиально распараллеливаются при программировании. Точность определения интеграла (7) ограничивается только нашими вычислительными возможностями.

2.3. Оценивание спектральной плотности, ковариационной функции и вариограммы

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

При оценивании спектральной плотности одномерного стационарного случайного поля Х(0 с нулевым средним, заданного на области В с Я2, традиционно используется периодограмма:

!п (1) =

1

I е

В

-і

X (ґ

где |В| - площадь области В. Известно, что если спектральная плотность распределения /(X) непрерывна в точке Х0, то тогда при неограниченном расширении области В

Е1В (10)  /(10).

Известно, что периодограмма 1В(Х) не является состоятельной оценкой спектральной плотности, однако линейные интегральные функционалы от нее состоятельно оценивают аналогичные линейные интегральные функционалы от спектральной плотности.

Непосредственное вычисление периодограммы в нашем случае невозможно, поэтому мы приближаем векторную периодограмму

2

2

Р

соответствующей интегральном суммой следующего вида:

ІВ (1) =

1 ( N

4 л2 ІВ

X

N

Е л( хк) Iе

к=1 Вк

—V

-і ( х,/

\

їх

X

Е Л(хк) Iе-(х,1)їх

к=1

где Вк, к = 1,...^ - суть некоторые подобласти области В (разбиение Вороного), с хк е Вк . Кроме того, нам требуется следующее усечение периодограммы:

ІВ (1) = ІВ (1), 1є[-10,1] ІВ (1) = о, 1ї[-1 о, 1 ]2,

1

- частота Найквиста, оп-

Л Ш1П хі - хл

і, 1 1 '

ределяемая из следующего соображения: наименьшая длина волны, которую мы можем себе позволить оценивать по имеющимся данным, вдвое больше наименьшего расстояния между скважинами.

Наконец, для того, чтобы частично компенсировать погрешности, возникающие в (8), (9), перенормируем оценку (9), исходя из следующего свойства спектральной плотности:

Г/(1)Л1 = Ел(х)л(х)Т »

Я 2

1 N

XЛ(хкЖхк) \Вк\ = Cо,

В к=1

так что в качестве оценки спектральной плотности можно взять величину:

Л 1 -1 -1 1

/(1)=С02Я 2/В(1)Я2С02,

где

р = | іВ (1)ї1.

Поэтому, очевидно, оценка спектральной плотности будет пронормирована следующим образом:

Л

Г / (1)Л1 =С0.

я 2

Далее мы можем вычислить ковариационную функцию и полувариограмму:

Л Л

С(и) = Г е”г(и,1) / (1)Л1,

Я 2

Л Л Л

у(и) = С(0) - С(и) = С0 - С(и).

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

2.4. Моделирование остатка Л

Как уже было отмечено, остаток А(х, к) в формуле (3) представляет собой некоторую высокочастотную составляющую кривой каротажа (остаток ГИС'а).

Ввиду ортогональности функции А

многочлену Р0 (х) ° 1 имеем:

1 1

Г А( х, к)Лк = 2 Г А( х, к) Р0 (к)Лк =0,

-1 -1

то есть среднее значение функции А(х, к) равно нулю. При этом расчеты показывают, что остаток А(х, к) не коррелирует со значениями 0(0 -А(х, к), в той же скважине и в других скважинах, а также не коррелирует со значениями остатков А(х, к) в других скважинах. Тем самым, мы можем моделировать значения А(х, к) в точках, удаленных друг от друга на расстояния не менее минимального расстояния между скважинами, как независимые случайные процессы. Принимая для остатков модель стационарного гауссовского процесса с нулевым средним, оценим в каждой к-й скважине спектральную плотность сначала периодограммой:

Ік <*•> = ^7

2л/

Ее'^ А к

1=1

2

, И = -1 + 21,

а затем сглаженной периодограммой:

л

/к (1) = |*(/ - г)Ік (1)1г,

где ^(Х) - некоторое сглаживающее финитное ядро, I - число замеров каротажа на скважине. Оценку (10) далее используем для моделирования остатка ГИС'а в выбранной для этого точке.

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

В таких условиях восстановление остатков ГИС'ов между скважинами можно производить

Л

В

к

Л

2

Р

по методу взвешенных расстояний. Опишем схематично соответствующий алгоритм.

С помощью датчика случайных чисел «построим» дополнительные гипотетические

скважины в точках хь / > Ы, заполняющие «куб» из того расчета, что минимальное расстояние между новыми скважинами и старыми скважинами, а также между собственно новыми скважинами будет равно примерно г. Во всех новых скважинах моделируются остатки ГИС'ов Д(хг, И), / = N + 1,...М, как независимые гауссовские стационарные процессы с нулевыми средними и спектральными плотностями

g fk

fr (l) =

k=1

x - x

5

1

x

5

k

и далее остаток ГИС'ов, соответствующий произвольной точке x, моделируется по анало-

M D

5

k=1 x - xk

гичной формуле

A( x, h) = M 1 .

-----F

k=1 |x - Xk\

3. ПРИЛОЖЕНИЕ РЕЗУЛЬТАТОВ

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

Рис. 1. Трехмерная модель пористости

Используемая техника моделирования в будущем даст нам возможность сравнивать различные реализации и по входным данным предсказывать количественную неоднородность пласта, его гидродинамические характеристики. Расчеты в примерах показали: а) согласованность полученных изображений с су-

ществующими геологическими представлениями о строении пачек горных пород;

б) согласованность смоделированных и реальных каротажных диаграмм в «свежих» скважинах, не участвовавших в расчете реализации геофизического поля (метод кросс-валидации);

в) хорошую сходимость средней пористости в пределах скважин и в пределах всей смоделированной реализации.

Рис. 2. Г еологические принципы

На рис. 2 (верхние два и нижние два рисунка) представлены две модели одного и того же геофизического поля с одними входными данными. Верхняя модель построена с использованием предлагаемой технологии, нижняя с использованием известной стандартной технологии. Нетрудно видеть, что реализации новых методов построения удовлетворяют основным принципам геологии (Короновский, [8]) и, по сравнению с известными методами, более реалистичны.

4. ВЫВОДЫ

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

1. Автоматизация. Снижение степени и защита от исследовательского субъективизма при построении геологической модели. Ввиду существенной зависимости результатов расчетов от субъективного выбора интерпретатора (например, при работе с вариограммами) желательно автоматизировать расчеты в соответствующих сегментах общей концептуальной модели. Если полуавтоматические системы обработки данных ориентированы на опытных

А

Л

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

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

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

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

СПИСОК ЛИТЕРАТУРЫ

1. Matheron G. The theory of regionalized variables and its applications. Fontainebelau: Center of Geostatistics, 1971. 212 c.

2. Matheron G. The intrinsic random functions and their applications // Adv. Appl. Probability. 1973. Вып. 5. С. 439-468.

3. Matheron G. The intrinsic random functions and their applications. Fontainebelau: Center of Geostatistics, 1978. 175 c.

4. Krige D. G. A statistical approach to some mine valuations and allied problems at the Witwatersrand // Master's thesis: University of Witwatersrand, 1951.

5. Deutsch C. V., Journel A. G. GSLIB, Geosta-tistical software library and User's guide // New York: Oxford university press, 1992. 340 c.

6. Хеннан Э. Многомерные временные ряды. М.: Мир, 1974. 575 с.

7. Дюбрул О. Использование геостатистики для включения в геологическую модель данных. Изд-во SEG, EAGE, 2002. 296 c.

8. Короновский Н. В. Общая геология. М.: МГУ, 2002. 448 с.

ОБ АВТОРАХ

Байков Виталий Анварович,

зав. каф. матем., зам. ген. дир. по разр. «РН-УфаНИПИнефть». Д-р физ.-мат. наук по матем. физике (ИПМ им. М.В. Келдыша, 1991). Иссл. в обл. физ. и матем. моделирования, группового анализа дифференциальных уравнений.

Бакиров Наиль Кутлужано-вич, внс ИМВЦ УНЦ РАН; снс РН-УфаНИПИнефть. Д-р физ.-мат. наук по теории вероятности и матем. статистике (СПб-ГУ, 1995). Иссл. в обл. теории вероятности и матем. статистики.

Яковлев Андрей Александрович, нс РН УфаНИПИнефть; доц. каф. матем. Канд. физ.-мат. наук по диф. уравнениям (ИМВЦ УНЦ РАН, 2008). Иссл. в обл. операторов на многообразиях и в случайных средах.

й:к