Корреляционная обработка кубоида инфракрасных изображений, получаемых с беспилотных летательных аппаратов. Часть 2. Метод обработки инфракрасных сигнатур эталонных объектов на основе численного решения нелинейной задачи теплообмена Текст научной статьи по специальности «Математика»

Journal of Siberian Federal University. Engineering & Technologies, 2016, 9(3), 376-384

УДК 518.6, 536.6

Correlation Processing Cuboid Infrared Images from Unmanned Aerial Vehicle. Part 2. Processing Method of Infrared Signatures of Reference Objects on the Basis of the Numerical Solution of the Nonlinear Heat Transfer Problems

Igor N. Ischukab, Andrey M. Filimonova, Konstantin V. Postnova, Evgeniy A. Stepanova and Dmitry D. Dmitriev*

aMilitary Training and Research Center of the Air Force «Air Force Academy ft. Professor N.E. Zhukovsky and Y.A. Gagarin» 54a Starykh Bol'shevikov Str., Voronezh, 394064, Russia

bSiberian Federal University 79 Svobodny, Krasnoyarsk, 660041, Russia

Received 01.12.2015, received in revised form 24.12.2015, accepted 17.01.2016

A version of a mathematical model of the thermal field of multilayer isotropic material, based on the solution of two-dimensional heat conduction problem of direct numerical method using implicit difference schemes, taking into account the mixed boundary conditions and interface thermal contact resistance. The results of the correlation processing algorithm of dynamic infrared signatures obtained by the daily thermographing "reference zones".

Keywords: numerical methods, thermal parameters, direct heat conduction problem.

Citation: Ischuk I.N., Filimonov A.M., Postnov K.V., Stepanov E.A., Dmitriev D.D. Correlation processing cuboid infrared images from unmanned aerial vehicle. Part 2. Processing method of infrared signatures of reference objects on the basis of the numerical solution of the nonlinear heat transfer problems, J. Sib. Fed. Univ. Eng. technol., 2016, 9(3), 376-384. DOI: 10.17516/1999-494X-2016-9-3-376-384.

© Siberian Federal University. All rights reserved Corresponding author E-mail address: flyfil87@mail.ru

Корреляционная обработка кубоида инфракрасных изображений, получаемых с беспилотных летательных аппаратов. Часть 2. Метод обработки

инфракрасных сигнатур эталонных объектов на основе численного решения нелинейной задачи теплообмена

И.Н. Ищука' б, А.М. Филимонов", К.В. Постнова, Е.А. Степанов3, Д.Д. Дмитриев6

аВУНЦ ВВС «ВВА им. проф. Н.Е. Жуковского и Ю.А. Гагарина» Россия, 371600, Воронеж, ул. Старых Большевиков, 54а бСибирский федеральный университет Россия, 660041, Красноярск, пр. Свободный, 79

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

Ключевые слов: численные методы, теплофизические параметры, прямая задача теплопроводности.

Введение

В последние годы как в России, так и в других странах мира в практику дистанционного зондирования Земли всё более ускоряющимися темпами внедряются технологии многоспектральной съемки. Данный вид съемки представляет собой новое перспективное направление исследования объектов поверхности Земли посредством изучения свойств объектов на основе анализа информации о распределении отраженного от них излучения в зависимости от длины волны. Многоспектральная съемка реализует для заданного района местности одновременное получение десятков и сотен изображений одной и той же сцены, зафиксированных в видимом и инфракрасном диапазонах длин волн в различные промежутки времени. При многоспектральной съемке наиболее значимыми являются данные, полученные в инфракрасном (ИК) диапазоне длин волн, так как информация о температурах объектов, расположенных на поверхности Земли, позволяет дистанционно определять их теплофизические свойства. Однако данная задача требует построения сложной математической модели процессов теплообмена, а также учета всего многообразия факторов, влияющих на процесс формирования радиационных температур и определяющих климат на поверхности Земли. Следовательно, можно сделать вывод об актуальности и целесообразности разработки принципиально новых программно-математических методов об- 377 -

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

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

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

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

Математическая модель

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

р>с> 1' г=1>2> от oz\ oz J ох \ ох ,

7 дв ~c*pkTz

OZ

oz

-Lpk —

г=0 dz

z=o dz

-B*=-q(T);

1

:=-hy

RT

-TL-

(1)

I z=-\

z=-hx OZ

z=+h

ТI =Т

где р* - плотность материала, кг-м-3; с* - удельная теплоемкость материала, Джкг-1К-1; Т - избыточная температура, К; т - время, с; X - теплопроводность, Втм-1К-1; ср - удельная изобарная теплоемкость воздуха, Джкг-1К-1; р - плотность воздуха, кг-м-3; к - коэффициент турбулентности в приземном слое атмосферы, м2-с; 9 - потенциальная температура, К, з - массовая доля водяного пара, %, 2 - вертикальная координата в глубину, м; х -координата в ширину, м; Ь - удельная теплота парообразования, Дж-кг;

