Результаты расчета всех гармоник

Номер гармоники k	Длина волны L / k	Коэффициент	Ампли-туда A	Дисперсия отклонений D (o, k)	Дисперсия гармоник
аk	bk	D (k)	D (k),%
	–	3,140	0,000	0,000	2,790	0,000	0,0
	24,0	0,335	–0,058	0,340	2,734	0,056	2,0
	12,0	–0,376	0,054	0,380	2,665	0,070	2,5
	8,0	–1,828	0,535	1,905	0,912	1,753	62,8
	5,0	–0,154	–0,070	0,169	0,898	0,014	0,5
	4,8	0,084	0,040	0,093	0,894	0,004	0,2
	4,0	0,597	–0,883	1,066	0,317	0,577	20,7
	3,4	0,469	0,042	0,471	0,210	0,107	3,8
	3,0	0,241	0,178	0,300	0,166	0,044	1,6
	2,7	–0,041	0,128	0,135	0,156	0,009	0,3
	2,4	–0,146	0,263	0,301	0,110	0,046	1,7
	2,2	–0,243	0,326	0,406	0,027	0,084	3,0
	2,0	–0,088	0,000	0,088	0,027	0,000	0,0
Сумма	–	–	–	–	–	2,764	99,1

Коэффициенты a_k и b_k позволяют рассчитать амплитуды всех гармоник по формуле (5.16), а отклонения – их дисперсии по формуле (5.22).

Результаты расчетов по всем 12 гармоникам суммированы в табл.5.8. Отметим, что отношение L / k – это длина волны гармоники k. В таблице приведены длина волны, коэффициенты a_k и b_k, амплитуда А, дисперсии отклонений D (o, k), их разности – дисперсии гармоник D (k) и они же, выраженные в процентах от дисперсии пространственной переменной.

Вычисления показали, что главное значение имеют третья гармоника с длиной волны 8 м (рис.5.12, б), ее доля составляет 62,8 %, и шестая гармоника с длиной волны 4 м (рис.5.12, в), ее доля 20,7 %. Остальные гармоники играют малую роль.

Следовательно, можно ограничиться двумя важнейшими гармониками, которые в сумме дают 83,5 % от дисперсии пространственной переменной (рис.5.12, г).

Данные табл.5.8 позволяют построить спектр амплитуд (рис.5.13).7

 

 

5.4.3. Периодограммный анализ

 

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

Таблица 5.9 Расчет периодической изменчивости и отклонений при длине волны L = 5 м Номер столбца           Исходные данные 2,36 2,08 1,95 7,68 5,72 3,08 2,20 2,76 1,32 1,02 3,24 4,90 4,48 2,30 1,97 2,27 2,08 1,45 3,85 5,18 6,34 3,93 2,33 2,11 1,90 Средние по столбцам (периодическая изменчивость) 3,46 3,04 2,59 3,45 3,16 Отклонения от средних –1,10 –0,96 –0,64 4,23 2,56 –0,38 –0,84 0,17 –2,13 –2,14 –0,22 1,86 1,89 –1,15 –1,19 –1,19 –0,96 –1,14 0,40 2,02 2,88 0,89 –0,26 –1,34 –1,26 ________________ Дисперсия исходных данных 2,790 Дисперсия отклонений 2,689 Дисперсия периодической изменчивости 0,101

Если имеется ряд значений пространственной переменной f (x), измеренных с шагом h, имеющих дисперсию D, то, задавая длину волны, можно выявить периодическую изменчивость w(х) при данной длине волны и вычесть ее из исходных данных, что дает отклонение d(х) = f (x) – w(х). Разность дисперсии исходных данных и дисперсии отклонений дает дисперсию, поглощенную периодической изменчивостью. Чем больше поглощенная дисперсия, тем сильнее проявлена периодическая изменчивость. Задавая различную длину волны, можно найти наибольшую дисперсию, что и лежит в основе периодограммного анализа.

