способ верификации системы наземных измерений состояния атмосферы мегаполисов

Формула изобретения

Способ верификации системы наземных измерений состояния атмосферы мегаполисов включает дистанционное получение генерализованного спектрозонального снимка региона, содержащего контрольные промышленные площадки в виде цифровых значений спектральной яркости изображения I(x, y) в красной полосе видимого диапазона (570 670 нм), калибровку пикселей яркости в значениях ПДК по относительной яркости пикселей контрольных площадок, вычисление пространственного спектра Фурье F(l/m) прокалиброванной матрицы сигнала и интервала r (км) между дискретными измерениями на местности: , программный расчет (N) локальных максимумов на изображении и отождествление их с максимумами загрязнения атмосферы, определение среднегеометрического расстояния r (км) между точками забора проб, отсчитываемого от локальных максимумов,

и объема измерений системы V из соотношений V=S/r²,

где F_max - максимальная частота пространственного спектра Фурье прокалиброванной матрицы изображения;

S - площадь изображения мегаполиса, км ²;

- среднее расстояние между локальными максимумами на изображении.

Описание изобретения к патенту

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

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

В соответствии с положением «О создании единой государственной системы экологического мониторинга» (ЕГСЭМ), Минприроды РФ, региональные органы обязаны осуществлять плановый контроль и проведение измерений показателей объектов экологического контроля. Мегаполис содержит сотни площадных (промышленных) и линейных (транспортных) источников выбросов в атмосферу. Согласно общесоюзному нормативному документа [ОНД-86 «Методика расчета концентраций в атмосферном воздухе вредных веществ, содержащихся в выбросах предприятий», СССР, Гидрометеоиздат, Лен, 1987 г., стр.32, пункт 5.12] расчет концентрации выбросов площадных целей проводится в сетке, содержащей порядка 1500 узлов.

Объем измерений достигает нескольких тысяч проб при стоимости химического анализа одной пробы ~10 тыс.руб. и времени анализа 8-12 ч. Таким образом, существует методологическая проблема научного обоснования как объема измерений, так и местоположения забора проб по площади мегаполиса.

Известен способ оценки состояния атмосферы путем расчета суммарного индекса q . Обычно q рассчитывают для пяти составляющих, определяющих основной вклад в загрязнение атмосферы при условии, что q ₁>q₂>q₃>q₄>q ₅ [см., например, Ежегодник состояния загрязнения атмосферы в городах на территории России, под редакцией Э.Ю.Безугловой, ГГО им. А.И.Воейкова, Санкт-Петербург, 1994 1996 гг. - Аналог]:

где m_i - [мг/м³] - средняя за год концентрация i-го вещества в атмосфере, определяемая по локальным измерениям на местности;

СН_i - предельно допустимая санитарная норма концентрации i-го вещества в атмосферном воздухе, согласно ГОСТ;

j - показатель степени изоэффективности вредного вещества, равный 0,85; 1; 1,3; 1,5 для веществ IV, III, II и I классов опасности.

Недостатками известных аналогов являются:

- статистическая неустойчивость метода единичных локальных измерений на местности в контрольных точках как таковых;

- неопределенность выбора самих контрольных точек забора проб и существенное влияние на результат измерений розы ветров и случайных завихрений атмосферы в точках забора;

- неоперативность методов аналитической химии путем забора проб, составляющая от 6 до 12 ч по одной пробе, высокая стоимость химического анализа пробы;

- невозможность одновременного контроля всего множества ингредиентов (до 600) из-за неоперативности и высокой стоимости;

- невозможность документальной проверки субъективно составляемых протоколов по результатам локальных измерений.

Ближайшим аналогом к заявляемому техническому решению является «Способ оценки загрязнения атмосферы», патент RU № 2117286, 1998 г.

Способ ближайшего аналога включает дистанционное получение спектрозонального снимка региона, включающего контрольные промышленные площадки, в виде цифровых значений спектральной яркости I(х, y) изображений G, R видимого диапазона, размерностью |m×n| элементов каждая, с линейным разрешением каждого элемента больше фрактального участка подстилающей поверхности, поэлементную логическую сортировку пикселей в обоих матрицах в соответствии с алгоритмом, если R G, то R, если R<G, то R=R_max-|k|-G, где k - коэффициент корреляции хроматических коэффициентов r, g, получают результирующую матрицу тех же размеров, вычисляют числовые характеристики результирующей матрицы - математическое ожидание, дисперсию, огибающую пространственного спектра, автокорреляционную функцию, рассчитывают гистограмму распределения пикселей по яркости, осуществляют привязку полученного относительного закона распределения к абсолютным значениям индекса состояния атмосферы региона по его значениям и площади контрольных площадок.