Рис. 1. Схема распространения тепла в цилиндрических телах

В0 = 5 оТ04 - 8а<тТд - эффективное излучение поверхности первого слоя, Вт-м 2; д - по-глощательная способность поверхности первого слоя; о - постоянная Стефана - Боль-цмана, Втм-2К-4; д - плотность теплового потока, Вт-м-2; Ят - коэффициент термического сопротивления, м2-К-Вт-1; Н - координата нижней границы многослойного материала; Тсрс _ среднесуточная температура воздуха, К.

Разностный аналог выражения (1) представляется в следующем виде [3, 4]:

tk+1 _ tk tk+1 — 2tk+1 4- tk+1

n n -L n n rmk п «+ 1 n p n n—1

= al[TpJ

_ tk+1 — 2tk+1 + tk+1

p,n ^ p,n _ \Tk 1 1 P,n ^ 1 P~\,n ^ |J P,n+l P,n ' 1 p,n-\ ;

At

Tk _ Tk

ККЛ p^p+u-a-T1

Az2

Ax2

i = 1,2;

p,n

Win]-Ь 2a№„] Дг '

J1 к _ J1 к

rpk _

Tk _ 'Tk

// + 1 И //V

P,n J rjik

О гтк /ли 1/i+\,n 1 ju+l,n 1fi-\,n .

- —:--

Tt

k

1 H,n - Ср.с '

(2)

где Ткп - сеточная функция; к - номер отсчета по времени, к е [0, К]; р - номер отсчета в глубину, р е [0, Н]; п - номер отсчета по горизонтали, п е[0, Я]; а - коэффициент температуропроводности слоя, м2-с-1; а - коэффи,циент теплоотдачи поверхности, Втм-2К-1; / - номер отсчета, соответствующий границе сопряжения в многослойной среде; Ах - величина шага дискретизации по времени, с; А2 - величина шага дискретизации по пространству в глубину слоя, м; Ах - величина шага дискретизации по ширине, м; V - коэффициент тепловых потерь.

1-2/ / ехр

_ V

2/

ехр

к-Ат (к-К*)-Ат

+ ехр

-ехр

к-Ат

v A j j j (к - К*)-At

Pi

Е, к<К*;

Е, к>К*;

(3)

где К - номер отсчетов по времени, соответствующий времени действия теплового потока; рь в2 - параметры релаксации теплового потока на стадии нагрева и остывания соответственно, с1; Е - энергетическая светимость ИК-излучателя, Вт-м-2.

v;

k

q

Метод решения и исходные данные

Рассмотрим алгоритм расчета нестационарного теплового поля согласно (2) с использованием продольно-поперечной схемы и метода прогонки [5, 6]:

Шаг 1. Задание исходных данных, а также начальных и граничных условий, необходимых для расчета: к, ah a2, a, h¡, h2, Ах, Timp, Az, x, в\, Pi, E, X2, R, v, Тсрс, H = h+ h2. Задание граничных и начальных условий: T0(z,x) и Tr(zr,xr, т).

Шаг 2. Введение нелинейности теплофизических параметров. Аппроксимация теплофи-зических параметров исследуемых материалов кусочно-линейными функциями

P.+V-PnT1'^

/гР\ <Ткрп </2 Je 1,2; Pl-tj-Pl,T^>T2,

Р,г.

" Р,п

где Т1, Т2 - пороговые значения температур, при которых происходят изменения теплофизических параметров материалов; п - параметры зависимости теплофизических свойств от изменения температуры; р - теплофизические параметры исследуемых материалов. Шаг 3. Вычисление плотности теплового потока С согласно (3).

Шаг 4. Последовательное прохождение временных отсчетов по к от 0 до К с шагом

Ах.

Шаг 4.1. Прогонка пространственных отсчетов по п от -Я до Я с шагом Ах. Шаг 4.2. Прогонка по глубине р от 0 до Н с шагом А2 при целом временном шаге.

Шаг 4.2.1. Вычисление вспомогательных коэффициентов акр, Ькр, скр, dкр для р е [0,к ]:

а\тк 1-Аг 2• ал \тк 1-Ат ал\тк 1-Ат

_"1|_ р>пл . ь _ ^ Р*А 1.Р; _ ^ Р>п1 . н - тк-1

Чр А? ' кр А? ' к'Р р~1>п'

Шаг 4.2.2. Вычисление вспомогательных коэффициентов акр, Ькр, скр, dк>р на границе сопряжения двух сред (р = к):

а2[г;„]-Дг 2-а2[т^}Ат а2[т^]-Ат

ЧР= - . ,—; ьк,Р=----1; ск„= , ,—;

