ALGORITHM OF RESTORATION OF THE GROUND SURFACE REFLECTION COEFFICIENT IN VISIBLE AND NEAR IR- RANGES USING MODIS SATELLITE DATA

Abstract

An algorithm of restoration of the ground surface reflection coefficient is proposed. The basis of the algorithm is an algorithm for atmospheric correction of images, based on a rigorous solution of the radiation transfer equation for an inhomogeneous Earth surface.

Keywords

• atmospheric correction
• Monte Carlo method
• spectral ground surface reflection coefficient

-

Тарасенков М.В.( 1 ), Белов В.В.( 1 )( 2 ), Кирнос И.В. ( 1 )( 2 ) Предлагается алгоритм восстановления коэффициента отражения земной поверхности. Основу алгоритма составляет алгоритм атмосферной коррекции изображений, основанный на строгом решении уравнения переноса излучения для неоднородной земной поверхности.

Ключевые слова: атмосферная коррекция, метод Монте-Карло, спектральный коэффициент отражения земной поверхности.

Введение. Проблема восстановления спектральных коэффициентов отражения земной поверхности по спутниковым измерениям является предметом многочисленных исследований на протяжении нескольких десятилетий (начиная с таких работ как [1] и заканчивая современными работами, такими как [2-3]). Эта проблема остается актуальной до сих пор. Получение качественной спутниковой информации для многих условий наблюдения невозможно без устранения искажающего влияния атмосферы на принимаемый сигнал - атмосферной коррекции. Все методы атмосферной коррекции можно разбить на две группы. В первую группу можно отнести методы, выполняющие атмосферную коррекцию, исходя из информации, содержащейся в исходном снимке без привлечения информации о состоянии атмосферы и решения уравнения переноса. Это, например, такие алгоритмы как [4]. Для выполнения атмосферной коррекции в этих алгоритмах применяются методы статистической обработки, привлекается информация о границах водоемов и др. Однако вопрос о погрешности и границах применимости таких алгоритмов остается открытым. К второй группе методов можно отнести методы, привлекающие для атмосферной коррекцию информации об оптических параметрах атмосферы и построенные на решении уравнения переноса излучения. К этому типу алгоритмов можно отнести алгоритмы [2-3,5]. Эти методы атмосферной коррекции, вследствие учета процесса переноса излучения, способны давать более качественную информацию об отражательных свойствах земной поверхности, но при условии качественной информации об оптических параметрах атмосферы. Помимо этого строгое решение уравнения переноса излучения требует значительных временных затрат. Для уменьшения машинного времени в большинстве алгоритмов вводятся упрощения в процесс переноса излучения. Так в алгоритме [5] коррекция на первом этапе выполняется в предположении об однородности земной поверхности для каждого отдельно взятого пикселя. Далее, исходя из полученных приближенных значений коэффициентов отражения, для каждого отдельно взятого пикселя значения корректируются с учетом бокового подсвета. Это значительно ускоряет расчет, но вносит дополнительную погрешность в результат коррекции. Эта погрешность незначительна для условий низкой мутности атмосферы, но для высокой мутности атмосферы эта погрешность может оказаться значительной. По этой причине в работах [6-7] нами разрабатывается алгоритм атмосферной коррекции, в котором с одной стороны выполняется более точное решение поставленной задачи, а с другой введено ряд приемов для ускорения процесса получения результата. Далее рассматривается предлагаемый нами алгоритм коррекции, и приводятся примеры сравнений с алгоритмом MOD09 [11].

Алгоритм атмосферной коррекции. Если рассматривать атмосферу как рассеивающую и поглощающую аэрозольно-газовую среду, а земную поверхность считать неоднородной ламбертовской поверхностью с неизвестным распределением коэффициента отражения, то интенсивность принимаемого спутниковой системой излучения определится как: Isum,i  Isun,i  rsurf,i Esum,i exp i    rsurf xw, yw Esumxw, yw hxw,i , yw,i , xw, yw dxwdyw  S где

Esum,i  E0  E0  rsurf xw , yw h1xw,i  xw , yw,i  yw dxwdyw 

S    E0   rsurf xw , yw h1xw  xw , yw  yw dxwdyw   rsurf xw , yw h1xw  xw , yw  yw dxwdyw  ...

S  S  где i – наблюдаемый пиксель, Isun – интенсивность излучения Солнца, невзаимодействовавшая с земной поверхностью, rsurf – коэффициент отражения земной поверхности, Esum – суммарная освещенность земной поверхности с учетом переотражений в системе атмосфера-земная поверхность, h – ФРТ канала формирования бокового подсвета, (xw,yw) – поверхностные координаты точки на земной поверхности, E0 – освещенность земной поверхности без учета переотражения излучения в системе атмосфера-земная поверхность, h1 – ФРТ канала формирования дополнительной освещенности переотраженным излучением, S – вся земная поверхность.

