Восстановление альбедо поверхности по данным высокого пространственного разрешения Текст научной статьи по специальности «Физика»

Восстановление альбедо поверхности по данным высокого пространственного разрешения

О.В. Николаева1 1 Институт прикладной математики РАН, Москва, Россия

Аннотация

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

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

Представлены результаты тестирования точности всех формул при решении задач атмосферной коррекции данных высокого (30 м) пространственного разрешения. Тестирование выполнено на задачах со слоем аэрозоля оптической толщины от 0,2 до 2 для длины волны X = 0,55 мкм и всевозможными значениями альбедо в целевом и окружающем пикселях (от 0,1 до 0,9).

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

Ключевые слова: альбедо поверхности, атмосферная коррекция, высокое пространственное разрешение.

Цитирование: Николаева, О.В. Восстановление альбедо поверхности по данным высокого пространственного разрешения / О.В. Николаева // Компьютерная оптика. - 2022. - Т. 46, № 3. -С. 406-414. - DOI: 10.18287/2412-6179-C0-1046.

Citation: Nikolaeva OV. Surface albedo retrieval based on high spatial resolution data. Computer Optics 2022; 46(3): 406-414. DOI: 10.18287/2412-6179-C0-1046.

Введение

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

Эта формула используется для атмосферной коррекции данных дистанционного зондирования, полученных в низком пространственном разрешении [2, 3], или же данных, зарегистрированных над пространственно-однородным океаном вдали от береговой линии [4].

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

При выполнении атмосферной коррекции по данным высокого пространственного разрешения требуется учитывать, что в коэффициенте яркости, зарегистрированном над целевым пикселем поверхности, присутствует вклад излучения, отраженного соседними пикселями поверхности (так называемый боковой подсвет). Для учета этого эффекта используется модифицированная формула, в которой два из трех альбедо в целевом пикселе заменены на среднее альбедо в окружающих пикселях [6 - 9]. В этой полуэмпирической формуле альбедо в целевом пикселе используется только в том слагаемом, которое описывает прямой переход отраженного от поверхности излучения на верхнюю границу атмосферы. Среднее аль-

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

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

Для детального моделирования переотражения света между различными пикселями поверхности было предложено использовать систему алгебраических уравнений [10]. Однако решение данной системы уравнений требует большого времени.

Поэтому в [11] была предложена формула, связывающая альбедо земных пикселей с коэффициентами яркости пикселей на верхней границе атмосферы. Формула содержит коэффициенты отражения и пропускания для атмосферы с черной поверхностью и с черно-белыми поверхностями. Черно-белой называется поверхность, отражающая все упавшее на нее излучение по закону Ламберта в одном пикселе и поглощающая все излучение в остальных пикселях. Задача о переносе излучения в атмосфере с такой поверхностью является трехмерной, даже если атмосфера представляется горизонтально-однородной.

Формула позволяет явным образом найти как альбедо по значениям коэффициентов яркости, так и коэффициенты яркости по значениям альбедо. Кроме того, в предложенной формуле не задается форма пикселей и не требуется, чтобы все они имели одну и ту же площадь. В частности, для нахождения альбедо в целевом пикселе можно использовать систему вложенных окружающих пикселей [12].

Была предложена аналогичная формула, содержащая коэффициенты отражения и пропускания для задач с черной и черно-красными поверхностями [13]. Черно-красной называется поверхность, изотропно излучающая в одном пикселе и поглощающая все излучение в остальных пикселях. Показано, что восстановление альбедо поверхности по данным низкого пространственного разрешения (250 м, что отвечает прибору MODIS) с помощью «черно-белой» и «черно-красной» формул имеет высокую точность.

В настоящей работе исследуется точность «черно-белой», «черно-красной» и полуэмпирической [6 - 9] формул при определении коэффициента яркости по значениям альбедо и при восстановлении альбедо по значениям коэффициентов яркости на модельных данных высокого пространственного разрешения (30 м, что отвечает прибору HYPERION). «Черно-белая» и «черно-красная» формулы рассматриваются в приближении двух пикселей, целевого и окружающего, для корректности сравнения их точности с точностью полуэмпирической формулы.

1. Задача о переносе излучения в атмосфере

Рассмотрим уравнение переноса монохроматического излучения в трехмерной области

GI = LI - SI = 0, (1)

- dI dI 9I ,, ... , ...

LI = Ц—+^—+ц—+к (х, y, z) I (x, y, z, ц, ф), (2) oz ox oy

