CC BY
60
Вестник Кузбасского государственного технического университета. 2015. № 6 97
ЭЛЕКТРОТЕХНИКА
УДК 622.822.2: 622.271: 519.6
ВОССТАНОВЛЕНИЕ ФОРМЫ ИСТОЧНИКА ЭЛЕКТРИЧЕСКОГО ПОЛЯ КАК РЕШЕНИЕ ОБРАТНОЙ ЗАДАЧИ ГЕОЭЛЕКТРИКИ
Сирота Дмитрий Юрьевич,
кандидат техн. наук, доцент, e-mail: [email protected]
Кузбасский государственный технический университет имени Т.Ф. Горбачева, 650000, Россия, г. Кемерово, ул. Весенняя, 28
Аннотация
Актуальность работы. Производство горных работ проходит обычно в сложных горно-геологических условиях нарушенности вмещающих пород. Дистанционное (с дневной поверхности), предварительное определение местоположения таких нарушенностей (сместителей) является весьма важной задачей для обеспечения безопасности ведения горных работ. Кроме того, при производстве подземных горных работ появляются области перепадов механических напряжений, которые являются потенциальными очагами горно-тектонических ударов. Определение местоположения этих областей также ведёт к повышению безопасности ведения горных работ.
Цель работы. Разработка численного алгоритма и его программной реализации для решения обратной задачи воостановления формы источника электрического поля по измеренным значениям потенциала на дневной поверхности.
Методы исследования. Методы математической физики - для формулировки задачи восстановления формы источника, как обратной, некорректной задачи решения интегрального уравнения 1-го рода.
Методы математического программирования - для формирования численной расчётной схемы приближённого последовательного восстановления формы источника поля.
Результаты. Произведён расчёт величины потенциала электрического поля, порождённого тестовой областью эллиптической формы, расположенного в трёхслойной среде. Используя полученные значения, решена обратная некорректная задачи восстановления формы плоского источника электрического поля, как решение интегрального уравнения Фредгольма - Урысона Нго рода. Проведён сравнительный анализ формы и результатов расчёта прямой задачи от точно заданной и приближённо найденной области.
Ключевые слова. Нелинейная некорректная обратная задача, интегральное уравнение Фредгольма-Урысона 1-го рода, методы регуляризации.
Введение. Одной из многочисленных задач интерпретации геофизических данных является задача определения формы и глубины залегания источника поля того или иного вида (электрического, гравитационного, теплового, магнитного или другого). В дальнейшем мы будем рассматривать источник поля электрической природы, порождённое перепадами механических напряжений [1]. Особенностью данного вида источников является то, что они не имеют чётко очерченных границ в отличие от, например, рудных включений, которые порождают магнитные аномалии. Отсюда следует отсутствие необходимости в предельно точном (до сантиметров) определении размеров источника поля.
Решению задач интерпретации геофизических данных посвящено достаточно много работ, в которых рассматриваются как натурные, экспериментальные методы интерпретации [2-4], так и теоретические, в которых рассматриваются численные, компьютерные методы [4-9].
В представленной статье рассматривается реализация численного подхода к вопросу интерпретации геофизических измерений, основанного на решении интегральных уравнений Фредгольма -Урысона 1-го рода. Уравнения такого сорта относится к классу некорректных задач, решению которых посвящены работы [10 - 12]. Отдельно следует упомянуть работы [13 - 15], так как в них рассматривается по сути та же задача, но без учёта геологического строения вмещающих пород.
1. Постановка и решения прямой задачи. При интерпретации геофизических данных методами обратной задачи придерживаются следующего общего подхода:
1) Формируют оператор прямого моделирования - А, который позволяет получить модельные геофизические данные для тестовых параметров геологической или геофизической модели.
В нашем случае в качестве оператора А будем рассматривать функционал расчёта потенциала естественного электрического поля (ЕЭП), порож-
98
Д. Ю. Сирота
денного плоским источником тока произвольной формы, расположенного на известной глубине в трёхслойном изотропном и однородном в пределах каждого слоя пространстве.
2) На базе полученного оператора А, решают обратную задачу восстановления параметров модели по имеющимся экспериментальным геофизическим данным.
В нашем случае в качестве параметров модели будет рассматриваться геометрическая форма источника ЕЭП.
Рис 1. Схема вмещающего пространства и плоского источника поля
Перейдём к формированию оператора прямого моделирования. Пусть ЕЭП порождается плоским горизонтальным источником произвольной формы, который расположен на некоторой глубине в нижнем слое трёхслойного однородного и изотропного вмещающего пространства (рис. 1).
В этом случае потенциал точечного источника будет определяться известной формулой [4]
U"=krJ0(m-r)eM-m-zM)dnu (1)
* 1 + К ехр(-т • И)
где к - коэффициент, характеризующий силу электрического тока источника; -
функция Бесселя нулевого
порядка;
г ~ т]хм + Ум ~ расстояние от точечного ис-
точника до точки измерения поля, м; h - мощ-
ность второго слоя, м; К = ——— - коэффици-
Pi+Pi
ент отражения второго слоя; р12~ удельное электросопротивление первого и второго слоёв (Омм).
После интегрирования (1) по произвольной плоской области Sp получим формулу для потенциала простого слоя
U м = J'uixM ~хр->Ум ~ Уp->zм >C)dSр ,(2)
где вместо г используется величина D = ^xM-xl,)2 + (yM-yl,)2 - расстояние от
точки области интегрирования Р до точки измерения поля М, м; С - параметры, которые характе-
ризуют вмещающее пространство.
Следуя [10, 16], будем, ограничиваясь только звёздным контуром, определять границу области в виде функции р{(р) в полярной системе координат. Переход к последней производится стандартным способом: xp=rcos(p, у р — г sin (р, где
(р е [0;2тс\ и г е [0; /?(^>)]. Далее представим выражение (2) в виде произведения размерного и безразмерного сомножителя, сделав две замены
переменных: во внешнем интеграле = w,
где под w подразумеваются координаты точки измерения поля хм, ум, zM\ расстояние г; граница области интегрирования р((р), и, как следствие,
m-h = m - во внутреннем. Тогда интеграл (2) примет вид
См =
00
= h-kjr J
J0(m-D)exp(-m-zM) d_
drd(p
(3)
1 + К exp(-m)
где D = ,J(xM -rcos<p)2 + (yM -r sirup)2 .
Для расчёта внутреннего несобственного интеграла произведём аппроксимацию выражения вида
1
— = Z *?('к) ехР[“ m (к ~ Н >
1 + Кехр(-т) ы
где коэффициенты q(k) определим методом наименьших квадратов.
После применения интеграла Вебера - Лиф-шица получим выражение вида
2л-Г и
о
р{ср)
uM=k-h\ ) 5=
d(p, (4)
где R = г2+ В-Г + А,
А = Хм+Ум + (zM + к-1)2,
В = -2{хм cos(р-ум sinср).
Здесь и далее будем считать, что известна глубина залегания источника ЕЭП ZM = Н и используются только данные измерений вдоль оси Ох, то есть ум= 0 .
Внутренний интеграл в (4) лёгко вычисляется:
р(<р)
rdr
\ PR
2 2 J
Р(<Р)
= F(xu,p)
Таким образом, для решения задачи определения формы источника поля по измерениям потенциала на дневной поверхности получаем интегральные уравнения (ИУ) Фредгольма - Урысона 1-го рода относительно функции р{ср) :
2 л
\xU{xM,yM,zM,p(cp))d<p = U', (5)
О
Вестник Кузбасского государственного технического университета. 2015. № 6 99
п
где KU = ^q(k)F(xM,p((p)) - ядро нелиней-*=]
ного вида; U* - измеренные на дневной поверхности значения потенциала ЕЭП.
2. Решение обратной задачи.
В серии работ [5, 13, 14], уравнения (5) после дискретизации рассматривают как систему нелинейных уравнений и решают её методом Левен-берга-Марквардта. В данной статье будем следовать базовым работам [11, 12, 16], посвящённым исследованию некорректных задач, где решение (5) ищется как результат минимизации сглаживающего регуляризованного функционала невязки
Ф (p)==\[AxM’P)-U*{xMtf dxM +>(6)
С
+ а • Cl(p) —» min
2к
где Л(хм,р)= | KU(xM ,p)d(p - нелинейный
о
интегральный оператор прямого моделирования
2 л
(4); Q(p) = J*p2d(p - стабилизирующий функ-
0
ционал; KU - ядро оператора, U*(xM} - правая часть ИУ, которая на практике задаётся экспериментально, а для тестовых задач определяется путём решения прямой задачи по формулам (4); а -параметр регуляризации, который определяется по формуле геометрической прогрессии [4, 6]
«*=«„■/“Л (7)
где q - номер итерации, / 6 [О; l] - знаменатель
прогрессии, а0 - стартовое число, которое определяется из соотношения (6) [4, 6]
и
\[А(хм,р0)~ f{xM)}2 dxM
ап = ■
Фо)
(8)
Минимум функционала (8) будем искать методом сопряжённых градиентов [], общая итерационная схема которого имеет вид
р(ч+]) =p^+k{q)-I{q), (9)
где = —1q - номер итера-
ции; к - шаг минимизации; S - линейная комбинация производных на предыдущих шагах; коэффициент определяется по формуле Поллака -
Рис 2. Пошаговые приближения от эллипса с полуосями а = 400, b = 250 м, I - № шага.
Рис. 3. Пошаговые приближения от круга радиуса а = 400,1-№ шага.
Рис. 4. Пошаговые приближения от эллипса с полуосями а = 100, b = 65 м, I -№ шага.
Рис. 5. Пошаговые приближения от круга радиуса а = 100,1 - № шага.
100
Д. Ю. Сирота
Рис. 6. Пошаговые приближения от эллипса с полуосями а = 100, b = 65м, I-№ шага.
Рис. 7. Пошаговые приближения от круга радиуса а = 100,1-№ шага
Рис. 8. Пошаговые приближения от эллипса с полуосями а = 10, b = 6 м, I - № шага.
Ребьера [17], Фд - производная функционала (6),
которая определяется по формулам [11]
2 п
4>/p=ap(t)+ \T[(p,t,p((p),p(t)\l<p- >(10)
о
~F[t.p(t)\
d
где T = J К[хм, р(ф)\ К![хм, p{t)\BM ,
С
d
Р = \к[хм>Р({)У\хм)с1хм, 0<t<2x -
С
фиксированный угол для переменной р{(р), по
которой вычисляется производная Фд .
Производные от ядра KU определяется по формуле
KUL = Yq(k)
k=1
Для вычисления интегралов в (6) используется метод квадратур, что преобразует интегральное представление (6) к дискретному на сетке по пе-
ременным хм , (р, t. Количество узлов в сетках по
переменным (р и t должно совпадать; шаг сетки
Рис. 9. Пошаговые приближения от круга радиуса а = 100,1-№ шага.
выберем равным постоянному значению h{ = ,
п
при этом стартовое значение (pQ = tQ = 0; узлы сетки по переменной X м будут также идти с по-. d-с
стоянным шагом п2 =-----. Применяя квадра-
т
турную формулу Симпсона к каждому интегралу из выражений (6), получим дискретную задачу минимизации
фа(А)->т1П 02)
для переменных рк, к — 1,..., п .
3. Реализация численного алгоритма.
Прежде всего, проверим точность вычисления интегралов выражения (10) по формуле Симпсона. Рассмотрим тестовую область в виде эллипса, безразмерное уравнение которого в полярной системе координат имеет вид
_, ч ab
Р(<Р) = —П-----;-----;---— • 03)
hyja sin <p + b cos (p
Рассмотрим два вида эллипсов: 1) <2=300 , Ь=96 м и 2) <2=50 , Ь=\5 м; размерный коэффициент в обоих случаях равен h=25 м. Расчёты показали, что точность определения производной
Вестник Кузбасского государственного технического университета. 2015. № 6 101
Рис. 10. Пошаговый расчёт потенциала от круга радиуса а = 400 м
для 1-й тестовой области; 1-№ шага (10) в первом случае составляет величину
А = 10-6, а во втором - А = 10-3.
Далее рассмотрим реализацию итерационной схемы (9) на тестовых областях.
В качестве первой тестовой области возьмём эллипс с полуосями а = 300, b = 96 м, положив в качестве стартовых 4 области: внешний эллипс с полуосями а = 400, b = 250 м (рис. 2); внешний круг радиуса а = 400 м (рис. 3); внутренний эллипс с полуосями а = 100, b = 65 м (рис. 4); внутренний круг радиуса а = 100 м (рис. 5).
В качестве второй тестовой области возьмём эллипс с полуосями а = 50, b = 16 м, положив в качестве стартовых 4 области: внешний эллипс с полуосями а = 100, Ъ = 65 м (рис. 6); внешний круг радиуса а = 100 м (рис. 7); внутренний эллипс с полуосями а = 10, b = 6 м (рис. 8); внутренний круг радиуса а = 10 м (рис. 9).
На основании рис. 2-9 можно заключить, что минимизация одномерного интегрального
Рис. 11. Пошаговый расчёт потенциала от эллипса с полуосями а = 10, b = 6 м для 2-й тестовой области; 1-№ шага
функционала (6) приводит, вообще говоря, к ложному «квазирешению» и определяется в первую очередь стартовой геометрией области: эллипс переходит в эллипс, окружность переходит в окружность. Это можно объяснить фактом наличия эквипотенциальных поверхностей: для одного и того же распределения величины потенциала можно подобрать как минимум два разных по форме источника поля. Однако, несмотря на несовпадение формы, вычисленные по формуле (4) приближённые и точные значения потенциала на дневной поверхности будет весьма близки друг другу (рис. 10, 11).
Можно сделать предположение о том, что если в (6) рассматривать не одномерный по переменной Хм , а двумерный функционал по плоской
области Хмум, то это рассмотрение приведёт к более точному воспроизведению двумерной формы тестовой области.
СПИСОК ЛИТЕРАТУРЫ
1. Тарасов, Б.Г. Физический контроль массивов горных пород / Б. Г. Тарасов, В. В. Дырдин, В. В. Иванов, А. Н. Фокин // М.: Недра, 1994, - 240 с.
2. Дашевский, Ю. А. Обратные задачи электрических зондирований в сейсмоактивных зонах / Ю.А. Дашевский, А. А. Мартынов // - Новосибирск: НГУ, - 2002 год, - 52 с.
3. Дашевский, Ю. А. Прямые и обратные задачи геоэлектрики в неразрушающих методах контроля свайных фундаментов / Ю. А. Дашевский, И. В. Суродина, О. Ю. Дашевский, В. П. Соколов // Сиб. журн. индустр. матем., 8.2 -2005 г., - с. 57-69.
4. Жданов, М.С. Геофизическая электромагнитная теория и методы / М. С. Жданов // М.: Научный мир, 2012, - 680 с.
5. Трушникова, Н. В. Об определении формы объёмного теплового источника в массиве горных пород по измерениям электрического потенциала на поверхности/ Н. В. Трушникова // Геодинамика и напряжённое состояние недр Земли. Труды международной конференции. - Новосибирск, изд-во ИГД СО РАН, 2001, с 165 -167
6. Жданов, М. С. Теория обратных задач и регуляризация в геофизике/ М. С. Жданов //-М.: Научный мир, -2007 710 с.
7. Preto, German A. Preto Geophysical inverse theory Universidad de los Andres / German A. Preto // El-siver -2009, - 48 p.
8. Tarantola, A. Inverse problem theory and methods for model parameter estimation / A. Tarantola // - SIAM, -2005, - 358 p
102
Д. Ю. Сирота
9. Кабанихин, С.И. Обратные и некорректные задачи/ С.И. Кабанихин// - Сибирское научное издательство, Новосибирск, -2009, 458 с.
10. Бакушинский, А.Б. Итеративные методы решения некорректных задач / А. Б. Бакушинский, А. В. Гончарский // -М: Наука, 1989, 127 с.
11. Верлань, А. Ф. Интегральные уравнения: методы, алгоритмы, программы / А. Ф. Верлань, В. С. Сизиков / /-Киев: Наукова думка, -1986, 544 с.
12. Васильева, А. Б. Интегральные уравнения / А. Б. Васильева, Н. А. Тихонов // -М.: ФИЗМАТЛИТ, 2002.- 160 с.
13. Трушникова, Н. В. Обратная задача локации очага самовозгорания угля по измерениям потенциала термоэлектрического поля на поверхности угольного разреза/ Н. В. Трушникова //Вестник КузГТУ, 2010, №3 .-с. 3-6.
14. Трушникова, Н. В. Определение некоторых параметров объёмных тепловых источников в массиве горных пород по распределению градиента электрического потенциала на поверхности / И. В. Трушникова // Вестник КузГТУ, 2006, № 6. -С. 3-7.
15. Kress, R Nonlinear integral equations and the iterative solution for an inverse boundary value problem/ R. Kress, W. Rundell// Inverse problems. - 2005. - № 21. - p. 1207 - 1223.
16. Бакушинский А. Б. Некорректные задачи. Численные методы и приложения / А. Б. Бакушинский, А. В. Гончарский //-М.: МГУ, 1989, -199 с.
17. Пантелеев, А. В. Методы оптимизации в примерах и задачах / А. В. Пантелеев, Т. А. Летова // -М.: Высшая школа, -2005 год. - 544 с.
Поступило в редакцию 6.10.2015
UDC 622.822.2: 622.271: 519.6
THE DETERMINATION FORM OF SOURCE OF ELECTRIC FIELD AS THE INVERSE PROBLEM OF GEOELECTRICS
Sirota Dmitriy Yu.,
C. Sc. (Engineering), Associate Professor, e-mail: [email protected]
T.F. Gorbachev Kuzbass State Technical University, 28 street Vesennyaya, Kemerovo, 650000, Russian Federation
Abstract
The urgency of the discussed issue. Mining operations usually takes place in difficult geological conditions of disturbance of the host rock. Remote (from the surface), the preliminary determination of the location of such forced is a very important task to ensure the safety of mining operations. In addition, in the manufacture of underground mining area there are changes of the mechanical stresses which are the potential hotbed of mining and tectonic shocks. The location of these areas also increases the safety of mining operations.
The main aim of study. Development of numerical algorithm and its software implementation to solve the inverse problem of restoration of shape of the source of electric field on the measured value of potential on the surface.
The methods used in the study Methods of mathematical physics is to formulate the problem of recovering the shape of the source as an inverse, ill-posed problems the solution of inte-integral equations of the 1st kind. Methods of mathematical programming - to create a numerical calculation scheme of approximate sequential restoring the shape of the field source.
The results. The calculation of the magnitude of the electric field potential, generate daily test area is elliptical in shape, located in a three-layer fry de. Using the obtained values, we solved the inverse incorrect problems res tore-ing form flat electric fields, solution of integral equations of Fredholm - Uryson 1-St kind. A comparative analysis of the form and calculation results of the direct problem from the exact and approximated found the field.
Keywords. Nonlinear inverse problem, integral equation of Fredholm-Uryson of the first kind, regularization methods.
REFERENCES
1. Tarasov B. G., Dyrdin, V. V., Ivanov V. V., Fokin A. N. Physical controls MASSY Islands of rocks. -M.: Nedra, 1994, 240 p.
Вестник Кузбасского государственного технического университета. 2015. № 6 ЮЗ
2. Dashevsky, Y. A., Martynov A. A. Inverse problem of electrical sounding in seismic zones. - Novosibirsk: NSU, 2002,, 52 p.
3. Yu. A. Dashevskii, V. I. Surodina, O. Yu. Dashevsky , Sokolov V. P., “Direct and inverse problems of geoelectrics in nondestructive inspection of pile funda-ments”, Sib. Phys. industry. Mat., 8.2 -2005, -p. 57-69.
4. Zhdanov M. S. Geophysical electromagnetic theory and methods. - M.: Scientific world, 2012, - 680 p.
5. Trushnikova N. V. definition of volumetric heat source in a rock mass from measurements of the electric potential on the surface //Geodynamics and stress state of the Earth's interior. Proceedings of the international conference. -Novosibirsk, publishing house of the Institute of mining of SB RAS, 2001, - p. 165 - 167.
6. M. S. Zhdanov, Theory of inverse problems and regularization in Geophysics -M.: Scientific world, pp. 710-2007
7. German A. Preto Geophysical inverse theory Universidad de los Andres, -2009, - p 48.
8. Albert Tarantola Inverse problem theory and methods for model parameter estimation, SIAM, -2005, - p 358
9. Kabanikhin S. I. Inverse and ill-posed problems, Siberian science statelist, Novosibirsk, -2009, p. 458
10. Bakushinsky A. B., Goncharsky A.V., Iterative methods for solving ill-posed problems. -M.: Nauka, 1989, p. 127
11. Verlan A. F., Sizikov V. S. Integral equations: methods, algorithms, programme-Kiev: Naukova Dumka, -1986, 544 p.
12. Vasileva A. B., Tikhonov N. A. Integral equations. -M.: FIZMATLIT, 2002. - 160 p.
13. Trushnikova, N. V. The inverse problem of the location of foci of spontaneous combustion of coal by measure-ments of the thermoelectric potential field on the surface coal mine //Vestnik KuzSTU, 2010, № 3 . -p. 3-6.
14. Trushnikova, N. V. Determination of some parameters of volumetric heat sources in the array of rocks on the distribution of the gradient of the electric potential-on the surface of La / N. Trushnikova, V.// Bulletin Of KuzSTU, 2006, № 6. -P. 3-7.
15. Kress, R Nonlinear integral equations and the iterative solution for an inverse boundary value problem/ R. Kress, W. Rundell// Inverse problems. - 2005. - No. 21. - p. 1207 - 1223.
16. Bakushinsky A. B., Goncharsky A. V. Incorrect problems. Numerical methods and applications.-M.: MGU, 1989, p.-199
17. Panteleev A. V., T. A. Letova The optimization methods in examples and problems. - M.: Higher school, 2005. - 544 p.
Received 6 October 2015
