О влиянии шумов в спектрах на ошибки восстановления концентраций веществ при их обнаружении в атмосфере по данным фурье-спектрорадиометра Текст научной статьи по специальности «Математика»

электронное научно-техническое издание

НАУКА и ОБРАЗОВАНИЕ

Эл № ФС 77 - 30569. Государственная регистрация №0421100025. ISSN 1994-040S

О влиянии шумов в спектрах на ошибки восстановления концентраций веществ при их обнаружении в атмосфере по данным фурье-спектрорадиометра

77-30569/255578

# 11, ноябрь 2011

Морозов А. Н., Фуфурин И. Л.

УДК 517.984.54

МГТУ им. Н.Э. Баумана [email protected]

Введение

Обработка результатов измерений является одной из ключевых задач построения фурье-спектрорадиометрических систем (ФСРС) и подразумевает разработку алгоритмов, позволяющих как выполнять преобразование регистрируемой информации в форму спектральных данных, так и интерпретировать спектральную информацию. Последнее требование подразумевает выполнение требования устойчивости получаемого приближенного решения для концентраций веществ в смысле численных методов [1]. Зарекомендовавшим себя методом анализа и интерпретации спектров, является регрессионный анализ нормированного спектра целевых вещества в атмосфере [2, 3]. Для расчета нормированного спектра регистрируются два спектра: суммарный спектр трассы наблюдения при ожидаемом присутствии целевого вещества и спектр подстилающей поверхности без наличия целевых веществ на линии обзора. Однако в ряде случаев принципиально невозможно регистрировать фоновый спектр трассы наблюдения, и необходимо проведение газоанализа в открытой атмосфере по единичному зарегистрированному спектру [3]. В [4] описана методика идентификации и восстановления концентраций веществ по единичному спектру, полученному по пассивной методике регистрации на ФСРС.

В основе интерпретации единичного спектра лежит линейная слоевая модель переноса излучения, позволяющая свести задачу восстановления концентраций к системе линейных алгебраических уравнений (СЛАУ). Данная задача является некорректно поставленной, а матрица системы в ряде случаев может быть плохо обусловленной. В данной статье предлагаются численные процедуры, позволяющие на основе предварительного моделирования спектра атмосферной трассы исключить из

рассмотрения атмосферные составляющие. Затем на основе подходов корреляционного анализа многомерная СЛАУ заменяется одной или несколькими одномерными задачами, для которых выполнены все условия корректности.

В настоящей работе исследовано влияние малых шумов во входной информации на решение задачи, определены уровни ошибок в решении обратной задачи для характерных шумов во входных спектрах.

Постановка задачи

В настоящей работе предлагается подход к решению задачи идентификации веществ и определения их концентраций (не являющихся непосредственно наблюдаемыми величинами) по свойствам их спектров, которые наблюдаются непосредственно и напрямую зависят как от искомых величин, так и от внешних факторов.

Рис. 1. Упрощенная физическая модель методики беспробоотборного мониторинга

веществ в открытой атмосфере

Фурье-спектрорадиометр (ФСР) регистрирует собственное тепловое излучение всех объектов (атмосфера, целевые вещества, подстилающая поверхности) на линии обзора. Считаем, что температуры газовых компонентов на протяженности трассы наблюдения постоянны, что оправдано для приземных трасс наблюдений, где температурные перепады вдоль трассы не превышают нескольких градусов. Если данное требование не выполнено, то считаем, что газовая составляющая излучает при некоторой эффективной температуре. Считаем, что подстилающая поверхность излучает как абсолютно черное тело, что оправдано для большинства как естественных, так и искусственных поверхностей. Рассеянием излучения пренебрегаем. В [4] показано, что в данном случае нормированный спектр газовой составляющей запишется в виде:

(V) = е

-О (V) = В(V) - Р (Г19 V) (1)

Р(Г0,V) - Р(Г19V)'

где т, (V) - нормированный на единицу коэффициент пропускания, содержащий в себе

сигнатуры всех объектов и веществ на линии наблюдения ФСР; В(у) - спектр,

регистрируемый в эксперименте; Р(Т, V) - функция Планка; Т1 - эффективная

температура газовой составляющей трассы наблюдения; Т0 - эффективная температура