К недостаткам ближайшего аналога следует отнести:

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

- неочевидность калибровки относительного закона распределения пикселей по яркости измерениями контрольных площадок.

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

Технический результат достигается тем, что способ верификации системы наземных измерений состояния атмосферы мегаполисов включает дистанционное получение генерализованного спектрозонального снимка региона, содержащего контрольные промышленные площадки в виде цифровых значений спектральной яркости изображения I(х, y) в красной полосе видимого диапазона (570 670 нм), калибровку пикселей яркости в значениях ПДК по относительной яркости пикселей контрольных площадок, вычисление пространственного спектра Фурье F(1/m) прокалиброванной матрицы сигнала и интервала r (км) между дискретными измерениями на местности: , программный расчет (N) локальных максимумов на изображении и отождествление их с максимумами загрязнения атмосферы, определение среднегеометрического расстояния r (км) между точками забора проб, отсчитываемого от локальных максимумов и объема измерений системы V из соотношений: , W=S/r²,

где F_max - максимальная частота пространственного спектра Фурье прокалиброванной матрицы изображения;

S - площадь изображения мегаполиса, км²;

- среднее расстояние между локальными максимумами на изображении.

Изобретение поясняется чертежами, где

фиг.1 - фрагмент исходного спектрозонального изображения мегаполиса;

фиг.2 - гистограмма распределения пикселей изображения по яркости;

фиг.3 - функция распределения загрязнения атмосферы по площади мегаполиса, прокалиброванная значениями ПДК контрольных площадок;

фиг.4 - огибающая пространственного спектра Фурье функции сигнала изображения;

фиг.5 - нахождение локальных максимумов в окне |3×3| элемента методами конечных разностей;

фиг.6 - функциональная схема устройства, реализующего способ.

Техническая сущность изобретения состоит в следующем.

Согласно аналогам и «Положению о создании ЕГСЭМ» наземные измерения концентрации вредных выбросов в атмосферу проводят на выходе источников. Приземная концентрация вредных веществ в атмосфере зависит от множества факторов: количества источников, объема выбросов, розы ветров, стратификации градиента температур, времени года, суток, рельефа, инфраструктуры региона, азимутов направлений улиц и транспортных потоков. Поэтому результат любого локального измерения на местности представляется случайной величиной. В пределах мегаполисов, над ложбинами, могут наблюдаться «висячие» смоги или хорошо «проветриваемые» участки. Из математики известно, что положение кривой полностью определяется точками максимума, минимума и точкой перегиба [см. Н.С.Пискунов «Дифференциальное и интегральное исчисления для ВТУЗов», том 1, 5-е изд., М.: Наука, 1964 г., стр.168-172, §9. Выпуклость и вогнутость кривой. Точки перегиба]. С другой стороны, согласно теореме отсчетов Котельникова-Шеннона, непрерывная функция, заданная на конечном интервале [0 R], однозначно определяется совокупностью дискретных отсчетов, проводимых через интервал , где F_max - максимальная частота спектра измеряемой функции [см., например, «Теоретические основы радиолокации», под ред. В.Е.Дулевича, учебник, М., Сов. Радио, 1964 г., стр.211]. Данные математические положения использованы при верификации системы наземных измерений. В качестве измеряемой функции в заявляемом способе используют космический снимок мегаполиса в виде коэффициента спектральной яркости функции I(х, y) системы «атмосфера-подстилающая поверхность».

Дистанционно регистрируемый уровень отраженной солнечной радиации является суммой отражений от подстилающей поверхности и атмосферной дымки. При замутненной атмосфере происходит «замывание» спектрально-энергетических и пространственных характеристик подстилающей поверхности. Существуют методы параметрического разделения эффектов подстилающей поверхности в результирующей яркости [см., например, «Итоговый отчет по исследованию параметров Атмосфера-Поверхность дистанционными методами», Эксперимент МКС-М-МКФ-6 на станции Салют-7, 1983 1985 гг., М.: ИКИ, АН СССР, стр.23-31]. Одним из методов параметрического разделения является генерализация изображения, т.е. использование снимков невысокого разрешения в красной полосе (570 670 нм) видимого диапазона.

Взаимодействие солнечного излучения с антропогенными частицами происходит на молекулярном уровне. При сталкивании фотонов светового потока с молекулами газов происходит передача квантов энергии (h *) молекулам, которые переходят в возбужденное состояние. При всех видах возможного взаимодействия светового потока с молекулами смогов над мегаполисами, как то: поглощение, рассеяние, флуоресцентное переизлучение - интегральный эффект состоит в смещении спектра видимого диапазона в его длинноволновую часть (красная область) [см., например, Р.Межерис, Лазерное дистанционное зондирование, перевод с англ., М.: Мир, 1987 г., стр.124, табл.3.4 Волновые числа комбинационного сдвига на длине волны 337,1 нм] Ниже представлены некоторые извлечения из данной Таблицы для некоторых «парниковых» молекул смогов.