SI=k,(r)J%(r, у,)I(r,Ц,Ф')dЦйф, (3)

Q

| = sin9cosф, "n = sin9sinф, ^ = cos8,

-X<x<X, -Y<y<Y, 0<z<H, -1<ц<1, 0<ф<2л .

Здесь решение I (r, ц, ф) - интенсивность излучения в пространственной точке r = (х, y, z) в направлении Q (ц = cos 8, ф); углы 8 и ф задают единичный вектор Q в сферической системе координат; ц -координаты вектора Q в прямоугольной системе. Коэффициент экстинкции k (r), коэффициент рассеяния k, (r) и индикатриса рассеяния % (r, у,) зависят как от высоты z, так и от горизонтальных координат x, y. Индикатриса зависит от косинуса угла между направлениями (ц, ф) и (ц', ф'), величины

у, (ц, ц', ф, ф')=|V 1-ц2 .^Д-ц'2 cos (ф-ф') .

Оператор L описывает движение излучения в среде и его поглощение средой. Оператор S описывает переизлучение средой поглощенного ею излучения, то есть процесс рассеяния. Оператор G = L - S описывает взаимодействие излучения со средой.

На верхней границе z = 0, -X < х <X, - Y<y < Y зададим условие параллельного потока излучения, падающего в направлении Q (м-0 = cos 80, ф0)

I(х,y,0, ц,ф) = f (ц,ф) приц< 0. (4)

Здесь / (ц, ф) = I0 5 (ц - ^0) 5 (ф - ф0), величина I0 -интенсивность внеатмосферного излучения, 5 - дельта-функция Дирака.

На нижней границе z = H, -X< х < X, - Y < y < Y зададим отражение по закону Ламберта с альбедо a (х, y)

I(х,y,H,ц,ф) = a(х,y)ЭТ приц < 0, (5)

где оператор отражения определяется равенством

1 2л

ЭТ = (1/л)J dц' ц'J dф' I(х, y, H, ц', ф'). (6)

0 0

На боковых границах зададим нулевые краевые условия

I(-X, y,z,!, ф) = 0 при ^>0, - Y < y < Y, 0 < z < H ,(7)

I(X,y,z, ц,ф) = 0 при ^<0, -Y < y < Y, 0 < z < H, (8)

I(х,-Y,z,ц,ф) = 0 при 0,-X<х<X, 0<z <H, (9)

I^Y,^!,ф) = 0 при^<0, -X<х<X, 0<z<H.(10)

Задача (1 - 10) имеет единственное решение.

2. Задача атмосферной коррекции

Нижнюю границу разобьем на два пикселя: целевой и, = {- хо < х < х0, -у0 <у <уо} и окружающий

ио ={- X < х< X,-У < у <У }\{-х0 < х< х0, - у0 < у < у0}.

Также выделим целевой и окружающий пиксели на верхней границе.

Введем средние значения альбедо по пикселям нижней границы

а/ = Р/а(х,у) = Цйхйу а(х,у)/|П'/1,/ = ,,о, (11)

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

Rj = ZjI = Pj I (X, y, 0, Д, p) л/ (д 0 Io), j = i, o.

(12)

Сформулируем задачу атмосферной коррекции: по средним значениям коэффициентов отраженной яркости Я/, / = ,,о, найти среднее значение альбедо поверхности в целевом пикселе а,-.

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

Полоса съемки

Целевой пиксель Окружающий пиксель

Рис. 1. Целевой и окружающий пиксели на модельном снимке

3. Основная формула для коэффициента яркости

Введем базовые задачи.

Задача с черной поверхностью

01ь = 0,1Ь(х,у,0, д,р) = /(д,р) при д>0, (13)

1ь(х,у, Н, д, р) = 0 при д < 0.

Задачи с черно-белыми поверхностями

01] = 0, 1] (х, у,0, д, р) = /(д, р) при д > 0, (14)

17 (х, у, Н, д, р) = А/ЭТ17 при д < 0, А/ = 1 при (х,у)е V/ -А/ = 0 при (х, у) £ V/.

Решение задачи (1 - 10) представим в виде суперпозиции решений базовых задач (13) и (14)

I = uIW + vIw0 + (1 - u- v) Ib

(15)

где и, V - некоторые константы.

Эта функция подчиняется уравнению (1) и краевым условиям (4), (7 - 10). Подставим представление (15) в краевое условие (5), (6) и применим к получившемуся соотношению оператор Р/, см. (11). Получаем систему из двух уравнений

/TW = aj\uTw + vTW + (1 -u -v)Tb], t Too = aB [u ToW + vToW + (1 - u - v) Tb ].

(16) (17)

Величины / = ЭТ 1] , Ть = Р/ ЭТ1Ь суть коэффициенты пропускания для решений базовых задач. Равенство (16) обеспечивает выполнение краевого условия (5) в целевом пикселе, равенство (17) - в окружающем.

