Результат интеллектуальной деятельности: Способ получения и обработки изображений, искажённых турбулентной атмосферой

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

Рассмотрим особенности получения и обработки изображений, искаженных турбулентной атмосферой.

Наличие турбулентной атмосферы Земли между зондируемым участком земной поверхности и космическим аппаратом дистанционного зондирования земли (ДЗЗ) существенно ограничивает информационные возможности систем ДЗЗ. Возникают две проблемы: проблема «видения» через турбулентную атмосферу и проблема «изопланатичности» зондируемого участка земной поверхности.

Суть этих проблем состоит в том, что проблема «видения» накладывает ограничения на минимальный размер деталей, разрешаемых системой атмосфера–телескоп на зондируемом участке земной поверхности, а проблема «изопланатичности» ограничивает максимальный размер зондируемого участка земной поверхности, который ещё является пространственно–инвариантным к атмосферным искажениям, то есть эта проблема ограничивает поле зрения системы ДЗЗ.

Эти проблемы существенно зависят от условий получения, и, в частности, от условий регистрации изображений ДЗЗ. Так, если время регистрации (экспозиции) τ_Э превышает интервал временной корреляции атмосферных флуктуаций τ_A (так называемого времени «замороженности» турбулентностей атмосферы, τ_A ≈ 1 мсек) говорят о длинно–экспозиционной регистрации, а, если время регистрации τ_Э меньше τ_A, то говорят о коротко-экспозиционной регистрации. Эти два крайних случая отличаются характером атмосферных искажений. Так, если длинно-экспозиционное изображение, усредненное по атмосферным искажениям за время τ_Э > τ_A, обладает худшим разрешением (резкостью), чем мгновенное коротко–экспозиционное изображение, регистрируемое за время τ_Э < τ_A, то оно является пространственно–инвариантным к атмосферным искажениям во всём поле зрения системы атмосфера–телескоп ДЗЗ, в отличии от коротко–экспозиционного изображения, состоящего в этом поле из ряда мгновенных областей изопланатичности, пространственно-неинвариантных к атмосферным искажениям (Свиридов К.Н. «Технологии достижения высокого углового разрешения оптических систем атмосферного видения», изд. «Знание», 2005 [1]).

В связи с отмеченными особенностями атмосферных искажений, желание работать в широком поле зрения стимулировало получение длинно-экспозиционных изображений ДЗЗ, как в отечественных космических аппаратах ДЗЗ: «Ресурс-ДК1» (Петри Г. Российский спутник «Ресурс-ДК1»: альтернативный источник данных сверхвысокого разрешения. Геоматика, № 4, с. 38-42. 2010) [2] и «Ресурс-П» (Кирилин А.Н. и др. Космический аппарат «Ресурс-П». Геоматика, № 4, с. 23-26. 2010) [3], так и в американских космических аппаратах ДЗЗ: «QuickBird», «WorldView» и «GeoEye» (Лавров В.В. Космические съёмочные системы сверхвысокого разрешения. Геоинформационный портал ГИС Ассоциации, № 2, 2010) [4].

Эти публикации [2, 3, 4] могут быть рассмотрены, как аналоги предлагаемого здесь способа в части получения изображений. Используемая в них при детектировании технология ВЗН (временной задержки и накопления) приводит к регистрации длинно–экспозиционного изображения ДЗЗ, усредненного по атмосферным искажениям.

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

Характеристикой этого изображения является средняя длинно-экспозиционная оптическая передаточная функция, определяемая соотношением (FriedD.L. «Optical Resolution through a Randomly Inhomogeneous Medium for Very Long and Very Short Exposures», JOSA, v.56, № 10, p. 1372-1379, 1966) [5]

<ф д-э ()>= <ф()> = ф0() exp{-3,44[(F0 (, H)]5/3} (1)

Здесь =_{1–р2/F – пространственно-частотный вектор в апертуре телескопа ДЗЗ, длина волны солнечного излучения подсвета земной поверхности, преимущественно = 0,5 мкм, F – фокусное расстояние телескопа ДЗЗ, ф0() – оптическая передаточная функция телескопа ДЗЗ, определяемая для круглой апертуры диаметра D соотношением}

τ0() =  , (2)

а r₀(H) – пространственный радиус корреляции атмосферных флуктуаций светового излучения на высоте H [км] космического аппарата ДЗЗ, определяемый соотношением (Fried D.L. «Limiting Resolution Looking DownThrough the Atmosphere», J. Opt. Soc. Am., v.56, №10, p.1380, 1966) [6]

где I(H/L) = Г(7/3; H/L)/Г(7/3), Г (а, в) = ^{a-1e-x– неполная гамма функция, а L – высота турбулентного слоя атмосферы, преимущественно L ≈ 10 км.}

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

Регистрируемое искаженное атмосферой изображение можно представить интегралом свёртки вида

Iи() = ] + (), (4)

