Метод наименьшей кривизны

⇐ ПредыдущаяСтр 9 из 10Следующая ⇒

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

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

В качестве сплайн-модели используют отрезки из кубических полиномов, которые плавно и гладко стыкуются в пунктах измерений. Каждый кубический полином имеет вид j(x) = a ₀ + a ₁ x + a ₂ x ² + a ₃ x ³и содержит четыре коэффициента, которые нужно найти. Если имеется n пунктов измерений вдоль линии, то между ними должно располагаться n -1 кубических полиномов. Следовательно, всего нужно найти 4(n -1) неизвестных коэффициентов путем решения системы линейных уравнений. При большом количестве точек система уравнений получается громоздкой. Для сокращения объема вычислений применяется так называемый скользящий сплайн, когда для составления системы уравнений и нахождения коэффициентов сплайна используют четыре соседних пункта измерений. Между ними получается три отрезка, и система уравнений содержит всего 12 неизвестных коэффициентов. Вычисленные коэффициенты полинома применяются к среднему отрезку. Далее вычисления перемещаются (скользят) на одно измерение, и расчет производится снова. Перемещение производится до тех пор, пока не кончатся пункты измерений.

Для составления системы уравнений сплайна используются следующие соотношения:

1. Каждый полином должен проходить через левую точку отрезка, что дает четыре уравнения.

2. Каждый полином должен проходить через правую точку отрезка, что дает еще четыре уравнения.

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

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

5. На конечных пунктах измерений (а их два) необходимо задать граничные условия. Чаще всего задается условие, что кривизна в крайних точках (вторая производная) должна быть нулевая. Иногда задается наклон в крайних точках (т.е. первая производная). Граничные условия дают два уравнения.

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

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

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

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

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

Точечный кригинг

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

Таблица 9

Значения свойств j

в пунктах наблюдений

N	j	N	j	N	j	N	j
	7.14		13.84		11.63		10.44
	13.51		13.38		12.94		1.99
	8.80		13.97		11.92		8.18
	1.59		6.53		12.20		6.87
	6.66		6.47		14.48		11.25

Требуется определить прогнозное значение свойства j в пункте Z. Для расчетов используем данные рис.12, где практически одинаковые дисперсии свойств D и радиус автокорреляции R =29 м. Там же приведена формула теоретической вариограммы:

g(h) = 1,183 h – 0,03285 h ²+ 0,000286 h ³(34)

Рис.20. Пример точечного кригинга

Рис.20. Схема расположения пунктов наблюдений

Таблица 10

Расстояния между пунктами в метрах

N							Z






Z

Из рисунка 20 видно, что в пределах радиуса автокорреляции, проведенного вокруг пункта Z, расположены 6 пунктов с номерами 7, 12, 13, 15, 18 и 19. Все они будут участвовать в расчете прогнозного значения свойства в пункте Z. Вначале необходимо определить все расстояния между всеми пунктами, включая пункт Z (табл.10).

Следующий пункт заключается в вычислении вариограммы между всеми пунктами наблюдений по формуле (2), используя расстояния, приведенные в таблице 10. Напомним, что за пределами радиуса автокорреляции R = 29 м, вариограмма равна дисперсии D = 13.655.

Составим систему уравнений кригинга, используя формулу 18:

13.655 1.655 0.483 0.000 0.000 0.000 0.142

1.655 13.655 0.650 1.965 0.000 0.032 6.886

0.483 0.650 13.655 0.000 0.483 0.000 3.199 (35)

0.000 1.965 0.000 13.655 0.000 1.075 0.650

0.000 0.000 0.483 0.000 13.655 0.231 0.483

0.000 0.032 0.000 1.075 0.231 13.655 0.483

Решение системы дает следующие коэффициенты:

-0.583 -.0.505 0.211 -0.028 0.028 0.036 сумма -0.841

Разделив коэффициенты на их сумму, получим весовые коэффициенты кригинга:

0.693 0.600 -0.251 0.033 -0.033 -0.043

Далее по формуле (17) вычисляем прогнозное значение свойства в пункте Z:

j_z = 0.693*13.38 + 0.600*11.92 – 0.251*12.94 + 0.033*14.48 – 0.033*8.18 – 0.043*6.87= 13.09