Далее применим к представлению (15) оператор 2//, см. равенство (12); получим соотношения

Rj = uRW + vRW + (1 -u -v)Rb .

(18)

Здесь величины Я^ 1 = 1], ЯЬь = 1Ь суть коэффициенты отражения для решений базовых задач.

Найдем из уравнений (16), (17) коэффициенты и, V и подставим их в равенство (18) при / = ,. Получим выражение для коэффициента яркости в целевом пикселе Я, через значения альбедо и функционалов решения базовых задач а, и ао

R = Rb +-

aj C + ao C2 + a ao C3

(19)

D0 - at D1 - ao D2 + at ao D3

C = Tb ToW (RW -Rb), C2 =Tb TW R -Rb), (20)

C3 = C3 + C3', (21)

C3 =(RW -Rb) \Tb (ToW-Tb )-Tob (TW -Tb)],

C3'=(RW - Rb) \Tob (TW -Tb) - Tb (ToW - Tb) ],

D0 =TW0TW, D = ToW (TW -Tb), (22)

D2 =TW (ToW -Tb), (23)

D3 = (TW-Tb) (ToW -Tob) - (ToW -Tb )(TW -Tb). (24)

В другом подходе вместо базовых задач с черно-белыми поверхностями (14) используются задачи с черно-красными поверхностями

GIj = 0, I)(x,y,0, д,p) = 0 при д>0, 1г (x,y, H, д, p) = Aj при д< 0.

(25)

Если в задачах (14) один пиксель поверхности отражает падающее на него излучение по закону Ламберта, то в задачах (25) пиксель не отражает излучение, а сам испускает постоянный поток излучения независимо от падающего на него излучения. Если в задачах (25) исходящее от поверхности излучение пространственно-однородно в пределах каждого пикселя, то в задачах (14) - нет.

Решение задачи (1 - 10) представим в виде суперпозиции решений базовых задач (13) и (25).

I = Ib + йП + vir

(26)

и

где и , V - некоторые константы. Применим к этому равенству оператор . Получим

R. = üRjj + vRjB + Rb.

(27)

Здесь величины = IО суть коэффициенты отражения для решения базовых задач (25). Подставим представление (26) в краевое условие (5), (6) и применим к получившемуся соотношению оператор Р . Получаем систему из двух уравнений

й=a¡ [üTij + vT¡j + Tj ], v = aB \üTj + vTO + Tj ].

(28) (29)

Здесь величины ] = Р ^ 10 суть коэффициенты пропускания для решений базовых задач (25).

Найдем из уравнений (28), (29) коэффициенты и , V и подставим их в равенство (27) при] = ,. Получим для коэффициента яркости в целевом пикселе выражение вида (19). Коэффициенты С1, С2, Сз, Бо, А, А, Б3 в этом случае имеют вид

С =ТьЯ0, С2 =Т0ьЩ0, (30)

Сз = щ [Т4 ТО - Ть ТО ]+яо [Ть ТГ - ТЬ Т- ], (31)

Бо =1, б=Т0, А =ТОО, А=Т?Т;0 -ТТ. (32)

Таким образом, равенство (19) выражает значение коэффициента яркости Я, в целевом пикселе через значения альбедо поверхности а, аО и коэффициентов отражения и пропускания базовых задач, определяемых равенствами

RW,t = Zj I?, Rrj,i = Zj Ij, Rj = Zj Ib, Tw, = Pj Ü I?, j = P ÜIj, Tb = p Ü Ib.

(33)

(34)

В зависимости от выбранного приближения используются коэффициенты (20 - 24) или (30 - 32).

Известная полуэмпирическая формула может быть получена из равенств (19), (30 - 32). Предполагаем атмосферу пространственно-однородной и имеем Tob = Tj = Tb, Rj = Rj = Rb. Далее в формуле (19) пренебрегаем величинами второго порядка малости ai aB. Получаем

Ri = Rb + Tb

ai Rj + ao RjB 1_ a.Tj - a Tj '

í <Ui ±ii t^B j- BB

(35)

Коэффициент ЯГ находим, пренебрегая рассеянием излучения в атмосфере:

н

ЯГ « ехр(- т/Д), где т = | й2к(0,0,2) -

о

оптическая толщина слоя атмосферы. Здесь учитывается, что площадь целевого пикселя и, много меньше площади окружающего пикселя иО. В знаменателе формулы (35) пренебрегаем величиной а1Т? (многократным отражением от целевого пикселя). Окончательно, имеем полуэмпирическую формулу

R = Rb + ai ехр(-т/Д) + aBRjB Tb ' 1_a Tj

í <ub j- bb

(36)