Отыскание наилучшей длины волны начинается с наименьшей длины волны L = h. Затем длина волны последовательно увеличивается до достижения половины ряда наблюдений N = n /2.

Таблица 5.10   Дисперсия различных длин волн   Длина волны, м Дисперсия отклонений Дисперсия периодической изменчивости D D, %   2,790 0,000 0,0   2,783 0,007 0,2   2,759 0,031 1,1   2,279 0,511 18,3   2,689 0,101 3,6   2,702 0,088 3,2   1,724 1,066 38,2   0,445 2,345 85,1   1,436 1,354 48,5   2,367 0,423 15,2   2,148 0,642 23,0   2,073 0,717 25,7

8 Пример 5.10. Измерения пространственной переменной приведены в табл.5.7. Необходимо рассчитать периодическую изменчивость ряда значений путем перебора различных длин волн.

Покажем порядок расчета при длине волны L = 5 м. Исходные данные «разрежем» на части по пять измерений и запишем один под другим в пять столбцов (табл.5.9). Средние по столбцам – это периодическая изменчивость при длине волны L = 5 м. После вычитания из исходных данных периодической изменчивости получим отклонения. Далее рассчитаем их дисперсию и вычтем из дисперсии исходных данных.

Подобные расчеты выполнены для всех длин волн от L = 2 до L = 12 (табл.5.10). Наибольшая дисперсия периодической изменчивости 85,1 % установлена для длины волны 8 м.

Наиболее сильно выраженная волна периодической изменчивости, показанная на рис.5.14, имеет следующий вид:

{1,91 1,52 3,01 5,92 5,51 3,10 2,17 2,38}.

Можно совместить график исходных данных с графиком периодической изменчивости и убедиться в удовлетворительном их совпадении (рис.5.15). Следует отметить, что дисперсия периодической изменчивости в данном случае несколько больше, чем при гармоническом анализе с использованием двух важнейших гармоник.7

 

 

ОСНОВЫ ГЕОСТАТИСТИКИ

5.5.1. Вариограмма и ее аппроксимации

 

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

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

γ(h) + K (h) = D. (5.25)

График вариограммы зависит от характера дискретности пространственной переменной. Для непрерывных пространственных переменных (например, мощности пласта) вариограмма начинается с нулевой отметки и возрастает до дисперсии исходных данных (рис.5.16). Для дискретных пространственных переменных вариограмма начинается с некоторой величины С (рис.5.17), называемой эффектом самородков, потому что он был вначале установлен на месторождениях золота.

Вариограмма имеет радиус влияния R, который идентичен радиусу автокорреляции в случайной функции. За пределами радиуса влияния вариограмма постоянная и равна дисперсии D, т.е. ее влияние отсутствует. Как упоминалось, радиус R является векторной величиной. В изотропных геологических телах по всем направлениям радиус влияния описывает окружность, за пределами которой влияние вариограммы отсутствует. В анизотропных геологических телах значения радиуса R по разным направлениям столь различны, что напоминают восьмерку.

В пункте с координатой R вариограмма может иметь горизонтальную касательную и плавно переходить в линию дисперсии.

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

Первая модель линейная. В ней вариограмма аппроксимируется прямой линией (рис.5.19), что нередко близко к действительности. Во второй модели вариограмма аппроксимируется кубическим полиномом, который не имеет общей касательной с линией дисперсии (рис.5.20). Эта модель лучше всего соответствует действительности. Третья модель (сферическая) также аппроксимируется кубическим полиномом, но он имеет горизонтальную касательную при соприкосновении с линией дисперсии (см. рис.5.16 или рис.5.17). Четвертая модель (модель Де-Вийса) характеризуется тем, что шаг вариограммы откладывается в логарифмическом масштабе, а сама вариограмма представляется в виде прямой линии (рис.5.21).

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

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

