CC BY
25
-------------------- © И.Ю. Рассказов, В.И Мирошников, В.А. Луговой,
Ю.И Болотин, А.В. Гладырь, А.Ю. Искра, 2011
УДК 622.831.32:534.213
И.Ю. Рассказов, В.И. Мирошников, В.А. Луговой,
Ю.И. БолотинА.В. Гладырь, А.Ю. Искра
ОПРЕДЕЛЕНИЕ ЭНЕРГИИ ВОЛН ГЕОДИНАМИЧЕСКОГО СОБЫТИЯ АКСЕЛЕРОМЕТРИЧЕСКИМ ДАТЧИКОМ
Приведена математическая модель и алгоритм преобразования электрического сигнала акселерометрического датчика в плотность потока энергии от точечного источника геодинамического явления. Алгоритм включат в себя определение спектра сигнала ускорения с помощью Фурье преобразования с учетом неравномерной частотной и направленной чувствительности датчика, преобразование его в спектр скоростей, восстановление временной зависимости скорости с помощью обратного преобразования Фурье, численном интегрировании последней зависимости.
Ключевые слова: удароопасные условия, горные работы, рудные месторождения, безопасность горных работ.
ш ж ри изучении геодинамических процессов явлений важную роль играет точ-
-1А. ность количественной оценки параметров регистрируемых геодинамических событий, в первую очередь их энергии, которая необходима для определения масштабов образования дефектов и разрушений массивов горных пород. Согласно приведенным в работе [1] обобщениям существующие представления о методах расчета энергии землетрясений не вполне устойчивы и иногда дают основания сомневаться даже в порядке получаемых оценок. Учитывая сложность геосреды, проблема повышения точности определения энергии геодинамических явлений до сих пор сохраняет свою актуальность.
Задача определения энергии геодинамического явления имеет две составляющие:
1) корректное преобразование данных измерений сейсмического или акустического сигнала в энергетическую величину; 2) оценку выделения энергии в очаге геодинамического явления, достаточно полно изложенное, например, в [2]. Рассмотрим первую составляющую данной задачи.
Акселерометрические датчики, применяемые в автоматизированных системах контроля горного давления и иных измерительных комплексах, измеряют ускорение в определенном частотном диапазоне и имеют неравномерно направленную чувствительность. По характеристикам сигнала требуется определить плотность потока упругой энергии (количество энергии на единицу поверхности как результат единичного акустического события в определенном промежутке времени). В силу неравномерно направленной чувствительности датчика необходимо знать направление движения фронта волны. Этому должно предшествовать локация координат источника сейсмоакустического события. Поскольку характеристика направленности датчика (зависимость чувствительности от направления) неравномерная, то в случае интерференции волн от разных источников нет возможности отделить волны разных направлений. Так же в случае прохождения прямых и отраженных волн через одну точку их энергии учитываются дважды, как сумма. В работе [3] в расчете энергии
колебаний предлагается ограничиться в основном первым периодом колебаний с добавлением части энергии последующих колебаний. В качестве обоснования данного подхода приводится аналогия с колебаниями маятника, когда его энергия с каждым периодом не возрастает (а наоборот, затухает из-за трения). Аналогичная картина может возникнуть из-за образования стоячих волн в результате пе-реотражения упругих колебаний в ограниченной области массива горных пород или впадения в резонанс самого датчика. Такие особенности процесса учесть нам не представляется возможным. Мы считаем, что одиночный импульс в виде со-литона при распространении на большие расстояния распадается на несколько периодов с сохранением суммарной величины энергии, поэтому в расчете энергии необходимо учитывать все периоды колебаний. Вместе с тем не исключено, что прием датчиком волн от других источников может привести к некоторому завышению энергии исследуемого сигнала.
Рассмотрим процесс распространения от точечного источника с координатами
К = [г\ ГЛГ3 ] упругой сферической волны, в общем случае неравномерной по направлениям. Волна проходит через контрольную поверхность, в которой установлен датчик с координатами гп = [г1,г2,гП] , где п - идентификатор датчика. Радиус контрольной сферы с центром в источнике R = \гп - г |. Вектор направленности фронта потока энергии п/ = (гп -гг)/R . Датчик колеблется в унисон с колебаниями среды и воспринимает эти колебания как плоский фронт волны.
Энергия (поток энергии) Еу переносимая упругой волной равна кинетической энергии в области объемом V , превышающем объем возбужденной части среды:
.р(¥-V)
ЕV«) = 1“^^^ , (1)
где V(г,t) - вектор скорости, являющийся функцией координат и времени; р -плотность среды. С течением времени происходит поглощение кинетической энергии, поэтому функция Еу (Г) - убывающая. С другой стороны за промежуток времени At волна проходит через контрольную поверхность (сфера радиуса R) по нормали (по радиусу):
rр\v\2 1Т) р\v\2 1у) р VI2 _
FCpdt = | FCpdt = FCp | ^2- А = Е, (R), (2)
At 2 t1(R) 2 ) 2
где Л^) - время, когда вся волна находится внутри сферы радиуса R; t2(R) - время, когда вся волна вышла из сферы; Ср - скорость продольных волн (не имея возможности разделить продольные и поперечные волны мы считаем, что весь поток упругой энергии переносится продольной волной. Это допущение вносит погрешность при определении направленной чувствительности датчика); F - площадь контрольной поверхности.
Плотность потока энергии в точке установки датчика при неравномерном распределении потока энергии по направлениям определим из выражения как функцию координат:
В случае измерения скорости (велосиметром) плотность потока энергии определяется по формуле (3) численным интегрированием. При измерениях ускорения требуется преобразование.
Гармонический Спектральный анализ
Упругие механические колебания можно характеризовать величиной перемещения, скорости или ускорения. Так, как нет возможности отделить волны разных направлений, дальнейшее изложение приводим в скалярных величинах. При гармонических колебаниях эти характеристики запишем как:
W(t) = -Y Mrn2 sin(ot + ф) = -YSa2 sin(ot) - YCrn2 cos(ot) =
2 2 (6) = -YS (2 л f) sin (2 л ft)- YC (2лf) cos (2л ft) = WS sin (2л ft) + WC cos (2л ft);
где t - время; ф и со - сдвиг по фазе и угловая скорость; f - частота гармоники; YM ,YS ,YC ,VM ,VS ,VC ,WM ,WS ,WC - коэффициенты разложения (амплитуды гармоник) для перемещения YM, YS, YC ; скорости VM ,VS VC и ускорения WM ,WS ,WC ; соответственно при синусе, при косинусе (могут иметь разные знаки) и по модулю -Ym=VYT+Yt и т.п.; YS =YM cos(^); YC =YM sin(^). Если заданы коэффициенты для
ускорения, то разложение для скорости определяется через эти коэффициенты и имеет вид:
V(t) = VC cos (2л ft) + VS sin (2л ft) + const =
W W W W (7)
cos(ot) +—- sin(ot) + const =-----— cos(2л ft) +--— sin(2 л ft) + const,
со со 2л f 2л f
В этом случае выражение (7) определено с точностью до константы, что делает, в какой-то степени, неопределенной задачу оценку энергии. В выражениях (4-6) постоянная составляющая отсутствует, однако в реальных сигналах она присутствует в силу произвольности в выборе анализируемого интервала времени.
Использование комплексных переменных (в анализе рядов Фурье) Гармонические колебания легко описываются комплексными функциями. Экспонента комплексного числа есть комплексное число
exp(a + i ф) = AM (cos(<p) + i sin(<p)) = Ar + iAi , (8)
где AM = exp(a) - модуль комплексного числа; exp(i ф) = cos(<p) + i sin(<p) - описывает повороты (вращения) на комплексной плоскости; i - мнимая единица; exp(i 0) = 1, Re(exp(i ф)) = cos^), Im(exp(i ф))=+ sin(ф). Для удобства дальнейшего анализа представим разложение вида (5) в комплексной форме.
Y (t) = YM sin(ot + ф) = YS sin(ot) + YC cos(at) = YS sin (2л ft) + YC cos (2л ft); (4)
V(t) = YM a cos(at + ф) = YSa cos(at) - YCa sin(ot) = VC cos (2л ft) + VS sin (2л ft); (5)
Рассмотрим комплексную функцию
Vz (t) = V(ф)ехр( i at) = (VRe + i VIm )(cos(ot) + i sin(ot))
= (VRe cos(ot) - Vm sin(ot)) + i (VRe sin(at) + Vta cos(ot))
где Vz - комплексная амплитуда гармоники; VRe = VC, Vta +VS;
VM = VVS2 +VC2 =-yJVZVZ - модуль амплитуды.
Действительная часть напоминает разложение вида (5) и (7), но имеет отличие в знаке. При замене знака в (9) получим:
V (t) = VzO)eXP( - iat) = (VRe + i VIm )(c0s(-at) + i sin(-fflf))= (10)
= (VRe + i Vm ) (c0s(®t) - i sin(fflt)) = (VRe C0s(®t) + Vta sin(®t)) + i (-VRe sin(®t) + Vta C0s(®t)).
Эта функция описывает вращение вектора длиной VM в комплексной плоскости с угловой скоростью а начиная от положения, сдвинутого на угол q>. Для описания колебаний используется только действительная часть функции V (t) = Re (Vz (t)). Выражение (10)
используется в обратном дискретном преобразовании Фурье.
При возбуждении датчика механическим воздействием вырабатывается электрический сигнал, нелинейно зависящий от параметров волны в точке от соотношения направлений нормали фронта и оси датчика, электромеханических свойств преобразователя:
U = CCn (nd • nf )Cf (f )X , (11)
где X - измеряемая величина, X = {V, W} - скорость или ускорение; Cem - коэффициент электромеханического преобразователя (осевая характеристика чувствительности преобразователя по паспорту, В-с2/м); Cn (nd • nf) - поправочный коэффициент направленной чувствительности (при совпадении направления оси датчика (нормаль nd) и направления распространения фронта волны (нормаль nf), когда модуль скалярного произведения векторов нормалей равен единице, в других случаях меньше единицы,
nd • nf = cos(nd nf)); Cf (f) - поправочный коэффициент нелинейности частотной (f)
характеристики преобразователя (на участке стабильности равен единице, за пределами его - отличен от единицы, применяется при измерениях, выходящих за пределы паспортного частотного диапазона или с приборами с заведомо нелинейной характеристикой).
Таким образом, по измеренному электрическому сигналу находится определяемая величина:
X = Ud/Co, Co = Co (nd, nf, f) (12)
Разложение сигнала в ряд Фурье
Периодический сигнал (колебания скорости) произвольной формы можно представить в виде ряда Фурье
V (t) = XVC,h C0S (2nhfot)+ VSh sin (2”hfot) = Re f'ZVZ,h eXP O' 2nhfot)"] , (13)
h V h )
где fo - частота первой гармоники fo • tm = 1; tm - период или время (интервал) сигнала; h - номер гармоники; Vz h - комплексная амплитуда гармоники с номером
h . Число гармоник должно быть достаточным для представления сигнала с требуемой точностью.
Аналогично для ускорения
W(t) = X Wch cos (2nfht) + wsh sin (T-nfht) = Re (xwz,h exP 0 2nfht)] (14)
h V h )
Между амплитудами гармоник в выражениях (13) и (14) имеются соотношения
Wz,h = i2nfhVZM, (15)
—1
Vz,h = ^~т Wz,h , (16)
2^fh
которые позволяют делать преобразования спектральных характеристик сигналов скорости в ускорения и наоборот.
Дискретное преобразование Фурье
Оцифрованный сигнал датчика представляет собой набор цифр в квантах АЦП. Будем считать коэффициент перевода ускорения в величину электрического напряжения известной. Оцифрованный с шагом St сигнал функции представляется набором чисел
Uj = U (tj) , (17)
где tj = j •St; j = 0..Kj — 1, Kj - количество точек. Эти значения преобразуются в
ускорения согласно выражению (12):
Wj = Uj/Co . (18)
Использование допущения Co = const равносильно неучету значений ускорения за пределами или на краях частотного диапазона чувствительности датчика. Процедура дискретного преобразования Фурье позволяет получить соответственное количество точек комплексных коэффициентов разложения в ряд Фурье. Формула дискретного преобразования:
Wz,h = Wch + iWSh = Kj TWj exp [i j . (19)
Аналогично для скорости
Vz,h = Vc,h + iVSh = К IV exp [i j . (20)
где h = 0..Kh — 1; Kh - количество гармоник.
Обратное преобразование осуществляется аналогичным образом: по известным значениям комплексных коэффициентов разложения Фурье Vzh вычисляются значения в узлах.
Wj = Re
SWz,h exP[ , (21)
Kj
J
Vj = Re
SVz,h exP [ ]]. (22)
h z,h 4 Kj
Причем, значения моментов времени как при прямом, так и при обратном преобразовании не участвуют. Для основной (первой) гармоники ^т = 1, где tm - период основной гармоники. Поэтому роль времени играет соотношение
1 t.
1 = -^. (23)
К t
■) т
Существуют быстрые алгоритмы для дискретного преобразования Фурье, но за скорость приходится платить ограничениями: К] = 2п, где п - целое. Соответственно
^ = 2"-1. При этом первый элемент всегда действительное число.
Определение плотности потока энергии
Численное интегрирование сигналограммы приводит к выражению
Алгоритм определения плотности потока энергии следующий
- Исходные данные: набор точек и = и), где tJ = ]-8^ ] = 0..К] -1, К] -
количество точек.
- Определяем ускорения по (18) Wj = и]/Со.
Предлагаемая методика позволяет более точно определять плотность потока энергии упругих волн, а при достоверных данных о поглощающих свойствах среды спрогнозировать выделение энергии из источника сейсмоакустического события.
1. Ризниченко Ю.В. Проблема величины землетрясения. //Сб. Магнитуда и энергетическая классификация землетрясений. М.: ИФЗ АН СССР, 1974. - Т. 1. С. 43-78.
2. Коган С.Я. Сейсмическая энергия и методы ее определения. М. Наука, 1975. - 152 с.
3. Сторчеус А.В. Заметки к методике расчета сейсмической энергии взрывов и землетрясений. // Материалы конференции, посвященной Дню вулканолога, 27-29 марта 2008 г. - Петропавловск-Камчатский: ИВиС ДВО РАН, 2008. С. 274-281.
КОРОТКО ОБ АВТОРАХ ---------------------------------------------------------------------------
Рассказов И.Ю. - доктор технических наук, директор, ИГД ДВО РАН,
Мирошников В.И. - кандидат технических наук, старший научный сотрудник,
Луговой В.А. - доктор физико-математических наук, гл. научный сотрудник,
Болотин Ю.И. - доктор физико-математических наук, вед. научный сотрудник,
Гладырь А.В. - научный сотрудник,
Искра А.Ю. - старший научный сотрудник,
Института горного дела ДВО РАН, г. Хабаровск, [email protected]
(24)
Определяем комплексные амплитуды ускорения прямого преобразования
Определяем комплексные амплитуды скорости прямого преобразования
Определяем значения скорости обратным преобразованием Фурье по (22)
Определяем значение плотности потока энергии по (24)
СПИСОК ЛИТЕРА ТУРЫ