Здесь ТОГО можно рассматривать как сферическое альбедо атмосферы, Я[О - как диффузную часть коэффициента пропускания от поверхности к верхней границе атмосферы. Фактически, здесь многократное отражение излучения от поверхности учитывается в одномерной геометрии. Формула (36) может быть записана в форме (19) при следующих коэффициентах

Q = ехр(-т/Д)Tb, С2 = TbRjo, С3 = 0, Do =1, Di = 0, D2 =TojB, D3 = 0.

(37)

(38)

4. Основная формула для альбедо поверхности

Не менее, чем формула для коэффициента яркости, важна формула для альбедо поверхности в целевом пикселе а,. Реально эта формула позволяет выполнять атмосферную коррекцию. Величину ai в «черно-белом» приближении находим из формулы (16), используя коэффициенты и, V, найденные из уравнений (18) при] = i, о. Имеем

a' = ■

(r¡ - rb) e1 - (rb - rb) e2

(Ri - Rb) E3 - (Rb - Rb) E4 + TbE,

T7 _T^pw T7 _T^DW T?b\

E1 = Tii (Roo - RB Л E2 = Tii \Ríb - Ri Л

E3 = (T?-Tb )(RB - Rb) - (TW-Tb)(RW - Rb), E4=(TW -Tb)(R? - Rb) - (TW-Tb )(RW - Rb),

T7 _/ Dw J?b\fr>w Db \ í T>w J?b\ÍJ?w J?b \

E5 =( R¡i - R )( Rbb - Rb )-(Roi - Rb )(R¡B - R¡ ) .

(39)

(40)

(41)

(42)

(43)

Выражение для величины аг в «черно-красном» приближении находится из формулы (28) с использованием коэффициентов и , V, найденных из уравнений (27) при] = г, о. Имеем равенство (39) со следующими коэффициентами Е\, Е2, Е3, Е4, Е5

Е\ = яОО , Е2 = ЯО , Е3 =ТЩо -ТьЩг, (44)

Т7 _ Г>Г 17 _ Т>гт>т Г>Г Г>Г /'/(^Л

Е4 = ЯюТц - Яц1ю , Е5 = яЦяОО -ЯoiЯio . (45)

Чтобы найти величину аг с использованием полуэмпирической формулы (36), найдем сначала альбедо аО в окружающем пикселе. Для этого используем уравнения (29) и (27) при ] = о и и = 0 ; исключая из этих равенств величину V, получаем

яО - Яь

(Яо - Яь )ТОО + ТьЯ0о

Подставляя это выражение в равенство (36), приходим к соотношению (39) со следующими коэффициентами

Е = ЯОО , Е2 = ЯО, Е3 = 0, (46)

Е4 =- ехр(-Т Д) ТО, Е5 = ехр(-Т Д) Я0. (47)

Таким образом, для нахождения альбедо поверхности в целевом пикселе а,- нужно выполнить следующие действия:

1) решить базовые задачи (13), (14) или (13), (25);

2) найти функционалы их решений (33), (34);

3) найти коэффициенты (40 - 43), или (44 - 45),

или (46 - 47);

4) найти искомое альбедо по формуле (39).

5. Численные результаты

Численное исследование точности полученных формул (19) для коэффициента яркости в целевом пикселе на верхней границе атмосферы и формулы (39) для альбедо поверхности в целевом пикселе выполним для модельного снимка с параметрами прибора HYPERION [14] (пространственное разрешение d = 30 м, ширина полосы съемки D = 7680 м).

Пусть на полосе съемки выбран квадратный фрагмент ширины D. Центральный пиксель этого фрагмента назначим целевым, тогда параметры х0 = y0 = d / 2 = 15 м. Остальные пиксели фрагмента объединим в окружающий пиксель; тогда параметры X = Y = D / 2 = 3840 м, см. рис. 1. Такой выбор окружающего пикселя позволяет максимально учесть влияние альбедо поверхности на значение коэффициента яркости в целевом пикселе.

Рассмотрим набор альбедо в целевом и в окружающем пикселях ai, ao = 0,1; 0,3; 0,5; 0,7; 0,9.

Атмосферу будем считать горизонтально-однородной высотой 100 км. Используем модель аэрозоля, построенную для Белоруссии и сопредельных регионов [15]. Коэффициенты экстинкции k, рассеяния ks и индикатрису рассеяния % (r, ys) найдем по теории Ми [16] для длины волны X = 0,55 мкм.

Используем оптические толщины аэрозоля т = 0,2; 0,4; 0,6; 0,8; 1; 1,2; 1,4; 1,6; 1,8; 2 (такой же интервал значений используется при тестировании в [17]). Зададим направление падения солнечного излучения углами 60 = 40°, ф0 = 0.