Таким образом, на результат прогноза влияют все значения в пределах радиуса влияния вокруг пункта Z. Наиболее сильно влияют близлежащие пункты 13 и 12, у них самые большие весовые коэффициенты.

Кстати, применение метода обратных расстояний с использованием всех пунктов в пределах радиуса автокорреляции дает прогнозное значение в пункте Z = 11.0.

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

Площадной кригинг

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

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

Площадной кригинг обычно осуществляется путем скольжения (сканирования) по всем пунктам внутри области измерений. Для каждой прогнозной точки выбирают все пункты измерений не далее радиуса R. Далее находят все расстояния между точками, как указано в пункте (3). Для каждого расстояния вычисляют выражение D - g(h). Из этих выражение составляют систему уравнений кригинга (18), приведенную выше. Решение системы уравнений дает коэффициенты кригинга, сумма которых приводится к единице. После этого рассчитывают прогнозное значение вначале в первой точке, а потом, скользя по графику (рис.20), и в остальных точках. В тех точках, где прогнозный пункт совпадает с измеренным пунктом, прогнозное значение принимается равным измеренному значению.

Рассмотрим двухмерный кригинг на примере главного рудного тела апатитового месторождения Коашва. В качестве исследования взяты координаты кровли рудного тела. Положение пунктов пересечения скважин с кровлей рудного тела в плане показано на рис. 21. Рудное тело падает на северо-запад к центру Хибинского массива под углом около 45^о.

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

соответствует 5м на местности.

Рис.21. План расположения пунктов пересечения скважин с кровлей рудного тела

Прежде чем рассчитывать вариограмму, из координат кровли рудного тела вычтены значения тренда, получен остаток от тренда. По остаткам вычислена эмпирическая вариограмма в предположении изотропной среды по простиранию и падению рудного тела (рис.23). Эмпирическая вариограмма аппроксимирована трендом в виде кубической параболы. На вариограмме показана дисперсия D = 657,0 и радиус влияния каждой скважины R = 450 м. Имеется незначительный эффект самородков, равный 55 м.

Имея вариограмму и радиус влияния, можно выполнить интерполяцию значений свойств (координат кровли рудного тела) между скважинами. Порядок вычисления следующий. Вокруг каждой из точек в пределах контура рудного тела (см. рис.22) проводится радиус влияния R. В эту окружность в данном примере может располагаться различное количество скважин – от нуля до 15. В зависимости от количества скважин в окружности принимаются различные решения.

Рис.23. График эмпирической (ломаная линия) и теоретической (плавная линия) вариограмм остатков от тренда кровли рудного тела

Рис.22. Тренд поверхности рудного тела, падающего на северо-запад. Глубина тренда дана в абсолютных отметках координаты z

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

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

Если в окружности находится одна скважина, то можно рекомендовать найти расстояние от прогнозной точки до скважины, вычислить “вес” этой скважины по вариограмме и определить прогнозное значение как сумму тренда и значения свойства j в скважине с весом p:

(36)

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

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

Подобные задачи можно решать также с помощью пакета Micromine и других подобных пакетов.

Рис.24. Пример разделения площади рис.22 на равные квадраты

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

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

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

Используя разделение пощади на квадраты (рис.24), учитывая вариограмму (рис.23), тренд и кригинг остатков, для центра каждого квадрата рассчитано прогнозное значение глубины кровли рудного тела (рис.25). Квадраты можно брать различного размера, здесь приняты квадраты 20*20 пикселей. На рис.25 также видно погружеиие кровли рудного тела на северо-запад, но имеются осложнения, вызванные кригингом. Отклонения результатов кригинга от тренда достигает иногда 100 м, тем самым учитывается неровность рельефа кровли рудного тела.

Рис.25. Изображение кровли рудного тела с применением тренда и кригинга.

Справа шкала глубин в абсолютных отметках (метры)

Дата добавления: 2016-01-06; просмотров: 41; Мы поможем в написании вашей работы!

Поделиться с друзьями:

⇐ Предыдущая 1 2 3 4 5 6 7 8910 Следующая ⇒

Мы поможем в написании ваших работ!