Если предположить, что в пределах пикселя изображения поверхность однородна, то ( 1 ) преобразуется к виду:

Isum,i  Isun,i  rsurf,i Esum,i exp i    rsurf, j Esum, j Aij  rsurf,i Esum,i Aout,i ,

N  j 1

Aij   hxw,i , yw,i , xw , yw dxwdyw ,

S j Aout,i   hxw,i , yw,i , xw, yw dxwdyw ,

N

S\ j1Sj где rsurf ,i - средний коэффициент отражения земной поверхности в окрестности i-го пикселя, полученный в однородном приближении, Esum,i - суммарная освещенность i-го пикселя в од

N нородном приближении при rsurf ,i  rsurf ,i ; S \  S j - область на земной поверхности вне облаj1 сти восстановления.

Если учитывать неоднородность отражения земной поверхности только для однократного переотражения, то в предположении об однородности поверхности в пределах пикселя ( 2 ) преобразуется к виду:

N Esum,i  E0  E0  rsurf, jCij  E0rsurf,iCout,i  E0  (rsurf,i 1)2  E0  (rsurf,i 1)3  ..., j1

Cij   h1xw,i  xw , yw,i  yw dxwdyw ,

S j Cout,i   h1xw,i  xw , yw,i  yw dxwdyw ,

N

S\ j1S j  1   h1xw  xw , yw  yw dxwdyw

S Тогда решение задачи сводится к решению системы линейных уравнений для поиска Q=rsurfEsum и нелинейных уравнений для поиска rsurf:

Q N Isum,i  Isun,i  i exp i    Ai, jQj  Aout,iQi ( 5 )  j1 ( 2 ) (3) (4) QE0i  rsurf,i 1 jN1 Ci, jrsurf, j  Cout,irsurf,i  1rsurrsf,uirf1,i21  (6) где Qi  rsurf,i Esum,i .

Для ускорения решения ( 5 )-(6) предлагается ряд приемов: 1) величина I sun,i рассчитывается не для всех пикселей по отдельности, а только для 35 узловых ситуаций и используется аппроксимационная формула, описанная в [6]; 2) функция h как видно из (3) зависит от положения наблюдаемой точки на земной поверхности, но приближенно можно разделить поверхность Земли на подобласти, где ее можно считать с заданной точностью постоянной и использовать одну функцию для всей области. Границы этих областей задается критерием, описанным в [6]; 3) расчет функций h и h1 входящих в (3) и (4) выполняется не по всей поверхности Земли S, а по поверхностям ограниченных радиусами R и R1 от центра наблюдаемого пикселя соответственно. Условия для задания радиусов описаны в [6]; 4) если на участке, где выполняется восстановление коэффициентов отражения, располагается плотная облачность, то предполагается, что для таких участков Qj  Qi и rsurf , j  rsurf ,i .

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

Сравнение результатов с алгоритмом MOD09 для безоблачных ситуаций. Для апробации предлагаемого алгоритма выполнялось сравнение с результатами, полученными алгоритмом MOD09 для участка Юга Томской области 55.95–56.850 с.ш. и 84.05-84.950 в.д. с 13 по 17 июля 2013 г. и участка пустыни Такла-Макан с координатами 38.55-39.449°с.ш. и 84.69185.396 °в.д. с 12 по 21 июля 2013 г. для 5 каналов MODIS (0.65, 0.47, 0.55, 1.24 и 0.41 мкм). Данные участки были выбраны потому, что доля облачных пикселей на рассматриваемых снимках составляла меньше 30%. Первый участок был выбран потому, что на его территории располагается станция Aeronet, что позволяет получать более качественную информацию об оптическом состоянии атмосферы. Второй участок был выбран по причине того, что это однородная территория без городов, рек, дорог и других неоднородных поверхностей.

Оптические параметры атмосферы для первого участка восстанавливались, исходя из данных Aeronet об аэрозольной оптической толщине (АОТ). Для снимков первого участка АОТ атмосферы при λ=0.55 мкм лежала в пределах 0.092-0.211. Исходя из значений АОТ, среди моделей LOWTRAN-7 [8] для лета средних широт выбиралась наиболее близкая по АОТ. Молекулярное рассеяние для первого участка определялось с использованием спутниковых данных о температуре и давлении и данных из работы [9]. Молекулярное поглощение бралось для модели средних широт LOWTRAN-7 и соответствующих длин волн. Пример восстановленных коэффициентов отражения, полученных предлагаемым алгоритмом, алгоритмом MOD09 и без коррекции приведен на рисунке 1. Из рисунка видно, что результаты двух алгоритмов коррекции полностью согласуются.