Введем по переменной х пространственную сетку из 58 ячеек (из них в двух ячейках помещается целевой пиксель). Такую же сетку введем по переменной y. Тогда целевому пикселю на поверхности отвечает 4 ячейки сетки. По высоте z введем пространственную сетку из 10 ячеек.

Введем по углам 6, ф равномерную квадратуру Карлсона ES12 из 168 узлов [18]. Рассматриваем направление отражения излучения в зенит ( Д = -1, ф = 0) и включаем в квадратуру соответствующий узел с нулевым весом.

Решение основной задачи (1 - 10) и базовых задач (13), (14), (25) находим сеточным методом дискретных ординат [19] в трехмерной (х, y, z) геометрии. Значение альбедо поверхности в каждой ячейке сетки на нижней границе z = 0 задаем в зависимости от того, в какой пиксель, целевой или окружающий, попадает эта ячейка.

После выполнения расчетов и получения сеточных величин в ячейках сетки на верхней границе z = H области находим:

1) точные значения коэффициентов яркости в целевом Я, (а,, а0, т) и окружающем Яо (а,, ао, т) пикселях по формуле (12);

2) коэффициенты отражения и пропускания по формулам (33), (34);

3) коэффициенты С1, С2, С3, Б0, А, Б2, 03 по формулам (20 - 24), или (30 - 32), или (37 - 38);

4) приближенные значения коэффициента яркости в целевом пикселе в «черно-красном» приближении ЯЬг (а,, а0, т), в «черно-белом» приближении Я'7 (а,, а0, т) и по полуэмпирической формуле Яет (а,, а0, т) с помощью равенства (19);

5) коэффициенты Е\, Е2, Е3, £4, Е5 по формулам (40 - 43), или (44 - 45), или (46 - 47);

6) приближенные значения альбедо поверхности в целевом пикселе в «черно-красном» приближении аЬг (а,, а0, т), в «черно-белом» приближении аь7 (а ¡, а0, т) и по полуэмпирической формуле а"т (а,, а0, т) с помощью равенства (39).

Далее находим относительные погрешности (%) определения коэффициента яркости и альбедо поверхности в целевом пикселе по формулам

еЬЯг (а,, а0, т) = 100% (1 - ЯЬГ (а,, а0, т)/Я, (а,, а0, т)), еЯ7 (а,, а0, т)=100% (1-Я'7 (а,, а0, т)/Я (а,, а0, т)), еТ (а1, а0, т)=100% (1-Я'т (а,, а0, т)/Я (а1, а0, т)), еьаг (аI, а0, т) = 100% (1 - аЬг (а,, а0, т)/а,), е'7 (а1, а0, т)=100% (1 - аЬ7 (а1, а0, т)/а1), ееат (аI, а0, т)=100% (1 - ает (аI, а0, т)/а,).

Находим минимальные, максимальные и средние по оптической толщине т значения погрешностей. Эти величины для всех рассмотренных значений альбедо а,, ао приведены на рис. 2 - 7 (средние значения обозначены штриховой линией, область всех возможных значений закрашена серым). Погрешности определения альбедо для самого яркого окружающего пикселя ао = 0,9 для разных оптических толщин аэрозоля т приведены на рис. 8 - 10. Погрешности определения альбедо для однородной поверхности (а, = ао) для всех значений т приведены на рис. 11.

Представленные результаты позволяют сделать следующие выводы.

1) Погрешность определения коэффициента яркости Я слабо зависит от оптической толщины т, см. рис. 2 - 4.

2) Большая зависимость от т наблюдается для погрешности определения альбедо а , полученного в «черно-красном» приближении и по полуэмпирической формуле, см. рис. 5, 7.

3) Зависимость от т погрешности определения альбедо а, в «черно-белом» приближении значительно меньше, см. рис. 6.

4) Модуль погрешности определения коэффициента яркости Я, в «черно-красном» при-

ближении растет с ростом альбедо целевого и окружающего пикселей, см. рис. 2, и может достигать 30 %.

5) Модуль погрешности определения коэффициента яркости Я, по полуэмпирической формуле растет с ростом альбедо окружающего пикселя и может достигать 16 %, см. рис. 3.

6) Коэффициент яркости Я,, полученный в «черно-красном» приближении и по полуэмпирической формуле, завышен в сравнении с точным значением, см. рис. 2, 3.

7) В «черно-белом» приближении коэффициент яркости Я, находится с погрешностью менее 0,3 %, см. рис. 4.

8) Модуль погрешности определения альбедо а, в «черно-красном» приближении в целом растет с ростом альбедо как целевого, так и окружающего пикселя, см. рис. 5, и может достигать 50 %.