где – истинное распределение интенсивности объекта, – импульсный отклик пространственно-инвариантной системы атмосфера-телескоп (функция рассеяния точки (ФРТ), является Фурье преобразованием от ОПФ) и () – некоторый аддитивный шум в детектируемом изображении.

Пространственная фильтрация изображений вида (4), впервые предложенная и исследованная в 50-х годах 20 века Аббе, Портером, Цернике и Марешалем, описанная в работе (TsujinchiJ. «Correction of optical images by compensation of aberrations and by spatial frequency filtering» in «Progress in Optics». /ed. E. Wolf, North-Holland, Amsterdam, v.2, p.130, 1963) [7],требует выполнения традиционной последовательности операций, представленной на Фиг. 1., а именно, преобразования зарегистрированного изображения по Фурье в область пространственного спектра, пространственной фильтрации спектра умножением его на фильтрующую функцию и восстановления улучшенного фильтрацией изображения при обратном Фурье преобразовании от фильтрованного пространственного спектра.

В зависимости от вида фильтрующей функции различают: инверсную, винеровскую, гомоморфную и некоторые другие виды пространственной фильтрации (Оппенгейм А.В., Шафер Р.В. «Цифровая обработка сигналов», изд. «Связь», М., 1979) [8].

Для сравнения возможностей использования этих фильтрующих функций применительно к проблеме атмосферного видения был проведен эксперимент по машинному моделированию искаженного атмосферой и аддитивным шумом изображения участка лунной поверхности и его последующей пространственной фильтрации (BreedloveJ.R., Jr. «Digital Image Processing of Simulated Turbulence and Photon Noise Degraded Images of Extended Objects» in Proc. of the SPIE «Imaging through the Atmosphere», v.75, p.155, 1976) [9].

В эксперименте были использованы следующие фильтрующие функции:

1. Фильтрующая функция инверсного фильтра

() = (5)

2. Фильтрующая функция гомоморфного фильтра

(6)

3. Фильтрующая функция линейного винеровского фильтра

(7),

где и - величины спектральных плотностей объекта и шума, авеличина пространственного радиуса корреляции атмосферных флуктуаций светового излучения в пределах турбулентного слоя L у земной поверхности, определяемая соотношением [5]

r0(L) = [0,4232(h)dh]-3/5 (8)

Результаты этого эксперимента для различных отношений сигнал/шум представлены на Фиг. 2 и Фиг. 3. Отсюда видно, что при точном знании импульсного отклика системы атмосфера-телескоп (ФРТ) и при большом отношении сигнал/шум возможно восстановление дифракционно-ограниченных изображений. При этом наблюдается малое отличие в результатах фильтрации различными фильтрующими функциями, так как при большом отношении сигнал/шум все они стремятся к инверсной фильтрации. Гомоморфный фильтр визуально дает более удовлетворительное качество изображения, чем линейный винеровский фильтр. Это можно объяснить лучшим согласованием нелинейного гомоморфного фильтра с нелинейными характеристиками зрительной системы человека.

Эту публикацию [9] можно рассматривать в качестве аналога предлагаемого нами способа в части обработки изображения. Недостатком этого аналога является то, что в нем реализуется пространственная фильтрация с априорно точно известными характеристиками турбулентных искажений ФРТ и аддитивного шума.

В идеальном случае для точного определения ОПФ системы атмосфера-телескоп необходимо иметь некоторый точечный источник в изопланатичной с объектом области. Реально при наблюдениях ДЗЗ подобная возможность практически отсутствует и возникает необходимость апостериорного определения ОПФ системы ДЗЗ, сформировавшей изображение, по самому зарегистрированному изображению. Из представленных выше выражений для фильтрующих функций (4), (5) и (6) видно, что неизвестным параметром, необходимым для их реализации, является величина пространственного радиуса корреляции атмосферных флуктуаций а также величины спектральных плотностей объекта и шума . Заметим, что в последних нет необходимости при инверсной пространственной фильтрации, предлагаемой нами для изображений ДЗЗ.

Наиболее близким техническим решением к предлагаемому нами способу является способ фильтрации длинно-экспозиционного изображения астрономического объекта, описанный в работе (Дудинов В.Н., Цветкова B.C., Кришталь В.А., Хованский Н.А. «Обработка астрономических изображений методами когерентной оптики», Вестник ХГУ, № 117, стр.19, 1974) [10].

