Ця тема оформилась у мене цілком спонтанно. Причина написання статті у тому, що поступово зростає кількість халтури, яка прикривається гордим ім’ям "геостатистичний аналіз". Халтурять усі: геологи, які виконують первинне документування; геологи, які проектують та заповнюють бази даних; геологи, які за допомогою модних на сьогодні методів, намагаються підрахувати запаси родовищ.
Біда в тому, що більшість споживачів подібної халтури (читай - керівників та інженерних співробітників підприємств гірничодобувної промисловості) навіть не розуміє особливостей застосування подібних методик та їхніх обмежень. В результаті ми отримуємо замкнене коло - через декілька років на підприємстві стає зрозумілим, що побудована модель зовсім не відповідає реальним показникам геологічного об’єкту. Підприємство звертається за допомогою до іншого виконавця. Інший виконавець розказує - як все погано і неправильно, і в решті-решт підсовує замовнику таку-саму халтуру, тільки по іншому розфарбовану. Якщо Вам неправильно інтерпретували результати аналізу сечі або крові, і декілька років лікували від хибного діабету або гастриту - якими будуть Ваші дії? А якщо Вам неправильно інтерпретували результати опробування і Ви декілька років намагались видобувати руду там, де її взагалі немає - якими будуть Ваші дії у цьому випадку?
Біда в тому, що більшість споживачів подібної халтури (читай - керівників та інженерних співробітників підприємств гірничодобувної промисловості) навіть не розуміє особливостей застосування подібних методик та їхніх обмежень. В результаті ми отримуємо замкнене коло - через декілька років на підприємстві стає зрозумілим, що побудована модель зовсім не відповідає реальним показникам геологічного об’єкту. Підприємство звертається за допомогою до іншого виконавця. Інший виконавець розказує - як все погано і неправильно, і в решті-решт підсовує замовнику таку-саму халтуру, тільки по іншому розфарбовану. Якщо Вам неправильно інтерпретували результати аналізу сечі або крові, і декілька років лікували від хибного діабету або гастриту - якими будуть Ваші дії? А якщо Вам неправильно інтерпретували результати опробування і Ви декілька років намагались видобувати руду там, де її взагалі немає - якими будуть Ваші дії у цьому випадку?
У цій статті ми спробуємо розібратись із декількома простими, але доволі дієвими методами "виведення на чисту воду" виконавців. Слід зауважити на тому, що деякі з виконавців навіть не можуть внятно роз’яснити - які методи обчислення вони використовували і якими методами перевіряли отримані моделі. Звіти, які вони подають замовнику - здебільшого являють собою стандартний набір не пов’язаних між собою картинок і таблиць.
Зазвичай, виконавці уповають на порівняння центральних моментів (середньої арифметичної, медіани, дисперсії) двох вибірок: даних опробування і результатів моделювання. Але, я вважаю, не треба пояснювати читачам те, що центральні моменти можуть співпадати навіть для даних, отриманих за допомогою генератору псевдовипадкових послідовностей. Більш показовими будуть функції розподілу значень для даних опробування і результатів моделювання, але і вони не дадуть вам гарантії просторової відповідності між первинними даними і кінцевими результатами. Нам треба дослідити - як у тривимірному просторі співпадають дані опробування і результати моделювання. Наскільки зміщені показники кінцевої блокової моделі родовища або рудного покладу.
Це можна доволі швидко перевірити не витрачаючи сили та час на мудрування із інтерфейсами геоінформаційних систем. Все, що нам треба - це експортовані у формат CSV блокові моделі із обчисленими параметрами вмістів компонентів та дані просторового розташування інтервалів опробування у просторі (також у форматі CSV). До речі - коли на запитання "в якому вигляді Вам надати дані", Ви відповідаєте: "дайте Вашу блочку в форматі CSV" - на Вас починають дивитись як на чаклуна, який щойно зламав сервер прибульців та отримав доступ до потаємніших знань! Для порівняння функцій розподілу ми використаємо базові можливості R, а для дослідження просторової сумісності результатів - пакет spatstat, який не входить до базового набору. Встановлювати пакет spatstat краще за допомогою прямого завантаження з CRAN - цей пакет має доволі багато залежностей, які також не входять у базовий набір R.
Для аналізу ми використаємо результати опробування на компонент "Ікс" та результати інтерполяції вмістів компоненту у блоковій моделі. Файл з первинними даними (результати опробування) - proby.csv, файл блокової моделі - bl_model.csv.
Тож почнемо! У першу чергу перевіримо - чи співпадають функції розподілу - побудуємо графіки кумулятивних функцій розподілу:
>proby<- read.table("proby.csv", header=T, sep=";", dec=",")#завантажуємо файл із даними опробування
>bl_model<- read.table("bl_model.csv", header=T, sep=";", dec=",")#завантажуємо файл із блоковою моделлю
>summary(proby)#дивимось загальну статистику для даних опробування
NUM XX YY ZZ
Min. : 1.00 Min. :6837 Min. :20147 Min. :-267.77
1st Qu.: 37.25 1st Qu.:6846 1st Qu.:20223 1st Qu.:-169.40
Median : 73.50 Median :7046 Median :20298 Median :-115.06
Mean : 73.50 Mean :7023 Mean :20302 Mean :-118.00
3rd Qu.:109.75 3rd Qu.:7194 3rd Qu.:20379 3rd Qu.: -59.87
Max. :146.00 Max. :7344 Max. :20449 Max. : -5.28
Iks
Min. :0.000001
1st Qu.:0.013957
Median :0.052824
Mean :0.101881
3rd Qu.:0.144562
Max. :0.827947
>summary(bl_model)#дивимось загальну статистику для блокової моделі
NUM NUM_BL XX YY
Min. : 1 Min. : 0 Min. :6695 Min. :19965
1st Qu.: 3524 1st Qu.: 3523 1st Qu.:6815 1st Qu.:20195
Median : 7048 Median : 7046 Median :6935 Median :20285
Mean : 7048 Mean : 7046 Mean :6959 Mean :20272
3rd Qu.:10571 3rd Qu.:10570 3rd Qu.:7085 3rd Qu.:20365
Max. :14094 Max. :14093 Max. :7435 Max. :20475
ZZ RAZM_X RAZM_Y RAZM_Z Iks
Min. :-395.0 Min. :10 Min. :10 Min. :10 Min. :0.0001564
1st Qu.:-205.0 1st Qu.:10 1st Qu.:10 1st Qu.:10 1st Qu.:0.0182835
Median :-125.0 Median :10 Median :10 Median :10 Median :0.0457895
Mean :-137.7 Mean :10 Mean :10 Mean :10 Mean :0.0914940
3rd Qu.: -65.0 3rd Qu.:10 3rd Qu.:10 3rd Qu.:10 3rd Qu.:0.1297130
Max. : -5.0 Max. :10 Max. :10 Max. :10 Max. :0.6591007
#нас цікавлять мінімальні та максимальні значення вмістів компоненту Ікс в обох вибірках - нам треба накласти кумулятивні криві розподілу на один графік
>plot(ecdf(proby$Iks), col="red", xlim=c(0,0.85), main="Кумулятивні криві розподілу вмісту компоненту Ікс")
>plot(ecdf(bl_model$Iks), col="blue", add=T)
>legend("bottomright", legend=c("Дані опробування","Блокова модель"),pch=20, col=c("red","blue"), cex=1.2)
>proby<- read.table("proby.csv", header=T, sep=";", dec=",")#завантажуємо файл із даними опробування
>bl_model<- read.table("bl_model.csv", header=T, sep=";", dec=",")#завантажуємо файл із блоковою моделлю
>summary(proby)#дивимось загальну статистику для даних опробування
NUM XX YY ZZ
Min. : 1.00 Min. :6837 Min. :20147 Min. :-267.77
1st Qu.: 37.25 1st Qu.:6846 1st Qu.:20223 1st Qu.:-169.40
Median : 73.50 Median :7046 Median :20298 Median :-115.06
Mean : 73.50 Mean :7023 Mean :20302 Mean :-118.00
3rd Qu.:109.75 3rd Qu.:7194 3rd Qu.:20379 3rd Qu.: -59.87
Max. :146.00 Max. :7344 Max. :20449 Max. : -5.28
Iks
Min. :0.000001
1st Qu.:0.013957
Median :0.052824
Mean :0.101881
3rd Qu.:0.144562
Max. :0.827947
>summary(bl_model)#дивимось загальну статистику для блокової моделі
NUM NUM_BL XX YY
Min. : 1 Min. : 0 Min. :6695 Min. :19965
1st Qu.: 3524 1st Qu.: 3523 1st Qu.:6815 1st Qu.:20195
Median : 7048 Median : 7046 Median :6935 Median :20285
Mean : 7048 Mean : 7046 Mean :6959 Mean :20272
3rd Qu.:10571 3rd Qu.:10570 3rd Qu.:7085 3rd Qu.:20365
Max. :14094 Max. :14093 Max. :7435 Max. :20475
ZZ RAZM_X RAZM_Y RAZM_Z Iks
Min. :-395.0 Min. :10 Min. :10 Min. :10 Min. :0.0001564
1st Qu.:-205.0 1st Qu.:10 1st Qu.:10 1st Qu.:10 1st Qu.:0.0182835
Median :-125.0 Median :10 Median :10 Median :10 Median :0.0457895
Mean :-137.7 Mean :10 Mean :10 Mean :10 Mean :0.0914940
3rd Qu.: -65.0 3rd Qu.:10 3rd Qu.:10 3rd Qu.:10 3rd Qu.:0.1297130
Max. : -5.0 Max. :10 Max. :10 Max. :10 Max. :0.6591007
#нас цікавлять мінімальні та максимальні значення вмістів компоненту Ікс в обох вибірках - нам треба накласти кумулятивні криві розподілу на один графік
>plot(ecdf(proby$Iks), col="red", xlim=c(0,0.85), main="Кумулятивні криві розподілу вмісту компоненту Ікс")
>plot(ecdf(bl_model$Iks), col="blue", add=T)
>legend("bottomright", legend=c("Дані опробування","Блокова модель"),pch=20, col=c("red","blue"), cex=1.2)
 |