подстилающей поверхности; О, (V) - оптическая плотность трассы наблюдения.

Для смеси невзаимодействующих газов (отсутствие ассоциации, диссоциации и иного химического взаимодействия) оптическая плотность - аддитивная величина [2], поэтому уравнение (1) можно переписать в виде СЛАУ:

Л (2)

X К (V) ■ С1 = О, (V)'

1=\

где к¡(V) - спектральный коэффициент экстинкции данного вещества; С1 - интегральные

концентрации веществ находящихся на линии обзора. В матричном виде уравнение (2) можно записать следующим образом:

КС = О,, (3)

С е[0,СП'

где К = [кг (V])}- матрица референтных спектров веществ (I - количество референтных спектров веществ, входящих в газовую смесь; J - размер спектра по волновым числам); С = [С±}=1 1 - массив (вектор) искомых интегральных концентраций (которые из физических соображений должны быть положительны и не превышать концентрации насыщения Сшах); О, = {О,(V])}^~l"J - вектор значений оптической

плотности трассы наблюдения, содержащий ошибки измерения. Значения ку, составляющие матрицу К, как правило, определяются в лабораторных экспериментах, и поэтому также содержат ошибки измерения.

Задачи подобного рода относятся к классу так называемых обратных задач [1], т.е. задач, где необходимо определение интересующих нас количественных характеристик явления по результатам измерений их косвенных проявлений. Поскольку матрица К системы и правая часть О, содержат погрешности измерений, то система вида (3) почти

всегда несовместна. В этом случае рассматривается решение методом наименьших квадратов:

С = (КТК )-1 Ктвъ, (4)

которое, во-первых, может не принадлежать множеству определения искомых переменных С1 е [0,С™ах], а, во-вторых, в ряде случаев, может быть неустойчиво в асимптотическом смысле. Более того, если даже формально решение устойчиво, но норма

матрицы

КТК

может быть очень мала. Это означает, что небольшие ошибки во

входных данных приводят к огромным погрешностям в решении, что делает невозможной физическую интерпретацию поученного решения. Так в [5] показано, что для решения вида (4) 5%-ый шум во входном спектре ОЕ может привести для ряда веществ к 100 %

ошибке, а 10%-ый шум (например, для смеси веществ этилен, аммиак, метанол, этанол) уже к ошибкам порядка 1000%. Исходя из выше сказанного, обратная задача (3) относится к классу некорректно поставленных задач [1], для решения которой необходимо применение специальных методов решения, о которых будет рассказано в следующем разделе. Нужно отметить, что в случае критически низкой обусловленности матрицы

системы Кт К задачу следует сводить к идентификации не отдельных веществ, а классов веществ, для которых обусловленность системы будет гораздо выше. В этом случае можно говорить об идентификации в данном эксперименте какого-то вещества из данного класса и определения количественного вклада данного класса веществ.

Процедуры численного решения обратной задачи

В настоящей работе предложены численные процедуры, позволяющие решать систему вида (3) за строго ограниченное время, что позволяет использовать их при построении систем дистанционного газоанализа, работающих в режиме реального времени.

Правая часть системы (3) содержит сигнатуры как целевых веществ, так и атмосферных составляющих. На предварительном этапе обработки экспериментального спектра проводится математическое моделирование спектра атмосферных

составляющих, при этом не восстанавливаются точные значения концентраций атмосферных поглотителей. Полученный в результате математического моделирования спектр служит для выделения нормированного спектра целевых веществ, концентрации которых необходимо восстановить, при этом значительно сокращается размерность исходной линейной задачи, а, значит, и время необходимое на ее решение, что важно для систем, работающих в режиме реального времени.

Принимая во внимание аддитивность оптической плотности, можно записать разложение суммарного нормированного спектра на атмосферную составляющую и спектр целевых веществ:

т(v) т(v) т(У),Ф) е [0,1],

где Та(т (у) - смоделированный спектр атмосферной трассы, тшЬ (у) - спектр целевых веществ. Ограничения на спектры запишутся в виде:

Т (V) — Т ат (V). (5)

Далее покажем, что моделирование атмосферного спектра можно проводить не только на основе спектров газовых составляющих, но и на основе спектров атмосферных

трасс. Для этого рассмотрим М атмосферных спектров В1, у = 1.М , каждый из которых

определяется вкладом поглотителей Н20, С02, Оз, N02, N20 и некоторых других. Для смеси невзаимодействующих газов верно соотношение:

N

в = у , ¿=1

где 0] -у'-й спектр (оптическая плотность) атмосферной трассы; г - номер атмосферного газа; у - номер экспериментального спектра атмосферы; N - число атмосферных газов;