9) Модуль погрешности определения альбедо а, по полуэмпирической формуле в целом растет с ростом альбедо как целевого, так и окружающего пикселя, и может достигать 140 %, см. рис. 7.

10) В «черно-белом» приближении альбедо а, находится с погрешностью менее 6 %, см. рис. 5.

11) Наибольшая по модулю погрешность вычисления альбедо получается в случае темного целевого пикселя в ярком окружении, см. рис. 5 - 7.

12) Модули погрешностей определения альбедо а, в ярком окружении (ао = 0,9) растут с ростом оптической толщины т и контраста ао - а,, см. рис. 8 - 10.

13) Для однородных поверхностей (ао = а,) погрешности определения альбедо малы (< 10 %) только для «черно-белого» приближения, а для «черно-красного» приближения и полуэмпирической формулы - только при малых оптических толщинах (т < 0,5), см. рис. 11.

0 —

=0,1

-10 -

-20 -

-30

0,1 0,5 0,9 0,1 0,5 0,9 0,1 0,5 0,9 0,1 0,5 0,9 0,1 0,5 0,9

Рис. 2. Погрешности определения коэффициента яркости в целевом пикселе Я1 в «черно-красном» приближении

-12 -

-16 -

0,1 0,5 0,9 0,1 0,5 0,9 0,1 0,5 0,9 0,1 0,5 0,9 0,1 0,5 0,9 а,

Рис. 3. Погрешности определения коэффициента яркости в целевом пикселе Я1 по полуэмпирической формуле

0.3 -1

0.2 -

0.1 -

0 0,1 0,5 0,9 0,1 0,5 0,9 0,1 0,5 0,9 0,1 0,5 0,9 0,1 0,5 0,9 а,

Рис. 4. Погрешности определения коэффициента яркости в целевом пикселе Я1 в «черно-белом» приближении

20 -1

10 -

-10 -

-20 -

-30 -

-40 -

-50 -

0,1 0,5 0,9 0,1 0,5 0,9 0,1 0,5 0,9 0,1 0,5 0,9 0,1 0,5 0,9 а,

Рис. 5. Погрешности определения альбедо в целевом пикселе а, в «черно-красном» приближении

0

ао=0,1 ао=0,3 ао=0,5 ао=0,7 ао=0,9

ао=0,7 ао=0,9

ао=0,3

ао=0,5

а

0

а

а =0,1

0 — 0

-4 -

0,1 0,5 0,9 0,1 0,5 0,9 0,1 0,5 0,9 0,1 0,5 0,9 0,1 0,5 0,9 а,

Рис. 6. Погрешности определения альбедо в целевом пикселе а, в «черно-белом» приближении

г=0,3 ао=0,5

0 --20 ■ -40 --60 -80 --100 -120 --140 -

0,1 0,5 0,9 0,1 0,5 0,9 0,1 0,5 0,9 0,1 0,5 0,9 0,1 0,5 0,9 а,

Рис. 7. Погрешности определения альбедо в целевом пикселе а, по полуэмпирической формуле

20 ■

-10 •

-30 •

-50 •

\

\

\

\

\

а,=0,3 а,=0,5

\

\

\

\

\

I 1 I 1 I 1 I 1 I 1 I 1 I 1 I 1 I

0,2 0,4 0,6 0,8 1 1,2 1,4 1,6 1,8 2 т

Рис. 8. Погрешности определения альбедо в целевом пикселе а, в «черно-красном» приближении для ао = 0,9 и разных оптических толщинах аэрозоля т

-40 ■ -60 •

-120 -

-140 -

\

\

\

\

\

а,=0,1

\

\

а,=0,5

\

\

• а,=0,9

-I—I—I—I—I—I—г

1 1 Г

I 1 I 1 I 1 I

0,2 0,4 0,6 0,8

1 1,2 1,4 1,6 1,8 2

Рис. 9. Погрешности определения альбедо в целевом пикселе а, по полуэмпирической формуле для ао = 0,9 и разных оптических толщинах аэрозоля т

I 1 I 1 I 1 I 1 I 1 I 1 I 1 I

0,2 0,4 0,6 0,8 1 1,2 1,4 1,6 1,8 2 Т

Рис. 10. Погрешности определения альбедо в целевом пикселе а, в «черно-белом» приближении для ао = 0,9 и разных оптических толщинах аэрозоля т

20 _ а,=ао=0,1. а,=ао=0,3 а,=ао=0,5 а,=ао=0,7 а,=ао=0,9

10 -

0 ■ -10 ■ -20 ■ -30 ■ -40 ■ -50 ■ -60 ■ -70 ■ -80 ■ -90 ■ -100 ■

"Черно-красное" приближение "" "Черно-белое" приближение Полуэмпирическая формула