| Кумулятивні криві розподілу вмісту компоненту Ікс для первинних даних опробування та блокової моделі |
Як видно на рисунку - обидві вибірки характеризуються схожими функціями розподілу вмістів компоненту. Але для більш точного порівняння краще використати не власний глазомір, а який-небудь статистичний критерій. Спробуємо застосувати критерій Ансарі-Бредлі:
>ansari.test(proby$Iks, bl_model$Iks)
Ansari-Bradley test
data: proby$Iks and bl_model$Iks
AB = 482610.5, p-value = 0.1319
alternative hypothesis: true ratio of scales is not equal to 1
Як бачимо - досягнутий рівень значущості не дозволяє нам спростувати нульову гіпотезу, і ми погоджуємось із тезою про подібність функцій розподілу (не забувайте, що правильно звучить не "подібність", а "відсутність відмін"). До речі: критерій Смірнова (функція ks.test) показує досягнутий рівень значущості p-value=0.303, що також не дозволяє нам спростувати нульову гіпотезу.
>ansari.test(proby$Iks, bl_model$Iks)
Ansari-Bradley test
data: proby$Iks and bl_model$Iks
AB = 482610.5, p-value = 0.1319
alternative hypothesis: true ratio of scales is not equal to 1
Як бачимо - досягнутий рівень значущості не дозволяє нам спростувати нульову гіпотезу, і ми погоджуємось із тезою про подібність функцій розподілу (не забувайте, що правильно звучить не "подібність", а "відсутність відмін"). До речі: критерій Смірнова (функція ks.test) показує досягнутий рівень значущості p-value=0.303, що також не дозволяє нам спростувати нульову гіпотезу.
Тепер перейдемо до аналізу просторової сумісності даних опробування та результатів обчислення якості в блоковій моделі. Для цього ми зробимо із блокової моделі вибірку найближчих до первинних проб блоків, і порівняємо значення у цих пробах та блоках. Зараз деякі "геостатистики" почнуть кидати у мене цеглинами, запевняючи що ці методи інтерполяції - "згладжувальні", і результати повинні відрізнятись від первинних даних. Звісно - вони повинні відрізнятись. Але коли вони відрізняються майже у два рази - є над чим замислитись. Для вибірки найближчих блоків ми використаємо можливості пакету spatstat і невеличкий цикл, який дозволить нам знайти позиції цих блоків у первинній таблиці. Зверніть увагу на те, що ми оперуємо координатною системою із прямокутним розташуванням осей, і всі вісі мають однаковий масштаб. У принципі - ми повинні поміняти місцями вісі X та Y, оскільки в геодезичних системах координат вони міняються. Але на даний момент це для нас не критично.
>install.packages("spatstat")#встановлюємо пакет spatstat
>library(spatstat)#завантажуємо пакет spatstat
>proby.pp3<- pp3(proby$XX, proby$YY, proby$ZZ, box3(c(6800, 7350), c(20100,20500), c(-300,0)))#будуємо тривимірний "бокс" із даними опробування
>bl_model.pp3<- pp3(bl_model$XX, bl_model$YY, bl_model$ZZ, box3(c(6650, 7500), c(19950,20500), c(-400,0)))#будуємо тривимірний "бокс" із даними блокового моделювання
>vidst<- crossdist(proby.pp3, bl_model.pp3)#обчислюємо відстані від кожної проби до кожного блоку
Результатом виконання функції crossdist є матриця, яка складається із 146 рядків (стільки у нас первинних проб, від яких ми вимірювали відстані) та 14094 стовпців (стільки блоків у нашій блоковій моделі). За допомогою функцій min, max та median можна подивитись основні параметри отриманої матриці:
>min(vidst)
[1] 0.6394529
>max(vidst)
[1] 857.6206
>median(vidst)
[1] 254.9157
Для більшої наглядності ми використаємо функцію hist, щоб відобразити графічно розподіл відстаней між блоками моделі та первинними пробами:
>hist(vidst, breaks=100, col="darkorange", main="Гістограма відстаней між первинними\nпробами та центрами мас блоків", xlab="відстані", ylab="частість")
>install.packages("spatstat")#встановлюємо пакет spatstat
>library(spatstat)#завантажуємо пакет spatstat
>proby.pp3<- pp3(proby$XX, proby$YY, proby$ZZ, box3(c(6800, 7350), c(20100,20500), c(-300,0)))#будуємо тривимірний "бокс" із даними опробування
>bl_model.pp3<- pp3(bl_model$XX, bl_model$YY, bl_model$ZZ, box3(c(6650, 7500), c(19950,20500), c(-400,0)))#будуємо тривимірний "бокс" із даними блокового моделювання
>vidst<- crossdist(proby.pp3, bl_model.pp3)#обчислюємо відстані від кожної проби до кожного блоку
Результатом виконання функції crossdist є матриця, яка складається із 146 рядків (стільки у нас первинних проб, від яких ми вимірювали відстані) та 14094 стовпців (стільки блоків у нашій блоковій моделі). За допомогою функцій min, max та median можна подивитись основні параметри отриманої матриці:
>min(vidst)
[1] 0.6394529
>max(vidst)
[1] 857.6206
>median(vidst)
[1] 254.9157
Для більшої наглядності ми використаємо функцію hist, щоб відобразити графічно розподіл відстаней між блоками моделі та первинними пробами:
>hist(vidst, breaks=100, col="darkorange", main="Гістограма відстаней між первинними\nпробами та центрами мас блоків", xlab="відстані", ylab="частість")
 |