Таблица 5.11 Содержание золота в пробах, г/т Пункт измерений x, м Содержание золота в пробах f (x), г/т Пункт измерений x, м Содержание золота в пробах f (x), г/т   5,2   4,6   3,5   4,7   4,6   6,7   5,2   6,8   6,8   4,9   6,1   3,0   7,0   3,5   7,4   3,8   7,8   5,4   7,6   4,1   6,2   4,0   5,2   5,5   2,6   5,1

 

8 Пример 5.11. В горной выработке определено содержание зо­лота через 1 м (табл.5.11, рис.5.22). Требуется построить вариограмму содержаний и аппроксимировать ее какой-либо моделью.

Эмпирическая вариограмма, вычисленная по формуле (5.4), показана на рис.5.18. Она содержит 13 точек, т.е. рассчитана до 12-го шага (половина длины ряда наблюдений).

Начальная часть ва­риограммы (первые четыре точки) расположена практически на одной линии, поэтому для построения теоретической вариограммы используем линейную модель, как на рис.5.19. Уравнение теоретической вариограммы, рассчитанное по четырем точкам, γ(h) = 0,8 h. Радиус автокорреляции R = 2,5 м получен при пересечении линии теоретической вариограммы с линией дисперсии исходных данных D = 2,03.

 

5.5.2. Влияние на вариограмму геометрической

базы измерений

 

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

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

Пусть в геологическом поле имеется множество точек, в каждой из которых пространственная переменная имеет значение f (x). Если геологическое поле разделить на области объемом v и в каждом из них найти среднее из N точечных значений, то получим новую пространственную переменную f_v (x):

при , (5.26)

где p (m) – весовая функция, зависящая от взаимного расположения точек в объеме v.

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

при . (5.27)

Нахождение средних значений f_v (x) по точечным значениям f (x) называется регуляризацией (сглаживанием) пространственной переменной в объеме v. Меняя объем v, будем получать различные регуляризованные пространственные переменные, соответственно, будут меняться и их характеристики.

Изменение дисперсии по мере увеличения объема v описывается формулой

, (5.28)

где D – дисперсия пространственной переменной на точечной базе; N – число точек в объеме v.

Вычитаемое в формуле (5.28) представляет собою среднее значение вариограммы, получаемое при всех возможных положениях точек i и j в объеме v. При переходе к бесконечному множеству точек формула (5.28) приобретает интегральный вид:

. (5.29)

Из формул (5.28) и (5.29) следует, что по мере увеличения объема v дисперсия пространственной переменной уменьшается.

Подобные формулы существуют и для вариограммы:

(5.30)

или в интегральной форме

(5.31)

Первое слагаемое в этих формулах характеризует среднюю ковариограмму, которая получается как средняя из вариограмм между всеми точками в объеме v, а вычитаемое – в таком же объеме, смещенном на расстояние h.

Формулы (5.28)-(5.31) выведены в предположении, что объем v мал по сравнению с объемом изучаемого геологического объекта, иначе сказывается граничный эффект, который можно учесть введением специальной весовой функции – геометрической ковариограммы p (h):

, (5.32)

где k (х) = 1, когда точка х находится внутри геологического объекта, k (х) = 0, когда точка х находится за его пределами.

Геометрический смысл ковариограммы заключается в степени перекрытия двух объемов, смещенных относительно друг друга на расстояние h. С учетом функции p (h) формулы (5.29) и (5.31) приобретают следующий вид:

; (5.33)

. (5.34)

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

 

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

Будем объединять и усреднять исходные данные по два, три и четыре соседних измерения путем сглаживания методом скользящего окна, что аналогично увеличению размеров проб исходных данных v. Количество исходных данных будет уменьшаться. Рассчитаем дисперсию D усредненных данных по формуле (2.2) и вариограмму g(h) по формуле (5.4) по каждому варианту усреднения. Результаты расчета изображены на рис.5.23 и 5.24. На рисунках видно, что при увеличении размера проб дисперсия уменьшается, по сравнению с рис.5.18 убывает также и вариограмма.7

 

 