0,2 1 2 0,2 1 2 0,2 1 2 0,2 1 2 0,2 1 2 т

Рис. 11. Погрешности определения альбедо в целевом пикселе а, для однородной поверхности (при а, = ао) и разных оптических толщинах аэрозоля т

ао=0,9

ао=0,3

ао=0,5

ао=0,7

0

а =0,3

а==0,7

т

ао=0,9

ао=0,1

ао=0,7

10

0 —

а=0,1

а =0,7

а=0,9

Заключение

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

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

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

Широко используемая полуэмпирическая формула для коэффициента яркости в целевом пикселе получена из формулы, найденной в «черно-красном» приближении. Предложена явная полуэмпирическая формула для альбедо поверхности.

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

Приведены результаты исследования точности представленных формул при определении как коэффициента яркости, так и альбедо поверхности. Тестирование выполнено на имитированных данных прибора высокого пространственного разрешения HYPERION (ширина полосы съемки 7680 м и пространственное разрешение 30 м) для оптических толщин аэрозоля т из интервала (0,2; 2) для длины волны X = 0,55 мкм.

Получено, что в «черно-белом» приближении погрешность определения коэффициента яркости не превышает 0,3 %, а альбедо поверхности -6 %. В «черно-красном» приближении погрешности определения коэффициента яркости и альбедо растут с ростом альбедо целевого и окружающего пикселей и достигают 30 % и 50 % соответственно. При исполь-

зовании полуэмпирических формул максимальные погрешности достигают значений 16 % и 140 %.

Во всех случаях наибольшая погрешность восстановления альбедо достигается при темном целевом пикселе в ярком окружении при наибольшей оптической толщине аэрозоля.

При однородной поверхности (одинаковых средних альбедо целевого и окружающего пикселей) погрешность определения альбедо для всех толщин т мала только в «черно-белом» приближении; в остальных случаях эта погрешность мала (меньше 10 %) только при малых значениях т (т < 0,5).

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

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

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

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

References

[1] Kaufman YJ, Tanre D, Gordon HR, Nakajima T, Lenoble J, Frouin R, Grassl H, Herman BM, King MD, Teillet PM. Passive remote sensing of tropospheric aerosol and atmospheric correction for the aerosol effect. J Geophys Res 1997; 102(D14): 16815-16830. DOI: 10.1029/97JD01496.

[2] Vermote EF, Kotchenova S. Atmospheric correction for the monitoring of land surfaces. J Geophys Res 2008; 113: D23S90. DOI: 10.1029/2007JD009662.

[3] Belyaev MY, Belyaev BI, Ivanov DA, Katkovsky LV, Martinov AO, Riazantsev VV, Sarmin EE, Siliuk VA, Shukajlo VG. Atmospheric correction of data registered on board the ISS. Part I. Methodology for spectra [In Russian]. Modern problems of remote sensing the Earth from Space 2018; 15(6): 213-222. DOI: 10.21046/2070-74012018-15-6-213-222.

[4] Ibrahim A, Franz B, Ahmad Z, Healy R, Knobelspiesse K, Gao BC, Proctor Ch, Zhai PW. Atmospheric correction for

hyperspectral ocean color retrieval with application to the Hyperspectral Imager for the Coastal Ocean (HICO). Remote Sensing of Environment 2018; 204: 60-75. DOI: 10.1016/j.rse.2017.10.041.

[5] Kokhanovsky A, Box JE, Vandecrux B, Manko K, Lamare M, Smirnov A, Kern M. The determination of snow albedo from satellite measurements using fast atmospheric correction technique. Remote Sens 2020; 12(2): 234-251. DOI: 10.3390/rs12020234.

[6] Katkovsky LV, Martinov AO, Siliuk VA, Ivanov DA, Kokhanovsky AA. Fast atmospheric correction method for hyperspectral data. Remote Sens 2018; 10(11): 1698-1715. DOI: 10.3390/rs10111698.

[7] Lisenko SA. Atmospheric correction of multispectral satellite images based on the approximate model of the solar radiation transfer [In Russian]. Atmospheric Ocean Opt 2017; 30(9): 775-788. DOI: 10.15372/AOO20170906.

[8] Belov AM, Myasnikov VV. Atmospheric correction of hyper-spectral images using approximate solution of MODTRAN transmittance equation. Computer Optics 2014; 38(3): 489493. DOI: 10.18287/0134-2452-2014-38-3-489-493.

[9] Denisova AY, Juravel YN, Myasnikov VV. Estimation of parameters of a linear spectral mixture for hyperspectral images with atmospheric distortions. Computer Optics 2016; 40(3): 380-387. DOI: 10.18287/2412-6179-2016-403-380-387.

