Моделирование инфракрасного излучения лимба Земли

Введение

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

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

Поглощающими газами в земной атмо­сфере являются в основном полярные моле­кулы, так как с электромагнитным полем эф­фективно взаимодействуют только молекулы с постоянным или индуцированным дипольным моментом. При этом существуют спектральные области, в которых поглощение невелико, так на­зываемые окна прозрачности атмосферы: корот­коволновое в диапазоне длин волн 3,5...4,1 мкм и длинноволновое - 8.13 мкм (рис. 2).

Компоненты атмосферы можно разделить на постоянные, относительное содержание кото­рых постоянно вплоть до высот порядка 80 км, и переменные, содержание которых меняется в зависимости от высоты, температуры и геогра­фического положения. Из постоянных компо­нент углекислый газ CO₂ наиболее интенсивно поглощает ИК-излучение [4]. Главной перемен­ной компонентой атмосферы, поглощающей ИК- излучение, является водяной пар, количество которого быстро уменьшается с высотой. Дру­гая переменная компонента атмосферы - озон - на уровне моря присутствует в очень неболь­ших количествах. Количество озона возрастает с высотой (при достижении максимума на высо­те примерно 20-26 км) и вновь уменьшается на б0льших высотах.

Применение в инженерных расчетах упрощения, когда излучение системы поверх­ность Земли - атмосфера принимается рав­ным излучению абсолютно черного тела, для рассматриваемого типа приборов не является достаточным, поскольку имеет значительную спектральную неравномерность [5], которая связана в основном с большим количеством спектральных линий поглощения воды, угле­кислого газа, озона и др. При этом в общем случае также следует учитывать различия в интенсивности излучения в зависимости от ге­ографических точек наблюдения, условий ос­вещенности, сезонности, климатических и по­годных условий. На рис. 3 приведены данные для ночного времени суток, апрель 2012 г. [6].

Целью исследования является разработка математической модели и проведение расчетов излучения системы поверхность Земли - ат­мосфера в среднем и дальнем ИК-диапазонах электромагнитного спектра для формирования исходных данных в заданном угловом направ­лении на лимб Земли.

Описание математической модели

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