В этом способе, апостериорной пространственной фильтрации при фильтрации длинно-экспозиционного изображения астрономического объекта, планеты Сатурн, осуществляли перебор величин неизвестного параметра в диапазоне его возможных значений. При этом пространственная фильтрация оказывается неточной. Примеры недокорректированного и перекорректированного изображений, полученных при переборе параметра в инверсном фильтре (4), представлены соответственно на Фиг. 4^б и Фиг. 4^в.

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

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

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

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

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

K()= exp{3,44[(F (, H)]} ,

[где длина волны солнечного излучения подсвета [мкм], преимущественно = 0,5 мкм, F – фокусное расстояние телескопа ДЗЗ, формирующего изображение [м], =_{1– 2F – пространственно-частотный вектор в плоскости апертуры телескопа [цикл/мм] или [пар. лин/мм], (, H) – адаптивно изменяемым в процессе фильтрации пространственный радиус корреляции атмосферных флуктуаций светового излучения на высоте H [км] космического аппарата ДЗЗ [м], а К – номер итерации адаптивного поиска глобального максимума функции резкости фильтруемого изображения] и итерационно изменяют в ней величину параметра(, H), определяемого для первой итерации по эмпирической формуле}

=

[где L – высота турбулентного слоя атмосферы над земной поверхностью [м], преимущественно L≈10 км, а – величина пространственного радиуса корреляции атмосферных флуктуаций на границе турбулентного слоя [м], преимущественно≈0,1 м] при этом (К=1) фильтруют пространственный спектр длинно-экспозиционного изображения первой итерацией фильтрующей функции ^{1() и при обратном Фурье преобразовании восстанавливают первую итерацию адаптивно фильтруемого длинно-экспозиционного изображения, а для организации последующих шагов итерационного процесса адаптивной фильтрации, формируют функцию резкости фильтруемого изображения в виде}

S^K=^2d2d,

[где – область изображения и – область пространственного спектра, а – оценка пространственного спектра К-той итерации отфильтрованного изображения ] и максимизируют ее путем изменения, а точнее, уменьшения величины параметра в выражении для фильтрующей функции ^{K() и последующего повторения итерационных шагов описанной последовательности операций адаптивной фильтрации, завершая процесс адаптивной фильтрации при достижении функцией резкости SK глобального максимума.}

Дадим краткое математическое обоснование предлагаемого способа и рассмотрим блок-схему его практической реализации (Фиг.5).

Начнем рассмотрение с выбора фильтрующей функции. При анализе задачи дистанционного зондирования Земли необходимо отметить достаточно яркую подсветку Земли изучением Солнца, обеспечивающую возможность пренебрежения аддитивным шумом () в регистрируемом изображении I_и() (4). В связи с этим при фильтрации предлагается использовать инверсный фильтр, аналогичный (5), но для задач ДЗЗ отличающийся тем, что вместо (параметра Фрида), используемого для фильтрации астрономических изображений, здесь в соответствии с (1) необходимо использовать в выражении для фильтрующей функции параметр определяемый (3). Отсюда выражение для фильтрующей функции инверсного фильтра длинно-экспозиционного изображения ДЗЗ представимо в виде

() = (9)

Что касается величины по формуле (3), то для практических целей наиболее удобно использование эмпирической формулы, полученной в работе (Свиридов К.Н. «О предельном разрешении аэрокосмических систем дистанционного зондирования Земли (ДЗЗ)». Ракетно-космическое приборостроение и информационные системы, т. 1, вып. 1, с. 34-40, 2014) [11].

В связи с тем, что задачи ДЗЗ решаются в ближней зоне (зоне Френеля) (D>, от каждой точки земной поверхности распространяется сферическая волна. С увеличением высоты Н над турбулентным слоем L для этой расходящийся сферической волны величина увеличивается практически линейно, как

≈  (10)

Эта эмпирическая формула позволяет оценить максимальное и первое для итераций значение величины в зависимости от высоты H космического аппарата ДЗЗ.

Заметим, что высота H орбиты космического аппарата ДЗЗ является, как правило, известной и определяется траекторией его движения в Космосе.

Для автоматизации процесса адаптивной фильтрации введена функция резкости изображения S.

Здесь в качестве функции резкости была использована функция резкости вида [1]

SК =2d2d, (11)

где – область изображения, – область пространственного спектра, – оценка пространственного спектра для оценки адаптивно отфильтрованного изображения , а К – номер итерации адаптивного поиска глобального максимума функции резкости.

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

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

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

1) регистрируют усредненное по атмосферным искажениям длинно-экспозиционное изображение зондируемого участка земной поверхности

<Iи()> = , (12)

где – координаты в плоскости изображения на борту космического аппарата ДЗЗ, – истинное распределение интенсивности изображения зондируемого участка земной поверхности, – координаты в картинной плоскости объекта, – усредненный по атмосферным искажениям импульсный отклик системы атмосфера-телескоп ДЗЗ.

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

2) преобразуют зарегистрированное длинно-экспозиционное изображение ДЗЗ по Фурье в область его пространственного спектра