- спектр ¿-го атмосферного газа; 2]1 - относительная интегральная концентрация ¿-го

составляющего атмосферы, соответствующая у'-му атмосферному спектру 0], равная отношению текущей интегральной концентрации к концентрации, для которой получен референтный спектр. Для расчета атмосферного спектра используем линейную комбинацию спектров атмосферных трасс:

М , , (6) В = У В1 у1

ат ¿.^ I ' 1=1

где ВаШ = ВаШ (у) - спектр атмосферной трассы для текущего измерения, у] -коэффициенты пропорциональности. Покажем возможность получения суммы следующего вида:

N

Вт =У В2\ ,

¿=1

где 2 '¿ - относительные интегральные концентрации атмосферных составляющих, соответствующие экспериментальному спектру суммарной оптической плотности В2 (у),

M N M N N M N

Datm = IyjIdz =цDZ =ID= IDZ' .

j=1 г =1 j=1 ¿=1 ¿=1 j=1 ¿=1

Коэффициенты пропорциональности при оптических плотностях будут иметь вид:

M j ■ (7) Iy j

J=1

'jzj = z V.

Таким образом, без шумов в спектрах и при выполнении условий линейности атмосферного спектра относительно спектров поглотителей, для разрешения системы (7) относительно (у 1} (т.е. для того, чтобы была теоретическая возможность смоделировать любой атмосферный спектр в рамках линейной модели) необходимо, чтобы число спектром атмосферных трасс было не меньше числа основных атмосферных поглотителей

М > N, а также, чтобы эти спектры были линейно независимы гапк^^ } > N.

Для реализации предложенного алгоритма создана база данных нормированных спектров атмосферных трасс, зарегистрированных при различных углах возвышения и различных погодных условиях (влажность, время года), что позволило получить наиболее независимые спектры. Размер базы определялся числом основных поглотителей атмосферы, существенных для приземных трасс в окне прозрачности атмосферы 700-1500 см-1.

М

-Е уа (V)

Согласно (6) и с учетом (5) рассчитывается атмосферный спектр е к= , невязка которого с экспериментальным спектром т, (V) минимальна. Для этого строится функционал вида:

max

Ф(Z) = J (р2(Г,тБ(v)) + aQ(r,TS(v)))wXe(v)dv^mm

(8)

-Е ук°к (V) р(Г, тБ (V)) = е к= -тБ (V),