Для второго участка по причине отсутствия данных Aeronet об АОТ использовались спутниковые данные MODIS. АОТ для рассматриваемых снимков лежала в пределах 0.1-0.46. Аналогичным образом среди тропических моделей LOWTRAN-7 выбиралась наиболее близкая по АОТ, но в отличие от предыдущего случая была доступна только АОТ при λ=0.55 мкм, поэтому модели подбирались экстраполяцией. Молекулярное рассеяние и поглощение в силу отсутствия информации брались, исходя из модели LOWTRAN-7 для соответствующей длины волны. Как показывает сравнение результатов, предлагаемый алгоритм дает значения выше, чем алгоритм MOD09. Причиной, по всей видимости, является то, что используемая оптическая модель плохо соответствует рассматриваемым условиям. 56,8 56,7 56,6 56,5 ад56,4 гр56,3 56,2 56,1 56,0 56,8 56,7 56,6 56,5 ад56,4 гр56,3 56,2 56,1 56,0 84,1 84,2 84,3 84,4 84,5 84,6 84,7 84,8 84,9  , град а Алгоритм MOD09 84,1 84,2 84,3 84,4 8,4г,р5ад84,6 84,7 84,8 84,9 0,000,00 0rs,u0rf,61 0,12

в г Рис. 1. Сравнение коэффициентов отражения, полученных без коррекции (а), алгоритмом MOD09 (б) и предлагаемым алгоритмом (в); сопоставление результатов двух алгоритмов (г). Юг Томской области, дата измерений 14.07.2013 г. Длина волны λ=0.55 мкм.

Возникает вопрос об оценке погрешности восстановления коэффициентов отражения. Для рассматриваемых диапазонов дат можно считать, что от снимка к снимку не происходит значительных изменений коэффициента отражения и считать это изменение линейным. Тогда погрешность восстановления для безоблачных участков можно описывать соотношением [2]:  

N 1 1  i2 ti1  ti1

N 1 1  i2 ti1  ti1  * 2 rsurf,i  rsurf,i 

, rs*urf,i  ti  ti1 rsurf,i1  ti1  ti rsurf,i1 , ti1  ti1 (7) ( 8 ) где N – количество анализируемых снимков; ti – время i-го измерения; rs*urf,i – приближенная оценка коэффициента отражения по значениям днем раньше и позже; rsurf,i – коэффициент отражения земной поверхности, восстановленный алгоритмом коррекции для i -го дня.

Используя (7)-( 8 ), была выполнена оценка средних погрешностей восстановления коэффициентов отражения земной поверхности для двух рассматриваемых участков. На рисунке 2 приведены средние значения коэффициентов отражения для рассматриваемых участков и средние значения погрешностей для алгоритма MOD09 и предлагаемого алгоритма.

Сравнение средних значений коэффициентов отражения земной поверхности и оценок погрешностей для Юга Томской области показывает, что данные MOD09 и предлагаемого алгоритма полностью согласуются между собой. Отличие результатов находится в пределах погрешности расчетов. Погрешности расчетов двух алгоритмов также почти совпадают. 0,1 0,6 0,8 1,0 1,2 0,4 0,6 0,8 1,0 1,2 , мкм , мкм а б Рис. 2. Сравнение средних коэффициентов отражения rsurf и средних погрешностей σ для участка Томской области (а) и участка пустыни Такла-Макан (б).

Аналогичное сравнение для участка пустыни Такла-Макан показывает, что восстановленные предлагаемым алгоритмом значения коэффициенты отражения значительно отличаются от алгоритма MOD09. Сравнение с коэффициентами отражения песка из [10] показывает, что результаты MOD09 полностью согласуются с [10]. Причиной некачественной работы предлагаемого алгоритма является недостаток информации об оптическом состоянии атмосферы.

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

Работа выполнена при финансовой поддержке РФФИ (гранты №15-01-00783-А, №1507-06811-А, №16-31-00033-мол_а).

References

[1] Otterman J. , Fraser R.S. Adjacency effects on imaging by surface reflection and atmospheric scattering: cross radiance to zenith // Appl. Opt. 1979 . V. 18, N 16. P. 2852 - 2860 .

[2] Breon F.-M. , Vermote E. Correction of MODIS surface reflectance time series for BRDF effects // Remote Sensing of Environment. 2012 . V. 125. P. 1- 9 .

[5] Vermote E.F. , Vermeulen A . Atmospheric correction algorithm: spectral reflectances (MOD09) . Algorithm Theoretical Background document, version 4.0 . 1999 . 107 p.

[8] Kneizys F.X. , Shettle E.P. , Anderson G.P. , Abreu L.W. , Chetwynd J.H. , Selby J.E.A. , Clough S.A. , Gallery W.O. User Guide to LOWTRAN-7 . ARGL-TR- 86 - 0177 . ERP 1010. MA. Hansom AFB. 1988 . 137 p.

[9] Bucholtz A . Rayleigh-scattering calculations for the terrestrial atmosphere // Appl. Opt. 1995 . V. 34 , iss. 15. P. 2765 - 2773 .