[10] Belov VV, Tarasenkov MV. On the accuracy and operation speed of RTM algorithms for atmospheric correction of satellite images in the visible and UV ranges [In Russian]. Atmospheric Ocean Opt 2013; 26(7): 564-571.

[11] Nikolaeva OV. A new algorithm of retrieving the surface albedo by satellite remote sensing data. Atmospheric Ocean Opt 2016; 29(4): 342-347. DOI: 10.1134/S1024856016040102.

[12] Nikolaeva OV. Studying the accuracy of the algorithm for retrieving the surface albedo with high spatial resolution from a fragment of a satellite Image. Atmospheric Ocean Opt 2016; 29(6): 526-532. DOI: 10.1134/S1024856016060105.

[13] Nikolaeva OV. On variants of the main atmospheric correction formula. Light Eng 2020; 28(4): 39-46. DOI: 10.33383/2019-061.

[14] EO-1 User Guide v.2.3. Source: (https://www.usgs.gov/centers/eros/earth-observing-1-hyperion-documentation).

[15] Katsev IL, Zege EP, Prikhach AS. Atmosphere aerosol microphysical model for Belarus and adjacent regions [In Russian]. Atmospheric Ocean Opt 2016; 29(7): 572-578. DOI: 10.15372/AOO20160707.

[16] Mishchenko MI, Dlugach JM, Yanovitskij EG, Zakharova NT. Bidirectional reflectance of flat optically thick particu-late layers: an efficient radiative transfer solution and applications to snow and soil surfaces. J Quant Spectrosc Radiat Transf 1999; 64: 409-432. DOI: 10.1016/S0022-4073(99)00028-X.

[17] Bassani C, Cavalli RM, Antonelli P. Influence of aerosol and surface reflectance variability on hyperspectral observed radiance. Atmos Meas Tech 2012; 5: 1193-1203. DOI: 10.5194/amt-5-1193-2012.

[18] Carlson BG. A method of characteristics and other improvements in solution methods for the transport equation. Nucl Sci Eng 1976; 61(3): 408-425. DOI: 10.13182/NSE76-A26927.

[19] Nikolaeva OV, Bass LP, Germogenova TA, Kuznetsov VS. Algorithms to calculation of radiative fields from localized sources via the code RADUGA-5.1(P). Transp Theory Stat Phys 2007; 36(4); 439-474. DOI: 10.1080/00411450701468308.

Сведения об авторе

Николаева Ольга Васильевна, 1971 года рождения, в 1993 году окончила факультет вычислительной математики и кибернетики Московского государственного университета им. М.В. Ломоносова. Старший научный сотрудник Института прикладной математики РАН. Область научных интересов: прямые и обратные задачи для уравнения переноса излучения.

ГРНТИ: 27.35.47

Поступила в редакцию 13 сентября 2021 г. Окончательный вариант - 8 ноября 2021 г.

Surface albedo retrieval based on high spatial resolution data

O.V. Nikolaeva1 1 Keldysh Institute of Applied Mathematics RAS, Moscow

Abstract

The paper aims to compare the accuracy of three methods of solving the atmospheric correction problem for a Lambertian surface using high spatial resolution remotely sensed data.

Three couples of formulas are presented. Each couple contains a formula for expressing the reflectance in a target pixel at the upper boundary of the atmosphere via albedo in target and adjacent pixels of the surface and a formula for the surface albedo in a target pixel via reflectances. The derivation of each couple of formulas is given. Derivation conditions are presented. Formulas of only one couple are found by solving the radiation transport equation in 1D geometry. Formulas of two other couples include values obtained by solving the transport equation in 3D geometry.

Results of testing the accuracy of all formulas when solving an atmospheric correction problem on data of high (30 m) spatial resolution are given. Problems with aerosol optical depths from 0.2 to 2 for a wavelength of X = 0.55 |im for all possible albedo values (from 0.1 to 0.9) in target and adjacent pixels are considered.

It is shown that only one couple of formulas out of the three gives high accuracy under any condition. Formulas of the two other couples give sufficient accuracy (with less than 10 % error) only for a small value of the aerosol optical depth and a small difference of the albedo of the target and adjacent pixels.

Keywords: surface albedo, atmospheric correction, high spatial resolution.

Citation: Nikolaeva OV. Surface albedo retrieval based on high spatial resolution data. Computer Optics 2022; 46(3): 406-414. DOI: 10.18287/2412-6179-C0-1046.

Author's information

Olga V. Nikolaeva (b. 1971) graduated from the Computational Mathematic and Cybernetics department of the Moscow State University. She works as senior researcher in the Keldish Institute of Applied Mathematics RAS. Research interests are direct and inverse problems for radiation transport equation.

Received September 13, 2021. The final version - November 8, 2021.