П((Г, тБ (V)) = е"р( Г,т*

где М - число спектров базы данных, используемых в оптимизации функционала; Г = (у-1}-1 л'м - коэффициенты пропорциональности; м?т (у) - весовая функция;

z

M

Ч™ -I YkDk(v)

J p (Г, Tz (v))w (v)dv - квадрат невязки расчетного e k= и экспериментального

max

Tz (v) спектров; J Q((r, Tz (v))wx^ (v)dv - штрафная функция, позволяющая

v

V

проводить минимизацию на множестве, описываемом неравенством (5); а - параметр штрафной функции. Весовая функция нормирована на единицу и отражает спектральную ошибку:

(V) =

Е О)

где Сх 2(у) - дисперсия спектральной ошибки, которая определяется как погрешность косвенных измерений:

Г л \2 (

2

_I

v

дВ

2

P (T0, v) - P(Tl9 v)

Здесь а2в = const - дисперсия спектральной ошибки в измеренном спектре B(v). Поскольку функционал (8) дифференцируем по Г, то для указанной оптимизации применяется метод градиентного спуска.

После моделирования спектра атмосферной трассы (точные значения концентраций атмосферных составляющих при этом не представляют интерес) вычисляют спектр целевых веществ:

Т-(v) = е-(v) е [0,1]. (9)

"zyA (v)

е

В результате система (3) преобразуется к виду:

'KC = Dmt, (10)

е [0,C,"x]'

Спектр (9) может содержать сигнатуры одного или нескольких целевых веществ, причем неизвестно, каких именно. Решение СЛАУ (10) затрудняется некорректностью постановки задачи и слабой обусловленностью матрицы системы, о чем сказано в разделе постановка задачи. В настоящей работе предлагается следующий подход: вместо решения многомерной СЛАУ используется корреляционный анализ и определяется целевое вещество, спектр которого статистически наиболее близок к экспериментальному Tsub (v). Для этого предварительно создана база данных референтных спектров целевых веществ, где для известных концентраций C'db хранятся референтные спектры т' (v) = е Ddb(v).

При проведении корреляционного анализа для каждого вещества из базы данных и экспериментального спектра рассчитывается коэффициент корреляции:

<т b-Т b,Т -Т1 > (11)

r (т Т1 ) = su b su b' db db

гА1шЫ1л) I— 2 1 —2 '

V< (Тsub Tsub ) ><Tdb -Tdb ) >

V

max

где <T1(v),T2(v) >= J T1(v)T2(v)w(v)dv - скалярные произведения,

V

min

V

max

Т = J T(v)w(v)dV - средние значения, w(v) - весовая функция. Вещество, для

V

min

которого коэффициент корреляции оказался наибольшим и превысил заранее заданное пороговое значение, считается идентифицированным для задачи в постановке (10). Пороговые значения определяются эмпирически. Для восстановления концентраций минимизируется функционал вида:

V/ , ч2 (12)

Ф(Z) = J(e"Ddb(V)Zl -Tsub (V)) w^ (v)dv^min,

V •

min

C

где z i = —j- - относительные интегральные концентрации, Wsub (v) - весовая функция,

Cdb

отражающая спектральную ошибку, а также учитывающая предварительные этапы обработки спектра. С учетом линейной зависимости спектральной ошибки о и нормированного спектра т, а также принимая во внимание (9), весовую функцию можно записать как

СТтЕ

CT =-TT-*-,

Wt„ (V)

-£ zkDk (v) 1

e k=1

Tsub

T„

< (V)

Поскольку функционал (12) дифференцируем по 2, его оптимизация проводится методом градиентного спуска. После идентификации одного вещества и восстановления его концентрации из спектра (9) вычитается вклад этого вещества, и процедуры (11) и (12) повторяются заново. Таким образом, вместо решения многомерной СЛАУ предлагается решать одну или несколько одномерных задач, для которых выполнены условия корректности, и решения которой в обычных на практике условиях позволяют идентифицировать целевые вещества и их интегральные концентрации

Исследование влияния малых изменений входной информации на восстановленные концентрации

Покажем, что для разработанной методики решения обратной задачи шумы во входном спектре, сравнимые по значению с шумами в реальном эксперименте, приводят к ошибкам в решении не больше допустимых значений для данного класса аппаратуры и методики регистрации спектральных данных. В противном случае, если малые изменения во входной информации приводят к значительным изменениям решений задачи, невозможно говорить о физической интерпретации полученного решения. Среднеквадратическое значение шумовой составляющей в спектре В(у),

полученного на фурье-спектрорадиометре, постоянно во всем спектральном диапазоне. В этом случае можно говорить об аддитивной гауссовской модели шума в исходном спектре:

В (у) = В(у) + 8В(у), (13)

где В (у) - спектр с учетом шумовой составляющей; В (у) - спектр без шума, 8В(у) -аддитивный шум в спектре, с математическим ожиданием Е(8В) = 0 и дисперсией

0(8В) = оВ2. Применяя зависимость (13), можно плавно изменять уровень входного

шума оВ в спектре В (у) и, зная точное значение концентрации целевого вещества в

данном численном эксперименте, исследовать возможность идентификации целевого вещества, а также зависимость восстановленных интегральных концентраций от среднеквадратического значения шумовой составляющей спектра.

В ходе вычислительных экспериментов проводилось математическое моделирование спектров атмосферной трассы при наличии целевого вещества (аммиака, как типичного вещества с узкополосным спектром, и изопропанола, как типичного вещества с широкополосным спектром). Для обоих случаев варьировались следующие параметры: исходная концентрация целевого вещества, температурный контраст, уровень шума в спектре трассы. Для заданных исходных концентраций целевого вещества и температурного контраста исследовалась зависимость относительной ошибки восстановления целевого вещества от уровня шума в спектре. Относительная ошибка восстановления концентрации веществ в конкретном эксперименте является случайной величиной, зависящей от реализации шума, как случайного процесса. На рис. 2 показана зависимость относительной ошибки восстановления концентрации аммиака от среднеквадратического значения входного шума, рабочее значение шума в исходном

спектре В(V) составляет порядка ~ 0,1 - отношение сигнал/шум),

температурный контраст «5°. Для каждого значения среднеквадратического шума рассчитано десять реализаций шумовой составляющей, и для каждой реализации получено значение ошибки восстановления. На рис. 2 8с - относительная ошибка

восстановленной концентрации. На рис. 3 а8 - среднеквадратическое значение

относительной ошибки восстановленной концентрации.

Рис. 2. Зависимость относительной ошибки восстановления концентрации аммиака от среднеквадратического значения входного шума на рабочем диапазоне шума. Жирной линией показана линейная регрессия.

\ISNR

Рис. 3. Зависимость среднеквадратического значения относительной ошибки восстановления концентрации аммиака от среднеквадратического значения шума в спектре Б(у). Жирной линией показана линейная регрессия. Из рис. 2 видно, что относительная ошибка имеет систематическую составляющую. Это объясняется тем, что для решения исходной задачи проводится математическое моделирование атмосферного спектра и, в итоге, вместо исходной СЛАУ решается одномерная система. В результате, появляется систематическая ошибка в решении, которой для больших значений шумов является определяющей в ошибке решения задачи. Применяя подобные численных эксперименты, можно выяснить при каких пороговых значения входных помех ошибка не превышает установленное допустимое значение и тем самым определить пределы применимости предлагаемой методики для разного типа веществ.

Далее проанализируем случайную составляющую относительной ошибки. Для восстановления концентрации идентифицированного вещества решается одномерная линейная задача:

ВЛъ (V) 2 = ВшЪ (V). (14)

Запишем оптическую плотность экспериментального спектра целевых веществ, выделив аддитивную шумовую составляющую:

Dsub (v) = Dsub (v) + 5Ds (v),

где Dsub (v) - зарегистрированное значение, Dsub (v) - точное значение, 8Dsub (v) -шумы в экспериментальном спектре. Т.к. шумы на разных частотах входного спектра, вообще говоря, независимые случайные величины, то

cov(SDsubk, Wsu\) = а2 8И,

где а = a2 (vk) - дисперсия шума на частоте vk, Skl - символ Кронекера. Не нарушая общности, нормируем исходное уравнение на а2 :

1 ^ 1

z—Dkz = — D

l ak ak

sub

k •

Для нормированного уравнения верно равенство СОУ(8О™ к, 1) = 8 Е(8ОетЪк ) = 0, т.е. шум равномерный. Решение уравнения (14) будем искать в виде:

и

k

\sub

DZ - Ds

^ min.

ZZ

(15)

где О - нормированная матрица уравнения. Решение Т уравнения (15) запишется в виде:

Т = (отО)-1 ВтОшЪ. (16)

Поскольку БзиЪк содержит шумы, то решение Т. является приближенным и запишется в

виде Т. = Т + 8Т . Тогда решение вида (16) можно переписать отдельно для точного значения входной информации и для шумов:

Т = (ОтО)-1 ВтВшЪ, 8Т = (ОтО)-1 От,

а поскольку ОаШ(у) - величина расчетная, то очевидно что 8О™Ъ = 8О,. Принимая во внимание (1) и экспоненциальную зависимость оптической плотности от нормированного коэффициента пропускания т^), очевидно, что среднеквадратичное отклонение