5.5.3. Понятие о кригинге

 

На основе геостатистических моделей создан новый метод интерполяции результатов между пунктами измерений, который получил название кригинг (крайгинг)в честь автора метода Д.П.Крига [7]. Как показано Д.Матероном [12], кригинг дает минимальную сумму квадратов отклонений прогнозных значений от фактических значений пространственной переменной. Существует несколько видов кригинга. Вначале рассмотрим простейший – точечный кригинг.

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

Необходимо рассчитать прогнозное значение пространственной переменной в точке Z. Предположим, что геологический объект изотропный, тогда вокруг точки Z можно провести окружность с радиусом влияния R. Все измерения внутри окружности влияют на значение пространственной переменной в точке Z. За пределами окружности результаты измерений в прогнозировании не участвуют.

Прогнозное значение находится в результате выполнения ряда последовательных операций:

1. Рассчитывается эмпирическая вариограмма, как на рис.5.18.

2. Осуществляется аппроксимация эмпирической вариограммы теоретической вариограммой (см. рис.5.19-5.21), т.е. выбирается геостатистическая модель.

3. В области с радиусом R рассчитываются все расстояния между пунктами измерений, а также между пунктами измерений и пунктом прогноза.

4. Составляется симметричная система уравнений кригинга следующего вида:

. (5.35)

Здесь g₁₂ – значение теоретической вариограммы между пунктами 1 и 2 (учитывается расстояние между этими пунктами); m – количество пунктов, участвующих в расчете (в данном случае m = 5).

5. Решение системы дает весовые коэффициенты p ₁, p ₂, p ₃, …, p_m.

6. Сумма весовых коэффициентов должна быть равна единице, для чего все коэффициенты делятся на их сумму.

7. С помощью весовых коэффициентов вычисляется прогнозное значение свойства в пункте Z:

, (5.36)

где z_i – значения свойства в пунктах c 1-го по m -й.

Если понадобится вычислить прогнозное значение в следующем пункте, то операции 3-6 придется выполнить снова. Если в пределах радиуса влияния не окажется ни одного пункта измерения, то кригинг нельзя применять. Приходится обращаться к другим методам интерполяции данных, чаще всего используется метод обратных расстояний.

Иногда при составлении системы уравнений (5.35) в матрице коэффициентов оказывается много нулей (расстояния между точками больше радиуса автокорреляции). В этом случае система уравнений является несовместимой, т.е. не имеет решения (или имеет много решений). Тогда весовые коэффициенты лучше всего находить путем деления значений D – γ_{1 z} , D – γ_{2 z} , D – γ_{3 z} , ¼, D – γ _mz на дисперсию D и потом приводить сумму коэффициентов к единице.

У точечного кригинга имеется один недостаток – в окрестностях пунктов измерений касательная к прогнозным значениям ориентирована горизонтально, как в полиномиальной модели (см. рис.5.2). Чтобы устранить этот недостаток, применяется универсальный кригинг – комбинация тренда с кригингом. Вначале вычисляется тренд f (x), потом остаток от тренда и по остатку осуществляется кригинг. Прогнозное значение находится по формуле

(5.37)

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

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

 

8 Пример 5.13. Нужно провести точечный кригинг, используя схему рис.5.25 и дополнив недостающие данные.

Примем изотропную линейную модель, радиус влияния R = 14 м, дисперсию исходных данных D = 0,42. Исходные данные приведены в табл.5.12. Вариограмма для такой модели имеет следующий вид: до 14 м и после 14 м соответственно

γ(h) = D / R = 0,03 h;	(5.3)
γ(h) = D = 0,420.

Таблица 5.12

Данные для кригинга, м

Номер пункта	х	у	z
			3,3
			3,1
			3,0
			2,6
			2,1
			3,0
			2,8
			2,5
			2,6
			2,0
Z			?