CC BY
10
УДК 543.4, 54.08, 53.088
ОСОБЕННОСТИ РЕШЕНИЯ ОБРАТНОЙ ЗАДАЧИ ОПТИЧЕСКОГО СМЕШЕНИЯ МЕТОДОМ СИНГУЛЯРНОГО РАЗЛОЖЕНИЯ МАТРИЦЫ ДЛЯ АНАЛИЗА СОСТАВА
ПОЛИДИСПЕРСНЫХ ВЗВЕСЕЙ НАНОРАЗМЕРНЫХ ЧАСТИЦ
1 1 12 2 Г.М. Янковский , Д.В. Кузнецов , С.Э. Кондаков ' , М.Я. Мельников
(1НИТУ МИСиС; 2Химический факультет МГУ им. М.В. Ломоносова; email: gyankovskii@mail.ru)
Представлены результаты численного решения обратной задачи спектроскопии оптического смешения методом сингулярного разложения матрицы (SVD). Используемый алгоритм учитывает положительность решения и постоянную фоновую составляющую. Приведены результаты восстановления модельных моно-, би- и полимодальных распределений. На примере численно сгенерированных (измеряемых с учетом шума) спектров изучены влияние шума измеренного спектра, величины фоновой составляющей и выбора частотного диапазона измерения сигнала на характеристики распределения, восстановленного методом SVD.
Ключевые слова: сингулярное разложение матрицы, спектроскопия оптического смешения, наночастицы в растворе, распределение по размерам.
Ранее нами были рассмотрены результаты численного решения обратной задачи оптического смешения методом регуляризации Тихонова для полидисперсных систем [1].
В данной работе рассматривается решение обратной задачи спектроскопии оптического смешения методом сингулярного разложения матрицы СЛАУ, известного в англоязычной литературе как метод SVD или TSVD (truncated singular value decomposition method или TSVD method) [2].
Решаемая задача представляет собой систему уравнений
относительно функции Z(r), характеризующей распределение рассеивающих частиц взвеси по размерам. Величина Г связана с радиусом R частиц:
г П16п2 . 2 в F = D——sin
X2 2
D =
h т
67ГГ|7?
где D - коэффициент диффузии частицы радиусом R в жидкости, X - длина волны монохроматического света в жидкости, 9 - угол рассеяния анализируемого излучения, Т - температура жидкости, п - динамиче-
ская вязкость жидкости. По виду функции Z(Г) можно установить наличие в составе взвеси частиц тех или иных размеров [3], а зная оптические свойства рассеивающих частиц, установить распределение концентраций частиц в зависимости от R.
Для решения системы уравнений (1) сведем ее к СЛАУ:
(2)
или
где
A z = g,
(2а)
Как было показано нами ранее [1], при применении метода регуляризации Тихонова для приближенного решения СЛАУ (2а) чувствительность решения г к погрешности Sg правой части (2) определяется свойствами матрицы £ в сингулярном разложении матрицы A: А = и£Ут, где и и V- ортогональные матрицы, £ - диагональная матрица с неотрицательными диагональными компонентами. Выберем £ так, чтобы ее компоненты на диагонали располагались в невоз-
растающем порядке. Точное решение г системы (2а), метода 8УБ, а также свойства получаемых решений
формально определяемое как
z = A-'g = (USVT) 1 g = VST1 UTg, при наличии погрешности Sg равно
z + Sz = VS_1 UT (g + Sg).
При этом
N
где
Такой выбор позволяет повысить устойчивость решения г к погрешности правой части (2а):
При этом невязка, даваемая решением z, равна
в зависимости от параметров алгоритма и параметров обрабатываемого сигнала. Описываемый алгоритм предназначен для решения задачи (2а), обобщенной с учетом фонового сигнала Ф, не зависящего от частоты:
Az + Ф = g,
(3)
где
отношение максимальном и минимальном
где
компонент матрицы которое может быть очень велико.
Метод 8УБ состоит в замене матрицы 5 матрицей нолучающейся из £ обнулением всех ее диагональных компонентов, кроме первых (наибольших)
АТ-компонентов: £ = diag(^s'1, 52, ..., % 0 ..., 0).
~ т ~ ~ ~
Перепишем новую СЛАУ LS'Fz = gв виде = Ь,
где р = Утъ и Ь = ит%. Приближенное решение этой
СЛАУ можно выбрать, например, в виде
При решении системы (3) фоновый сигнал Ф моделировался дополнительным лоренцианом большой ширины:
Фоновая составляющая Ф в спектрах оптического смешения может являться следствием шума фотоприемника, тракта регистрации и усиления сигнала, а также может включать достаточно широкие по сравнению с полосой обзора экспериментальной установки спектральные составляющие, обусловленные рассеянием лазерного излучения растворителем.
Будем искать такое приближенное решение г задачи (3), когда 2] > 0 при у = 1, ..., N + 1. Этого можно добиться, если вспомнить, что выбор решения (2)-(3) в форме
где
в значительной мере произволен. Воспользуемся вместо него в качестве решения (3) вектором
min(|piPi -"щ р
2II Н2^
■У Рп ,
(4)
Последняя сумма дает вклад в невязку ||г|| , который не может быть уменьшен и не зависит от применяемого приближенного метода решения (2а). Меняя параметр К, можно осуществлять баланс между малым значением ||5г ||/||г || (при малом К) и малым значением ||г|| (при большом К).
Рассмотрим практическую реализация алгоритма решения задачи оптического смешения с помощью
при условии х > 0.
Здесь р: = (рр р^ рп = р+р PN), Ь1 = (Ь^ ЬК), у - параметр метода, - матрица размером К*К
вида = diag(s1, ..., зК). При решении (3)-(4) будем предполагать, что величина погрешности |^|| правой части (3) известна. Описываемый алгоритм включает: 1) задание начального (достаточно малого) значения К < N
2) задание начального (достаточно большого) значения у;
3) решение (4);
4) вычисление невязки ||Аг - g|| по найденному в предыдущем пункте решению г;
5) уменьшение у (в том случае, если ||Аг - g|| >
5) при условии у > умин возврат к решению (4); в ином случае увеличение К на 1 и возврат к пункту 2.
Таким образом, для найденного решения
> 0 и ||Аг - < М.
Продемонстрируем работу алгоритма на примере модельных распределений в зависимости от вида распределения г, величины фона Ф и величины погрешности правой части |^||. Шум компоненты gl вектора g (с учетом фона) моделировался нормально распределенной случайной величиной с нулевым математическим ожиданием и среднеквадратичным отклонением о^.
Определение характеристик мономодального распределения размеров
Приведем результаты работы алгоритма на примере мономодальных распределений в зависимости от величины фона Ф, величины о, характеризующей относительную погрешность правой части и положения моды Г задаваемого распределения. Далее везде предполагаем ю. = I, I = 1, ..., N Г. = 1, 1 = 1, ..., N N = 50. При решении (3)-(4) для модельного g фоновый сигнал моделировался дополнительным лоренцианом с полушириной Г№1 = 20ю^ При этом на всем диапазоне рассматриваемых частот от ю1 до юN «фоновый» лоренциан меняется на 0,25%.
На рис. 1.2-1.8 показано влияние величины шума о в правой части (2а) на вид восстановленного распределения Z(Г). При достаточно малом значении о с большой точностью определяется положение моды Г заданного распределения, а также значение фона Ф. Определялась также величина ошибки восстановления фона ДФ/Ф. Эта величина представлена в подписях к рисункам в тех случаях, когда вычисленное значение фона не равняется 0. Ошибка в восстановлении фона определяется как
где Фвыч - вычисленное значение фона, Ф - действительное (заданное в распределении) значение фона.
При величине шума о в диапазоне от 0,001 до 0,01 уширение распределения все еще невелико, его центр совпадает с Г заданного модельного распределения. При дальнейшем росте шума погрешность определения фона резко возрастает, восстановленное распределение оказывается либо размытым, либо сдвинутым вправо относительно Г модельного распределения.
Рассмотрим влияние величины фона Ф на характеристики восстановленного распределения Z(Г). Рост фона Ф приводит к росту правой части (2а) при неизменном модельном распределении Z(Г). Вследствие этого при неизменной относительной величине шума о абсолютная погрешность правой части (2а) растет. В результате при очень больших значениях Ф существенный вклад в общую погрешность дает величина оФ. Это приводит к размытию восстановленного распределения и к сдвигу центра распределения в сторону больших значений Г. Данные закономерности иллюстрируют рис. 1.5-1.8.
Рассмотрим влияние положения моды Г модельного распределения (рис. 2, 3) во всем рассматриваемом диапазоне значений Г. от 1 до 50 (данный диапазон в нашем случае также совпадает с диапазоном рассматриваемых частот ю.). На рис. 2.2, 2.3, 3.2 и 3.3 показаны восстановленные распределения при модельных распределениях, изображенных на рис. 2.1 и 3.1. Сравнение этих рисунков показывает, что сдвиг Г модельного распределения в сторону больших значений приводит к уширению восстановленного распределения Z(Г). При этом точность определения фона также падает. Такое ухудшение восстановленного распределения вызвано тем, что для достаточно широких относительно диапазона рассматриваемых частот (до ю^ лоренцевых кривых абсолютная погрешность задания высокочастотной части спектра оказывается выше, чем для узких кривых. Поэтому следует выбирать юN и ^ так, чтобы для предполагаемых значений Г искомого распределения выполнялось условие
Г < 0^ (или Г < 0,4 Г^.
При этом значение Г должно быть достаточно велико по сравнению с Г1 (несколько Г1).
Определение характеристик бимодального распределения размеров частиц
Рассмотрим влияние шума и расстояния между модами на характеристики восстановленного распределения Z(Г) в случае бимодального модельного распределения (рис. 4, 5).
Рис. 1. Восстановление мономодального распределения. Влияние параметров а и Ф
1,2
1,0
л
£ 0,8
ч с 0,6
< 0,4
0,2
0,0
1,2 1,0
л
< 0,4 0,2 0,0
^
к
ч
с
<
10
20
30
10
20
30
1,2 1,0 0,8 0,6 0,4 0,2 0,0
10
20
30
Рис. 2.1. Модельное распределение Ъ (20) = 1
Рис. 2.2. а = 0,001
Ф = 0,01 АФ/Ф = 2,1
Рис. 2.3. а = 0,01 Ф = 0,01 АФ/Ф = 25%
Рис. 2. Влияние положения моды Г на восстановление распределения
1,2-.
§ 0,8 н к
0,4-
0,0
20
30
40
Рис. 3.1. Модельное распределение Ъ (30)= 1
л
1,2-,
5 0,8 Н £
| 0,4-
0,0'
■20
Ж
30
40
Рис. 3.2. а = 0,001 Ф = 0,01 Дф/ф = 17%
1,2
I 0,8
к
I олА
0,0
20
30
40
Г
Рис. 3.3. а = 0,01 Ф = 0,01
Вычисленное значение фона Ф = 0,00171
Рис. 3. Влияние положения моды Г на восстановление распределения
Рис. 4. Влияние шума а на восстановление бимодального распределения
На рис. 4.2-4.4 видно, что с ростом шума а ухудшается разрешение мод бимодального распределения. При достаточно большом значении а восстановленное распределение (рис. 4.4) уже не разрешает двух мод и имеет один максимум. При увеличении расстояния между модами последние начинают разрешаться при шумах а = 0,01 (рис. 5.2). При уменьшении расстояния между модами модельного распределения (рис. 5.1) моды перестают разрешаться даже при малых шумах а (рис. 5.2-5.3). Так, для близких к реальным значениий а, равных 0,001 и 0,01, минимальные значения относительного расстояния между двумя модами составляют соответственно
Д Г/ Гср. - 0,4 и А Г/ Гср.
0,6;
где Г1 и Г2 - полуширины двух мод заданного распределения.
Анализ полимодальных распределений размеров
частиц
Рассмотрим примеры восстановления распределений 2(Г) в случае полимодальных распределений (рис. 7, 8). При восстановлении тримодально-го распределения (рис. 7.1) удается получить три моды (рис. 7.2) только при низких уровнях шума (а = 0,0001), что соответствует идеализированным условиям. Средняя мода и крайняя правая мода (с
Рис. 5. Влияние расстояния между модами Г1 и Г2 на восстановление бимодального распределения
Рис. 6. Влияние расстояния между модами Г1 и Г2, и шума о на восстановление бимодального распределения
наибольшей частотой) сдвинуты вправо, и только крайняя низкочастотная мода восстанавливается с правильной частотой. При увеличении уровня шума до более реалистичного уровня (о = 0,001) в тримо-дальном распределении удается восстановить только две моды: крайнюю низкочастотную и среднюю, имеющую сдвиг вправо.
При восстановлении пентамодальных распределений при уровне шума измеряемого сигнала о = 0,001 полученное распределение имеет три моды. Точно восстановить удается только крайнюю моду, имеющую наименьшую частоту. Одна из средних мод распределения и крайняя мода, имеющая наи-
большую частоту, восстанавливаются со сдвигом вправо.
Заключение
В работе рассмотрено решение обратной задачи оптического смешения методом сингулярного разложения матрицы (SVD) для разных распределений размеров динамически рассеивающих частиц. Описан алгоритм поиска приближенного решения задачи оптического смешения с учетом положительности компонент вектора решения, а также с учетом постоянной фоновой составляющей сигнала. Учет положительности компонент решения СЛАУ (к
Рис. 7. Влияние шума а на восстановление тримодального распределения
Рис. 8. Влияние шума а на восстановление пентамодального распределения
которой сводится задача оптического смешения) позволяет повысить устойчивость приближенного решения обратной задачи. Для определения момента остановки алгоритма требуются априорные сведения о величине шума входного сигнала. Приведены результаты восстановления модельных мономодальных, бимодальных и полимодальных распределений. При этом правая часть решаемой СЛАУ генерировалась численно, шум в правой части СЛАУ моделировался случайной величиной, распределенной по Гауссу с нулевым математиче-
ским ожиданием и дисперсией, пропорциональной сигналу на данной частоте.
Рассмотрен вопрос о влиянии шума измеряемого сигнала (правой части СЛАУ) на характеристики решения обратной задачи методом 8УБ для моно- и бимодальных распределений размеров частиц в анализируемой взвеси. При достаточно малом значении шума а программа с большой точностью определяет положение Г моды заданного мономодального распределения, а также значение фона. При величинах шума а от 0,1 до 1,0% уширение распределения
все еще невелико, его центр совпадает с положением моды заданного модельного распределения. При дальнейшем росте шума погрешность определения фона резко возрастает, восстановленное распределение оказывается либо сильно размытым, либо сдвинутым относительно моды модельного распределения. Анализ влияния шума на разрешимость двух мод бимодального распределения показал, что для близких к реальным значениям о (0,1 и 1,0%) минимальные значения относительного расстояния между полуширинами Г двух мод равны соответственно
А Г/ ГСр. « 0,4 и А Г/ ГСр. « 0,6.
Рассмотрено влияние величины фона Ф на характеристики восстановленного распределения. Рост фона Ф приводит к росту измеряемого сигнала при неизменном модельном распределении, вследствие чего при неизменной относительной величине шума о абсолютная погрешность правой части СЛАУ растет. В результате при очень больших значениях Ф существенный вклад в общую погрешность дает величина оФ. Это приводит к размытию восстановленного распределения и к сдвигу центра распределения в сторону больших Г.
Изучение вопроса о том, как следует выбирать диапазон рассматриваемых частот юмакс и ширины ло-ренцианов Гмакс, основываясь на априорных знаниях об искомом распределении рассеивающих частиц по размерам, показало, что юмакс и Гмакс следует выбирать так, чтобы для предполагаемых значений Г искомого распределения выполнялось условие
Г < 0,4 «макс (или Г < 0,4ГМаКс).
При этом значение Г должно быть достаточно велико по сравнению с полушириной Г1 самого узкого рассматриваемого лоренциана (несколько Г1). Выбор меньших юмакс и Гмакс приводит к росту абсолютной погрешности задания высокочастотной части спектра (так как «хвосты» лоренцианов вычисляемого распределения оказываются еще слишком большими на границе юмакс рассматриваемого частотного диапазона) и ухудшению качества решения.
При восстановлении полимодальных распределений (три и более мод в распределении частиц по размерам) при уровне шума выше 0,1% не удается идентифицировать моды исходного модельного распределения. Однако при достаточно малых значениях шума для тримодальных распределений удается получить достаточно близкое к исходному распределение. При восстановлении пентамодальных распределений удается достаточно близко восстановить крайние моды распределений и одну из средних мод. При увеличении расстояния между модами, возможно, удастся более точно восстанавливать полимодальные распределения. Результаты этой работы позволяют надеяться, что при хорошем качестве приборной составляющей аппаратно-программного комплекса, обеспечивающей низкий шум и фон, и при усовершенствовании алгоритма 8УЭ удастся восстанавливать с приемлемой точностью полимодальные распределения. Сравнение результатов, полученных методом 8УЭ и методом регуляризации Тихонова, рассмотренных нами в работе [1], показывает, что использование метода 8УЭ более предпочтительно для восстановления полимодальных распределений.
При одинаковых параметрах численного моделирования метод 8УЭ позволяет получить вдвое лучшее разрешение мод. Для бимодального распределения при о = 0,1% минимальные значения относительного расстояния между полуширинами двух мод равны (для метода 8УЭ А Г/ Гср ~ 0,4, а для регуляризации Тихонова А Г/ Гср ~ 0,8).
После доработки метод 8УЭ предполагается использовать для анализа данных при исследовании биологических растворов, содержащих большое количество компонент, в частности сыворотки крови. Помимо увеличения соотношения сигнал/шум улучшению качества восстанавливаемых распределений будут способствовать правильный подбор диапазона рассматриваемых частот и увеличение набора базисных лоренцианов, а также априорная информация о составе мод восстанавливаемого распределения.
Работа выполнена при поддержке Министерства образования и науки РФ в рамках Государственного контракта № 14.518.11.7041, соглашения № 14.А18.21.0842 и программы фундаментальных исследований Президиума РАН № 28, подпрограмма 1, раздел «Физика, химия и биология воды»
СПИСОК ЛИТЕРАТУРЫ
1. Янковский Г.М., Кузнецов Д.В., Кондаков С.Э., Мельников М.Я. // Вестн. Моск. ун-та. Сер. 2. Химия. 2013. 54. С. 278.
2. Hansen P.C. The truncated SVD as a method for regulariza-tion // BIT. 1987. 27. Р. 543.
3. Тихонов А.Н., Гончарский А.В., Степанов В.В., Ягола А.Г. Численные методы решения некорректных задач. М., 1990.
Поступила в редакцию 20.06.13
SOLUTION OF THE INVERSE PROBLEM OF LIGHT BEATING SPECTROSCOPY USING SINGULAR VALUE DECOMPOSITION (SVD) METHOD FOR THE ANALYSIS OF POLYDISPERSE SUSPENSIONS OF NANOPARTICLES
О.М. Yankovskii, D.V. Kuznetsov, S.E. Kondakov, М.^ Мelnikov
(National University of Science and Technology, Faculty of Chemistry M.V. Lomonosov Moscow State University)
The results of numerical solution of the inverse problem of light beating spectroscopy using using singular value decomposition (SVD) method are presented. The developed algorithm takes into acount the positivity of the solution and the constant background signal. Reconstructed uni and multimodal distributions are presented. The effect of the noise of the input signal, the background signal and the frequency bandwidth of the measured signal on the characteristics of the reconstructed distributions is discussed.
Key words: singular value decomposition, SVD, light beating spectroscopy, nanoparticles, solution, size distribution.
Сведения об авторах: Янковский Георгий Маркович - науч. сотр. кафедры ФНС и ВТМ НИТУ МИСиС, канд. физ.-мат. наук (gyankovskii@mail.ru); Кузнецов Денис Валерьевич - зав. кафедрой ФНС и ВТМ НИТУ МИСИС, канд. техн. наук (dk@misis.ru); Мельников Михаил Яковлевич - зав. лабораторией химической кинетики химического факультета МГУ, докт. хим. наук (melnikov@excite.chem.msu.ru); Кондаков Сергей Эмильевич - вед. науч. сотр. кафедры химической кинетики химического факультета МГУ, докт. фарм. наук (kse@excite.chem.msu.ru).