восстановленной концентрации а линейно зависит от среднеквадратичного значения

шума во входном спектре ав, что видно на рис. 3. Далее можно показать количественную связь среднеквадратического отклонения восстановленной концентрации и среднеквадратического значения шума в исходном спектре В^) :

Оя

ОС

1

^ (17)

г 1

2 О (Р - В)

Ч^1

где О^ = D(Vj) - оптическая плотность для искомого вещества, соответствующая референтному спектру; Р. = Р(Т1, Vj); В^ = B(Vj) . В серии численных экспериментов

показано, что количественные значения, описываемые зависимостью (17), и получаемые экспериментально являются величинами одного порядка. Расхождение в значениях коэффициентов пропорциональности можно объяснить тем, что для каждого нового численного эксперимента заново моделируется атмосферный спектр (7), что искажает спектр (9).

Заключение

Продемонстрированы численные процедуры решения обратной задачи спектроскопии, позволяющие на основе анализа единичного спектра собственного теплового изучения веществ в открытой атмосфере проводить их идентификацию и восстанавливать значения интегральных концентрации. Данная методика основывается на промежуточном моделировании атмосферного спектра; идентификации одного или нескольких целевых веществ корреляционными методами; решении вместо исходной многомерной СЛАУ одной или нескольких одномерных задач, для которых выполнены условия корректности, и решения которых в обычных на практике условиях позволяют восстанавливать интегральные концентрации целевых веществ.