Возьмем луч направления s и рассмо­трим элемент среды в виде цилиндра единич­ного сечения, ось которого совпадает с направ­лением луча (рис. 4). Пусть луч пересекает перпендикулярные ему основания в точках M и M′, находящиеся друг от друга на расстоянии ds. Интенсивность излучения в этих точках рав­на I_λ  (M, s) и I_λ  (M', s) соответственно.

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

где β_λ - коэффициент ослабления;

J_λ - функция источника;

da - количество вещества в цилиндре.

Количество вещества da в общем случае может быть определено разными способами. Для того чтобы произведение β_λda являлось безразмерным, в описываемой математической модели используется массовый коэффициент ослабления β_λ. При этом величина а опреде­ляется как масса, приходящаяся на единицу площади:

a = ρs,

где ρ - плотность вещества.

Если бесконечно малое приращение мо­нохроматического оптического пути опреде­лить как

du_λ = -β_λ da,

то, подставив du_λ в выражение (1), получим уравнение переноса излучения в среде

Как было отмечено выше, процессы из­лучения и рассеяния с учетом введенных при­ближений являются линейными, поэтому ко­эффициент ослабления β_λ можно представить суммой коэффициентов поглощения к_х и рас­сеяния σ_λ:

β_λ = к_λ + σ_λ.

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

Положим, что рассматриваемый участок среды (атмосферы) находится в состоянии ло­кального термодинамического равновесия [7].

Интенсивность излучения такого участка за­висит только от длины волны и температуры. Тогда функция источника J_λ (M, r) будет за­даваться функцией Планка:

где B_λ (T) - спектральная плотность энерге­тической яркости в единичном телесном угле;

- постоянная Планка; c - скорость света в вакууме; k - постоянная Больцмана.

Введем плоскопараллельную модель ат­мосферы (рис. 5), для которой:

где z - вертикальная координата; θ - зенитный угол.

Представим уравнение (2) в следующем виде:

С учетом введенных приближений далее достаточно рассматривать только поле уходяще­го вверх излучения . Направления лучей s при этом образуют область направлений 

Если входящие в уравнение (3) перемен­ные задать в виде функции от z, то можно по­лучить обыкновенное дифференциальное уравнение первого порядка, решение которого будет иметь вид:

где C - постоянная интегрирования;

v , w - переменные интегрирования по оси z.

Для задания граничного условия для функ­ции  на уровне z = 0 учитывается собствен­ное тепловое излучение земной поверхности (без учета отраженного потока нисходящего те­плового излучения), имеющей температуру T₀:

Здесь ε_λ (θ) - коэффициент излучения земной поверхности.

Введем монохроматическую функцию пропускания τ_λ (θ, Z₁, z₂) , характеризующую долю излучения, пропущенную под углом θ к вертикали слоем атмосферы между уровнями z₁ и z₂:

где u_λ(z₁, z₂) - оптическая толщина при θ = 0.

При этом

где ρ( z) - средняя плотность поглощающего вещества в слое z₁ - z₂.

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

Для расчета интенсивности излучения для заданного положения наблюдателя, на­ходящегося за пределами лимба Земли, необ­ходимо учитывать сферичность атмосферы Земли (рис. 6):

R - радиус Земли;

h₀ - расстояние от поверхности Земли до рассматриваемого слоя атмосферы, отсчиты­ваемое по направлению в зенит.

Тогда функция пропускания примет сле­дующий вид:

где h₁, h₂ - нижняя и верхняя границы рас­сматриваемого слоя.

Дальнейшее упрощение выражения (4) состоит в переходе от интеграла к алгебраиче­ской сумме интенсивностей каждого из слоев стратифицированной атмосферы [9] и преоб­разованию подынтегрального выражения [10]. Тогда запишем конечную форму выражения для интенсивности в узком спектральном ин­тервале:

Из выражения (5) видно, что интенсив­ность ИК-излучения в верхнюю полусферу направлений складывается из интенсивности излучения поверхности Земли, ослабленного поглощением в атмосфере, и собственного из­лучения каждого слоя атмосферы, также осла­бленного поглощением вышележащих слоев.

Если линия наблюдения не пересекается с поверхностью Земли, то в формуле (5) ис­пользуется только правая часть, при этом слой n = 1 соответствует высоте h₀ (см. рис. 6).

Исходные данные для моделирования

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

В расчетах использованы 4 атмосферные модели: среднеширотная, арктические зима и лето, тропики. Исходные данные распре­делений по высоте атмосферы парциальных концентраций паров воды, озона, углекислого газа и других веществ, давления и темпера­туры взяты из модели RFM [11]. Данные по углекислому газу обновлены в связи с тем, что его концентрация в атмосфере по сравнению с 2001 г. выросла. Содержание углекислого газа на апрель 2017 г. составило 410 ppmV.

Значения спектральных коэффициентов поглощения атмосферных газов взяты из ин­формационной системы SPECTRA [12] для за­данных значений температуры, давления, типа уширения линий, спектрального диапазона с заданным разрешением по волновому числу.

Базы данных информационной системы SPECTRA на сегодняшний день являются наи­более полным и надежным признаваемым ми­ровым научным сообществом открытым источником информации о параметрах атмосферных газов. Эти параметры получены из решений уравнений квантовой механики и подтвержде­ны результатами многократных прямых изме­рений с помощью специального Фурье-спек- трометраMcMath-Pierce (Аризона, США).

Верификация математической модели и результаты расчетов

Существует ряд теоретических моделей по­верхность Земли - атмосфера, реализован­ных в виде специализированного программ­ного обеспечения [13]. Однако, некоторые из них, как, например, MODTRAN (Moderate Resolution Atmospheric Transmission) [14] военно-воздушных сил США или SciTran (Scientific Transparency, США и Герма­ния) [15], обладают высокими требования­ми к производительности вычислительных средств и при этом являются фактически не­доступными для предприятий оборонно-про­мышленного комплекса России. В других используются слишком грубые приближения в части учета спектральных параметров погло­щения или геометрической модели атмосферы как, например, CRTM (Community Radiative Transfer Model, США) [16], CRM (Column Radiation Model, США) [17], COART (Coupled Ocean-Atmospheric Radiative Transfer code, США) [18]. По этим причинам рассматривать на практике эти теоретические модели как источники данных для верификации разраба­тываемого программно-алгоритмического ап­парата в целом не представляется возможным.

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

Верификация модели проводилась в два этапа.

На первом этапе проводился расчет ве­личины спектрального пропускания атмосфе­ры по высотам в направлении лимба (рис. 7), результат сравнивался со значениями обработанных экспериментальных данных, полу­ченных специализированным космическим аппаратом MetOp (прибор IASI). При этом в качестве исходных данных о состоянии атмос­феры были приняты соответствующие условия космической съемки [19]. Отличия графиков в областях, выделенных красными и голубы­ми овалами (см. рис. 7), в основном связаны с присутствием в спектрах пропускания, по­лученных экспериментальным путем, спек­тральных линий дополнительных веществ, что легко устранить в разработанной математиче­ской модели, введя дополнительные исходные данные по содержанию веществ в атмосфере.

На втором этапе для заданного состоя­ния атмосферы проводился расчет величины спектральной плотности потока излучения от системы поверхность Земли - атмосфера в направлении в надир, который сравнивался с экс­периментальными данными со спутника MetOp (прибор IASI) [2] и спутника Nimbus 4 [20] (рис. 8, 9).

Хорошо видно, что излучение системы поверхность Земли - атмосфера в целом по­вторяет форму графика функции Планка (пун­ктирные линии на рис. 8, 9) при определенной температуре, однако для исследований опти­ко-электронных систем необходимо учесть полосы поглощения атмосферных газов.

При анализе графиков видно, что опи­сываемая модель не в полной мере учитывает все спектральные линии поглощения, которые присутствуют в экспериментальных данных, что ограничивает возможности ее практиче­ского применения при рассмотрении узких спектральных диапазонов (в диапазоне вол­новых чисел менее ~100 см^-1). Как было отмечено выше, данный недостаток модели устра­няется введением дополнительных исходных данных по содержанию веществ в атмосфере.

Полученные расчетным путем абсолют­ные значения спектральной плотности энер­гетической яркости незначительно отличают­ся от экспериментальных данных только для результатов со спутника Nimbus 4 (см. рис. 9), что объяснимо отсутствием точных исходных температурных данных по слоям атмосферы для данных географического региона и кли­матических условий.

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

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

Преимуществом данной разработки по сравнению с другими, заложенными в ресурсо­емкое и дорогостоящее программное обеспече­ние (см. например, [14, 21]), является возмож­ность использовать математическую модель в инженерных расчетах и программном обеспе­чении стендов полунатурного моделирования в режиме реального времени.