Тип молекулы вещества	NO2	SO2	СO 2	NH3	С2Н 2	H2 S	CO	NO	Н2 O
Длина волны рассеянного излучения, нм	345,7	350,8	352,5	378,8	369	380,3	363,9	365,9	384,4
Абсолютная величина смещения , нм	8,6	13,7	15,4	42,7	43,2	32	16,8	18,8	47,3

В результате комбинационного рассеяния солнечного света происходит перераспределение энергии между спектральными составляющими видимого диапазона, а регистрируемое спектральное изображение антропогенно загрязненных участков приобретает преимущественно красноватый или темно-вишневый оттенок. Кроме того, в случаях мощных антропогенных дымок органического и фотохимического происхождения имеет место дополнительное интенсивное поглощение (красного и ближнего ИК диапазонов), сопровождаемое «замыванием» пикселей изображения подстилающей поверхности [см., например, «Коэффициенты спектральной яркости воздушной дымки», Л.И.Чапурский. «Отражательные свойства природных объектов в диапазоне 400 2500 нм», ч.1, Мин. Обороны СССР, 1986 г., стр.103-107, 134-137, табл.П.10].

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

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

Функция распределения значений ПДК по площади мегаполиса, прокалиброванная значениями ПДК контрольных площадок, иллюстрируется фиг.3. Прокалиброванную матрицу изображения мегаполиса |m×n| элементов подвергают программной обработке, для чего определяют скорость флуктуации искомой функции сигнала по пространственным координатам вычислением ее пространственного спектра.

По определению, двумерный пространственный спектр сигнала матрицы G(F_x, F_y), с учетом линейных размеров (m, n) анализируемого участка, рассчитывают по формуле:

Интегрируют двумерный пространственный спектр по кольцевым сегментам и получают одномерную огибающую пространственного спектра сигнала в функции полярного радиуса Ф(1/R). Данная операция реализуется на основе алгоритмов быстрого преобразования Фурье (БПФ) [см. Марпл С.А. «Цифровой спектральный анализ», перев. с англ., М.: Мир, 1990 г., стр.77-79]. Программы БПФ входят в комплект специализированного программного обеспечения [см., например, «Специализированное программное обеспечение MATH САД.6.0 PLVS», издание 2-е, стереотипное, М.: Информ. издат. дом «Филинъ», 1977 г., стр.411]. Расчетное значение огибающей пространственного спектра анализируемого изображения иллюстрируется фиг.4. Максимальная пространственная частота функции (фиг.4) соответствует . Интервал дискретизации функции сигнала по теореме Котельникова-Шеннона составляет . Следующей задачей анализа является нахождение локальных максимумов загрязнений атмосферы и отождествление их с точками забора проб наземной системы измерений. Поиск локальных максимумов двумерной функции сигнала цифровой матрицы реализуется следующим алгоритмом. Как известно, в точке экстремума производная функции равна нулю. Для нахождения производной, последовательно от начала массива, разбивают матрицу на окна размером [3×3] смежных элемента и рассчитывают конечные результаты по координатам x, y относительно центрального пиксела.

;

;

Процедура нахождения локального максимума в окне |3×3| элемента иллюстрируется фиг.5. Текст программы нахождения локальных максимумов приведен в примере реализации. После нахождения (N) точек локальных максимумов рассчитывают среднее расстояние между ними. Учитывая, что между локальными максимумами существуют точки перегиба (изменения знака производной), можно рассчитать среднюю площадь S/N, приходящуюся на один максимум, откуда среднее расстояние между максимумами . Поскольку расстояние между точками забора проб рассчитывают двумя методами ( r и b), то в качестве оптимального расстояния принимают среднегеометрическое .

Пример реализации способа

