Просмотр кода страницы Ортокоррекция космических снимков с использованием RPC

== Необходимые данные ==

Для проведения ортокоррекции космического снимка необходимо:
* Космический снимок в одном из форматов поставки (обычно GeoTIFF)
* Данные RPC к снимку 
* Информация о рельефе в виде DEM (Digital Elevation Model)

=== Космические снимки ===

Рассмотрим подробнее комплекты поставки основных компаний рынка детальной космической съемки: GeoEye и DigitalGlobe (образы продукции собраны на [[Образцы данных съемки|отдельной странице]]).
==== GeoEye ====
Продуктовая линейка GeoEye включает в себя<ref name="GeoEye_products">[http://geoeye.com/CorpSite/products-and-services/imagery-collection/satellite-imagery-products/Default.aspx Satellite Imagery Products], Copyright © 2012 GeoEye, Inc.</ref>:
* Geo
* GeoProfessional
* GeoStereo
Подробнее о продуктах компании можно почитать скачав брошюру [http://geoeye.com/CorpSite/products-and-services/imagery-collection/satellite-imagery-products/Product_Guide.aspx Product Guide] (требует регистрацию, но потом дает прямую ссылку).

Нас в первую очередь интересует продукт '''Geo''' - который представляет набор данных подготовленный к ортокоррекции (он же самый дешевый).

Комплект поставки продукции GeoEye включает в себя (на примере [http://geoeyemediaportal.s3.amazonaws.com/assets/images/ChannelPartnerPortal/Samples/GE-1_Transportation.zip Transportation GeoEye-1 Sample]):

 GeoEye_logo.tif
 po_344780_aoi.dbf
 po_344780_aoi.prj
 po_344780_aoi.shp
 po_344780_aoi.shx
 po_344780_component.dbf
 po_344780_component.prj
 po_344780_component.shp
 po_344780_component.shx
 po_344780_image.dbf
 po_344780_image.prj
 po_344780_image.shp
 po_344780_image.shx
 po_344780_metadata.txt
 po_344780_rgb_0000000.hdr
 po_344780_rgb_0000000.tfw
 po_344780_rgb_0000000.tif
 po_344780_rgb_0000000_ovr.jgw
 po_344780_rgb_0000000_ovr.jpg
 po_344780_rgb_0000000_rpc.txt
 SingleOrganization_license.txt

Как можно заметить все файлы начинаются на '''po_344780''' - это номер заказа по которому подготовлен это рабочий набор.

В комплекте поставки имеются три shape-файла (район заказа, проекции изображений), текстовый файл метаданных, собственно изображение в формате GeoTIFF, уменьшенное изображение в формате JPEG с привязкой и файл с RPC данными. Следует отметить, что программное обеспечение воспринимает изображение не как один файл в формате GeoTIFF, а как набор из файлов, в который входят, для нашего случая, еще и файлы с расширением hdr, tfw, _rpc.txt. 

Аналогичный состав имеет комплект поставки материалов съемки с КА Ikonos. 

А вот комплект поставки с третьего КА компании - OrbView-3, который стал свободно доступен (подробнее можно почитать на странице [[Каталог данных OrbView-3|описания каталога данных OrbView-3]]), имеет ряд нюансов.

Комплект поставки продукции GeoEye с КА OrbView-3 включает в себя (на примере [http://gis-lab.info/data/orbview-3/samples/rb/3v050909p0000897861a520004700712m_001631680.zip территории в Белоруссии]):
 3v050909p0000897861a520004700712m_001631680.att
 3v050909p0000897861a520004700712m_001631680.dbf
 3v050909p0000897861a520004700712m_001631680.eph
 3v050909p0000897861a520004700712m_001631680.jgw
 3v050909p0000897861a520004700712m_001631680.jpg
 3v050909p0000897861a520004700712m_001631680.prj
 3v050909p0000897861a520004700712m_001631680.pvl
 3v050909p0000897861a520004700712m_001631680.shp
 3v050909p0000897861a520004700712m_001631680.shx
 3v050909p0000897861a520004700712m_001631680.tif
 3v050909p0000897861a520004700712m_001631680_aoi.dbf
 3v050909p0000897861a520004700712m_001631680_aoi.prj
 3v050909p0000897861a520004700712m_001631680_aoi.shp
 3v050909p0000897861a520004700712m_001631680_aoi.shx
 3v050909p0000897861a520004700712m_001631680_rpc.txt
 3v050909p0000897861a520004700712m_001631680_src.dbf
 3v050909p0000897861a520004700712m_001631680_src.prj
 3v050909p0000897861a520004700712m_001631680_src.shp
 3v050909p0000897861a520004700712m_001631680_src.shx
 unrestricted_license.txt
Как видим названия файлов также имеют общую часть '''3v050909p0000897861a520004700712m_001631680''' - здесь мы видим дату съемки 050909 (2005-09-09T09:32:38.745600Z), номер в каталоге и уровень обработки 897861a (89786 и 1А), номер заказа 001631680 (OPS001631680). А вот состав файлов отличается. Но главное отличие, что поставляемые TIFF файлы не имеют географической привязки и описания СК, внедренные в файл.

Для ортокоррекции таких файлов необходима предварительная подготовка (см. например, [[Ортокоррекция космических снимков в wxGIS]] и [[Ортокоррекция данных OrbView-3 с помощью GDAL]])

==== DigitalGlobe ====
Продуктовая линейка DigitalGlobe включает в себя<ref name="DigitalGlobe_products">[http://digitalglobe.com/products#imagery Products], © 2012 DIGITALGLOBE, INC.</ref>:
* [http://www.digitalglobe.com/downloads/BasicImagery-DS-BASIC-Web.pdf Basic Imagery]
* [http://www.digitalglobe.com/downloads/StandardImagery-DS-STAND-Web.pdf Standard Imagery]
* [http://www.digitalglobe.com/downloads/Stereo-DS-STER-Web.pdf Stereo Imagery]
* [http://www.digitalglobe.com/downloads/AdvancedOrtho-DS-AOS-Web.pdf Advanced Ortho Series]

В продуктовой линейке нас интересует продукты Basic Imagery и Standard Imagery. Рассмотрим продукт Standard Imagery.
Комплект поставки продукции DigitalGlobe включает в себя (на примере [http://www.digitalglobe.com/downloads/product_samples/ortho/Ortho_Ready_Stnd_Bundle_16bit.zip QuickBird: Ortho ready standard Bundle 16bit]):

 |-005554608030_01
 | |-005554608030_01_P001_MUL
 | | |-05JUL04180116-M2AS-005554608030_01_P001-BROWSE.JPG
 | | |-05JUL04180116-M2AS-005554608030_01_P001.IMD
 | | |-05JUL04180116-M2AS-005554608030_01_P001.RPB
 | | |-05JUL04180116-M2AS-005554608030_01_P001.TIF
 | | |-05JUL04180116-M2AS-005554608030_01_P001.TIL
 | | |-05JUL04180116-M2AS-005554608030_01_P001.XML
 | | |-05JUL04180116-M2AS-005554608030_01_P001_README.TXT
 | | |-DEMONSTRATION.TXT
 | |-005554608030_01_P001_PAN
 | | |-05JUL04180116-P2AS-005554608030_01_P001-BROWSE.JPG
 | | |-05JUL04180116-P2AS-005554608030_01_P001.IMD
 | | |-05JUL04180116-P2AS-005554608030_01_P001.RPB
 | | |-05JUL04180116-P2AS-005554608030_01_P001.TIF
 | | |-05JUL04180116-P2AS-005554608030_01_P001.TIL
 | | |-05JUL04180116-P2AS-005554608030_01_P001.XML
 | | |-05JUL04180116-P2AS-005554608030_01_P001_README.TXT
 | | |-DEMONSTRATION.TXT
 | |-GIS_FILES
 | | |-005554608030_01_ORDER_SHAPE.dbf
 | | |-005554608030_01_ORDER_SHAPE.shp
 | | |-005554608030_01_ORDER_SHAPE.shx
 | | |-005554608030_01_PRODUCT_SHAPE.dbf
 | | |-005554608030_01_PRODUCT_SHAPE.shp
 | | |-005554608030_01_PRODUCT_SHAPE.shx
 | | |-005554608030_01_STRIP_SHAPE.dbf
 | | |-005554608030_01_STRIP_SHAPE.shp
 | | |-005554608030_01_STRIP_SHAPE.shx
 | | |-005554608030_01_TILE_SHAPE.dbf
 | | |-005554608030_01_TILE_SHAPE.shp
 | | |-005554608030_01_TILE_SHAPE.shx
 | |-005554608030_01_LAYOUT.JPG
 | |-005554608030_01_README.TXT
 | |-005554608030_01_README.XML
 |-005554608030_01.MAN
 |-005554608030_01_EOT.TXT
Как можно видеть структура папок и файлов в поставке с КА QuickBird отличаются от продуктов компании GeoEye.

В данном продукте имеются как панхроматический, так и мультиспектральный снимок. Панхроматический имеет размер пиксела 0,6 метра, а мультиспектральный - 2,4 метра. Каждый из снимков лежит в соответствующей папке 005554608030_01_P001_PAN и 005554608030_01_P001_MUL. 

Посмотрим, какие данные составляют снимок: 

* файл, который оканчивается на -BROWSE.JPG, представляет собой уменьшенную копию изображения
* файлы с расширениями IMD (image metadata), XML, _README.TXT содержат метаданные изображения
* файл с расширением RPB содержит RPC данные
* файл с расширением TIF - это само изображение
* файл с расширением TIL содержит в себе схему нарезки на тайлы.

Следует обратить внимание, что если изображение разделено на тайлы, то файл RPC данных будет один для всех растровых файлов. Например, рассмотрим такой набор данных.
 08JUL17070604-M2AS_R1C1-052053322020_01_P001.TIF
 08JUL17070604-M2AS_R1C2-052053322020_01_P001.TIF
 08JUL17070604-M2AS_R2C1-052053322020_01_P001.TIF
 08JUL17070604-M2AS-052053322020_01_P001.IMD
 08JUL17070604-M2AS-052053322020_01_P001.RPB
 08JUL17070604-M2AS-052053322020_01_P001.TIL
 08JUL17070604-M2AS-052053322020_01_P001.XML
 08JUL17070604-M2AS-052053322020_01_P001_README.TXT
 08JUL17070604-M2AS-052053322020_01_P001-BROWSE.JPG
 BASE.TXT

Здесь мы видим три растра и один TIL и RPB файл.
Рассмотрим, что из себя представляет TIL файл.

 bandId = "Multi";
 numTiles = 3;
 tileSizeX = 4096;
 tileSizeY = 4096;
 tileUnits = "Pixels";
 tileOverlap = 0;
 BEGIN_GROUP = TILE_1
 	filename = "08JUL17070604-M2AS_R1C1-052053322020_01_P001.TIF";
 	ULColOffset = 0;
 	ULRowOffset = 0;
 	URColOffset = 4095;
 	URRowOffset = 0;
 	LRColOffset = 4095;
 	LRRowOffset = 4095;
 	LLColOffset = 0;
 	LLRowOffset = 4095;
 	ULLon =     63.71122049;
 	ULLat =     56.27582320;
 	URLon =     63.86991680;
 	URLat =     56.27480990;
 	LRLon =     63.86791836;
 	LRLat =     56.18651934;
 	LLLon =     63.70958654;
 	LLLat =     56.18752929;
 	ULX = 544040.40000002;
 	ULY = 6237005.99999958;
 	URX = 553868.40000006;
 	URY = 6237005.99999954;
 	LRX = 553868.40000006;
 	LRY = 6227177.99999955;
 	LLX = 544040.40000002;
 	LLY = 6227177.99999959;
 END_GROUP = TILE_1
 BEGIN_GROUP = TILE_2
 	filename = "08JUL17070604-M2AS_R1C2-052053322020_01_P001.TIF";
 	ULColOffset = 4096;
 	ULRowOffset = 0;
 	URColOffset = 6907;
 	URRowOffset = 0;
 	LRColOffset = 6907;
 	LRRowOffset = 4095;
 	LLColOffset = 4096;
 	LLRowOffset = 4095;
 	ULLon =     63.86995555;
 	ULLat =     56.27480962;
 	URLon =     63.97888572;
 	URLat =     56.27399629;
 	LRLon =     63.97663705;
 	LRLat =     56.18570843;
 	LLLon =     63.86795703;
 	LLLat =     56.18651907;
 	ULX = 553870.80000006;
 	ULY = 6237005.99999954;
 	URX = 560617.20000015;
 	URY = 6237005.99999951;
 	LRX = 560617.20000014;
 	LRY = 6227177.99999952;
 	LLX = 553870.80000006;
 	LLY = 6227177.99999955;
 END_GROUP = TILE_2
 BEGIN_GROUP = TILE_3
 	filename = "08JUL17070604-M2AS_R2C1-052053322020_01_P001.TIF";
 	ULColOffset = 0;
 	ULRowOffset = 4096;
 	URColOffset = 4095;
 	URRowOffset = 4096;
 	LRColOffset = 4095;
 	LRRowOffset = 5325;
 	LLColOffset = 0;
 	LLRowOffset = 5325;
 	ULLon =     63.70958614;
 	ULLat =     56.18750773;
 	URLon =     63.86791788;
 	URLat =     56.18649778;
 	LRLon =     63.86732029;
 	LRLat =     56.15999958;
 	LLLon =     63.70909755;
 	LLLat =     56.16100852;
 	ULX = 544040.40000002;
 	ULY = 6227175.59999959;
 	URX = 553868.40000006;
 	URY = 6227175.59999955;
 	LRX = 553868.40000006;
 	LRY = 6224225.99999956;
 	LLX = 544040.40000002;
 	LLY = 6224225.99999959;
 END_GROUP = TILE_3
 END;

Как можно заметить, TIL файл содержит в себе перечень названий растров и их углы в географической, проекционной и пиксельной СК (некоторый аналог VRT формата из библиотеки GDAL). Программное обеспечение работает с таким файлом, как с целым изображением. 

Большинство программного обеспечения при работе с TIL выставляет СК WGS84 и использует градусные координаты, что приводит к тому, что растр "сплющивается" с потерей разрешения на местности, а значит и качества. В таком случае, чтобы корректно выполнить ортокоррекцию такого набора, необходимо выполнить мозаику тайлов в единое изображение с тем же именем, как и RPB файл. 

=== Данные RPC ===
Самый правильный способ обработки данных для ортокоррекции требует реконструкции геометрии датчика в момент регистрации для каждой строки изображения. Классический фотограмметрический способ основан на извлечении цифровой модели местности (ЦМР) с помощью обработки стереопары снимков и сопровождается ортотрансформированием одного из двух изображений; также возможно использовать существующие ранее ЦМР, качество которых соответствует масштабу конечного продукта, и использовать наземные точки привязки - (Ground Control Points – GCP). Это последнее звено обработки обычно следует за обработкой спутниковых данных, т.к. не все спутниковые датчики способны получать стерео пары, и в большинстве случаев требуется обработка единичных сцен. 

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

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

Модель камеры спутника в виде '''обобщенных аппроксимирующих функций (рациональных полиномов)''' можно считать одной из обобщающих моделей, используемых вместо строгих моделей уже более десяти лет (Greve и др. 1992 г .). Она выражает отношения между объектом и координатами изображения с помощью коэффициентов многочленов. Эти коэффициенты называются коэффициентами рационального многочлена ('''RPC – Rational Polynomial Coefficients'''). 

Надо полагать, что для получения неизвестных коэффициентов модели лучше ввести в вычисление строгие параметры датчика, даже в том случае, если обобщенная модель позволяет применять фотограмметрические процедуры без знания строгой физической модели датчика, типа датчика и процедуры получения изображения. Точность этих параметров влияет на точность конечных продуктов, даже если коэффициенты также могут быть вычислены, только начинаясь с соответствующего количества GCP, опознанных на изображении и измеренных по картам ( Tao и др. 2001 г .).<ref name="rpc">[http://www.sovzond.net/articles/article_geomobr.html Геометрическая обработка данных со спутника QuickBird]</ref>

Метод, основанный на RPC, был принят для обработки спутниковых данных в компаниях DigitalGlobe и GeoEye. Некоторые поставщики спутниковых данных, предоставляя коэффициенты рациональной функции, одновременно с этим частично скрывают информацию о датчике ради сохранения технологического «ноу-хау» . Для конечного же пользователя доступность RPC и программного обеспечения, основанного на них, облегчает точную геометрическую обработку изображений и позволяет обрабатывать данные, полученные от разнообразных датчиков, без необходимости добавлять новые модули к своему программному обеспечению. 

Для проведения ортокоррекции материалов космической съемки необходимы файлы с данными RPC, вычисленными с применением строгой модели сенсора. При использовании файла RPC можно вычислить нормализованную колонку и строку в изображении, как коэффициенты многочленов нормализованной геодезической широты, долготы и высоты.<ref name="rpc" />:

Перечень данных RPC приведен в таблице<ref name="rpc" />.
{| class="wikitable"
|-
! Наименование параметра !! Описание
|-
| HEIGHT_OFF || сдвиг высоты в метрах
|-
| HEIGHT_SCALE || нормализующий коэффициент для высоты
|-
| LAT_OFF || сдвиг широты в градусах в десятичном формате
|-
| LAT_SCALE || нормализующий коэффициент для широты
|-
| LINE_DEN_COEFF || список 20 коэффициентов для многочлена знаменателя строки
|-
| LINE_NUM_COEFF || список 20 коэффициентов для многочлена числителя строки
|-
| LINE_OFF || сдвиг строки в пикселах
|-
| LINE_SCALE || нормализующий коэффициент для строки
|-
| LONG_OFF || сдвиг долготы в градусах в десятичном формате
|-
| LONG_SCALE || нормализующий коэффициент для долготы
|-
| LONG_OFF || сдвиг долготы в градусах в десятичном формате
|-
| LONG_SCALE || нормализующий коэффициент для долготы
|-
| SAMP_DEN_COEFF || список 20 коэффициентов для многочлена знаменателя колонки
|-
| SAMP_NUM_COEFF || список 20 коэффициентов для многочлена числителя колонки
|-
| SAMP_OFF || сдвиг колонки в пикселах
|-
| SAMP_SCALE || нормализующий коэффициент для колонки
|}

Каждый многочлен используется до третьего порядка по P, L, H. Поэтому для каждой комбинации широты, долготы и высоты можно вычислить, какой пиксел (выборку, строку) следует рассматривать.

Для того чтобы использовать модель камеры спутника в виде обобщенных аппроксимирующих функций (рациональных полиномов), обязательно надо иметь DEM, иначе обработка должна выполняться с постоянным значением превышения, например для полностью плоских областей. Результаты могут быть улучшены, если использовать одну или более наземных точек привязки - GCP ( Ground Control Points )<ref name="rpc" />.

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

<math>~r_n = \frac{\sum^{20}_{i=1} LineNumCoef_i * \rho_i(P,L,H)}{\sum^{20}_{i=1} LineDenCoef_i * p_i(P,L,H)}</math>
<br /><br />
<math>~c_n = \frac{\sum^{20}_{i=1} SamNumCoef_i * \rho_i(P,L,H)}{\sum^{20}_{i=1} SamDenCoef_i * p_i(P,L,H)}</math>

где <math>~r_n, c_n, P, L, H</math> - нормализованные значения, выраженные как :<br />

<math>~P = \frac{(Latitude - LatOff)}{LatScale}</math>

<math>~L = \frac{(Longitude - LonOff)}{LonScale}</math>

<math>~H = \frac{(Height - HgtOff)}{HgtScale}</math>

<math>~r_n = \frac{(Row - LinOff)}{LinScale}</math>

<math>~c_n = \frac{(Column - SamOff)}{SamScale}</math>

и:

''LineNumCoef'' - список 20 коэффициентов для многочлена числителя строки;

''LineDenCoef'' - список 20 коэффициентов для многочлена знаменателя строки;

''SamNumCoef'' - список 20 коэффициентов для многочлена числителя колонки;

''SamDenCoef'' - список 20 коэффициентов для многочлена знаменателя колонки;

''LatOff'' - сдвиг широты в градусах в десятичном формате;

''LonOff'' - сдвиг долготы в градусах в десятичном формате;

''LinOff'' - сдвиг строки в пикселах;

''SamOff'' - сдвиг колонки в пикселах;

''HgtOff'' - сдвиг высоты в метрах;

''LinScale'' - нормализующий коэффициент для строки;

''SamScale'' - нормализующий коэффициент для колонки;

''LatScale'' - нормализующий коэффициент для широты;

''LonScale'' - нормализующий коэффициент для долготы;

''HgtScale'' - нормализующий коэффициент для высоты.

Числитель и знаменатель рациональной функции полиномиального уравнения представляет собой сумму 20-ти кубических полиномиальных функций и могут быть представлены в следующем виде:
<math>~\sum^{20}_{i=1} C_i * \rho_i(P,L,H) = C_1 + C_6*L*H + C_{11}*P*L*H + C_{16}*P^3</math>

<math>~+ C_2*L + C_7*P*H + C_{12}*L^3 + C_{17}*P*H^2</math>

<math>~+ C_3*P + C_8*L^2 + C_{13}*L*P^2 + C_{18}*L^2*H</math>

<math>~+ C_4*H + C_9*P^2 + C_{14}*L*H^2 + C_{19}*P^2*H</math>

<math>~+ C_5*L*P + C_{10}*H^2 * C_{15}*L^2*P + C_{20}*H^3</math><ref name="rpc_formula">[http://geotiff.maptools.org/rpc_prop.html RPCs in GeoTIFF]. Mathematical Model</ref>

Каждый многочлен используется до третьего порядка по P, L, H. Поэтому для каждой комбинации широты, долготы и высоты можно вычислить, какой пиксел (выборку, строку) следует рассматривать.<ref name="rpc" />

=== Информация о рельефе ===
Для выполнения ортокоррекции необходим файл рельефа в растровом виде (DEM, Digital elevation model). Высотные данные могут быть получены в результате наземных измерений, при помощи горизонталей с топографической карты, с помощью стереосъемки, по радарным данным или из общедоступных грубых ЦМР: [http://gis-lab.info/qa/srtm.html SRTM] (разрешение 30-90 м) и [http://gis-lab.info/qa/aster-gdem.html ASTER GDEM] (разрешение (15-90 м).

Наиболее удобным представляется использование данных ASTER GDEM, которые можно получить с [http://www.gdem.aster.ersdac.or.jp/index.jsp этого ресурса] (требуется предварительная регистрация). Для поиска необходимых участков ЦМР, очень удобно подставлять в качестве условий поиска шейп-файлы из поставки снимков.

[[Файл:Aster_gdem_search_interface.png|600px|thumb|center|Интерфейс поиска каталога ASTER GDEM]]