Космический мониторинг и математические модели лесных пожаров в районе аварий на газопроводах

Аварии на газопроводах часто сопровождаются взрывами, в результате которых могут возникать крупные лесные пожары, наносящие большой ущерб народному хозяйству и ОПС. Проблема обнаружения, регистрации и учета лесных пожаров (ЛП) имеет общепланетарное значение. ЛП являются значительным источником парниковых газов и аэрозолей и вносят существенный вклад в развитие парникового эффекта. Международная Конвенция 1997 года в Киото обязывает страны-участницы проводить оценки и всячески воздерживаться от увеличения выбросов парниковых газов в атмосферу, развивать системы национальных мониторингов, совершенствовать измерительные системы и наращивать банки данных по выбросам в атмосферу. Большинство мероприятий по контролю лесных пожаров проводится в рамках международного проекта по лесомониторингу GOFC (Global Observation of Forest Cover).

К настоящему времени разработаны методологические основы и ряд инструментальных средств для систем мониторинга лесных пожаров. Так, например, в соответствии с предложениями Европейского Космического Аген- ства (ESA - European Space Agency) для решения задачи обнаружения и регистрации послепожарных гарей может быть предложена система из шести скалярных индексов, называемых К-параметрами, суть которых состоит в следующем.

Известно, что неповрежденная растительность имеет большую отражательную способность в ближней инфракрасной области (NIR - Near Infra-Red, 0,7 - 1,3 мкм), по сравнению с другими видами подстилающей поверхности и с участками послепожарных гарей. Это приводит к тому, что участки гарей в этом спектральном канале будут менее яркими по сравнению с участками неповрежденного леса. Важной характеристикой послепожарных гарей является также большая вариация коэффициента отражательной способности.

Следующей важной физической характеристикой гарей является повышенная в дневные часы температура поверхности, что происходит за счет более интенсивного поглощения солнечного излучения частичками сажи и другими продуктами сгорания. Повышение температуры обусловлено также уменьшенной транспирацией растительного покрова на послепожарных гарях. Пепел и сажа образуют сухой слой на поверхности, который препятствует ее охлаждению. В результате температура гарей на 7-8 градусов больше, чем температура неповрежденных насаждений. Эти температурные изменения могут быть измерены с помощью "теплового" спектрального канала 11,0 мкм (TIR - Temperature Infra-Red, 10,4 - 12,5 мкм), который может быть проградуирован в единицах яркостной температуры поверхности.

Эти физические характеристики гарей используются ESA для получения главного интегрального индекса К1-параметра, для чего пиксели изображения наносятся на координатную плоскость с координатными осями NIR и TIR. Учитывая, что в NIR диапазоне пиксели смещаются влево, а в TIR диапазоне вверх, получаем в результате заметное смещение пикселей гарей по диагонали влево. Это смещение мало зависит, как от типа лесного покрова, так и от состояния атмосферы, что позволяет отделить кластер гарей от основного кластера неповрежденной территории. Для количественной оценки смещения в каналах NIR и TIR рассчитываются средние значения яркостей изображения - центральная точка, которая принимается в качестве начала отсчета. Параметр К1 определяется, как расстояние в координатах NIR-TIR от центральной точки до середины кластера пикселей гари. К1 -параметр является весьма значимой характеристикой, однако в предложениях ESA используется целая совокупность параметров.

Индексы К2 и КЗ рассчитываются аналогично индексу К1, но с использованием вместо NIR соответственно IR (Infra-Red, инфракрасный, 0,76 - 0,9 мкм) -канала и комбинации спектральных каналов, называемой глобальным индексом мониторинга ОПС и обозначаемой - GEMI (Global Environmental Monitoring Index).

Индексы K4, К5 и Кб определяются из анализа временных рядов космических снимков. Предварительно рассчитываются К4- и К5-изображения. К4- изображение представляет собой нормализованную разность между NIR- каналом и его средним значением для всего временного ряда. Аналогично в качестве К5-изображения используется нормализованная разность между TIR- каналом и его средним значением для того же временного ряда.

