CC BY
15
УДК 621.373.826
В.А. Трофимов, Ю.В. Трощиев2
О НЕКОТОРЫХ СПОСОБАХ ОПРЕДЕЛЕНИЯ СИГНАЛА ТЕРАГЕРЦОВОГО ЛАЗЕРНОГО ИЗЛУЧЕНИЯ И ДИНАМИКИ ЕГО ЧАСТОТ
В работе изучается возможность восстановления сигнала терагерцового диапазона, прошедшего через вещество и представленного в виде набора интегральных измерений. Такие задачи актуальны в связи с проблемами идентификации веществ. В работе предложены два алгоритма восстановления сигнала с помощью сплайна, позволяющие практически полностью уменьшить присутствие ложных частот. Принципиально, что сплайн с числом узлов, примерно равным числу измерений, не гарантирует высокого качества восстановления сигнала. Предложенный алгоритм с выбором оптимального числа измерений свободен от этого недостатка.
Рассмотрение проводилось как для модельного сигнала, так и для реального отклика среды.
Ключевые слова: терагерцовая спектроскопия, динамика спектральных линий.
Введение. Как известно, в настоящее время интенсивно развивается оптика малопериодных импульсов излучения терагерцового диапазона, что обусловлено его широким практическим приложением. Это излучение используется в дистанционном зондировании и спектроскопии химических, биологических, полупроводниковых объектов [1-4]. С его помощью можно проводить безвредную для человека диагностику [5-7], а также управлять химическими реакциями и электронными состояниями в квантовых ямах, разрабатывать системы безопасности и беспроводной связи компьютеров и периферии в помещениях [8, 9], определять наличие наркотиков и взрывчатых веществ [10-13]. Предложенное решение последних проблем в литературе основывалось на анализе спектра Фурье сигнала или отклика среды, полученного от исследуемого вещества. Однако вещества, которые представляют интерес для идентификации, часто имеют спектр, близкий ко многим веществам, используемым в повседневной практике [13]. Поэтому применение известных методов на основе фурье-анализа является малоэффективным в ряде случаев. Возможности идентификации веществ существенно расширяются при анализе динамики спектральных линий отклика вещества. Для ее определения предложены некоторые методы в [14-17], часть из которых апробированы в физических экспериментах [14, 18]. Ранее более простая модификация подобной методики применялась в [19] для динамики спектрального отклика среды при воздействии оптических фемтосекундных импульсов. Отметим, что многие применяемые сейчас в задачах пространственно-временной спектроскопии (time-domain spectroscopy) методы изложены в [20].
В настоящей работе предложены некоторые обобщения развиваемых алгоритмов, в частности сплайн-восстановления [17], обладающие рядом преимуществ по сравнению с вариантом, описанным в указанной работе.
Постановка задачи. Обозначим через f(t) исходный неизвестный сигнал (например, форму светового импульса), непрерывный на отрезке времени [tb,te]. В результате измерений в физическом эксперименте известен набор интегральных данных di, (¿2,..., (¿дг от функции f(t) (в случае лазерного импульса это энергия), определяемый следующим образом:
а для промежутков времени справедливы соотношения ^ — 0,5Тг ^ ¿5, + 0,5Тг ^ ¿е и ^ < ¿¿+1; г = 1,..., N. Причем первое измерение начинается в момент времени т. е. — = а последнее измерение заканчивается в момент времени ¿е, т. е. ¿дг + 0,5Тдг = Здесь N — число временных
хФакультет ВМиК МГУ, д.ф.-м.н., проф., e-maihvatroQcs.msu.su
2Факультет ВМиК МГУ, к.ф.-м.н., с.н.с., e-maihyuvtQcs.msu.su
* Работа выполнена при частичной финансовой поддержке гранта РФФИ 08-07-12007офи.
ti+Q,5Ti
(1)
ti-Q,5Ti
интервалов, — время начала регистрации сигнала, ¿е — время ее окончания, Т.\ — длина окна, т. е. длина интервала, на котором производится измерение. Отметим, что в данном случае термин "окно" не означает использования оконной функции (см., например, [14, 16, 18]). Для удобства определим также величину А^ = — ^ — величина сдвига окна, т.е. промежуток времени, на который сдвигается центр окна ^ в интеграле (1). Задача состоит в восстановлении сигнала /(¿) по набору этих интегральных измерений. Таким образом, эта задача представляет собой обратную задачу.
В качестве объектов, для которых исследуется эффективность используемого метода, рассматривается модельный сигнал
а также реальный сигнал — отклик молекулы Ь-сувЫпе при воздействии терагерцового импульса [14].
Восстановление сигнала без дополнительных условий на границах временного отрезка. В работе [17], в частности описан алгоритм восстановления сигнала по набору интегральных измерений, основанный на сплайн-интерполяции этого сигнала: применялся кубический сплайн. При всех положительных качествах этого алгоритма (высокая точность вычисления интегралов, гладкость сигнала на непрерывном множестве, быстродействие) остались и некоторые вопросы. Например, не исследована чувствительность алгоритма к способу расстановки точек на временном отрезке для сплайн-интерполяции или влияние задания нулевых граничных условий для записи дополнительных уравнений с целью однозначного определения коэффициентов сплайна. В этих работах число отрезков сплайна, восстанавливаемого по N измерениям, было равно N+1, а, узлы интерполяции расставлялись с помощью датчика случайных чисел исходя из разрешимости системы линейных уравнений. Это объясняется тем, что N измерениям соответствуют N+2 точки (середины измерений и границы отрезка). Поэтому узлы интерполяционного многочлена не совпадали с моментами времени, соответствующими измерениям. Учитывая, что при равномерном задании узлов и равномерном расположении измерений сплайн очень сильно отличался от исходного сигнала, а сами измерения перекрывались, система узлов интерполяции отличалась от системы измерений. Ниже выполнено дальнейшее развитие и улучшение алгоритма.
Отметим, что в случае осциллирующего около нулевого значения сигнала задание его нулевых значений (так же как и нулевых значений второй производной) на границах временного интервала не приводит к большой ошибке восстановления его формы. Однако желательно все же задавать более точные значения, например сигнала вблизи границ временного отрезка. Для этого может быть применен следующий подход. Пусть число измерений, как и прежде, равно N. Число отрезков кубического сплайна на восстанавливаемом отрезке времени уменьшим на 4 по сравнению со сплайном, предложенным в [17], т.е. число отрезков сплайна будет равно N — 3. При этом число уравнений, необходимых для построения этого сплайна, уменьшится на 16, так как каждый его элемент описывается четырьмя коэффициентами, а в системе алгебраических уравнений для определения его коэффициентов исчезнут 12 уравнений, записываемых из условий непрерывности сплайна и его производных. Таким образом, для построения такого сплайна требуется на 4 условия меньше, чем в случае сплайна с Ж + 1 отрезком. Поэтому для записи систем уравнений относительно его коэффициентов будет достаточно результатов измерений.
Следующим шагом к построению сплайна является расстановка точек для сплайн-интерполяции, так как от нее зависит качество восстановления. В [17] точки выбирались с использованием датчика случайных чисел. В приводимых ниже расчетах первая и последняя точки сплайна совпадают с границами временного отрезка, а внутренние точки сплайн-интерполяции совпадают с серединами отрезков измерений. Поскольку число внутренних точек сплайн-интерполяции на 4 меньше, чем число измерений, то по две точки с каждой стороны отрезка исключались. Соответственно на отрезке [¿ь,£е] исходный сигнал приближается кубическим сплайном /(£):
где вi — узлы интерполяционной сетки: 0О = 0\ = = 2, г = 2,..., N — 5, = ¿лг-ъ
0дг-з = ¿е- Коэффициенты сплайна Д, 7$, <5$, г = 0...../V — 1. находятся из условия равенства
, ( ¿-3° ! ¿-80
/(¿) = 8т(1,2£) + ) ЯП 1,5* + 2е"1—) 8т(М)
(2)
т = ЩЦ - вг)* + - вг)2 + - вг) + 8и I € [вьв^г)
(3)
интегралов от сплайна и измерений
и+ъ/ 2
¿г = ¿г, ¿г = J /(¿) (Й,
и-Тг/2
% = 1,
(4)
В точках интерполяционной сетки задаются условия непрерывности функции /(£), а также первой и второй производных, которые, как известно, позволяют дополнительно записать уравнения для построения сплайн-сигнала (3):
щАв\ + &Д0? + лАвг + 5г = 5г+ъ ЗщАв* + 2РгА0г +Л= Ъ+Ъ ЪщАвг + & = &+1,
Мг=и-и-Ъ / I.....Л— 1. (5)
Так как для любого временного окна //7; /2, trn + Тгп/2], тп (Е [1, /V^. могут быть найдены значения г и такие, что — наименьший покрывающий его отрезок:
Рг+Ъ — 1] [¿т Тт/2,*т
то значение измеренной интегральной величины с^ может быть представлено в виде
(6)
£т+Тт/2 01+1 j — 2 ^ + 1 ¿т+Тт/2
^т —
/(¿) Л =
№ &
ггп-тгп/2
ггп-тгп/2
Е
к=г+1 0
№ &
т <й =
= (-тш
+ 1ъАв* + 5г Авг
3 2
+ + + бгАвг
I •) 1
3-2
£
к=г+1
-лакАв1 + \ркАв1 + \1кАв2к + 6кАвк
+ ^ <*,-_! Д^ + + Ц-1 +
Ав = ¿то — Тто/2 — Д0 = ¿то + Тто/2 — ^--1. (7)
Следовательно, относительно коэффициентов сплайна (3) получается система линейных уравнений (5), (7). Нетрудно видеть, что число неизвестных коэффициентов 4(Ж^З) равно числу уравнений. Поэтому они могут быть найдены из решения этой системы уравнений с ленточной матрицей, ширина ленты которой зависит от того, сколько отрезков сплайна может быть покрыто окном одного измерения [1т 2,¿то +Тт/2], т.е. от значения тах(^(т) — г(т)) в формуле (6). Отметим, что значения
сплайна в граничных точках определяются из этой системы уравнений, а не задаются искусственно.
Прежде чем переходить к обсуждению точности восстановления формы сигнала, заметим, что время расчетов, проведенных на персональном компьютере с тактовой частотой 2 СНи с одинарной точностью (4 байта на число), не превышало 0,08 с (как правило, оно изменялось в интервале 0,015 ± 0,005). Это очень важно для практики.
Для исследования свойств описанного метода проводилось компьютерное моделирование как для модельного (2), так и для реального сигнала. Точность восстановления оценивалась двумя критериями: относительной ошибкой восстановления в норме 1. <
А Д2 =
/-/с
11/1
(8)
и относительной ошибкой восстановления в норме С
/-/с
А }с =
с
11/1
(9)
с
вычисленными по точкам, в которых задан сигнал, а также учитывалась максимальная разность между измерениями реального и восстановленного сигнала
Ай = тах
i=l,...,N
¿г — 4с,г
(10)
В формулах (8)—(10) /с и (¿с,г — соответственно численно найденный сплайн и интегралы от него (см. (4)). В общем случае А(I не равно нулю, например, из-за ошибок округления. Эта величина характеризует также применяемый алгоритм. В зависимости, например, от расстановки точек для сплайн-интерполяции значение величины Аё может существенно отличаться от нуля или система может оказаться вообще не решаемой. Поэтому малость этой величины (в соответствии с точностью представления чисел) может использоваться как необходимое условие успешного решения системы линейных уравнений при определении коэффициентов сплайна. Важно также подчеркнуть, что для вычисления величин Д/ь2 и А/с нужно знать исходный сигнал (очевидно, на практике это, как правило, не реализуется, особенно для реального сигнала), а для вычисления величины Ай достаточно значений измерений. Необходимо подчеркнуть, что, несмотря на малое значение величины Аё в точках измерений, отклонение формы восстановленного сигнала может быть существенным (рисунок, б). Поэтому предлагаются новые способы восстановления сигнала.
Некоторые результаты компьютерного моделирования приведены на рисунке, а, б — для модельного сигнала при N = 500 и реального сигнала при N = 426 в случае измерительной системы с окном Т = 1 и сдвигом окна А = 0,2. В случае модельного сигнала исходный и восстановленный сигналы практически совпадают, а в случае реального сигнала исходный сигнал в точках несовпадения виден как темная линия внутри осцилляций восстановленного сигнала. Как следует из табл. 1, 2, при увеличении сдвига окна А до 0,4; 0,6; 0,8 для обоих сигналов наблюдается немонотонность точности восстановления от числа измерений (см. варианты 1-4 табл. 1, 2): для А = 0,6 она выше, чем достигаемое ее значение при А = 0,4 и А = 0,8.
Табл и ца 1
Модельный сигнал
Вариант способа восстановления 1 2 3 4 5 6 7
(Д ,Т) (0,2; 1) (0,4; 1) (0,6; 1) (0,8; 1) (0,2; 1) (А) = 0,2, Т = 1 А = 0,2, (Т) = 1
ы 1,4 • Ю-5 1,2 • 10-® 5,7- Ю-7 1,4 • 10-® 0,196 0,196 0,268
А /ь2 6,05 • Ю-2 0,213 0,16 0,527 0,221 0,205 0,227
А/с 0,314 0,53 0,479 0,633 0,47 0,597 0,520
Л? 500 250 165 125 500 500 500
п — — — — 152 158 158
Табл и ца 2
Реальный сигнал
Вариант способа восстановления 1 2 3 4 5 6 7
(А, Т) (0,2; 1) (0,4; 1) (0,6; 1) (0,8; 1) (0,2; 1) (А) = 0,2, Т = 1 А = 0,2, (Т) = 1
Ы 1,4 • 10-® 2,8 • 10-® 1,3 • ю-7 6,0 • ю-8 0,0103 1,05 • Ю-2 2,549-Ю-2
А/ь2 0,766 2,591 0,369 0,677 0,20 0,398 0,150
А/с 0,134 1,419 0,263 0,586 0,12 0,919 0,102
Л? 426 213 142 106 425 425 425
п — — — — 144 163 151
Отметим также, что погрешность восстановления сигнала в непосредственной близости к границам временного отрезка обычно больше погрешности восстановления внутри него. Необходимо подчеркнуть, что, несмотря на малое значение величины Аё в точках измерений, отклонение формы восстановленного сигнала может быть существенным.
Восстановление сигнала по подмножеству измерений. Учитывая, что точность восстановления сигнала на основе сплайн-интерполяции зависит от числа интегральных измерений немоно-
-.-■-.-■-.-■-.-■-.-■ -1,0-^-1-■-1-■-.-■-1-
0 20 40 60 80 100 ь 0 20 40 60 Ь
Восстановление модельного и реального сигналов: динамика модельного (а) и реального (б) сигналов и их восстановленной зависимости измерительной системой с А = 0,2, Т = 1; точность совпадения всех измерений восстановленного модельного (в) и реального (г) сигналов в зависимости от числа выбранных измерений; форма восстановленных модельного сигнала, п = 152, (д) и реального сигнала, п = 144 (е)
тонно, предлагается следующий метод определения числа необходимых измерений. Пусть выполнено N измерений с длиной окна Т и шагом А. Как известно, интеграл от любой линейной комбинации гармонических функций с периодами, кратными Т, на любом отрезке длины Т равен нулю. Следовательно, сигнал может быть восстановлен с точностью до такой линейной комбинации и необходимо каким-то способом выбрать решение из всего множества функций, соответствующих набору измерений. Поэтому задача является некорректной.
Выбор функции по принципу минимизации частот, кратных Т, является неправильным, так как эти частоты могут реально присутствовать в сигнале. Очевидно, что для А > Т (отрезки измерений не пересекаются) увеличение числа измерений добавляет информацию об отдельных участках сигнала и обычно увеличивает точность восстановления. Поэтому при реализации случая А ^ Т целесообразно использовать для построения сплайна все N измерений. Если же А < Т, то последовательно, начиная с п, равного целой части числа N• ф (число измерений п = N• ^ соответствует А « Т), выбираем п из N
измерений, не изменяя их положений по временной оси, и производим по ним восстановление сигнала с помощью сплайна из п—3 отрезков по алгоритму, описанному в предыдущем пункте. Для построенного сплайна будем вычислять точность Ай восстановления всех N измерений. В этом случае близкое к нулю значение величины Аё свидетельствует не только об успешном решении системы линейных уравнений, но и об успешном восстановлении даже тех измерений, которые не использовались для восстановления сигнала. При этом учитывается, что для терагерцового излучения слабой мощности вероятность появления отклика среды с существенно более высокой частотой, чем воздействующее излучение, близка к 0. Отметим, что при выборе п из N равномерно расположенных измерений сдвиг окна перестает быть постоянным, так как некоторые из измерений пропускаются. Если длина окна достаточно велика, то сначала точность восстановления будет возрастать (речь идет о возрастании в среднем, так как зависимость может оказаться немонотонной). И таким способом найдем минимальное значение п, при котором точность восстановления измерений Ай максимальна. Сплайн, построенный при таком п, будем считать восстановленным сигналом. В этом сигнале ложные частоты практически отсутствуют.
Хотя описанный алгоритм не является строгим, при его практическом применении к реальному терагерцовому сигналу проблем не возникало. В качестве примера на рисунке, в, г приведены зависимости точности по всем N измерениям от числа выбранных измерений в случае модельного и реального сигналов. Для обоих сигналов хорошо виден участок постепенного возрастания точности, а затем участок ее убывания и немонотонного изменения. Для модельного сигнала вычисленное таким способом п оказалось равным 152 (Ж = 500). При этом точность восстановления по всем 500 измерениям равна Ай = 0,196 (рисунок, д; табл. 1, вариант 5). Для реального же сигнала соответствующее значение п оказалось равным 144 (Ж = 425). При этом точность восстановления по всем 425 измерениям равна Ай = 0,0103 (рисунок, е; табл. 2, вариант 5). Как видно из рисунка, е, реальный сигнал восстановился практически точно и ложные осцилляции, подобные рисунку, б, в нем отсутствуют. Резкое увеличение точности совпадения измерений при п = N соответствует восстановлению всех N измерений в соответствии с точностью вычислений.
Восстановление сигнала при случайных сдвигах и размерах окна. При случайном сдвиге и размере окна, очевидно, результаты измерений с^ есть значения функции Используя интерполяцию, их можно вычислить для точек (^,Т), где Щ = 1ь + А • (г — 1) +Т/2, соответствующих постоянному сдвигу и длине окна, и затем по ним восстановить сигнал. Следует, однако, подчеркнуть, что интерполяция при флуктуациях длины окна обычно приводит к потере точности. Поэтому целесообразно непосредственно применить описанный в предыдущем пункте метод восстановления. В этом случае для определения начального значения п используются средние значения сдвига и длины окна.
Эффективность предложенного способа восстановления сигнала иллюстрируется, в частности двумя случаями расчетов, представленных в табл. 1, 2. Так, варианты 6 табл. 1, 2 соответствуют флуктуации сдвигов окна ±49% при постоянном его размере. Средний размер окна совпадает с вариантами 5. Вариантам 7 соответствует флуктуация длины окна ±10% варианта 5 при постоянном сдвиге окна (поскольку окно в пять раз больше сдвига, то амплитуда флуктуации такая же, что и для сдвигов окна). Точность восстановления в обоих случаях равна точности, полученной для измерительной системы с регулярными параметрами.
Восстановление сигнала в случае последовательности измерений с общей границей.
Для практики представляет интерес измерительная система, в которой, например, момент окончания всех измерений совпадает, а начало каждого измерения сдвигается. Такая последовательность измерений может быть реализована проще, чем соответствующий сдвиг измерительных окон. При этом получаем последовательность измерений с уменьшающимся размером окна. В этом случае начало измерений будет производиться в следующие моменты времени:
и ~ у = + {г - 1)1е /,• ± у /, • г = 1,2,..., Ж,
т. е.
к = + Ъ = + / 1.2.....Л\ (И)
Тогда вычитанием — dj, г = 1, 2,..., iV — 1, получаем вместо (11) набор непересекающихся измерений:
¿, = (<-0,5)^^, ъ ; 1.2.....Лч (12)
Важно подчеркнуть, что точность восстановления при этом получается выше, чем при таком же числе пересекающихся измерений. Для примера приведем следующие значения, которые необходимо сопоставить с табл. 1, вариант 1: Ad = 5,96• 10"8; N = 500; А = 0,2; Т = 0,2; Д/с = 0,11; Д/ь2 = 1,346• 10"2. Отметим, что значения Д и Т соответствуют набору измерений (12), а соответствующее значение Д, приведенное в табл. 1, вариант 1, округлено.
Заключение. Предложены две модификации алгоритма сплайн-восстановления сигнала тера-герцового излучения, прошедшего через вещество, и его спектра по набору интегральных измерений от него. Показано, что практически точное восстановление всех измерений не означает такого же точного восстановления искомого сигнала. Предложенный алгоритм, состоящий в выборе подмножества из всего множества проведенных измерений, позволяет существенно повысить точность восстановления формы сигнала, несмотря на некоторое снижение точности восстановления измерений. Предложен легкореализуемый на практике критерий качества восстановления сигнала. Проведенные расчеты показали высокую эффективность рассматриваемого алгоритма. Показана возможность практически точно восстановить сигнал и его спектр.
СПИСОК ЛИТЕРАТУРЫ
1. Berry Е., Boyle R.D., Fitzgerald A.J. et al. Time frequency analysis in terahertz pulsed imaging // Computer Vision Beyond the Visible Spectrum. Springer-Verlag, 2005. P. 290-329.
2. Smye S. W., С ha mber la i n J. M., Fitzgerald A. J. et al. The interaction between Terahertz radiation and biological tissue // Phys. Med. Biol. 2001. 46. N 9. P. R101-R112.
3. Berry E. Risk perception and safety issues // Biol. Phys. 2003. 129. P. 263-267.
4. Прохоров A.C., Анзин В. Б., Витухновский Д. А. Терагерцовая спектроскопия спиновых стекол AuFe // ЖЭТФ. 2006. 130. Вып. 6(12). С. 1027-1034.
5. Fitzgerald A.J.,Berry Е.,Miles R. Е. et al. Evaluation of image quality in terahertz pulsed imaging using test objects // Phys. Med. Biol. 2002. 47. N 21. P. 3865-3873.
6. Stringer M.R., Lund D.N., Foulds A. P. et al. The analysis of human cortical bone by terahertz timedomain spectroscopy // Phys. Med. Biol. 2005. 50. N 14. P. 3211-3219.
7. Reed S., В e г г у E., F о u 1 d s A. P. et al. Evaluation of an explanation for boundary artefacts in terahertz pulsed images using tooth slices and step wedges // Proc. of Medical Image Understanding and Analysis. University of Bristol, 2005. P. 131-134.
8. Cole В. E., Williams J. В., King B.T. et al. Coherent manipulation of semiconductor quantum bits with terahertz radiation // Nature. 2001. 410. P. 60-63.
9. Zandonella C. Terahertz imaging T-ray specs // Nature. 2003. 424. P. 721-722.
10. Tribe W. R., Newnham D. A., Tad ay P. F. et al. Hidden object detection: security applications of terahertz technology // Terahetz and Gigahertz Electronics and Photonics III. Proc. of SPIE. 2004. 5354. P. 168-176.
11. Nagel M., Forst M., Kurz H. THz biosensing devices: fundamentals and technology // J. Phys. Condens. Matter. 2006. 16. S601-S618.
12. Khazan M., Meissner R., Wilke I. Convertible transmission-reflection time-domain terahertz spectrometer // Review of Scientific Instruments. 2001. 72. N 8. P. 3427-3430.
13. http://www.teraview.com.
14. Сафонов B.H., Трофимов В. А., Шкуринов А. П. О точности измерения мгновенных спектральных интенсивностей фемтосекундных импульсов // ЖТФ. 2006. 76. Вып. 4. С. 78-85.
15. Trofimov V. A., Varentsova S. A. New method for analysis of temporal dynamics of medium spectrum under the action of terahertz pulse // Laser Physics and Photonics, Spectroscopy and Molecular Modeling VII. Proc. of SPIE. 2007. 653703.
16. Варенцова С. А., Трофимов В. А. Восстановление сигнала и его мгновенных спектральных характеристик методом скользящих окон // ЖТФ. 2007. 77. Вып. 5. С. 58-64.
17. Trofimov V.A., Troshchiev Y.V., Varentsova S. A. High effective method for temporal terahertz spectroscopy under the condition of random probe signals // Nonlinear Laser Spectroscopy and High-Precision Measurements; and Fundamentals of Laser Chemistry and Biophotonics. Proc. of SPIE. 2007. 67271H.
18. Nazarov М.М., Mukina L.S., Shuvaev А. V. et al. Exitation and propagation of surface electromagnetic
waves studied by terahertz spectrochronography // Laser Phys. Lett. 2005. 2. N 10. P. 471-475.
19. Скрипов Д. К., T р о ф и м о в В. А. Гистерезисиая зависимость некоторых спектральных компонент фем-тосекундного импульса при его распространении в нелинейной среде // ЖТФ. 2004. 74. Вып. 2. С. 77-82.
20. Чуй К. Введение в вэйвлеты. М.: Мир, 2001.
Поступила в редакцию 30.04.2008
УДК 553.4
В.И. Дмитриев, А.А. Канцель2, Е.С. Куркина3
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССОВ РАСТВОРЕНИЯ И ОТЛОЖЕНИЯ ПРИ ФИЛЬТРАЦИИ РАСТВОРОВ В ПОРИСТОЙ СРЕДЕ
Предложена математическая модель, описывающая фильтрацию растворов в пористой среде, процессы растворения и вторичного отложения минерала. Найдены автомодельные решения, описывающие движение переднего фронта и заднего фронта, — границы зоны полного растворения. Исследованы основные закономерности рассматриваемых процессов, проведены численные расчеты. Показано, что данная модель хорошо описывает экспериментальные данные по выщелачиванию минералов сернокислыми растворами. Исследована динамика извлечения полезного ископаемого из продуктивных растворов в среде с неоднородным распределением кислотности. Показано, что в зонах с повышенным содержанием рН происходит растворение минерала, а в зонах с пониженным значением кислотности происходит его вторичное отложение. В этом случае после завершения работ в пласте могут остаться более или менее значительные запасы полезного минерала в зависимости от объема зоны с пониженным содержанием рН и близости ее к откачной скважине.
Ключевые слова: математическое моделирование, инфильтрационный метасоматоз, автомодельные решения, неоднородное распределение рН, процессы растворения и вторичного отложения минерала.
Введение. Целью работы является исследование математических моделей процессов растворения при фильтрации растворов в пористой среде. В процессе растворения происходит изменение химического состава породы за счет замещения одних минералов другими, так что в целом порода сохраняет твердое состояние. Такие процессы называются метасоматическими. Теория метасоматической зональности была разработана академиком Д.С. Коржинским в середине прошлого века [1].
Метасоматоз бывает диффузионным и инфильтрационным. При чисто диффузионном метасоматозе перенос вещества осуществляется медленно посредством диффузии через застойные поровые растворы. При инфильтрационном метасоматозе компоненты переносятся течением водных растворов, просачивающихся через поры горных пород. В природе инфильтрационные процессы всегда сочетаются с диффузионными. Растворы просачиваются вдоль стыков зерен или обтекают отдельные более плотные участки породы, тогда как замещение этих зерен и участков породы происходит в результате диффузии компонентов.
Метасоматоз используется при добыче полезных ископаемых способом подземного выщелачивания (СПВ), который является одним из наиболее экологически чистых и распространенных геотехнологических методов добычи [2, 3]. Он заключается в том, что в подземный рудоносный пласт из ряда скважин закачивается кислотный или щелочной раствор. Раствор, просачиваясь сквозь пористую среду, растворяет минералы полезных ископаемых и, дойдя до откачных скважин, откачивается.
Математические модели метасоматоза должны учитывать фильтрацию раствора между зернами минерала, диффузию раствора внутри пор зерен, растворение и осаждение минерала. Объемное растворение различно в разных участках пласта и происходит только тогда, когда концентрация в растворе меньше предельной концентрации насыщенного раствора. Понятно, что вблизи закачивающих
хФакультет ВМиК МГУ, д.ф.-м.н., проф., e-maihdmitrievQcs.msu.su
2Факультет ВМиК МГУ, асп., e-maihkantselQmail.ru
3Факультет ВМиК МГУ, д.ф.-м.н., в.н.с., e-maihelena.kurkinaQcs.msu.su
