CC BY
19
Вычислительные технологии
Том 19, № 3, 2014
Алгоритм атмосферной коррекции спутниковых изображений неоднородной земной поверхности в видимом и УФ-диапазонах длин волн*
М.В. Тарасенков1, В. В. Белов1'2 1 Институт оптики атмосферы им. В.Е. Зуева СО РАН, Томск, Россия 2 Томский государственный университет, e-mail: tmv@iao.ru, belov@iao.ru
Тарасенков М.В., Белов В.В. Алгоритм атмосферной коррекции спутниковых изображений неоднородной земной поверхности в видимом и УФ-диапазонах длин волн // Вычислительные технологии. 2014. Т. 19, № 3. С. 48-56.
Рассматривается алгоритм атмосферной коррекции спутниковых изображений для восстановления коэффициентов отражения земной поверхности в видимом и УФ-диапазонах. Алгоритм объединяет учёт основных факторов, влияющих на формирование изображения, и ряд приемов, позволяющих значительно сократить время расчёта.
Ключевые слова: метод Монте-Карло, атмосферная коррекция, уравнение переноса излучения.
Tarasenkov M.V., Belov V.V. Atmospheric correction algorithm for satellite images of nonuniform Earth's surfase in the visible and UV ranges // Computational Technologies. 2014. Vol. 19, No. 3. P. 48-56.
When solving problems for atmospheric correction of satellite images of the Earth surface in the visible and UV spectral ranges, it is important to consider such factors as multiple scattering, absorption of radiation, adjacency effect, and re-reflection. One of the approaches to solve these problems is connected with the solution of the stationary radiation transfer equation. It has no general analytical solution. Therefore, some of the above-listed factors are disregarded in widely used algorithms of atmospheric correction.
In the present paper, an algorithm for correction of the aerospace images due the distorting effect of the atmosphere is suggested. It takes into account all main processes important for the formation of images (based on the solution of the radiative transfer equation by the Monte Carlo method) and involves a number of computer procedures which considerably reduce the computational time (by a factor of 6-10).
Particular example for processing of the actual image of a fragment of the western coast of Africa with coordinates 27.5-27.90 N and 13.4-12.90 W (measurements were performed on August 19, 2011 at 11h 25 min) was considered. The image was recorded by MODIS in the third spectral channel (0.469 m). The correction for the distorting effect of the atmosphere was carried out using the suggested algorithm along with the algorithm of uniform correction and the MOD09 algorithm developed by NASA.
* Работа выполнена при поддержке Госконтракта 14.515.11.0030, Президентской программы "Ведущих научные школы РФ" (грант № НШ-4714.2014.5), РФФИ (грант № 12-05-33068 мол_а_вед) и Интеграционного проекта СО РАН № 131а.
A comparison demonstrates that in the considered situation, the application of the algorithm of uniform correction for the distorting effect of the atmosphere is sufficient, and the correlation coefficient between the reflection coefficients of the Earth surface computed by the suggested algorithm and the MOD09 algorithm makes 0.984. Results of the reconstruction of the reflection coefficients using the suggested algorithm are on average by 0.023 higher than the results obtained with the help of the MOD09 algorithm. This can be explained by the difference in the optical models of the atmosphere employed to solve the radiation transfer equation. On this basis, it is possible to assert that the suggested algorithm yields reliable results and can be used to correct aerospace images for the distorting effect of the atmosphere in the visible range of wavelengths.
Keywords: Monte Carlo method, atmospheric correction, radiative transfer equation.
Введение
Данные пассивного спутникового мониторинга об отражательных свойствах земной поверхности имеют повсеместное применение в самых разнообразных областях экологии, природопользования и др. Основное преимуществом спутниковых данных — их все-охватность и оперативность. Однако использование спутниковых данных имеет свою специфику. Важным здесь является вопрос пригодности спутниковых данных для решения конкретных задач. Один из факторов, который необходимо учитывать при восстановлении отражательных свойств земной поверхности, — процесс взаимодействия оптического излучения с атмосферой как мутной средой, вызывающий изменение восходящих от поверхности световых потоков. Для устранения указанного влияния осуществляют атмосферную коррекцию изображений. К решению этой задачи в настоящее время существуют несколько подходов, однако каждый из них имеет свои ограничения. Такие методы, как описанный в работе [1], требуют наличия априорной информации об участках земной поверхности с резким перепадом отражательных свойств. Кроме того, вызывает вопросы достоверность этих методов при высокой мутности атмосферы. Методы, подобные представленному в [2], могут давать результаты с высокими погрешностями в видимом и УФ-диапазонах в условиях высокой мутности атмосферы и при резких перепадах значений коэффициента отражения. Предложенный в [3] и ему подобные методы характеризуются наиболее точными результатами, но требуют существенных временных затрат, что приводит к необходимости поиска средств их уменьшения. На основе ранее полученных результатов [4-6] авторами настоящей работы получен один из способов решения задачи создания точного и одновременно достаточно быстрого метода атмосферной коррекции.
1. Алгоритм расчёта
Задача решается в следующей постановке (рис. 1). На высоте hd от сферической земной поверхности располагается пассивная спутниковая система, ориентированная в направлении Wd и ведущая наблюдение за участком земной поверхности. Земная поверхность ламбертова с неизвестным распределением коэффициента отражения. Пространственное разрешение оптического приёмника в пределах наблюдаемой области считается постоянным. На верхнюю границу атмосферы падает поток солнечного излучения в направлении Wsun. Требуется, зная оптические параметры атмосферы и значения интенсивности, измеренные спутниковой системой, восстановить коэффициент отражения Tsurf.
Рис. 1. Геометрическая схема постановки задачи
Решение задачи строится следующим образом [7]. Суммарная интенсивность принимаемого спутниковой системой излучения Isum состоит из интенсивностей солнечной дымки Isun — излучения Солнца, рассеянного в атмосфере и не отражённого земной поверхностью, нерассеянного излучения, отражённого наблюдаемым участком земной поверхности I0, и поверхностной дымки Isurf — рассеянного излучения, отражённого от земной поверхности. Если считать, что в пределах пикселя земная поверхность однородна, компоненты излучения меняются незначительно, и при учёте только дополнительной освещённости земной поверхности первой кратности для поиска rsurf необходимо решить нелинейную систему уравнений [7]. Решение системы делится на два этапа. На первом определяется величина Q = rsurf Esum из системы уравнений вида
где ^ — косинус угла между направлением на приёмную систему и вертикалью в наблюдаемой точке поверхности; т — оптическая длина трассы; Нг — функция размытия точки (ФРТ) канала формирования бокового подсвета при наблюдении 1-й точки; И^г — интеграл по поверхности к-го пикселя от ФРТ канала формирования бокового подсвета при наблюдении ¿-го пикселя. На втором этапе определяется коэффициент отражения из нелинейной системы вида
<
(1)
Qi
Eo
N
rsurf,ll 1 + ^ ] rsurf,kHk,l
k=1
N
Qn = r I 1 + V^ r H (1)
= rsurf,N | 1 + rsurf,kHk
k=1
H
(1) l,i
hi(xw xw,i} yw Vw,i)dxwdyw
Si
где rsurfi — искомое значение коэффициента отражения в i-м пикселе изображения; h1 — ФРТ канала формирования дополнительной освещённости; Hj]) — интеграл по поверхности 1-го пикселя от ФРТ канала формирования дополнительной освещённости.
Если необходим учёт больших кратностей переотражения, то его предлагается выполнять в однородном приближении. В этом случае коэффициент отражения определится по формуле
' surf,i (1 + ' surf, iYi)
rsurf,i
1 + r surf,iYl(1 + r surfiYl)
5)
где г^ — значение коэффициента отражения, полученное с учётом бесконечного числа кратностей переотражения; г8игЦ — результаты решения системы (3); 71 — доля дополнительной освещённости %-й кратности от дополнительной освещённости (% — 1)-й кратности при rsurf = 1.
Реализация данного подхода требует существенных затрат машинного времени. Для ускорения расчёта предлагается учитывать ряд следующих факторов.
1. Критерий изопланарности изображений.
Изображение можно разбить на зоны, в пределах каждой из которых ФРТ можно считать постоянной. Для их определения в [6] предложен критерий выделения изопла-нарных зон вида
dd>l+i = arcco/1 - ((moo(0°) - ) | exp(A)) ^
moo(ed,i) = moo(0°) - exp(A)(1 - cosed,i)N,
тоо(в^) = ^ ехр(—п) + у у Н(х'У],у'т; 6^)^Лу'^,
где А, N — аппроксимационные константы, определяемые методом наименьших квадратов по узловым значениям интеграла ФРТ канала формирования бокового подсвета, рассчитываемым с помощью алгоритма, описанного в [4]. 2. Радиус бокового подсвета.
Функция к во многих случаях является быстроубывающей, поэтому область интегрирования в (4) имеет смысл ограничить радиусом бокового подсвета Rsurf (радиус понимается в поверхностных координатах). Для задания Rsurf предлагается условие вида
кг(К8иг1,г) > 6 + (6 — 1)^ ехр(—'Тг)/П, (7)
где
ki(Rsurf,i; xw ■) yw )
J hi(xw xw ,yw yw )dxw dyw
S(Rsurf )
Hi (xw ■) yw )
Иг(хш,уш) J J Ь<1(хт х,ш,у,ш Уш)Лхи]<Лу,ш, (9)
здесь кг, Яяиг/,г, тг, Иг, ^г соответствуют 1-й изозоне; 8 задаёт точность определения Q (использовалось 8 = 0.99); Б(Rsurf) — область земной поверхности, ограниченная Rsurf, Б — вся земная поверхность.
3. Радиус области формирования переотражения.
В силу быстрого убывания функции Н1 область интегрирования в (4) можно ограничить радиусом формирования переотражения Rl. Для оценки Rl предлагается использовать условие
81 ( 81 тА1 - 71
где
I I Л-1(Хш — хш ,уш — уш )$х,ш
k,№) > -Ü- - 1), (10)
V 1 - Yi /
kl(Rl) = -----, (11)
/ / hi(xw — xw,yw — yw)dxwdyw
S
Yi — требуемая точность определения коэффициента отражения. 4. Формула для интенсивности солнечной дымки
В работе [5] для описания зависимости Isun от угла ориентации спутниковой системы предлагается формула вида
Isun(ed,^d)=- B+2fl- lACi3, (12)
A = C11 cos2 9d + C¿i(sin 9d cos <d)2 + Ci2 cos 9d sin 9d cos < - (sin 9d sin <d)2, (13)
B = C12 cos 9d + Ci3 sin dd cos <d, (14)
где 0d,<d — углы ориентации оси оптической системы; C11, C12, C13, C21, C22, C23, C31, C32, C33 — константы, получаемые путём аппроксимации узловых расчётов интенсивности солнечной дымки методом Монте-Карло для узловых направлений dd = 0,15,... , 60°, <d = 0, 30,... , 180° (35 узлов). Константы C21, C22, C23 используются при 0 < <d < 90°, а C31, C32, C33 — при 90 < < < 180°.
2. Сравнение алгоритмов по эффективности
Чтобы показать эффективность предлагаемого алгоритма, который далее будем называть базовым, по сравнению с упрощёнными алгоритмами, рассмотрим тестовую задачу восстановления коэффициента отражения участка земной поверхности, представленного на рис. 2. Для расчётов была выбрана длина волны Л = 0.469 мкм, что соответствует центру третьего спектрального канала прибора MODIS (спутники TERRA и AQUA).
Рис. 2. Взаимное расположение спутника и наблюдаемой области и распределение коэффициента отражения по земной поверхности
1.0-,
0.5
Ъ С
0.2
±3
А
10.8 11.0
Долгота, град
1.0
0.8 -
0.2
5
л ■ ¿а.
1
5а
11.2
10.8 11.0
Долгота, град
11.2
Рис. 3. Восстановленные различными способами значения коэффициента отражения для тестового примера (см. рис. 2) вдоль линии р = 0.45°: а — Л = 0.469 мкм, Бм = 1 км, б — Л = 0.469 мкм, Бм = 50 км
а
Результаты данных расчётов выбраны в качестве иллюстраций, поскольку они позволяют наглядно убедиться в необходимости учёта влияния процессов бокового подсвета и переотражения при решении задач восстановления отражательной способности наблюдаемой поверхности.
Результаты расчётов коэффициента отражения земной поверхности вдоль линии р = 0.45°, полученные с помощью базового алгоритма (символы 1) и производных от него, основанных на учёте переотражения в однородном приближении (символы 2) и на приближении однородности поверхности (символы 3), представлены на рис. 3 для различных метеоусловий. Проведённое сравнение показывает, что базовый алгоритм даёт результат с абсолютной погрешностью в пределах 0.044, восстановление с учётом переотражения в однородном приближении — 0.077, а атмосферная коррекция в однородном приближении при учёте переотражения и бокового подсвета — 0.308. Из этого следует, что в ряде случаев использование базового алгоритма позволяет определить
коэффициент отражения с заметно большей точностью, в сравнении с алгоритмами, учитывающими приближённо отдельные факторы влияния неоднородности поверхности на принимаемое излучение. Для получения результатов (см. рис. 3) базовому алгоритму требуется 90-100 мин работы ЭВМ производительностью 19.5 ГФлопс, алгоритму в однородном приближении при учёте переотражения — 40-60 мин, а алгоритму в однородном приближении — порядка 15 мин. Без использования предложенных приёмов время счёта составляет более чем в 6 раз большую величину.
3. Сравнение результатов расчётов алгоритмами на примере тестового участка земной поверхности
Следующим этапом тестирования стало сравнение результатов расчётов базовым и стандартным алгоритмами для данных М0В1Б алгоритма М0В09 с результатами расчёта алгоритмом однородной коррекции для тестового участка земной поверхности в ви-
27.90-
27.60-
-13.3 -13.2 -13.1
Долгота, град
surf, 1
surf, 1
Рис. 4. Коэффициенты отражения рассматриваемого участка земной поверхности, полученные базовым алгоритмом (rsurf,l) (а), алгоритмом МОБОЭ (гзиг/,2) (б), алгоритмом однородной коррекции (rsurf,3) (в); сравнение значений коэффициента отражения, вычисленных базовым алгоритмом rsurf,l, со значениями, вычисленными алгоритмом МОБОЭ rsurf,2 (слева) и алгоритмом однородной коррекции rsurf,3 (справа) (г)
в
г
де области западного побережья Африки с координатами 27.5-27.9° с.ш., 13.4-12.9° з.д., время измерений 19 августа 2011 г. в 11 ч 25 мин, выбранной в силу того, что примерно половину её занимает пустынная территория, а вторую половину — море. В качестве исходных данных для предлагаемого алгоритма и алгоритма однородной коррекции использовались данные MOD02 и MOD04 о распределении интенсивности в третьем канале прибора MODIS (0.47 мкм) с пространственным разрешением 500 м и об аэрозольной оптической толщине, которая в данном случае почти не менялась и составляла 0.307. Оптическая модель атмосферы для базового алгоритма и алгоритма однородной коррекции выбиралась из континентальных тропических моделей безоблачного неба, рассчитываемых генератором LOWTRAN-7 [8]. Наиболее близкой по аэрозольной оптической толщине оказалась модель, соответствующая метеорологической дальности видимости Sm = 24 км. Сравнение результатов (рис. 4) показывает, что коэффициент корреляции значений, вычисленных базовым алгоритмом и алгоритмом MOD09, равен 0.984 (рис. 4, г, левый график), а вычисленных базовым алгоритмом и алгоритмом однородной коррекции, — 0.9985 (рис. 4, г, правый график). На основе значений, полученных методом наименьших квадратов, были построены прямые, приведенные на рис. 4, г пунктиром. Из рисунка видно, что значения, найденные алгоритмом MOD09, в среднем меньше значений, вычисленных базовым алгоритмом, на величину 0.023. Это может быть связано с тем, что в указанных алгоритмах используются различные оптические модели. Что касается базового и однородного алгоритмов, то сравнение показывает, что в данном случае можно использовать алгоритм однородной коррекции.
Таким образом, разработан новый алгоритм восстановления коэффициентов отражения земной поверхности. Показано, что существуют ситуации, когда применение данного алгоритма необходимо. Сравнительный анализ результатов расчёта базовым алгоритмом и алгоритмом MOD09 для тестового участка показал, что они заметно различаются, но коэффициент корреляции, равный 0.984, косвенно указывает на то, что причина различий может быть связана с использованием разных оптических моделей атмосферы.
Авторы благодарят Н.В. Кабанову и Д.В. Соломатова за предоставленные спутниковые данные для тестового участка на земной поверхности.
Список литературы
[1] Протасов К.Т., Бусыгин Л.А., Белов В.В. Метод преобразования гистограмм яркостей и вейвлет-коррекция атмосферных искажений спутниковых изображений // Оптика атмосферы и океана. 2010. Т. 23, № 2. С. 136-142.
[2] Vermote E.F., Vermeülen A. Atmospheric Correction Algorithm: Spectral Reflectances (MOD09). Algorithm Theoretical Background Document: version 4.0, 1999.
[3] Reinersman P.N., Carder K.L. Monte Carlo simulation of the atmospheric point-spread function with an application to correction for the adjacency effect // Appl. Optics. 1995. Vol. 34, No. 21. P. 4453-4471.
[4] Белов В.В., Тарасенков М.В. Статистическое моделирование интенсивности световых потоков, отражённых сферической земной поверхностью // Оптика атмосферы и океана. 2010. Т. 23, № 1. С. 14-20.
[5] Белов В.В., Тарасенков М.В., Пискунов К.П. Параметрическая модель солнечной дымки в видимой и УФ-области спектра // Там же. 2010. Т. 23, № 4. С. 294-297.
[6] Белов В.В., Тарасенков М.В. Статистическое моделирование функции размытия точки в сферической атмосфере и критерий выделения зон изопланарности изображений // Там же. 2010. Т. 23, № 5. С. 371-377.
[7] Белов В.В., Тарасенков М.В. О точности и быстродействии RTM-алгоритмов атмосферной коррекции спутниковых изображений в видимом и УФ- диапазоне // Там же. 2013. Т. 26, № 7. С. 564-571.
[8] Kneizys F.X., Shettle E.P., Anderson G.P. et al. User Guide to LOWTRAN-7. ARGL-TR-86-0177. ERP 2010. Hansom AFB. MA 01731. 137 p.
Поступила в 'редакцию 28 февраля 2014 г.