Л fe] ^¡Л y.t-1

Дz

\

; Tl~n =d

м-

Шаг 4.2.3. Вычисление вспомогательных коэффициентов акр, Ькр, скр, dк>р для р е [к1, Н]:

*2[г;,„].Лг 2.д2[г;,„].Дг а2[г;,„].Дг

^ X? ' ^ Д? ' Кр д? ' •

Шаг 4.2.4. Вычисление прогоночных коэффициентов екр,

Л ск,р п dк,р „ п dк,р - ак, р-1 ' 4,р-1 в граничных точках при р = 0: е к, р = ——; I к, р = -—; при р = Н: I к, р = --;

Ь к, р Ь к, р Ь к, р - а к, р-1 ' е к, р-1 ск,р_ ш п = dк,р - ак,р-1 ' 4,р-1

в интервале р е (0, H) : ek р =----; fk,р =

bк,p aк,р-\ 'eк,p-\ bk,p aк,р-\ ' ek,р-\

Шаг 4.2.5. Определение значений температур на целом временном слое обратной прогонкой прир = Н: Ткр-- = Ак р; при р е (Н,0] с шагом Дг: Ткр_1л -е^-Ткр п. Шаг 4.3. Организация прогонки на временном полушаге вдоль оси х. Шаг 4.4. Прогонка по глубине по р от 0 до к2 с шагом Дг при временном полушаге с условиями:

ит {«1 [ткР,п\р^га .Г, 1 [ткР,п]р^га " К )р ^ЮТ ^ ^ 2 \р

Шаг 4.5. Прогонка по ширине по п от -Я до Я с шагом Дх при временном полушаге.

Шаг 4.5.1. Вычисление вспомогательных коэффициентов ak, n, bkn, ck, n, dk,n для n e [-R,R]:

a _ «К . b _ 2 • a К c _ a К . d _ T

ak, n---—2-, bk, n----—2--1 , ck, n---—2-, dk, n--T t

Ax Ax Ax

Шаг 4.5.2. Вычисление прогоночных коэффициентов ek, n fk, n:

• f,,„ ,

bk,n

при n e(-R, R): ek

b k, n

_ d^iL 1 bk,n ; при n = R: f Lk, n d k, n b k, n - a k, n-1 - a k, n-1

Ck, n • f d k, n - a k, n-1 • f Lk, n-1

a k, n-1 ' Lk, n ' ek,n-1 b k, n - ak,n-1 ' ek,n-1

Шаг 4.5.3. Определение значений температур на временном полушаге обратной прогонкой при n = R: Tp, = dk n; при n е (R,-R]c шагом Ах: тk= ,х - et,x • Tp\„.

Шаг 4.6. Организация целого временного шага для прогонки по оси z. Возвращение к шагу 4.1. Повторение расчетов до момента времени k = K.

Шаг 5. Формирование численных значений кубоида избыточных температур и построение соответствующих графиков.

Результаты исследования

Для оценки возможности применения корреляционного алгоритма, описанного в ч. 1 статьи, и построенного алгоритма численного решения прямой задачи теплопроводности были проведены дистанционные измерения теплофизических параметров эталонных объектов путем получения их ИК-сигнатур. Съемка экспозиции объектов производилась тепловизион-ным приемником FLIR Tau 2, размещённым на беспилотном летательном аппарате в период с 18:45 16.10.2015 до 07:45 19.10.2015 с высоты 10 м. Общее время съемки составило 61 ч. В качестве эталонов использовались материалы с известными теплофизическими параметрами: 1. Металлическая модель танка. 2. Бетонная плита (X = 0,7 Вт-м-1-К-1). 3. Мешок с песком (X = 0,35 Вт-м-1-К-1). 4. Плита из красного полнотелого кирпича (X = 0,67 Вт-м-1-К-1). 5. Деревянная плита (X = 0,18 Вт-м-1-К-1). 6. Мраморная плита (X = 2,9 Вт-м-1К-1). 7. Гранитная плита (X = 3,5 Вт-м-1-К-1). 8. Стальной лист (X = 58 Вт-м-1-К-1). 9. Плита из пенопласта (X = 0,04 Вт-м-1-К-1); 10. Газосиликатные блоки (X = 0,2 Вт-м-1-К-1). 11. Асфальт (X = 0,7 Вт-м-1-К-1) (рис. 2).

В результате применения алгоритма корреляционной обработки кубоида ИК-изображений были получены сегментированные изображения «эталонных» материалов. Принимая во внимание зависимость порога принятия решения об обнаружении согласно критерию Неймана -

а) б)

Рис. 2. Расположение эталонных материалов на грунте: а - в видимом диапазоне длин волн; б - в ИК-диапазоне длин волн

Таблица 1. Результаты оценки необходимого количества ИК-сигнатур