Заявляемый способ может быть реализован на базе устройства по схеме фиг.7. Функциональная схема устройства фиг.7 содержит орбитальный комплекс наблюдения 1, типа космического аппарата (КА) «Ресурс» с установленной на его борту цифровой фотокамерой высокого разрешения 2 (типа DCS760), для приема сигнала светового потока, отраженного от подстилающей поверхности 3 через светофильтр 4 в красной полосе видимого диапазона. Оперативную съемку запланированных объектов и включение фотокамеры 2 осуществляет бортовой комплекс управления (БКУ) 5 по командам, передаваемым из центра управления полетом (ЦУП) 6 по радиолинии командного управления 7. Отснятые кадры записывают в буферное запоминающее устройство 8 и в зонах видимости КА с наземных пунктов, по автономному высокоскоростному каналу передачи данных 9 сбрасывают на пункты приема информации (ППИ) 10, где осуществляют их запись на средствах 11. Информацию с ППИ перегоняют по наземным каналам связи 12 в экологический центр 13, где осуществляют выделение отснятых кадров по служебным признакам. Скомпонованные массивы изображений мегаполисов помещают в базу данных 14 сети Интернет 15 (сайт открытого доступа для пользователей «Ресурс»). Программную обработку изображений по операциям способа проводят на средствах ПЭВМ 16 в стандартном наборе элементов: процессор 17, винчестер 18, оперативное запоминающее устройство 19, дисплей 20, принтер 21, клавиатура 22. Программы вычисления локальных максимумов и спектра Фурье записывают на винчестер 18. В качестве примера обрабатывался спектрозональный снимок г.Мытищи, Московской обл. Снимок заказан в Госцентре «Природа», заказ № 11/93-42. Уровень квантования сигнала по амплитуде соответственно шкале (0 255). Фрагмент обрабатываемого изображения (фиг.1) соответствовал матрице |m×n|=|512×512| элементов.

Наибольшее загрязнение на изображении имеет атмосфера над территорией Мытищинского машиностроительного завода (ММЗ), выбранная в качестве эталонной промплощадки. Средний уровень ПДК по территории ММЗ на дату получения снимка составил 4,7 ПДК, I_min=92, I_max=248. Гистограмма распределения пикселей по яркости анализируемого изображения иллюстрируется фиг.2, где а) гистограмма абсолютных значений яркости, б) гистограмма распределения ПДК по площади мегаполиса. Минимальная яркость на изображении соответствует максимальным значениям ПДК (максимальное затухание сигнала). Для эталонной промплощадки (ПДК=4,7) относительное значение яркости пикселей соответствует: I_{пром.площ}/I_max 92/248=0,37, откуда калибровочный коэффициент обратно пропорциональной зависимости для любой точки относительной яркости гистограммы соответствует:

Для прокалиброванной (в значениях ПДК) матрицы изображения рассчитывался пространственный спектр. Огибающая пространственного спектра иллюстрируется фиг.4. Максимальное значение частоты пространственного спектра F_max=0,3·10 ^-3 1/м. Интервал пространственной дискретизации измерений искомой функции загрязнения атмосферы равен r=1/2 F_max=66 км.

Затем программным методом вычисляют количество (N) локальных максимумов анализируемой функции.

Текст программы определения количества локальных максимумов в матрице изображения

{Подсчет количества локальных максимумов}

program Extrem070206;

uses crt;

const maxx=500;

Kmem=100;

vaf S:real;

ssrstring;

Buf1, Buf2, Buf3:array[1 maxx] of real;

f, ff, T:text;

z:byte;

Pixelx, Pixely, I, j, k, kolx, koly, kolM,

zmax, zmin, kolExtr:integer;

mashtab, Lx, Ly, dx, dy, dh, dpx, dpy, Ax, Ay,

z1, mm, sigm, sigm2, SonTree, Distreal;

function Ex (i:integer):boolean;

var x:real;

begin

x:=Buf2(i);

Ex:=(x>=Buf1[i-1]) and (x>=Buf1[i]) and (x>=Buf1[i+1]) and

(x>=Buf2[i+1]) and (x>=Buf2[i+1])

and (x>=Buf3[i-1]) and (x>=Buf3[i]) and (x>=Buf1[i+1])

end;

begin (main)

clrscr;

readln (Pixelx); {Количество пикселей по x на снимке}

readln(Pixely); {Количество пикселей по у на снимке}

readln (Lx); {Размер снимка по х в cm}

readln (Ly); {Размер снимка по y в cm}

readln (mashtab); {Масштаб изображения в 1 см 550 м}

kolM:=0;

repeat

writeln ('BBeflMTe имя файла');

readln (ss);

assign (f, ss);

{SI-} reset(f); {SI+}

I:=IOresalt;

until j=0;

assign(ff, 'Rez.txt');

rewrite(ff);

assign (T, 'Maxim.txt');

rewrite (T);

koly:=0;

mm:=0; zmax:=0; zmin:=255;

while not eof (f) do

begin

koly:=koly+1;

kolx:=0;

while not eoln (f) do

begin

read (ff.z); kolx:=kolx+1; mm:=mm+z;

if z<zmin then zmin:=z;