Показано, что ошибка решения обратной задачи линейно возрастает с увеличением среднеквадратичного значения шума во входном спектре. В серии вычислительных экспериментов продемонстрировано, что значения относительной ошибки восстановления концентрации для рабочего уровня не превышают 30%, что является допустимым рабочим значением для беспробоотборных газоаналитических систем и методик. Методика реализована в программном комплексе и апробирована в серии как вычислительных, так и натурных испытаний и может быть использована для проведения дистанционного беспробоотборного газоанализа по единичному спектру.

Литература

[1] Тихонов А.Н., Арсенин В.Я. Методы решения некорректно поставленных задач. М: Наука, 1979. 285 с.

[2] Морозов А.Н., Светличный С.И. Основы фурье-спектрорадиометрии. М.: Наука, 2006. 275 с.

[3] Дворук С.К., Корниенко В.Н., Кочиков И.В. и др. Мониторинг загрязняющих веществ в атмосфере с помощью Фурье-спектрорадиометра // Опт. журн. 2004. Т. 71, № 5. С. 7-13.

[4] Кочиков И.В., Морозов А.Н., Светличный С.И., Фуфурин И.Л. Распознавание веществ в открытой атмосфере по единичной интерферограмме фурье-спектрорадиометра // Опт. и спектр. 2009. Т. 106, №5. С. 743-749.

[5] Козинцев В.И., Белов М.Л., Городичев В.А., Федотов Ю.В. Лазерный акусто-оптический анализ многокомпонентных газовых смесей. М.: Изд-во МГТУ им. Н.Э. Баумана, 2003. 352 с.

electronic scientific and technical periodical

SCIENCE and EDUCATION

_EL № KS 77 - 3Ü56'». .V;II421100025, ISSN 1994-jMOg_

Influence of spectral noise on errors of substance concentration recovery at substance detection in aero-space by Fourier spectroradiometer data

77-30569/255578

# 11, November 2011 Morozov A.N., Fufurin I.L.

Bauman Moscow State Technical University

[email protected]

The authors consider numerical solution algorythms for identification and concentrations calculation of substances in the atmosphere. Instead of solving the multidimensional system of linear equations it is proposed to solve one or several one-dimensional systems. Conversion to one-dimensional problem is performed by correlation analysis. The influence of noise in the initial spectrum on errors of inverse problem solution was studied.

Publications with keywords: inverse problem, FT1R spectroscopy, molecular spectroscopy, substance identification

Publications with words: inverse problem, FT1R spectroscopy, molecular spectroscopy, substance identification

Reference

1. Tikhonov A.N., Arsenin V.Ia., Methods of solution of incorrectly posed problems, Moscow, Nauka, 1979, 285 p.

2. Morozov A.N., Svetlichnyi S.I., Fundamentals of Fourier spectroradiometry, Moscow, Nauka, 2006, 275 p.

3. Dvoruk S.K., Kornienko V.N., Kochikov I.V., et al., Opticheskii zhurnal 71 (5) (2004) 713.

4. Kochikov I.V., Morozov A.N., Svetlichnyi S.I., Fufurin I.L., Optika i spektroskopiia 106 (5) (2009) 743-749.

5. Kozintsev V.I., Belov M.L., Gorodichev V.A., Fedotov Iu.V., Laser acousto-optical analysis of multicomponent gas mixtures, Moscow, Izd-vo MGTU im. N.E. Baumana -BMSTU Press, 2003, 352 p.