книги / Статистический анализ данных в геологии. Кн. 1

а

(L

Рис. 4.59. Произвольный вход со значениями [2 11 (а); требуемая форма от­ клика [3 0 01 (б); фильтр наименьших квадратов [1,43—0,57] (в); отклик филыра наименьших квадратов [2,9 0,3—0,6] (г)

Умножая слева обе части на [ j 3 j получаем

-о		-	2	, 0 '
	1	0 1	1	2	_ ' / о ’ _			* 2
									(О
0 2		1				1	I	О
			0	1 .			/ 1 ^
		.
				' 5	2	'	/о	■6 '
				2	5	_	А 1.		0

	1 0	' ■	3 '
	1		0
		.	0

Подстановка матрицы левой части дает

• 5	2 ‘ —1	‘ 0,2381	— 0 ,0 9 5 2 ‘
2	5 .	0 ,0 9 5 2	0 ,2 3 8 1 _

Искомые коэффициенты фильтра есть

/о	0,2 3 8 1	— 0 , 9 8 5 2 ]		1,4286'
А х .	- 0 .0 9 5 2	0,2 3 8 1 т	-	- 0 ,5 7 1 2

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

[С] = {2 1]* [1,4286 —0,5712].

Свертка имеет вид

302

г 9 1		- 1,4286	—0,5712'
С\ =	1 ! *	2,8572	— 1,1424
L	j	1,4286	— 0,5712.
= 12,8372	{1,4286— 1,1424»			—0,5712]
= ] 2,4572		0,2862		—0,5712]
Выход фильтра указан на		рис. 4.59; очевидно, мы не преуспели
в создании чистого острия		[3 0 0], но выход является очень
хорошей аппроксимацией.		Действительно,		это есть наилучшая

возможная аппроксимация для этого частного входа и фильтра, состоящего только из двух членов. Теперь мы можем измерить, насколько хорош фильтр, сравнив его выход с тем выходом, ко­ торый мы хотели получить. Это делается вычислением ошибки, определенной как сумма квадратов разностей между С и D:

в= (3—2,8572)2 + (0—0,2862)24- (0 + 0,5712) 2 = 0,1428.

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

Сглаживание и временной тренд-анализ

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

ряда. Это — произвольные фильтры, общее	назначение	кото­
рых— сгладить данную последовательность;	выход из	фильт­

р а —это аппроксимация, тесно связанная с входом. Обоснова­

ние необходимости такого процесса	фильтрации состоит в том,
что временной ряд состоит из двух	компонент — «длиннодейст­

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

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

3 0 3

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

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

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

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

Так как множество коэффициентов должно быть определено

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

оценки Yi в центральных точках сегментов. С этим ограничени­ ем можно переписать уравнения ортогональных многочленов так, чтобы получить альтернативное множество ортогональных чле­ нов, которые прямо обеспечат получение оценок в центральных точках. Детали этого метода приводятся Савицким и Голеем [48J.

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

Таблица 4.30 отражает изложенное выше; она дает ортого­ нальные полиномиальные весовые функции м для нечетною чис­ ла точек в последовательности для порядков кривых 2 и 3, а также для порядков 4 и 5.

Для выполнения сглаживания с использованием ортогональ­ ных многочленов мы просто свернем данную последовательность с множеством полиномиальных членов и затем разделим резуль­ тат на сумму членов, т. е.

	с = в *	[со]
		£<о

304

Так как множество членов симметрично относительно централь­ ного члена, результат одинаков, если мы применим метод сколь­ зящего среднего, которое есть просто скользящее кросспроизве­ дение между множеством членов и временным рядом, деленным на 2(о. Этот процесс, возможно, больше знаком геологам под на­ званием временного тренд-апалпза. Среди уравнений, использу­ емых для указанных выше целей, пятичленный фильтр Шеппар­ да— Тодда, обычно записываемый как скользящее среднее в виде

00 Если сравнить веса уравнения скользящего среднего Шеп­

парда с ортогональными полиномиальными членами, приведен­ ными в табл. 4,30, то мы увидим, что это просто квадратное уравнение, натянутое методом наименьших квадратов на пять точек. Другие сглаживающие уравнения, используемые во вре­ менном тренд-анализе, также либо идентичны, либо очень напо­ минают ортогональные многочлены. Эти сглаживающие функции описаны в классической работе Уиттекером и Робинсоном [60J. На рис. 4.60 представлены результаты сглаживания временного ряда каротажа скважин с помощью различных уравнений. Со­ всем, не очевидно, что эти функции скользящего среднего явля­ ются фильтрами, или что процесс скользящего среднего мате-

Т а б л и ц а 4.30

Ортогональные полиномиальные весовые функции для аппроксимации (сглаживания). Функция прямо дает оценку среднего значения

в последовательности из п точек, где п— нечетное число от 5 до 17 [48]

п

3 W

Квадратичные			и кубичные		члены			Члены четвертого и П Я ТО ГО порядков
5	7	9	11	13		15	17	л	7	9	11	13	15	17
- 3	— 2 —21		—36	- 1 1	—78 —21				5	15	18	п о	2145	195
12	3	14	9	0 —		13	—6		— 30 —55 —45 - 1 9 8 - -28S0 — 195
17	6	39	44	9		42	7		75	30 —	10 — 160--2937 —260
12	7	54	69	16		87	18		131	135	60	110	165	117
— 3	6	59	84	21		122	27		75	179	120	390	3755	135
	3	54	89	24		147	34		— 30	135	143	600	7500	415
	—2	39	84	25		162	39		5	30	120	677	10125	660
		14	69	24		167	42			— 55	60	600	11053	825
		— 21	44	21		162	43			15— 10		390	10125	883
			9	16		147	42			— 45		ПО	7500	825
			—36	9		122	39				18 —	160	3755	660
				0		87	34				—	198	— 185	415
				— 11		42	27					ПО--2937		135
					—	13	18						-2860 —	117
					— 78		7						2145 — 260
							—6						—	195
							— 21							195
Ю СО	21	231	429	143	1105		323	2<о	231	429	429	2431	46189	4199

20— 201

305

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

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

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

...0 0 0 2 1 0 0 ...,

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

вскользящем среднем, то результат будет такой:

. . . . 0	0	0		2		1	0	о . . . .
							-0,5712	1,4286
					-0	СЛ*2	1,4286	,------- '
								I
			.-0,5712		1,428o~~|
	4JJi712	1>ans	7		I
. J		1	i		i	j
		0	2,6572	0,2862		-0,5712		0 . . . .

Это в точности тот же результат, который мы получили приме­ нением операции свертки.

Производные

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

нения значений Y. Если Y увеличивается или уменьшается внут­ ри интервала с примерно постоянной скоростью, то производная

306

Рис. 4.60. Результаты сглаживания последовательности значений каротажной диаграммы в скважине е помощью различных уравнений [26, рис. 4.21:

а — исходные данные; б—д — число членов; 5 (б); 9 (S), 13 (г), 17 (с?)

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

Если подгоняемая кривая определена ортогональным много­ членом, то производная может быть выражена через ортогональ­ ные члены. Так как нас интересует только вычисление производ­ ной подгоняемой линии в центральной точке интервала, то поли­ номиальное уравнение может быть переписано так, чтобы пря­ мо дать значение производной. Этот процесс описан Савицким и Голеем [48J, которые дают также таблицы для нахождения пер­ вой п высших производных ортогональных многочленов. Табл. 4.31 содержит члены для нахождения первой производной, взя­ тые из работы [48J.

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

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

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

2 0 *

307

Т а б л и ц а 4.31

Ортогональные полиномиальные функции для первой производной. Функция прямо дает оценку первой производной в центральной точке полиномиальной кривой, подгоняемой в п точках методом наименьших квадратов, где п — нечетное число от 5 до 17

		Производная квадратичного многочлена
п	5	7	9	п	13	15	17
	- 2	—3	— 4	— 5	—6	— 7	—8
	— 1	— 2	—3	— 4	— 5	— 6	— 7
	0	— 1	—2	— 3	— 4	— 5	— 6
	1	0	— 1	— 2	—3	— 4	—5
	2	1	0	— 1	—2	—3	— 4
		2	1	0	— 1	— 2	—3
		3	2	1	0	— 1	—2
			3	2	1	0	— 1
			4	3	2	1	0
				4	3	2	1
				5	4	3	2
					5	4	3
					6	5	4
						6	5
						7	6
							7
							8
So)2	10	28	60	ПО	182	280	408
Производная кубичного многочлена и многочлена четвертой степени [481
п	5	7	9	1!	13	15	1Т
	1	22	86	300	1133	12922	748
	—8	— 67	— 142	— 294	—660	—4121	—98
	0	— 58	— 193	—532	— 1578	— 14150	—643
	8	0	— 126	—503	— 1796	— 18332	—930
	— 1	58	0	—296	— 1489	— 17842	— 1002
		67	126	0	—832	— 13843	—902
		—22	193	296	0	— 7506	—673
			142	503	832	0	-—Зь8
			- 8 6	532	1489	7506	0
				294	1796	13843	358
				—300	1578	17842	673
					660	18332	902
					— 1133	14150	1002
						4121	930
						— 12922	643
							98
							—748
Норма-	12	252	1188	5148	24024	334152	23256
лизую-
щий член

б

3 г

2 -

1 -

0L

Рис. 4.61. Ступенчатая функция, состоящая из последовательностей единиц и троек (а). Графическое изображение первой производной ступенчатой функ­ ции, вычисленной на основе пятичленной квадратичной полиномиальной ап­ проксимации (б)

ск

та

I

<

о

ю

п

X

о

са

1	5	10	15	20	2 5	30	33
		Номер	с к в а ж и н ы

Рис, 4.62. Измерения пористости, сделанные с интервалом 10 футов (3 м) в скважине, пробуренной сквозь нефтяной резервуар на арктическом побережье Аляски:

а — необработанные данные; б — первая производная необработанных данных

ной скважине. На рис. 4.62 также представлены их первые про­ изводные или, скорее, первые производные квадратичных функ­ ций, подогнанных к последовательным множествам из пяти то­ чек. Измерения пористости обладают большой изменчивостью, частично из-за изменчивости экспериментальных данных, но так­ же в силу того, что использование очень маленьких объемов проб приводит к непредставительной выборке из резервуара. Хо-

309

										Т а б л и ц а 4.32
Значения	пористости (в		%), измеренные в продуктивном слое нефтяной
скважины, пробуренной на арктическом побережье Аляски.
		Пробы взяты с интервалом в 3 м
(Верх)	1	13,2	9	6,7	17	18,9	25	21	2
	2	12,6	10	12,8	18	22,5	26	21,7
	3	14,3	И	14,2	19 21,6		27	23,0
	4	16,2	12	16,3	20	23,2	28	22,1
	5	15,2	13	15,4	21	22,8	од	22,4
	6	14,1	14	14,0	22	24,7	30	21,8
	7	11,2	15	15,0	23	20,0	31	19,8
	8	15,3	16	23,2	24	23,4	32	24,7		(Основание)

тя график первой производной ясно показывает максимальное изменение пористости в образце с номером 16, ошибки в исход­ ных данных также оказывают влияние на остатки.

На рис. 4.63 представлены измерения пористости после сгла­ живания пятнчленным квадратичным фильтром. Хотя различия в пористости между верхней и нижней частями резервуара со­ храняются, большие флуктуации от точки к точке погашаются. Если необработанные данные действительно составлены как смесь «истинных» пористостей со случайными ошибками, то гладкая аппроксимация может давать более близкую к дейст­ вительности картину изменения пористости, чем исходные дан­ ные.

Первая производная сглаженных данных также представле­ на на рис. 4.63. Явное изменение пористости между пробами с номерами 1 и 16 ясно показано. Как и следовало ожидать, все производные остатков сглаженной кривой близки к нулю.

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

Рис. 4.63. Сглаженные значения пористости из арктической скважины:

а-^данные, сглаженные с помощью пятичленного квадратичного многочлена; 6 — первая производная сглаженных данных

310

фильтрации с дисперсией после фильтрации. Необходимые сум­ мы квадратов для вычисления дисперсии таковы:

для исходного временного ряда

S*o =	----	П
=	-----L (s r r ,
		п*

для ряда, прошедшего через фильтр

SS, = 2 К2----—

п1

для отклонений

s s a. = - L z ( K - n a.

п*

Приближенное процентное отношение сумм квадратов имеет вид (SSj/SSc*) 100%. В этих равенствах все суммирования про­ водятся по п данным точкам исходной последовательности или по точкам сглаженной последовательности. В частности, отме­ тим, что вычисляются два значения для суммы квадратов пер­ воначальных данных. Первая из них учитывает все точки дан­ ной последовательности от i= 1 до п, вторая включает только те наблюдения, для которых вычисляются оценки Yi. Таким обра­ зом, из-за потери данных на концах сглаженной последователь­ ности мы имеем л* = п— (т-—1). Хотя в этом случае, как и в рег­ рессионном анализе, можно было бы ожидать выполнения ра­ венства SSo* = SS/-j-SS(/, тем не менее его нет в методе сколь­ зящего среднего. Это происходит по той причине, что вычисле­

ние оценок Yt вблизи концов последовательности данных отчас­ ти основано на использовании значений, не входящих в вычис­ ление SSo*. Поэтому значение процентного отношения сумм квадратов является приближенным, по его можно использовать как показатель эффективности процесса фильтрацнн.

АНАЛИЗВЗАИМОЗАМЕНЯЕМОСТИ

На рис. 4.64 представлен гипотетический стратиграфический разрез с закодированной последовательностью различных лито­ логических разновидностей пород. Если вы исследуете закоди­ рованную последовательность, то заметите, что в ней часто встречаются последовательности А-*-В-*С н A-+D-+C. Это на­ водит на мысль о том, что состояния В и D как-то связаны и одно может быть заменено другим в предлагаемой последова­ тельности. Это свойство двух или более состояний встречаться в одном и том же окружении называется взаимозаменяемостью,

311