CC BY
21
ГОРНО-ПРОМЫШЛЕННАЯ ГЕОЛОГИЯ
© В.А.: Ермолов, В.П.: Зервандова,
А.Н. Быховец, 2000
УДК 622:014.3:502.76
В.А. Ермолов, В.П. Зервандова, А.Н. Быховец
МАТЕМАТИЧЕСКОЕ ОБЕСПЕЧЕНИЕ ЭКОЛОГОТЕХНОЛОГИЧЕСКОГО РАЙОНИРОВАНИЯ ТЕХНОГЕННЫХ МЕСТОРОЖДЕНИЙ
Введение
Разработка и внедрение на горных предприятиях технологических приемов безопасного размещения и получения из твердых минеральных отходов дополнительных товарных продуктов, обладающих высокой потребительской стоимостью, в настоящее время является актуальной научной и производственной проблемой. Известно, что суммарная площадь, занимаемая отходами горнодобывающих предприятий в России, превышает 500 тыс. га, а территория, на которой явно ощущается их негативное влияние на окружающую человека среду, - в 10 - 15 раз больше [1, 2, 5].
Ежегодный выход хвостов в процессе обогащения руд составляет примерно 200 млн. м3, представленных, в основном, тонкоизмельченными минералами с относительно высоким содержанием в них полезных компонентов.
Оценка показывает, что их использование позволило бы не менее, чем на 15-20 % расширить минерально-сырьевую базу горно-металлурги-ческой, горно-химической и строительной отраслей промышленности. Однако, до настоящего времени уровень использования минерального вещества техногенных образований остается очень низким. Например, при потенциальной возможности прямого использования в качестве минерально-строительного сырья около 30 % извлекаемых из недр вскрышных, вмещающих пород и отходов обогащения сейчас утилизируется не более 4 %.
Проблема направленной утилизации отходов рудообо-гащения комплексная. Приоритетная роль в направленной утилизации принадлежит геологическому обеспечению. Из широкого спектра геологических задач, рассмотренных нами в работах [4-7], в настоящей статье освещаются вопросы математического обеспечения районирования экологотехнологических зон техногенных месторождений.
1. Методы и алгоритмы выделения экологотехнологических зон техногенных месторождений
Районирование техногенных массивов связано с выделением однородных эколого-технологи-ческих зон. Выделение однородных зон включает задачу классификации (распознавания) объекта и состоит из следующих основных этапов: нахождение оптимальной системы признаков разделения; обоснование решающего правила в системе признаков; выбор статистических критериев для проверки гипотезы об однородности.
Обоснование оптимальных признаков целесообразно осуществить с помощью методов факторного анализа, при котором исследуется матрица ковариаций или корреляций [3, 9]. Алгоритм решения задачи нахождения оптимальной системы признаков (сокращения пространства) с помощью факторного анализа предусматривает нормирование (стандартизацию) данных, вычисление элементов корреляционной матрицы, а также средних и дисперсий по признакам. Далее производится факторное преобразование и вычисляются: процент вкладов факторов в суммарную общность; коэффициенты при ортогональных факторах; корреляции признаков с косоугольными факторами; коэффициенты при косоугольных факторах; матрицы корреляций между косоугольными факторами; матрица коэффициентов при параметрах; значения факторов. После нахождения косоугольных факторов и их матрицы корреляции производится повторное факторное преобразование, причем косоугольные факторы принимаются за начальные данные. Так продолжается до тех пор пока не останется лишь один фактор ^го порядка.
Второй этап распознавания состоит в построении решающего правила в найденной оптимальной системе признаков, т.е. в проведении кластерного анализа (в данном случае речь идет о методах классификации без обучения). Для этого целесообразно использовать иерархическую кластеризацию, в которой множество мер сходства между всеми парами объектов представляется в виде симметричной матрицы [3], либо путем последовательного разбиения на классы с оценкой степени их однородности [10].
Иерархическая кластеризация на основе матриц. Для вычисления элементов симметричной матрицы порядка п X п следует провести транспонирование матрицы сходства исходных данных, порядок которой п X т, в матрицу т X п. В результате будет получена матрица сходства между переменными порядка т X т (в отличие от корреляционной матрицы сходства между наблюдениями п X п ). Элемент матрицы Су дает характеристику сходства между г-ым и у-ым объектами
Таблица 1
СТАТИСТИЧЕСКИЕ ХАРАКТЕРИСТИКИ ПОКАЗАТЕЛЕЙ ПО СЛОЮ +279 М И ГЕОЛОГОТЕХНОЛОГИЧЕСКИМ ЗОНАМ IИIIКОВДОРСКОГО ТЕХНОГЕННОГО МЕСТОРОЖДЕНИЯ
Показатели Среднее значение С, % Стандартное отклонение С, % Коэффициент вариации V, % Асиммет- рия А Эксцесс Е
Слой +279 м
Feобщ 3,49 0,47 13,3 1,29 2,31
Р2О5 11,1 1,58 14,2 0,68 3,31
СО2 12,0 2,19 18,2 -0,68 3,66
MgO 20,4 1,72 8,4 0,32 0,58
ZrO2 0,35 0,17 51,5 2,01 4,67
FR +3,0-0,14мм 40,0 16,6 47,5 0,06 -1,01
FR +0,074мм 30,9 16,5 53,4 1,43 2,51
FR -0,074мм 18,4 18,8 102,6 1,08 0,50
Р, мкр/ч 23,0 9,19 40,0 1,81 1,04
Слой +279 м - зона I
Feобщ 3,53 0,44 12,5 1,26 2,33
Р2О5 11,2 1,66 14,8 1,97 6,11
СО2 11,3 2,19 19,4 -0,94 7,73
MgO 20,8 1,57 7,52 0,06 1,32
ZrO2 0,34 0,17 50,7 1,84 3,75
FR +3,0-0,14мм 47,9 11,0 23,0 -0,05 -1,0
FR +0,074мм 23,6 8,13 34,5 0,17 -0,8
FR -0,074мм 15,5 7,06 45,6 0,33 -0,18
Р, мкр/ч 31,6 9,39 29,7 1,18 -1,3
Слой +279 м - зона II
Feобщ 3,44 0,32 9,31 0,48 -0,59
P2O5 10,0 1,51 15,1 -0,08 -0,91
CO2 12,0 1,81 14,8 -1,05 0,99
MgO 20,9 1,88 8,95 0,71 -0,31
ZrO2 0,24 0,05 20,2 0,98 0,03
FR +3,0-0,14мм 18,4 8,03 43,5 0,01 -1,35
FR +0,074мм 34,5 11,1 32,1 -0,27 -0,67
FR -0,074мм 44,3 14,7 33,4 0,49 1,09
Р, мкр/ч 17,9 3,32 18,6 1,39 -0,76
(скважинами или пунктами опробования с координатами х, .К 2).
Следующий этап - получение иерархических групп объектов, при которых объекты с наивысшими коэффициентами сходства размещаются вместе. Затем группы объектов соединяются в другие группы, с которыми они наиболее тесно связаны, и так продолжается до тех пор, пока не будет получена полная классификация в виде дендрограмм. Количество выделенных групп устанавливается на основании значимого коэффициента корреляции при соответствующем уровне разделения или технологическими схемами переработки отходов рудообогащения. Программная реализация данного алгоритма выполнена инженером-программистом кафедры геологии МГГУ Г.А. Воротниковой.
В соответствии с изложенным подходом для складированных отходов рудообогащения Ковдорского апатит-магнетитового месторождения на основании метода главных компонент определены нагрузки на собственные вектора (главные компоненты). Первые два собственных вектора представляют 94 % дисперсии изучаемого множества данных ^еобщ., Р205, ZrO2, MgO,CO2, выход фракций: +3,0 -
0.14 мм; 0,14 - 0,074 мм; -0,074 мм). Первый вектор, представляющий 58,1 % общей дисперсии, дает наибольшие вклады в переменные FR - 0,074 мм, а второй - представляющий 36,2 % - в переменные Р205, и ZrO2. Таким обра-
зом, первая главная компонента характеризует гранулярный состав техногенного месторождения, а вторая -его качественное состояние. Следует также отметить, что на третий собственный вектор приходится 4,1 % общей дисперсии, причем наибольший вклад составляет переменная FR +0,074 мм. Однако, по сравнению с первыми факторами, нагрузка данной компоненты несущественна. На рис. 1 дано графическое представление нагрузки на главные компоненты, которые показывают, что районирование техногенного массива может быть проведено на основе показателей: Р205, ZrO2, FR + 3,0-0,14 мм и FR - 0,074 мм.
Для выделения собственно эколого-технологи-ческих зон использован алгоритм группировки по коэффициентам корреляции с построением дендрограммы по значимым классификационным признакам. При критическом уровне сходства 0,4 (матрица 120 X 4) выделены две эколого-технологические зоны, статистические характеристики которых представлены в табл. 1.
Статистический анализ показателей зон показывает, что их выделение произошло в основном по параметрам гранулярного состава. Так, коэффициента вариации FR - 0,074 мм в зоне I составляет 45,6 %, в зоне II - 33,4 %, а в целом по слою без выделения зон 102,6 %. Аналогичная картина отмечается и для FR + 3,0-0,14 мм и FR + 0,074 мм. При этом в зоне I преобладает выход FR + 3,0-0,14 мм - 47,9 %, а зоне II - FR - 0,074 мм - 44,3 %.
«Контрастно-групповой анализ» статистических данных Характерной особенностью данного метода является последовательное разбиение исходной совокупности на все новые и новые подклассы. Сначала вся совокупность делится на два класса по значениям переменных (значимых факторов), каждый из которых далее разбивается на два подкласса и т.д.
Два подкласса, полученные при очередном разбиении, должны отвечать следующим условиям:
• внутренняя однородность каждого подкласса выше, чем однородность их объединения;
• взаимная однородность подклассов максимальная. Разбиения на подгруппы осуществляются на основании
оценок дисперсии (суммы квадратов),
Рис. 1. Нагрузка на главные компоненты:
А - 1-ая главная компонента; Б - 11-ая главная компонента; 1-8 - показатели компонент: 1 - Fe06щ; 2 - Р205$ 3 - ZrO2; 4 - С02; 5 - MgO; 6 -FR+0.14 мм; 7 - FR+0.074 мм; 8 - FR-0.074 мм
поскольку дисперсия является характеристикой однородности или неоднородности выборочной совокупности. Рассмотрим алгоритм «контрастно-группововго анализа».
Пусть
о 1 п _ ^
°2 =1Е (уг — У)2 (1Л)
п 1=1
- 1 п
и У = п Е Уг (12)
п г=1
2 1 п _
Величина п<5 = — Е ( У{ — У) (13)
п г=1
называется общей суммой квадратов (TSS).
При разбиении всей совокупности на два класса имеем,
а) в каждом классе будет ц элементов (у=1,2), причем 2
Е пу = п (1.4)
у=1 у
б) индекс г пробегает в каждом классе значения от 1 до
п .
Данные имеют вид Уу(г=1,...,щ- ]=1,2).
Пусть тогда
- 1 пу У.у = + Е У.
(1.5)
Т$>$> можно теперь преобразовать следующим образом:
2 2 пу - 2 2 пу __________2
= Д^/У/ — У)2 =Е1Е1((Угу—Угу)+(Уу—У))2 =
2 п] _2_________________2 __________________
=Д .Е ((Угу—УуГ + (Угу — У)2 +2( Угу — У у)( У у—У) =
2 пу _ 2 2__________2 2___________пу _
Е Е (Уц — Уу)2 +2 Е. п,(Уу-У)2 +2 Е (Уу — У> Е (Угг— Уу)
у=1г=1
Так как
пу\
Е (У г=1
у=1
п
у=1
г=1
П] -у
(1.6)
и
у = Д Уу — пуУу
■ 0.
2 2 пу _ 2
пст = Д (Е1(Уу—у +
2___________2
Е пу(Уу — У) у=1
Сумма
2 П у Ел Е(Уу J =1 г=1
У ■■)2
у
(1.7)
(18)
называется «внешней суммой квадратов» (№88). Это сумма на подклассах, она характеризует их однородность. Сумма
2 „ - 2 „ -2 . „ -2 „-2
у
1
пу(уч - У)2
■пуУ г1 + п2 У г2
пУ
(1.9)
называется «внешней суммой квадратов» (АХ?) и определяется средними на подклассах и средним по числу элементов в подклассах. Поэтому показатель BSS характеризует взаимную неоднородность подклассов.
Описание модели. Пусть даны: т независимых переменных Хг (г=1,...,т), причем каждая переменная Хг (г=1,...,т), причем каждая переменная Хг имеет ^ градаций Су у\=1,...,Ь); зависимая от этихХг переменная Y и выборка объема п. Пусть каждый из этих п элементов принадлежит в точности одному из kl, ^,..., ^ классов. Таким образом, элементы массива данных являются точками т-мерного пространства, при этом имеется взаимно однозначное соответствие между классами и наборами (1, у2,..., ут).
В каждом классе имеется пц j2,..,т элементов. Следовательно,
У п ,• ,• = п
(1,)
(1.10)
Вместо множества значений Т\уцр,...,т (^=1,...,пцр..., т) в некотором классе рассматривается их среднее арифметическое:
_ 1 Пі1...м
у =------------------------- У лу > ?■
Л ■’т І1...1 т Л=1 1
(1.11)
Среднее арифметическое по всему множеству равно
- 1 ^ -
У = — У п ,• у • У і і пґі ) п ) 11...]т ]1-]т
(.71 /-Л ]т'
(1.12)
Выполнив процедуру подготовки, осуществим теперь разбиение. Это означает разбиение множества всех значений независимой переменной х, где і выбирается произвольно, но затем фиксируется. На следующем шаге вычисляются значения BSS, WSS и ^ .
Для произвольно выбранного, но фиксированного і в качестве которого в последующих уравнениях фигурирует I, что не уменьшает общности, получим
(іЛі =^и ■*и ^Сі'ІІ С1,І } П {С1,ІІ } = 0 , (113)
где I С {і | ]1=1,...,к} = ЮІІ Далее найдем
У У п / ,• = пС т (1.14)
(і1 еіИіі 1) ]1і 11
и У У п, / = пСл „ (1.15)
(11ЄІІ)(1І Ф11) 11
Сг у £ Сг т или сщ соответственно на множестве Хг.
^ г ’
Это можно сделать лишь выполнив перебор. Однако следующее соображение облегчает задачу. Если упорядочить ^ классов Сг у независимой переменной XI по возрастанию
'Н
значений
У11-.-1г.
то легко показать, что оптимальным
Ч"' -^т
будет одно из следующих разбиений I = {1} , II = {2,...,кц}
I = {1, 2} , II = {3,...,к,}
I = {1,...,к-1}, II = {к}
Если оптимальных разбиений будет р, то величина
1 Р
Ж -1
будет характеризовать общую долю объ-
Отсюда
- _ 1 ^ -УС1.І = ”си(/е І)(1і У11 п^т-^Г1
_ (116)
У п . вычисляется аналогично. Отсюда имеем С1,І
Ш = Пс^ІІ(Ус1,ІІ -у)2 + пс1,І(Ус1,І -у)2
(1.17)
Назовем величину (BSS/TSS)^100 «процентом объясненного TSS».
Как уже отмечалось, BSS должно принимать максимальное значение после каждого расщепления. Таким образом, дело сводится к поиску оптимальной комбинации
ясненного ТХХ, причем р пробегает значения от 1 до
т
Е к — 1. Но поскольку целью данного метода является
г=1
получение определенного набора объясняющих значений, зададим некоторые правила разбиения и установки:
а) после разбиения Б(Б=1,...,р) следует анализировать подгруппу с максимальным значением ТЯЯ?, поскольку это необходимое условие максимизации значения ВХХХ; (ТЯЯ? /ТЯЯ)Х100 должно быть больше 1-2 %, в противном случае эту подгруппу практически можно считать однородной; число пя должно быть не меньше 10-20, в противном случае стандартная ошибка становится слишком большой;
б) (ВББХ /ТББ)Х100 должно быть 1-2 %, в противном случае разбиение не приведет к получению существенной объясняющей информации; значение р не должно быть слишком большим, так как это затрудняет интерпретацию.
Технологическими исследованиями [4] установлено, что при существующей технологии получения бадделеитового концентрата, принятой на Ковдорском ГОКе, потери при гравитационном обогащении ввиду переизмельчения зерен циркона во фракции -0,074 мм составляют 50 % и более в связи с чем получение данного концентрата из хвостов с преобладающим выходом FR -0,074 мм экономически нецелесообразно. Получение бадделеитового концентрата при существующей технологической схеме возможно лишь из лежалых хвостов зоны I, в которой суммарный выход фракций +3,0-0,014 мм и +0,074 мм в среднем составляет 70 %. Получение апатитового и фосфорно-магнезиального концентратов возможно как из зоны I, так и из зоны II. При этом лежалые хвосты зоны II могут быть направлены непосредственно на флотацию без дополнительного измельчения, в отличие от хвостов зоны I, требующих доизмельче-ния.
Мвбшовскии государственный горный
УнийГейсйтетКИБ ХАРАКТЕРИСТИКИ ПОКАЗАТЕЛЕЙ ПО ГЕОЛОГО-ТЕХНОЛОГИЧЕСКИМ У нпзОжМ;Ьь1деленным методом «КОНТРАСТНО-ГРУППОВОГО анализа»
Показатели Среднее значение С, % Стандартное отклонение С , % Коэффициент вариации v, % Асим- метрия А Эксцесс Е
Слой +279 м - зона I
Feoбщ З,61 0,З6 10,1 1,15 2,З7
P2O5 11,4 1,46 12,8 0,ЗЗ 1,15
CO2 11,З 1,86 16,6 -0,58 0,64
MgO 20,5 1,50 7,З 0,28 0,49
ZrO2 0,З7 0,1З З5,6 1,96 2,З1
FR +З,0-0,14мм 54,2 1З,1 24,1 0,15 -0,86
FR +0,074мм 17,1 5,З5 З1,З 1,12 1,96
FR -0,074мм 10,6 З,9З З7,1 0,88 0,56
Р, мкр/ч З4,4 8,77 25,5 1,З2 1,01
Слой +279 м - зона II
Feoбщ З,48 0,З4 9,8 1,1З 0,ЗЗ
P2O5 10,8 1,1 10,2 1,З5 1,52
CO2 10,9 1,77 16,2 1,З4 0,57
MgO 19,5 1,54 7,9 -0,58 0,64
ZrO2 0,З1 0,11 З4,4 -0,47 0,82
FR +З,0-0,14мм 44,2 10,19 2З,1 -0,66 1,02
FR +0,074мм 28,5 8,69 З0,5 0,15 1,26
FR -0,074мм 16,7 5,З6 З2,1 0,6З -0,51
Р, мкр/ч 28,6 6,98 24,4 0,28 0,86
Слой +279 м - зона III
Feoбщ З,З4 0,28 8,З6 0,85 0,67
P2O5 10,2 1,19 11,7 0,46 -0,80
CO2 11,8 1,6З 1З,8 -0,27 -0,86
MgO 20,2 1,60 7,9 0,59 0,6З
ZrO2 0,22 0,04 18,6 0,04 -1,25
FR +З,0-0,14мм 12,1 4,З1 З5,6 0,15 1,11
FR +0,074мм З7,5 10,9 29,1 0,ЗЗ 0,67
FR -0,074мм 50,2 15,1 З0,1 1,08 2,0З
Р, мкр/ч 15,З 2,26 14,8 0,61 0,6З
Выбор математических методов моделирования поверхностей (линий), ограничивающих эколого-технологи-ческие зоны техногенных массивов, определяется априорным характером зависимости между координатами характерных точек поверхности (линий), полученных в результате выделения однородных зон по данным опробования [8].
Рассмотрим случаи, возникающие при моделировании ограничивающих плоскостей. Для остроугольного треугольника (рис. 2) ближайшими районами каждой его вершины являются выпуклые четырехугольники S1, S2 и Sз с общей вершиной Q, лежащей внутри треугольника. Найдем значения площадей S1, S2, Sз и координаты точки Q.
Определим сначала длины сторон треугольника
а =М1М2 = ^(х1 -х2)2 + (у2 -у1)2.
ь=м 3 =
с = 2М 3 :
ХЗ— xi)2 + (уЗ— у1)2,
На основе контрастно-группового анализа показателей качества лежалых хвостов Ковдорского техногенного месторождения выделены три геолого-технологических зоны (табл. 2). Таким образом, зоны I и II, выделены методом «иерархической кластеризации», а зоны I и Ш - методом «контрастно-группового анализа», что,
R =
= д/(х3 — х2 )2 + (У3 — У2 )2 (2.1)
Радиус описанной окружности, центр которой находится в точке Q, найдем по формуле С
2 sin а
2л/ї
2
— cos^а
(2.2)
где а - угол между векторами а и b .
Используя координатную форму скалярного произведения, имеем
Ш (x2 — xi)(хЗ — xi) ,
cos а =---- —--------------і------—
lb
ab
(2.3)
где а и Ь определены формулами ( 2.1).
Найдем теперь расстояния от точки Q до сторон треугольника, как это показано на рис. 2:
Рис. 2. Схема расположения «ближайших районов» для треугольной блокировки
' ' 2 Г2
hi =
V4R 2 — а2
h2 =
^4R 2 — c 2
hз =
V4R2 —
(2.4)
как показывают исследования, позволяет дать более детальное описание строения техногенного массива, и, соответственно, обосновать наиболее эффективные технологии переработки отходов рудообогащения.
2. Алгоритм оконтуривания эколого-технологических зон при районировании техногенных массивов
Окончательно имеем:
а
2
2
2
ЬЫ
+ -
_ah1
Б1 —г
Б _ ahl , ь^ Б2 _ — + —
_ 0*2 + Ь1з
(2.5)
4 4
Для определения координат точки Q найдем параметрическое уравнение прямой, которой принадлежит отрезок P1Q. Орт вектора а имеет вид:
1а _
_ I х2 _ Х1 у2 _ у1
а
а
I
Легко проверяется, что орт вектора РО=р имеет вид:
Г *2 - *1 У2 - У11
а
прямой
а
а
и уравнение искомой
Х1 + х2 _ (У2 _ У)
2 а '
У1 + У2 _ (х2 _ Х1 >
(2.6)
x(t)■■
У(0 _:1Н^ _'"2 1 t
2а
Полагая в уравнениях (2.6) ґ=к1, где Н1 определяется одной из формул (2.4), находим координаты точки Q ХQ =
x(hl), УQ = У(*).
Для удобства вычислений построим два вектора Б и к. Первый из них имеет вид:
С _ |^1,5*2,53} , второй к _ |^1 к2 кз} (2 7)
где
(, 1, если проба в т. М1 сорта О, к = |
V 0, в противном случае. і = 1,2,3.
Тогда суммарная площадь S, включаемая в сорт О, легко находится по формуле:
Б = Бікі + Б2к2 + Бзкз. (2.8)
Определение координат начала и конца очередного сегмента контура, ограничивающего сорт О и принадлежащего рассматриваемому треугольнику в зависимости от числового значения переменной М = к1 + к2 + к3
Если треугольник не содержит вершин, принадлежащих многоугольной границе анализируемой области в плане, то искомые сегменты имеются лишь при М=1 и М=2. В каждом отдельном таком случае имеется ровно три способа кодирования вершин треугольника по их принадлежности к
типу О . Договоримся начало и конец сегмента контура, ограничивающего сорт О, выбирать так, чтобы оконтуривае-мая область оставалась слева при обходе контура. Учитывая вышесказанное, осуществим выбор сегментов для всех шести вариантов: М = 1
к _ 1,0,0}^ ро и 0Р3,
М = 2
р20 и 0р1
рзО и 0р2
р20 т 1с^ и
ро > Г4 1с^
рзО и 0Р1-
Если триангуляция плана содержит треугольники с тупыми углами, то вышеизложенный способ определения «ближайших районов» неприменим. Для правильного определения границ сортов сырья в пределах треугольного блока необходимо привлекать к расчетам остальные прилегающие к треугольнику пункты опробования.
Для построения универсального алгоритма используем непосредственное определение «ближай-шего района». Пусть на плоскости задано множество точек М={М1, М2,...,Мп}. «Ближайшим районом» V(Mi) точки М1 является ближайшей из точек множества М. Очевидно, что V(Mi) является выпуклым многоугольником (рис. 3).
При построении многоугольника можно поступить следующим образом. Соединим точку М1 со всеми отдельными точками множества {М} и проведем прямые через середины отрезков MMj перпендикулярно этим отрезкам. Тогда многоугольник V(Mi) совпадает с перечислением всех полученных полуплоскостей, содержащих точку М1 . Таким образом, задача заключается в нахождении координат вершин всех многоугольников V(Mi) .
Пусть теперь требуется найти координаты V(Mo) при наличии N соседних с координатами (х1,у1), (х2,у2), ..., (хм,ум).
Алгоритм определения вершин многоугольника V(M0) построим по итерационному принципу. Определим сначала вершины прямоугольника V(M0) с центром в точке Мо, содержащего внутри себя число точек из множества {М}, не менее чем К. Зададим параметр D, равный удвоенному среднему расстоянию между соседними точками М1 и Mj (М, , Mj €Е М). Выберем в качестве координат времени V(M0) следующие (рис. 4)
Wl(хo - D, уо - D); W2(хo + D, уо - D); Wз(хo + D, уо + +D); W4(хo - D, уо + D).
При обходе вершин V(M0) в порядке нумерации (против часовой стрелки) М0 всегда будет находиться слева. Если внутри прямоугольника V0 число точек из М оказывается меньше чем К, то размеры прямоугольника удвоим.
На этом этапе будем считать «ближайшим районом» к М0 прямоугольник V/). Искомый многоугольник V(M0) найдем последовательными итерациями за Р шагов, где Р -число точек внутри Vо ,
Рис. 3. Схема построения многоугольника У(М0)
исключая саму точку М0 . На каждом шаге будем отыскивать пересечение полуплоскости, ограниченной прямой Lk, проведенной перпендикулярно к середине отрезка М0Мк и содержащей точку М0 , и многоугольника Vк-1(Mo), полученного на предыдущем шаге. Результатом пересечения будет являться новый многоугольник Vк(M0) или сохранится прежний Vк-1(M0), если пересечение пусто.
Рассмотрим итерационный процесс построения V(M0), сначала в виде обобщенного алгоритма в предположении, что начальное приближение уже найдено в виде V(M0) и сформированы массивы координат точек {х1,х2,.,хр} и {у1,у2,.,уР}, при-
влекаемых к расчетам.
1. Положить к= 1, О = 4.
2. Найти параметры уравнения прямой Дк, проходящей через середину отрезка М0Мк.
3. Найти среди вершин Wi Є Vk-1(M0) слева от наиболее удаленную от Дк.
4. Если слева от отсутствуют вершины Wi, то положить (пока к < Р) к: = к+1 и перейти к п. 2.
5. Перенумеровать вершины Wi Є Vk-1(M0) так, чтобы
6. Положить і = 1.
7. Определить взаимное расположение Wi, Wi+1 и Дк.
8. Если Wi и Wi+1 слева от Дк, то положить і = і + 1 (пока і < О-1) и перейти к п. 7.
9. Если Wi слева от и Wi+1 справа от , то найти точку А пересечения и отрезка WiWi+1, Положить]': = і. Положить і =
і + 1 и перейти к п. 7.
10. Если Wi и Wi+1 справа от Дк, то положить і = і + 1 и перейти к п. 7.
11. Если Wi справа от и Wi+1 слева от , то найти точку В пересечения с отрезка WiWi+1, Положить q: = і.
12. Заменить вершину Wj точкой А.
13. Положить Т: = q -]'.
14. Если Т = 1, то между вершинами Wj+1 и Wj+2 вставить точку В. Положить О: = О + 1. Перейти к п. 16.
15. Если Т #1, то удалить из множества V вершины с номерами Wj+2, Wj+3, ,,,, Wq+2 и заменить вершину Wq точкой В. Положить О: = О - Т.
16. Если к<р , то положить к: = к + 1 и перейти к п. 2.
17. Конец алгоритма.
Перейдем к более подробному описанию отдельных предписаний алгоритма.
П. 2. Прямая проходит через середину Е(хЕ,уЕ) отрезка
М 0Мк
хЕ _
_х0 + хк
2 ' Уе 2 Направляющий вектор Б _ {/,т} прямой перпенди-
_у0 + ук
(2.9)
кулярен вектору М оМк = {а,Ь}: а = хк - хо ; Ь = ук - уо
Из условия перпендикулярности М0М и Б и учитывая ориентацию Lк (М0 должна лежать слева от Lk) имеем
Б = {?,т}={- Ь;а}, (2.10)
т.е. 1=уо - ук; т = хк - хо
Параметрические уравнения Lk имеют вид:
І х _ хе + ^ 1 у _ уЕ + mt
(2.11)
П. 3.4. Пусть массивы координат вершин Wi многоугольника Vk.1(M0) после к-1-ой итерации имеет следующий вид
х1 ,х2 ,...,хО} и |у1 у .....уо} (212)
Перед первой итерацией число вершин О = 4. После оче-
Рис. 4. Схема корректировки вершин многоугольника V(M0) на к-ом шаге итерационного процесса
редной итерации, как это видно из рис. 5, возможны три варианта корректировки массивов вершин (12) многоугольника V(M0), т.е. замены Vk.1(M0) на Vk(M0), в зависимости от числа вершин Wi, находящихся справа от прямой Дк.
Обозначим число таких вершин через Т. Очевидно, что при Т = 0 (все Mi слева от Дк). В этом случае многоугольник
Vk.1(M0) сохраняется и корректировка не производится.
W1 была наиболее удалена от прямой Дк.
Взаимное расположение точки Ші и прямой Ьк будем определять по знаку проекции Сг векторного произведения
С = 5 X М0Жі .
Сг = 1( уі - уо) - т( хі - хо )
Если Сг < 0 , то Ші расположена слева от Ьк. Если Сг > 0, то Ші - справа от Ьк .
Учитывая вышесказанное, рассмотрим алгоритм, реализующий п.п. 2 и 3 основного алгоритма.
1. Т = 0(число вершин справа от Ьк ).
2. d
1 _ V l2 + m2
+ m ; d = 0
3. і: = 1
4. Вычислить Сг
5. Если Сг > 0 то
6. а = (хі - ^0)т - (уі - ^0)ї\/,если а
> <і, то (іі=а; G=i), иначе Т=Т + 1
7. Если і < Q ,то і = і + 1, перейти к п. 4.
8. Конец.
П. 5. Если в результате вычислений значение Т переменной будет отлично от нуля, необходимо перенумеровать вершины по правилу:
х1 ,х2 }={ х'о ,хо+1,...,хё ,х1 ,х2 ,-,^-1
УрУ2’’Ув | _ {Уа>Уа+1--У0>У1>У2--У0-1
П.п. 6-11. Эти предписания алгоритма выполняются, если параметр Т Ф 0, т.е. пересечение
полуплоскости, ограниченной ориентированной прямой Lk, и многоугольника Vk-1(Mo) не пусто и необходима корректировка вершин многоугольника. Способ нумерации вершин
Vk-1(Mo) позволяет осуществить эту корректировку рациональным образом.
Действительно, при обходе вершин многоугольника в порядке следования вершин многоугольника и проверке на пересечения прямой Lk с очередной стороной многоугольника начало и конец нового ребра АВ появляются в необходимой последовательности. Для конфигурации, изображенной на рис. 4, новые вершины обнаруживаются на втором и четвертом шаге. Корректировка многоугольника сводится к замене вершины Wз на А и W4 на В , число вершин Q при этом не изменится (Т = 2).
Анализ взаимного расположения двух соседних вершин осуществляется также, как и в п.п. 3, 4 основного алгоритма. Для нахождения координат точки пересечения прямой Lk, параметрические уравнения которой имеют вид (11), со стороной WiWi+1 запишем уравнение прямой, которой принадлежит этот отрезок.
(х - х)(у+1 - у) - (у - у)(х+1 -х) = 0
Рис. 5. Варианты корректировки «ближайшего района» на к-ом шаге итерационного процесса
Решая это уравнение совместно с уравнением (11), находим координаты точки пересечения:
ф Х0 + Хк „ , ф
x
- +
(y0 - yk)t
У* _ У0 + Ук + (x0 - xk)t*
(2.13)
где
t* _ (xi+1 - xi)(yE - yi) - (yi+1 - yi)(xE - xi) (У0 - Ук)( Уі+1- Уі)+(xo - xk)( xi+l- xi)
1. і: = 1; G: = 0 (число пересечений Ьк с Vk-l(M0))
2. С!: = (уо - Ук) (Уі -уЕ) - (Хк- хо )(х,- - хе)
3. С2: = (уо - Ук) (Уі+і уЕ) - (Хк - Хо) (Х,+! - хе)
4. Если Сі > 0 и С2 > 0, то перейти к п. 5 4’ Если Сі > 0 и С2 <0, то вычислить х и у по формулам (13) положить хА = х ; уА = у ; ] = і ; G: = G + 1 Перейти к п. 5
4’’ Если Сі <0 и С2 > 0 , то вычислить х и у по формулам (13) положить хв = х ; ув = у ; q = і ; G: = G + 1
5. Если і < Q - і и G < 2 , то положить і: = і + 1 и перейти к п. 2
6. Конец.
В результате вычислений в соответствии с приведенным алгоритмом переменная принимает значение G = 0 либо G = q . В первом случае пересечений нет и многоугольник Vk(Mо) имеет те же вершины, что и многоугольник Vk.1(Mо), полученный на предыдущей итерации.
Площадь многоугольника, полученного на последнем шаге итерационного процесса, находится по формуле:
5 = 1 5 = 2
(x1 - x2)(y1 + y2) + (x2 - x3)(y2 + у3) +
• •• + (xn -xi)(Уп + yi)
(2.14)
Моделирование поверхностей, ограничивающих объекты, связано с построением интерполяционного бикубического сплайна S(x,y) и сводится к нахождению коэффициентов бикубического многочлена в каждой ячейке сети Di в виде [8, 9]
S(x,y) = Sif = | | akJ/x' -x)k(y-yJ)1 (8Л5)
k=0l=0 J
S(x,y) имеет непрерывные первые производные и вторые частные производные, S(xi,yj)=cij , т.е. принимают заданные
значения в узлах сети; d2S ^ = о, где V - внешняя нормаль dv2
к границе L области D.
СПИСОК ЛИТЕРАТУРЫ
1. Бедрина Г.П., Гончарук В.К., Месхи Н.Ж. Информационное обеспечение геоло-го-экологической оценки горно-
обогатительных предприятий // Горный информационно-аналитический бюллетень. - М.: МГГУ. - 1995. - Вып. 5.
2. Гальперин А.М., Ферстер В.Ф., Шеф Х.-Ю. Техногенные массивы и охрана окружающей среды. - М.: МГГУ, 1997.
3. Девис Дж. Статистический анализ данных в геологии. В 2-х кн. / под ред. Д.А. Родионова. М.: Недра, 1990.
4. Ермолов В.А., Быховец А.Н. Геоло-го-технологическая оценка техногенных месторождений в системе рационального природопользования // Экологические проблемы горного производства, переработка и размещение отходов. - М.: МГГУ, 1997.
5. Ермолов В.А., Мосейкин В.В. Научно-методические аспекты оценки и модели-
рования техногенных месторождений // Изв. Вузов. Геология и разведка. 1997, № 1.
6. Ермолов В.А., Быховец А.Н. Освоение Ковдорского техногенного месторождения. Горный журнал.
7. Ермолов В.А., Бедрина Г.П., Зер-вандова В.П., Мосейкин В.В. Теория и практика моделирования и ресурсной оценки техногенных месторождений. - Изв. вузов. Геология и разведка, 1998. № 6. С. 65-72.
8. Ершов В.В., Дремуха А.С., Трость ВМ. и др. Автоматизация геологомаркшейдерских графических работ. М.: Недра, 1991.
9. Ершов В.В. Геолого-маркшей-дерское обеспечение управления качеством руд. М.: Недра, 1986.
10. Математика в социологии: моделирование и обработка информации./ Перевод с англ. Под ред. А.Г. Аганбегяна, Ф.М. Бородкина. М.: «Мир». 1977.
шшшу'';'''''''''' /
1Р % / жж-М--■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ У Ермолов, Валерий Александрович докюр технических наук, профессор кафедры геолог ии МГГУ. у Зервандова Валентина Павловна ведущий инженер-программист кафедры геологии МГГУ. Быховец Александр Николаевич манный геолог Ковдорскою ГОКа.