| Гістограма відстаней між блоками блокової моделі та первинними пробами |
За бажанням - можна подивитись розташування первинних проб у просторі за допомогою команди plot(proby.pp3). В принципі - можна відобразити і просторове розташування блоків блокової моделі (точніше - їхніх центроїдів), але через надмірну концентрацію точок Ви зможете побачити лише зовнішні контури блокової моделі.
 |
| Тривимірне відображення розташування проб по свердловинах |
Тепер нам треба знайти найближчий блок для кожної із проб. Спочатку ми приведемо програмний код, а потім пояснимо - що і для чого ми робили. Кінцевим результатом приведеного нижче коду буде таблиця найближчих до первинних проб блоків блокової моделі. За своєю структурою (кількість та назви стовпців) ця таблиця буде подібна до таблиці повного переліку блоків bl_model.
>znach<- numeric(nrow(vidst))
>stovp<- integer(nrow(vidst))
>for (i in 1:nrow(vidst)){#починаємо цикл
+ znach[i]<- vidst[i,1]
+ for (j in 1:ncol(vidst)){#починаємо вкладений цикл
+ if (vidst[i,j]<znach[i]) {znach[i]<- vidst[i,j]; stovp[i]<- j}#заносимо логічну умову для порівняння відстаней від блоків до кожної проби
+ }}#закінчуємо обидва цикли
>vidst_table<- data.frame(znach=znach, stovp=stovp)#формуємо таблицю значень відстаней і номерів блоків
>names(bl_model)#перевіряємо імена стовпців у таблиці блокової моделі
[1] "NUM" "NUM_BL" "XX" "YY" "ZZ" "RAZM_X" "RAZM_Y" "RAZM_Z"
[9] "Iks"
>ind.bl<- match(vidst_table$stovp, bl_model$NUM)#створюємо вектор відповідності для порівняння між таблицями
>bl_poruch_pr<- bl_model[ind.bl,]#на основі вектору відповідності робемо вибірку з таблиці bl_model і заносимо її до таблиці bl_poruch_pr
>head(bl_poruch_pr)#перевіряємо результат (дивимось "шапку" таблиці)
NUM NUM_BL XX YY ZZ RAZM_X RAZM_Y RAZM_Z Iks
4511 4511 4510 6845 20145 -205 10 10 10 0.0008770
4514 4514 4513 6845 20155 -245 10 10 10 0.2451710
4509 4509 4508 6845 20145 -225 10 10 10 0.2102884
4510 4510 4509 6845 20145 -215 10 10 10 0.0136489
9527 9527 9526 7035 20215 -155 10 10 10 0.0041822
12347 12347 12346 7195 20385 -55 10 10 10 0.4323209
Давайте тепер більш детально розпишемо послідовність наших дій. По-перше: ми створили два тимчасових вектори: znach та stovp. Вони нам знадобляться для введення в дію циклу, який буде перебирати таблицю (матрицю) vidst і вибирати з неї найближчий блок для кожної із первинних проб. Далі ми починаємо цикл, і робимо вкладений цикл з логічною умовою "якщо елемент таблиці відстаней менше за первинний елемент вектору значень, тоді записати нове значення із таблиці відстаней, а в вектор номерів стовпців записати положення цього значення - номер стовпця матриці відстаней". По закінченні цикла ми робимо зведену таблицю vidst_table значень і номерів стовпців (ці номери взяті із матриці відстаней vidst і співпадають із номерами блоків у полі таблиці bl_model$NUM).
>znach<- numeric(nrow(vidst))
>stovp<- integer(nrow(vidst))
>for (i in 1:nrow(vidst)){#починаємо цикл
+ znach[i]<- vidst[i,1]
+ for (j in 1:ncol(vidst)){#починаємо вкладений цикл
+ if (vidst[i,j]<znach[i]) {znach[i]<- vidst[i,j]; stovp[i]<- j}#заносимо логічну умову для порівняння відстаней від блоків до кожної проби
+ }}#закінчуємо обидва цикли
>vidst_table<- data.frame(znach=znach, stovp=stovp)#формуємо таблицю значень відстаней і номерів блоків
>names(bl_model)#перевіряємо імена стовпців у таблиці блокової моделі
[1] "NUM" "NUM_BL" "XX" "YY" "ZZ" "RAZM_X" "RAZM_Y" "RAZM_Z"
[9] "Iks"
>ind.bl<- match(vidst_table$stovp, bl_model$NUM)#створюємо вектор відповідності для порівняння між таблицями
>bl_poruch_pr<- bl_model[ind.bl,]#на основі вектору відповідності робемо вибірку з таблиці bl_model і заносимо її до таблиці bl_poruch_pr
>head(bl_poruch_pr)#перевіряємо результат (дивимось "шапку" таблиці)
NUM NUM_BL XX YY ZZ RAZM_X RAZM_Y RAZM_Z Iks
4511 4511 4510 6845 20145 -205 10 10 10 0.0008770
4514 4514 4513 6845 20155 -245 10 10 10 0.2451710
4509 4509 4508 6845 20145 -225 10 10 10 0.2102884
4510 4510 4509 6845 20145 -215 10 10 10 0.0136489
9527 9527 9526 7035 20215 -155 10 10 10 0.0041822
12347 12347 12346 7195 20385 -55 10 10 10 0.4323209
Давайте тепер більш детально розпишемо послідовність наших дій. По-перше: ми створили два тимчасових вектори: znach та stovp. Вони нам знадобляться для введення в дію циклу, який буде перебирати таблицю (матрицю) vidst і вибирати з неї найближчий блок для кожної із первинних проб. Далі ми починаємо цикл, і робимо вкладений цикл з логічною умовою "якщо елемент таблиці відстаней менше за первинний елемент вектору значень, тоді записати нове значення із таблиці відстаней, а в вектор номерів стовпців записати положення цього значення - номер стовпця матриці відстаней". По закінченні цикла ми робимо зведену таблицю vidst_table значень і номерів стовпців (ці номери взяті із матриці відстаней vidst і співпадають із номерами блоків у полі таблиці bl_model$NUM).
Але, ми з Вами отримали дані номерів блоків та значення їхніх відстаней від первинних проб. А нам потрібні показники вмістів хімічних елементів та інші властивості, які містяться в цих блоках. Щоб отримати ці дані ми за допомогою функції match створюємо вектор відповідності ind.bl. У цей вектор будуть занесені номери блоків, властивості яких потрібні нам для подальшої роботи. Далі, на основі цього вектору ми "витягуємо" (з усіма полями) із таблиці bl_model лише ті блоки, які нам потрібні, і заносимо їх до таблиці bl_poruch_pr ("Блоки поруч із пробами"). Нам залишилось лише перевірити - чи перенесли ми дані в таблицю. Для цього ми використовуємо функцію head, яка виводить нам на екран заголовки полів та перші шість записів.
Отже, ми з Вами отримали таблицю, яка містить 146 блоків - найближчих до 146 первинних проб. Тепер ми можемо порівняти статистичні показники просторово суміщених даних, і подивитись - наскільки в реальності не співпадають дані моделювання із первинними значеннями (даними опробування). Ми не будемо наводити тут всі можливі варіанти порівняння даних і побудуємо лише графік кумулятивних функцій розподілу для значень вмісту компонента "Ікс" у первинних і змодельованих даних. Не забудьте перевірити граничні показники для обох таблиць, щоб задати обмеження для побудови красивого графіку:
>plot(ecdf(proby$Iks), col="red", xlim=c(0,0.85), main="Кумулятивні криві розподілу компоненту Ікс\nв пробах і найближчих блоках")
>plot(ecdf(bl_poruch_pr$Iks), col="blue", add=T)
>legend("bottomright", legend=c("Дані опробування", "Найближчі блоки"), fill=c("red", "blue"))
>plot(ecdf(proby$Iks), col="red", xlim=c(0,0.85), main="Кумулятивні криві розподілу компоненту Ікс\nв пробах і найближчих блоках")
>plot(ecdf(bl_poruch_pr$Iks), col="blue", add=T)
>legend("bottomright", legend=c("Дані опробування", "Найближчі блоки"), fill=c("red", "blue"))
 |
| Кумулятивні криві розподілу вмісту компоненту Ікс для первинних даних опробування та найближчих до них блоків блокової моделі |
Як видно з рисунку, розподіл вмістів у блоках в цілому відповідає первинним даним. Разом із тим, спостерігаються деякі невідповідності: максимальний вміст у найближчому блоці майже в два рази менший, ніж у пробі (кінцеві точки на обох функціях); на графіку значень по блоках спостерігаються "зриви" функції, і функція являє собою набір декількох відособлених ділянок. Подібні "зриви" можуть свідчити про погану (невідповідну, нерівномірну) інтерполяцію, або про ручне виправлення кінцевих даних. Але за значних розмірів блокової моделі ручне виправлення малоймовірне і ми зупиняємось на висновку про недостатню увагу моделювальників до особливостей будови рудного тіла та розподілу в ньому компонентів.
Немає коментарів:
Дописати коментар