Номер п/п Объекты, по которым производился расчет выборочного коэффициента корреляции Количество ИК-сигнатур, обеспечивающих устойчивое решение задачи сегментации

1 Плита из красного полнотелого кирпича 10

2 Стальной лист 10

3 Плита из пенопласта 12

4 Деревянная плита 14

5 Бетонная плита 9

Пирсона от заданной вероятности ложной тревоги, на едином сегментированном изображении в соответствии с заданным количеством «эталонов» можно выделить материалы с близкими к ним теплофизическими параметрами. Путем использования алгоритмов предварительной обработки изображений становится возможным получение четких контуров таких объектов. Однако не всегда есть возможность непрерывно и длительно производить съемку. Например, при ведении воздушной разведки в условиях работы противовоздушной обороны противника либо в условиях ограниченного запаса топлива. Возникает необходимость оценки количества ИК-сигнатур в ходе наблюдения одной и той же экспозиции объектов, обеспечивающих решение задачи сегментации изображения и последующее устойчивое решение коэффициентной обратной задачи теплопроводности. Результаты такой оценки представлены в табл. 1.

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

- наилучшее решение задачи сегментации достигается при построении кубоида ИК-изображения, содержащего 150 кадров с интервалом 10 мин в течение суток (рис. 3а);

- минимально возможным количеством ИК-изображений является три изображения, при этом интервал съемки не должен превышать 500 мин (рис. 3б);

а б

Рис. 3. Результат сегментации кубоида ИК-изображений: а - 150 кадров с интервалом 10 мин.; б - три кадра с интервалом 500 мин

Рис. 4. Экспериментальные и рассчитанные на ЭВМ графики суточного нагрева и остывания поверхности пенопласта

- оптимальным количеством обрабатываемых ИК-изображений является 14 изображений с интервалом получения каждого из них около 100 мин.

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

Построенное численное решение задачи (2) после восстановления функции источника тепла по данным метеорологических наблюдений погоды позволяет с приемлемой точностью рассчитывать термодинамические температуры в любой пространственно-временной точке кубоида ИК-изображений. Так, на рис. 4 представлен график зависимости изменения избыточных температур на поверхности эталонного материала «пенопласт», рассчитанных с помощью разработанного алгоритма и зарегистрированных в ходе суточного эксперимента.

Заключение

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

Работа выполнена при финансовой поддержке РФФИ (грант № 15-08-02611 А).

Список литературы

[1] Ishchuk, I.N., Parfiriev, A.V. The Reconstruction of a Cuboid of Infrared Images to Detect Hidden Objects. Part 1. A Solution Based on the Coefficient Inverse Problem of Heat Conduction (2014) Measurement Techniques, vol. 56, Issue 10, pp 1162-1166.

[2] Ishchuk, I.N., Parfiriev, A.V The Reconstruction of a Cuboid of Infrared Images to Detect Hidden Objects. Part 2. A Method and Apparatus for Remote Measurements of the Thermal Parameters of Isotropic Materials (2014) Measurement Techniques, vol. 57, Issue 1, pp 74-78.

[3] Громов Ю. Ю., Балюков А. М., Ищук И. Н., Ворсин И. В. Математическая модель автоматизированной системы испытаний инфракрасной заметности объектов в условиях неопределенности. Промышленные АСУ и кон- троллеры, 2014, 7, 12-15. [Gromov Y.Y., Biryukov A. M., Ischuk I.N., Vorsin I.V. A mathematical model of an automated test system infrared visibility of objects in condition of uncertainty. Industrial ACS and controllers, 2014, 7, 12-15. (in Russian)]

[4] Ищук И.Н, Обухов В.В., Парфирьев А.В., Филимонов А.М. Методика дистанционного контроля изотропных материалов путем редукции кубоида инфракрасных изображений Измерительная техника, 2015, 9, 41-45. [Ischuk I.N., Obukhov V.V., Parfir'ev A.V., Filimonov A.M. Method of remote control of isotropic materials by reduction cuboid infrared images. Measurement technology, 2015, 9, 41-45. (in Russian)]

[5] Дульнев Г.Н., Парфенов В.Г., Сигалов А.В. Применение ЭВМ для решения задач теплообмена. М.:Высш. шк., 1990. 207 с. [Dulnev G.N., Parfenov V.G., Sigalov A.V. // The use of computers for solving problems of heat transfer. Moscow: Graduate School, 1990. 207 (in Russian)]

[6] Ращиков В.И., Рошаль А.С. Численные методы решения физических задач: учеб. пособие. Техкнига, СПб.: Лань., 2005. 208 с. [Raschikov V.I., Roshal A.S. Z7 Numerical methods for solving physical problems: Textbook, SPb .: Lan., 2005. 208 (in Russian)]