Индекс Кб характеризует яркостную неоднородность гарей. Для расчета индекса используется предварительное преобразование красного канала (RED, 0,63 - 0,69 мкм) для каждого изображения временного ряда в дисперсию яркости пикселей, которая рассчитывается по скользящему окну размером 3x3 пикселя. Далее в качестве Кб-изображения используется нормализованная разность между результатом обработки и его средним значением для всего временного ряда.

Окончательно К4, К5 и Кб индексы получаются в результате обработки К4-, К5- и Кб-изображений пороговыми алгоритмами, которые являются адаптивными алгоритмами, коэффициенты которых зависят от параметров съемки, освещенности (высота Солнца), атмосферных условий и пр. Уменьшение индекса К4 и одновременное увеличение индексов К5 и Кб является важным критерием. В качестве более чуствительных показателей могут быть предложены следующие комбинации К4/К5 и К4/К6, при переходе которых через заданные пороги может приниматься решение об отнесении пикселя к классу гарей.

Задачи мониторинга ЛП можно условно разделить на три большие группы: А - задачи обнаружения ЛП; Б - задачи обнаружения и анализа послепо- жарных гарей; В - задачи анализа текущего состояния ЛП и прогнозирования его распространения.

Для решения задач группы (А) оптимальным является использование информации со спутника NOAA. Для решения задач группы (Б) требуется информация с многоканальных сканеров в спектральных каналах: GR (Green - зеленый, 0,52 - 0,6 мкм), RED (красный, 0,63 - 0,69 мкм), NIR, IR и TIR. Для решения задач группы (В) могут быть использованы космические изображения дымовых шлейфов ЛП с ИСЗ "Ресурс".

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

где Eq и Е - прямая и прошедшая через шлейф солнечная радиация, Р- коэффициенты отражения, поглощения и пропускания соответственно для дымового шлейфа (к = 1) и подстилающей поверхности (к - 2), -

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

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

где К . радиация от эталонных участков на входе в фотографическую аппаратуру. Формула (2.3.2) позволяет провести нормировку формулы (2.3.1):

или:

где Р2 - эффективный СКО подстилающей поверхности под шлейфом.

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

где С(х,у) - планарная, a q(x,y,z) - объемная концентрации.

В простейшем приближении прямой пропорциональной зависимости:

где к - коэффициент пропорциональности.

Для коэффициента пропускания в линейном приближении:

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

Подстановка (2.3.5) и (2.3.6) в (2.3.7) приводит к формуле для яркостных характеристик изображения в функции планарной концентрации дымового шлейфа:

где X - длина волны солнечного излучения.

Предположим также, что имеют место следующие соотношения:

где а - коэффициент пропорциональности, Ь, Ь2 - яркости пикселей изображения. При этом (2.3.8) примет следующий вид:

Для оценки параметров уравнения (2.3.10) рассмотрим задачу атмосферного переноса дымового шлейфа. Математическому моделированию ЛП посвящено множество научных работ.

Под пиролизом лесных горючих материалов (ЛГМ) будем понимать горение ЛГМ, в процессе которого при высокой температуре происходят химические реакции окисления с образованием углеродной пыли и газообразных горючих и негорючих продуктов: СН4, Н2, СО, СО?, Н20. Наиболее интенсивное горение происходит в части ЛП, называемой фронтом пожара. На космических изображениях (КИ) фронт ЛП представляет собой группу пикселей максимальной яркости. Иногда над фронтом ЛП образуются кучевые облака вследствие конденсации паров воды, образовавшихся в результате горения ЛГМ. Эти кучевые облака могут наблюдаться на КИ.

Тепловыделение ЛП на единицу длины фронта за единицу времени равно:

где к - коэффициент полноты сгорания ЛГМ, q - теплотворная способность топлива, w - скорость фронта ЛП, т - запас ЛГМ на единицу площади.

Для количественного описания ЛП используют следующие ярусы: 1 - подстилка; 2 - мох, лишайники и опад; 3 - кустарнички высотой до 20см; 4 - травянистые растения; 5 - кустарники высотой до 2м; 6 - кроны подроста высотой до 6м; 7 - кроны верхнего яруса высотой до 25м; 8 - приземный слой атмосферы высотой до 200м; 9 - пограничный слой атмосферы высотой до 2 км.