d = <() > = () (13)

Здесь () – пространственный спектр истинного распределения интенсивности объекта , а средняя ОПФ длинно-экспозиционного изображения, определяемая (1);

3) формируют фильтрующую функцию инверсного фильтра в виде

к() = (14)

4) для известной высоты H [км]космического аппарат ДЗЗ определяют максимальное значение параметра равное

= (15)

и используют его в качестве первой итерации параметра в фильтрующей функции (14);

5) осуществляют пространственную фильтрацию спектра длинно-экспозиционного изображения (13) первой итерацией фильтрующей функции

() >·1, (16)

6) получают первую итерацию отфильтрованного пространственного спектра в виде

= ()F); (17)

7) при обратном Фурье преобразовании отфильтрованного пространственного спектра (17) восстанавливают первую итерацию отфильтрованного длинно-экспозиционного изображения ДЗЗ

d = ; (18)

8) далее по восстановленному отфильтрованному изображению (18) формируют функцию резкости изображения для первой итерации процесса адаптивной фильтрации

S1 =2d ; (19)

9) изменяют, а точнее, уменьшают величину пространственного радиуса корреляции атмосферных флуктуаций

(20)

10) формируют вторую итерацию фильтрующей функции

() = ; (21)

11) повторяют последовательность операций адаптивной фильтрации 5) – 8) до формирования функции резкости отфильтрованного изображения второй итерации

S2 =2d (22)

12) сравнивают функцию резкости S^K+1 c функцией резкости S^K (в данном случае S² и S¹) и, если отказывается, что S^K+1<S^K, то поиск заканчивают, а в противном случае повторяют последовательность операций п.3 – п.8 до достижения функцией резкости S^K+1 глобального максимума S^K+1 = S^K = S_макс. При этом адаптивно отфильтрованное изображение п.7 обладает максимальной резкостью.

Сегодня отечественные средства ДЗЗ [2,3] обеспечивают пространственное разрешение не обработанного длинно-экспозиционного изображения 1 м. Зарубежные аналоги [4] обеспечивают при тех же условиях пространственное разрешение длинно-экспозиционного изображения 0,4м.

Пространственная фильтрация длинно-экспозиционных изображений астрономических объектов в прототипе [10] позволила улучшить пространственное разрешение в 2раза при отмеченной выше её не точной фильтрации.

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

Практическая реализация предлагаемого способа для обработки длинно-экспозиционных изображений ДЗЗ будет осуществлена в 2016 году.

Список цитированной литературы.

1. Свиридов К.Н. «Технологии достижения высокого углового разрешения оптических систем атмосферного видения», изд. «Знание»,2005 .

2. Петри Г. Российский спутник «Ресурс-ДК1»: альтернативный источник данных сверхвысокого разрешения. Геоматика, № 4, с. 38-42. 2010.

3. Кирилин А.Н. и др. Космический аппарат «Ресурс-П». Геоматика, № 4, с. 23-26. 2010.

4. Лавров В.В. Космические съёмочные системы сверхвысокого разрешения. Геоинформационный портал ГИС Ассоциации, № 2, 2010.

5. Fried D.L. «Optical Resolution through a Randomly Inhomogeneous Medium for Very Long and Very Short Exposures», JOSA, v.56, № 10, p.1372-1379, 1966.

6. Fried D.L. «Limiting Resolution Looking Down Through the Atmosphere», J. Opt. Soc. Am., v.56, №10, p.1380, 1966) [6]

7. Tsujinchi J. «Correction of optical images by compensation of aberrations and by spatial frequency filtering» in «Progress in Optics». /ed. E. Wolf, North-Holland, Amsterdam, v.2, p.130, 1963.

8. Оппенгейм А.В., Шафер Р.В. «Цифровая обработка сигналов», изд. «Связь», М., 1979.

9. Breedlove J.R., Jr. «Digital Image Processing of Simulated Turbulence and Photon Noise Degraded Images of Extended Objects» in Proc. of the SPIE «Imaging through the Atmosphere», v.75, p.155, 1976.

10. Дудинов В.Н., Цветкова B.C., Кришталь В.А., Хованский Н.А. «Обработка астрономических изображений методами когерентной оптики», Вестник ХГУ, № 117, стр.19, 1974.

11. Свиридов К.Н. «О предельном разрешении аэрокосмических систем дистанционного зондирования Земли (ДЗЗ)». Ракетно-космическое приборостроение и информационные системы, т. 1, вып. 1, с. 34-40, 2014.