CC BY
32
Оригинальная статья / Original article УДК: 550.83:519.6
СРАВНИТЕЛЬНЫЙ АНАЛИЗ ДВУХ МЕТОДОВ КОЛИЧЕСТВЕННОЙ ИНТЕРПРЕТАЦИИ АНОМАЛИЙ ПОТЕНЦИАЛЬНЫХ ГЕОФИЗИЧЕСКИХ ПОЛЕЙ
© В.С. Канайкин1, А.А. Субботин2, А.И. Булнаев3
1-3Иркутский национальный исследовательский технический университет, Российская Федерация, 664074, г. Иркутск, ул. Лермонтова, 83.
РЕЗЮМЕ. Цель данной статьи - анализ возможностей регрессионных методов количественной интерпретации аномалий потенциальных геофизических полей. Сравниваются два способа решения обратной задачи гравиразведки: методом наименьших квадратов с применением параметра регуляризации и методом сингулярного разложения матриц. Методы. Для успешного решения поставленной задачи была сформирована петрофизическая модель с известным распределением избыточной плотности аномалиеобразующих объектов (прямоугольных призм). На основе данной модели путем решения прямой задачи произведен расчет гравитационного поля, в результате чего получена схематическая физико-геологическая модель, которая соответствует массиву ортоамфиболитов, контролирующему размещение редкометалльных пегматитов. Результаты. Решена обратная задача гравиразведки методом наименьших квадратов и сингулярного разложения для двух ситуаций: первая - интерпретируемое гравитационное поле Ад, не осложненное ошибками наблюдений, вторая - гравитационное поле Ад1, осложнено помехой порядка ±3%. В итоге решения обратной задачи получены значения избыточной плотности аномалиеобразующих объектов. Решение задачи для поля, не осложненного помехой, восстанавливается достаточно точно, а для поля, осложненного помехой, решение «срывается» и становится непригодным для практического использования. Решение обратной задачи методом наименьших квадратов удается восстановить при использовании параметра регуляризации А.Н. Тихонова. Метод сингулярного разложения матриц также позволил восстановить решение, но для этого пришлось сформировать псевдообратную матрицу с учетом точности выполняемых работ. Выводы. Решение обратной задачи гравиразведки показывает, что для заданных модельных условий оба метода интерпретации дают схожие результаты вычислений избыточной плотности аномальных объектов. Ключевые слова: количественная интерпретация гравитационного поля, метод наименьших квадратов и сингулярное разложение, регуляризация решений обратной задачи.
Формат цитирования: Канайкин В.С., Субботин А.А., Булнаев А.И. Сравнительный анализ двух метдов количественной интерпретации аномалий потенциальных геофизических полей // Известия Сибирского отделения Секции наук о Земле Российской академии естественных наук. Геология, разведка и разработка месторождений полезных ископаемых. 2017. Т. 40. № 2. С. 95-100.
1 Канайкин Виктор Степанович, доцент кафедры технологии геологической разведки, е-mail: [email protected]
Viktor S. Kanaikin, Associate Professor of the Department of Geological Exploration Technology, е-mail: [email protected]
2Субботин Анатолий Александрович, студент, е-mail: [email protected]
Anatoly A. Subbotin, Student, е-mail: [email protected]
3Булнаев Андрей Иосифович, член-корреспондент Российской академии естественных наук, доктор геолого-минералогических наук, профессор кафедры технологии геологической разведки, e-mail: [email protected]
Andrey I. Bulnaev, Corresponding Member of the Russian Academy of Natural Sciences, Doctor of Geological and Mineralogical sciences, Professor of the Department of Geological Exploration Technology, e-mail: [email protected]
Известия Сибирского отделения Секции наук о Земле РАЕН. ISSN Геология, разведка и разработка месторождений полезных ископаемых Т. 40, № 2
2541-9455 Proceedings of the Siberian Department of the Section of Earth Sciences RANS.
Geology, Exploration and Development of Mineral Deposits Vol. 40, No. 2
COMPARATIVE ANALYSIS OF TWO QUANTITATIVE METHODS FOR POTENTIAL GEOPHYSICAL FIELD ANOMALY INTERPRETATION
V.S. Kanaikin, A.A. Subbotin, A.I. Bulnaev
Irkutsk National Research Technical University,
83, Lermontov St., Irkutsk, 664074, Russian Federation.
ABSTRACT. The purpose of this article is to analyze the capabilities of regression methods for quantitative interpretation of potential geophysical field anomalies. Two solution methods for the gravity inversion problem are compared: the method of least squares with the use of the regularization parameter and the method of singular matrix decomposition. Methods. A petrophysical model with the known distribution of excess density of anomaly-forming objects (rectangular prisms) has been formed in order to solve the set task. This model enabled to calculate a gravitational field through solving a direct problem. As a result, we obtained a schematic physico-geological model corresponding to a massif of ortho-amphibolites, which controls the location of rare metal pegmatites. Results. The inverse gravity problem is solved by the method of least squares and singular matrix decomposition for two situations: the first is the interpreted gravitational field Ag uncomplicated by observation errors, the second is the gravitational field Ag1 complicated by an interference of about ± 3%. Having solved the inverse problem, we obtained the values of the excess density of anomaly-forming objects. The problem solution for the inference uncomplicated field is restored quite accurately, whereas the solution for the field complicated by an interference fails and becomes unsuitable for practical use. The solution of the inverse problem by the method of least squares can be reconstructed using A.N. Tikhonov's regularization parameter. The method of singular matrix decomposition also allows to restore the solution, but requires the formation of a pseudo-inverse matrix taking into account the accuracy of the work performed. Conclusions. The solution of the inverse gravity survey problem shows that for the specified model conditions both interpretation methods give similar results in the calculation of anomalous object excess density that is proved by the values of mean square errors.
Keywords: quantitative interpretation of the gravitational field, least-squares method and singular value decomposition, regularization of inverse problem solution
For citation: Kanaikin V.S., Subbotin A.A., Bulnaev A.I. Comparative analysis of two quantitative methods for potential geophysical field anomaly interpretation. Proceedings of the Siberian Department of the Section of Earth Sciences of the Russian Academy of Natural Sciences. Geology, Exploration and Development of Mineral Deposits. 2017, vol. 40, no. 2, рр. 95-100. (In Russian).
Введение
Методов количественной интерпретации аномалий потенциальных полей достаточно много [3, 5, 6, 7]. В статье сравниваются два способа количественной интерпретации гравитационных аномалий, основанных на регрессионном анализе. Первый способ относится к методу наименьших квадратов с применением параметра регуляризации [8, 11], второй использует свойства сингулярного разложения матриц [9, 10, 11].
Материал и методы исследования
Решение обратной задачи рассматривается на модельном примере и нацелено на определение морфологии массива ортоамфиболитов, контролирующего размещение редкометалльных
пегматитов [1]. Массив вытянут в горизонтальном направлении (рис. 1), его пе-трофизическая модель состоит из двадцати прямоугольных призм (к = 20), вертикальные и горизонтальные мощности которых соответственно равны 1000 и 300 м. Избыточные плотности Дст призм известны. Координата точки верхнего левого угла набора призм относительно профиля измерения равна 2000 м, глубина h - 10 м.
Информация о петрофизической модели позволила сформировать физико-геологическую модель объекта путем решения прямой задачи гравираз-ведки, т.е. вычисления поля G по профилю наблюдения длиной 8000 м, что соответствует n = 80 пикетам при шаге съемки 100 м.
Известия Сибирского отделения Секции наук о Земле РАЕН. P« Геология, разведка и разработка месторождений полезных ископаемых Т. 40, № 2 ISSN
Proceedings of the Siberian Department of the Section of Earth Sciences RANS. 2541-9455
Geology, Exploration and Development of Mineral Deposits Vol. 40, No. 2
°0 20 40 60 35
v 0.35 v v 0.23 v V 0.30 v v 0.35 v
0.08 v 0.44 v „ 0.41 v 0.08
0.08 V 0.36 у V V 0.33 „ V 0.07
0.12 V 0.22 v V V 0-21 „ V 0.08
] 1 | VV | 2 | Q-35 | 3
Рис. 1. Схематический разрез массива ортоамфиболитов и гравитационное поле над ним:
1 - гравитационное поле G; 2 - ортоамфиболиты; 3 - избыточные плотности Ао призм, г/см3 Fig. 1. Schematic section of the ortho-amphibolite massif and the gravitational field above it:
1 - gravitational field G; 2 - ortho-amphibolites 3 - excess densities Ао of prisms, g/sm3
Решение обратной задачи выполнено в два этапа: на первом этапе использовано модельное гравитационное поле в без ошибок наблюдений, на втором - поле , измененное с учетом ошибок. Искомыми параметрами (или оценками) интерпретации при этом считались избыточные плотности Док прямоугольных призм, т.е. вектор, содержащий 20 значений избыточной плотности.
Данная задача сводится к решению следующих систем линейных уравнений: - в случае применения метода наименьших квадратов:
01 = а1 1 • Аст1 + а1 2 • Аст2 + • • • + а1 к + 81
02 = й2 .1 • Ао-! + аг 2 •Ас 2 + • • • + а г к • Ас к + 82
Оп = ап1 • Ао1 + ап. 2 • Ао2 + • • • + апк • Аок + 8;
- в случае - сингулярного разложения матриц:
01 = а1 1 • А01 + а1 2 • Ао2 + • • • + а1 к • Аок G2 = а21 • Ао1 + а2 2 • Ао2 + • • • + а2 к • Аок
6п = ап. 1 'А°1 + ап2 • Ао2 + ••• + апк • Аок ,
где Gi - значение гравитационного поля в /-ом пикете, равное сумме всех вкладов
от каждого / тела; эц - коэффициенты матрицы плана А, / = 1...п, п = 80; До? -значение избыточной плотности призм, / = 1...к, к = 20; 6/ - ошибка, связанная с неточностью выбора интерпретационной модели.
Для вычисления коэффициентов э/? матрицы плана А использовалась формула для расчета гравитационного поля от прямоугольной призмы [2], поэтому элементами этой матрицы являются значения базисных функций э/ которые соответствуют эффекту /-ой прямоугольной призмы (/ = 1...20) в /-ой точке профиля (/ = 1.80), причем избыточная плотность такой призмы принята равной 1.
Как известно, решение обратной задачи методом наименьших квадратов в матричной форме записывается в следующем виде [4, 11]:
До0 = (АгА)-1АгО, а для сингулярного разложения матриц эту запись можно представить в виде [9, 10]
Дор = АРв, где АР - псевдообратная матрица, формирование которой будет рассмотрено отдельно.
Известия Сибирского отделения Секции наук о Земле РАЕН. ISSN Геология, разведка и разработка месторождений полезных ископаемых Т. 40, № 2
2541-9455 Proceedings of the Siberian Department of the Section of Earth Sciences RANS.
Geology, Exploration and Development of Mineral Deposits Vol. 40, No. 2
Результаты исследования
Решение обратной задачи для гравитационного поля без учета ошибок наблюдений показало, что сумма квадратов невязки между исходными (модельными) значениями избыточной плотности призм А6к, и расчетными А60к и Аорк, определенная по формулам (1), (2), очень мала и для метода наименьших квадратов и для метода сингулярного разложения соответственно равна 1,9110-7 и 1,86■ 10-7 г/см3:
3 = ) - к У2/20 = (1)
= 1,9110 -7 2 / см 3,
(( А^к ) ~АаРк >у20 =
= 1,86-10 ~7 г / см 3.
(2)
Из этого следует вывод о том, что решение обратной задачи в случае использования при интерпретации неискаженного геофизического поля (модельные значения) позволяет получать приемлемые результаты решения.
Для того чтобы решение обратной задачи приблизить к более реальной ситуации и подчеркнуть ее неустойчивость к ошибкам измерений, выбора параметров петрофизической модели и вычислительных процедур, в исходное поле G была внесена помеха, равная ±(3-5)%. Поле с помехой, как уже было сказано выше, обозначено как G1.
Решение обратной задачи для поля с помехой методом как наименьших
квадратов, так и сингулярного разложения показало очень плохие результаты, оно «срывается» и становится непригодным для практического использования. Невязка резко возрастает и составляет 6 = 16,214 г/см3. Это объясняется тем, что решаемая система линейных уравнений вырожденная, определитель D = | АТА | = 3,1710-14.
Для восстановления решения при использовании метода наименьших квадратов мы использовали способ регуляризации, предложенный А.Н. Тихоновым [8]
АоRp = (АТА+арВ)-1(Ат^1+арАо1), где G1 - интерпретируемое поле; В -единичная диагональная матрица; ар = уар-1, у < 1, а1 - число, изначально принимающее достаточно большое значение, р = 1...пп, пп = 103...105 - возможное количество решений; Ао1 - вектор предполагаемых значений избыточной плотности .
В ходе вычислений рассчитывается функция цели SАоR как сумма квадратов невязки между исходными значениями избыточной плотности призм и расчетными (рис. 2). Результаты вычислений SАоR показали, что функция цели имеет значительный локальный минимум (SАоR = 0,028 г/см2) и эта область характеризуется стабилизацией решения. Отсюда следует вывод о том, что использование параметра регуляризации позволяет восстановить решение и получать приемлемые результаты.
Рис. 2. Сумма квадратов невязки между исходными и расчетными значениями избыточной плотности призм Fig. 2. Sum of squares of discrepancy between the initial and calculated values of excess prism density
20
Известия Сибирского отделения Секции наук о Земле РАЕН. р. Геология, разведка и разработка месторождений полезных ископаемых Т. 40, № 2 ISSN
Proceedings of the Siberian Department of the Section of Earth Sciences RANS. 2541-9455
Geology, Exploration and Development of Mineral Deposits Vol. 40, No. 2
Применение аппарата сингулярного разложения также позволило получить хорошие результаты. Для этого пришлось выполнить решения с учетом реальной точности исходных материалов (полевых работ), в нашем случае это 3-5%.
Для формирования оптимального варианта решения обратной задачи использовались процедура сингулярного разложения матрицы плана А. При этом в матрицу плана А был добавлен 21 -ый столбец R, все члены которого равны 1. Размерность R - (80*1). Данное действие позволило нам в дальнейшем получить информацию о значении фонового гравитационного поля по данному профилю.
Сингулярным разложением действительной матрицы А называется всякая ее факторизация вида А = V•S•U, где и, V - ортогональные матрицы; S - диагональная матрица, у которой диагональные члены уу называются сингулярными числами матрицы. S = (66,50; 45,19; 33,45; 26,78; 9,81; 7,58; 3,12; 2,46; 1,86; 1,34; 0,66; 0,34; 0,22; 0,16; 0,06; 0,02; 0,01; 0,008; 0,004; 0,001; 0,0007). Тогда число обусловленности матрицы А в = у1/у21 = 66,5/0,0007 = 9,5104 - это достаточно большая величина, что говорит о значительной вырожденности матрицы А и плохих результатах решения системы линейных уравнений обычными методами.
Если задать оптимальную точность решения 3%, то наименьшее значение сингулярного числа т = 0,03 /1 = 66,5 0,03 = 1,995, т.е. это соответствует использованию первых девяти сингулярных чисел, к = 9.
На следующем этапе следует
Библиографический список
1. Вахромеев Г.С., Давыденко А.Ю. Моделирование в разведочной геофизике. М.: Недра, 1987. 192 с.
сформировать три новых матрицы Sk, Uk и Vk, которые будут включать в себя определенные фрагменты соответствующих матриц S, U и V [9, 11].
Далее рассчитывается псевдообратная матрица АР как произведение соответствующих матриц:
AP = VkSkUTk, где UkT - транспонированная ортогональная матрица размерностью (9*80); Vk - ортогональная матрица размерностью (21*9); S - диагональная матрица размерностью (9*9).
Вектор искомых значений избыточной плотности 20 призм и уровень нормального гравитационного поля вычисляется по следующей формуле:
Дор = APG1, где АР - псевдообратная матрица размерностью (21*80), G1 - вектор столбец интерпретируемого гравитационного поля, содержащий помеху.
Результат решения обратной задачи для гравитационного поля с ошибкой наблюдений показал, что сумма квадратов невязки Дор между исходными значениями избыточной плотности призм Абк и расчетными Дорк мала: б = 0,032г/см2. Данный результат вполне приемлем для практических целей.
Заключение
Общий вывод проведенного исследования заключается в том, что метод наименьших квадратов с применением аппарата регуляризациии и метод сингулярного разложения матриц показывают вполне удовлетворительные, схожие результаты решения обратной задачи гра-виразведки для рассматриваемого модельного примера.
References
1. Vahromeev G.S., Davydenko A.Ju. Modelirovanie v razvedochnoj geof-izike [Modeling in geophysical prospecting]. Moscow, Nedra Publ., 1987.192 p.
Известия Сибирского отделения Секции наук о Земле РАЕН. ISSN Геология, разведка и разработка месторождений полезных ископаемых Т. 40, № 2 gg
2541-9455 Proceedings of the Siberian Department of the Section of Earth Sciences RANS.
Geology, Exploration and Development of Mineral Deposits Vol. 40, No. 2
2. Успенский Д.Г. Гравиразведка. Л.: Недра, 1968. 230 с.
3. Логачев А.А., Захаров В.П. Магниторазведка. Л.: Недра, 1979. 360 с.
4. Кудрявцев В.А., Демидович Б.П. Краткий курс высшей математики. М.: Наука, 1978. 620 с.
5. Никитин А.А. Теоретические основы обработки геофизической информации. М.: Недра, 1986. 337 с.
6. Никитин А.А. Статистические методы выделения геофизических аномалий. М.: Недра, 1979. 280 с.
7. Грознова А.А. Математические методы интерпретации магнитных аномалий. М.: Недра, 1985.150 с.
8. Старостенко В.Н. Устойчивые численные методы решения обратных задач. М.: Недра, 1988. 250 с.
9. Форсайт Дж. Машинные методы математических вычислений. М.: Мир, 1980. 280 с.
10. Яновская Т.Б., Прохорова Л.Н. Обратные задачи геофизики. СПб.: Изд-во СПбГУ, 2004. 214 с.
11. Блох Ю.И. Интерпретация гравитационных и магнитных аномалий. М.: Изд-во РГГРУ, 2009. 230 с.
2. Uspenskij D.G. Gravirazvedka [Gravity prospecting]. Leningrad, Nedra Publ., 1968. 230 p.
3. Logachjov A.A. Magnitorazvedka [Magnetic prospecting]. Leningrad, Nedra Publ., 1979. 369 p.
4. Kudrjavcev V.A. Kratkij kurs vyssej matematiki [A short course in higher mathematics]. Moscow, Nauka Publ., 1978. 620 p.
5. Nikitin A.A. Teoreticheskie osnovy obrabotki geofizicheskoj informacii [Theoretical bases of geophysical information processing]. Moscow, Nedra Publ., 1986. 337 p.
6. Nikitin A.A. Statisticheskie metody vydelenija geofizicheskih anomalij [Statistical methods of geophysical anomaly distinguishing]. Moscow, Nedra Publ., 1979. 280 p.
7. Groznova A.A. Matematicheskie metody interpretacii magnitnyh anomalij [Mathematical methods of magnetic anomaly interpretation]. Moscow, Nedra Publ., 1985. 150 p.
8. Starostenko V.N. Ustojchivye chislennye metody reshenija obratnyh zadach [Stable numerical methods of reverse problem solution]. Moscow, Nedra Publ., 1988. 250 p.
9. Forsajt Dzh. Mashinnye metody matematicheskih vychislenij [Computer methods of mathematical calculations]. Moscow, Mir Publ., 1980. 280 p.
10. Jasenovskaja T.B., Prohorova L.N. Obratnye zadachi geofiziki [Reverse problems of geophysics]. Saint Petersburg, Sankt-Peterburgskii gosudarstven-nyi universitet Publ., 2004. 214 p.
11. Bloh Ju.I. Interpretacija gravi-tacionnyh i magnitnyh anomalij [Interpretation of gravity and magnetic anomalies]. Moscow, Rossiiskii gosudarstvennyi ge-ologorazvedochnyi universitet imeni Sergo Ordzhonikidze Publ., 2009. 230 p.
Статья поступила 30.05.2017 г.
The article was received 30.05.2017.
Известия Сибирского отделения Секции наук о Земле РАЕН. Геология, разведка и разработка месторождений полезных ископаемых Т. 40, № 2 ISSN Proceedings of the Siberian Department of the Section of Earth Sciences RANS. 2541-9455
Geology, Exploration and Development of Mineral Deposits Vol. 40, No. 2