Если при ЛП сгорают (1) - (4) слои, то ЛП называется напочвенным. Если при ЛП повреждаются (6) и (7) слои, то ЛП называется верховым. ЛП называется огненным штормом, если он распространяется в условиях засухи и ураганного ветра.

Для математического моделирования ЛП первостепенное значение имеет схема параметризации конкретного горящего лесного выдела. Такая параметризация должна быть рассчитана в результате обработки таксационных баз данных по стандартным методикам. К сожалению в настоящее время не существует автоматизированных информационных систем (АИС), позволяющих производить подобные расчеты исходя из таксационных описаний. Важное значение имеет также модель кинетики дымообразования. Суммарная масса частиц дыма, возникающих при лесном пожаре составляет 1-5% первоначальной массы ЛГМ. Дым ЛП состоит из частиц сажи и золы.

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

где р и р0 - плотности соответственно частиц дыма и ЛГМ, А, В, к, Е - параметры модели, Т - температура. Можно использовать также другую математическую модель кинетики дымообразования, в которой полагается, что мощность ЛП прямо пропорциональна скорости пиролиза ЛГМ:

где а - 0,01-0,05, р и т - плотности соответственно частиц дыма и масса Л ГМ. Модель (2.3.13) должна быть дополнена моделью скорости пиролиза Л ГМ:

где т, тс - массы ЛГМ и коксика, Е - энергия активации пиролиза, к, ас - коэффициенты пропорциональности. Уравнения (2.3.14) справедливы лишь в тех случаях, когда термокинетические постоянные процесса пиролиза одинаковы для всех типов ЛГМ в лесном биогеоценозе.

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

где п, т, р - число атомов углерода, водорода и кислорода в молекуле ЛГМ. Оксид углерода в результате химической реакции горения превращается в углекислый газ:

Доля оксидов углерода можно оценить из уравнения химической реакции (2.3.15):

Для чистой целлюлозы: п= 1, т=2, р= 1 vСОСОг =0,875. Для хвои сосны

vco,co2 = 0» 7.

Исходя из массы коксика, которую можно рассчитать по уравнению (2.3.14), можно оценить с помощью (2.3.17) суммарную массу выделившихся при пиролизе ЛГМ газообразных оксидов углерода, что представляет интерес для биосферных оценок.

На рис. 2.3.1 представлено сканерное изображение части территории Хабаровского края, на котором виден дымовой шлейф лесного пожара в районе аварии на газопроводе. Географические координаты снимка: 53,502° С.Ш., 137,235° В.Д. - 53,176° С.Ш., 137,815° В.Д. Территория принадлежит Тахтин- скому лесничеству Тахтинского лесхоза Хабаровского края. Для анализа ЛП может быть использован план лесонасаждений, который необходимо совместить со снимком.

Сканерное изображение лесного пожара в районе аварии

Рис. 2.3.1. Сканерное изображение лесного пожара в районе аварии

газопровода

На первом этапе моделирования на космическом изображении выделяют дымовой шлейф ЛП. При этом для выделения однородных по яркости объектов могут быть использованы, как пороговые алгоритмы, так и алгоритмы кластеризации. Как было отмечено выше, в настоящее время разработано достаточно большое количество алгоритмов кластеризации с использованием, как четких, так и нечетких методов обработки. В общем, эти алгоритмы можно разделить на три больших класса: эвристические, иерархически и разделительные. На рис.2.3.2 представлены результаты кластеризации космического изображения ЛП (рис.2.3.1) с помощью метода нечетких с-средних. Для кластеризации может быть использована как однокластерная, так и многокластерная модель дымового шлейфа. В приближении однокластерной модели дымового шлейфа кластер пожара содержит п - 1913 пикселей, а его площадь S - 7173,75 га. В рамках этой модели предполагается, что шлейф имеет однородную непрозрачную оптическую плотность, задаваемую при помощи порогового значения концентрации дымовых частиц. В приближении многокластерной модели требуется задание нескольких порогов соответственно для каждого кластера, которые могут быть выбраны эквидистантными.

Результат кластеризации космического изображения ЛП с помощью метода нечетких с-средних для

Рис.2.3.2. Результат кластеризации космического изображения ЛП с помощью метода нечетких с-средних для: а - 30 и б - 40 итераций