if z>zmax then zmax:=z

end;

readln (f)

end; close (f);

mm:=mm/kolx/koly;

writeln ('Среднее значение z', mm: 12:7);

writeln ('Максимальное и минимальное значения', zmax:4, zmin:4);

writeln ('kolx=', kolx: 10);

writeln ('koly=', koly: 10);

writeln (ff, 'Среднее значение z', mm:12:7);

writeln (ff, 'Максимальное и минимальное значения', zmax:4, zmin:4);

writeln (ff, 'kolx=', kolx: 10);

writeln (ff, 'koly=', koly: 10);

if (kolx>maxx) or (koly>maxx) then

begin

writeln ('flaHHbie не умещаются в массиве'); halt

end;

dx:=Lx/(Pixelx-1) *mashtab;

dy:=Ly/(Pixely-1) *mashtab;

Ax:=Lx/(Picselx-1) *(kolx-1) *mashtab;

Ay:=Ly/(Picsely-1) *(koly-1) *mashtab;

reset (f);

sigm2:=0;

while not eof(f) do

begin

while not eoln(f) do

begin

read(f, z); sigm2:=sigm2+sqr(z-mm)/(kolx*koly)

end;

readln (f);

end;

close (f);

sigm:=sqrt(sigm2);

writeln(7tycnepcnfl', sigm:12:7);

writeln (ff, 'Дисперсия', sigm:12:7);

dh:=26.0/(2*sigm);

dh:=1;

reset (f);

kolExtr:=0;

for i:=1 to kolx do

begin

read (f, z); zl:=z*dh; Bufl(i):=zl;

end;

readln (f);

for i:=1 to kolx do

begin

read (f, z); zl:=z*dh; Bulf2(i):=z1;

end;

readin (f);

for i:=1 to kolx do

begin

read (f, z); zl:=z*dh; Bulf3(i):=zl;

end;

readln (f);

j:=2;

while not eof(f) do

begin

Buf1:=Buf2; Buf2:=Buf3;

for i:=1 to kolx do begin

read (f, z); zl:=z*dh; Buf3(i):=zl;

end;

readln (f);

for i:=1 to kolx-1 do

begin

if Ex(i) then begin

kolExtr:=kolExtr+1;

writeln (T, i*dx:15:2, **, j*dy:15:2);

end

end;

j:=j+1

end;

Son:=Ax*Ay/kolExtr;

close (f);

close (T);

writeln ('Имя файла-', ss);

writeln ('dx=', dx:10:5, 'dy=', dy:10:5, 'dh=', dh:10:5);

writeln ('Cтрок -', koly:5, 'Стобцов - ', kolx);

writeln ('Ширина области', Ax:15:7);

writeln ('Длина области', Ay:15:7);

writeln ('Проекция площади ', Ax*Ay:15:3);

writeln ('Количество локальных максимумов', kolExtr:10);

writeln ('Площадь под один максимум', SonTree:15:7);

writeln (ff, 'Имя файла -', ss);

writeln (ff, 'dx=', dx:10:5, 'dy=', dy:10:5, 'dh=', dh:10:5);

writeln (ff, 'Строк -', koly:5, 'Столбцов -', kolx);

writeln (ff, 'Ширина области', Ax:15:7);

writeln (ff, 'Длина области', Ay:15:7);

writeln (ff, 'Проекция площади', Ах*Ау:15:3);

writeln (ff, 'Количество локальных максимумов', kolExtr:10);

writeln (ff, 'Площадь под один максимум', Son:15:7);

close (ff);

end. (main)

С дискретизацией по амплитуде в один ПДК на обрабатываемом фрагменте выделено (N=3) локальных максимума. Масштаб спектрозонального снимка соответствовал: 1:55.000, площадь обрабатываемого фрагмента составила S=4 км. Выделены объекты локальных максимумов: ММЗ (ПДК=4,7), «Химволокно» (ПДК=3,5), «Стройпластмасс» (ПДК=2,5). Фоновый (минимальный) уровень ПДК по площади фрагмента 1,57, что соответствует данным наземных измерений. Кроме локальных максимумов, программа позволяет вычислить площадь под каждым локальным максимумом до точки перегиба. Среднее расстояние между локальными максимумами: Среднегеометрическое расстояние между точками забора проб: . Объем измерений загрязнений атмосферы анализируемого фрагмента: V=S/r²=4/1,93 2,1.

Согласно плану измерений СЭС г.Мытищи, измерения проводятся в 8 точках забора на выходе основных источников выбросов. Таким образом, избыточность информации, при установленном плановом контроле системы наземных измерений, достигает (8:2,1) 4 раз.

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