На следующем этапе анализа ЛП осуществляется моделирование дымового шлейфа. Рассмотрим ЛП в виде постоянно действующего источника дымовых частиц. Вследствие высокой температуры в зоне горения образуются мощные восходящие тепловые потоки, которые, увлекая за собой частицы пепла и золы, высоко поднимают их над лесом. Эта начальная высота подъема получила название высоты теплового подъема шлейфа. С этой высоты начинается обычный атмосферный перенос и рассеяние дымовых частиц ветровыми потоками. Краевая задача атмосферного переноса дымового шлейфа имеет следующий вид:

где Q - мощность ЛП, q - концентрация дымовых частиц, и - скорость ветра, Kv, К: - усредненные по высоте пограничного слоя коэффициенты турбулентной диффузии, vg - скорость оседания частиц, Я - высота теплового подъема, 5 - дельта-функция для моделирования точечного источника (приближение точечного источника допустимо, если начальные размеры дымового клуба много меньше размеров всего остального шлейфа). Задача (2.3.18) имеет аналитическое решение следующего вида:

где erfcQ - функция ошибок.

Более точные оценки можно получить, если в задаче (2.3.18) использовать полуэмпирические аппроксимации для профилей скорости ветра и коэффициентов турбулентной диффузии.

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

где q(x,y,z) - концентрация дымовых частиц, ? - удельный коэффициент поглощения, 1 - интенсивность излучения. Из (2.3.20) следует, что интенсивность монохроматического излучения уменьшается по экспоненциальному закону:

где т - коэффициент пропускания, Zmax - максимальная высота подъема дымового шлейфа. Коэффициент поглощения равен: р = 1 - т.

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

где С| и 02, d и d2 - соответственно поперечные сечения поглощения и диаметры частиц сажи и золы. Поперечные сечения поглощения рассчитывают по следующим формулам:

где cti и а2 - спектральные коэффициенты поглощения веществ, из которых состоят частицы сажи и золы. При проведении оценок можно положить для частиц сажи: а, = 1, dt = 5 i ei , а для частиц золы: а2 = 0,27, <7, = 2 / И .

Формулы (2.3.19) - (2.3.23) позволяют при наличии метеорологической информации провести весь необходимый расчет для получения планарного распределения относительного коэффициента оптического поглощения. При этом введение порога непрозрачности дымового шлейфа позволяет связать результаты моделирования атмосферного переноса и кластеризации. Примем порог непрозрачности на уровне коэффициента поглощения р = 75%.

Рассмотрим однопараметрическую обратную задачу атмосферного переноса дымового шлейфа. Для решения задачи требуется определить неизвестный параметр - мощность ЛП, для чего необходимо уравнять площадь шлейфа, полученную в результате расчетов по формулам (2.3.19) - (2.3.23), с площадью шлейфа, рассчитанной с помощью космического изображения. Из результатов кластеризации КИ следует, что площадь дымового шлейфа ЛП S = 7173,75 га. С другой стороны эта же площадь может быть выражена следующим образом:

[ 1, при х > О

где 0 - функция Хэвисайда: 9(х)= 1 „.

w [О, при х < О

Из уравнения (2.3.24) можно определить неизвестную величину мощности ЛП, а по формуле (2.3.17) оценку суммарной мощности ЛП по выбросам газообразного оксида углерода и углекислого газа. В результате моделирования могут быть получены следующие оценки: площадь дымового шлейфа S = 7,2-107 м2; расчетная скорость поступления сажи и золы в атмосферу (мощность ЛП) Qp = 332 кг/с; размеры зоны активного горения: = 3,4 га (180x180 м); расчетная скорость поступления в атмосферу углеродистых газов: Qg - 230 кг/с.

Описанная математическая модель может быть реализована с помощью инструментального пакета MathCad. Далее приведем распечатку программы решения обратной задачи о распространении дымового шлейфа ЛП с использованием инструментального пакета MathCad.

Программа решения обратной задачи о распространении дымового шлейфа (задача оценки мощности ЛП)

 
Посмотреть оригинал
< Пред   СОДЕРЖАНИЕ   ОРИГИНАЛ     След >