地球の画像

空間解像度が3メートルで、OpenStreetMap（PlanetScope衛星コンステレーション）の建物マスクがオーバーレイされた、偽色（緑、赤、近赤外線）の衛星画像



こんにちは、Habr！ 分析に使用するデータソースを常に拡張しているため、衛星画像も追加することにしました。 当社の衛星画像分析は、起業家精神と投資のための製品に役立ちます。 前者の場合、ジオデータ統計は、小売店を開く場所を理解するのに役立ちます;後者では、企業の活動を分析できます。 たとえば、建設会社の場合、1か月に建てられた床の数、農業会社の場合-栽培された作物のヘクタール数などを計算できます。



この記事では、地球の宇宙画像の大まかな概念を説明しようとします。衛星画像の操作を開始するときに遭遇する可能性のある困難について説明します：前処理、分析のアルゴリズム、および衛星画像と地理データを操作するためのPythonライブラリです。 だから、コンピュータビジョンの分野に興味のある方は、カットしてください。



それでは、衛星が撮影しているスペクトル領域から始めて、いくつかの衛星の撮影機器の特性を考えてみましょう。

スペクトルシグネチャ、大気の窓、衛星

異なる地上表面には、異なるスペクトルシグネチャがあります。 たとえば、新鮮な玄武岩質溶岩とアスファルトは、異なる量の赤外線を反射しますが、可視光では似ています。





乾いた草、アスファルト、新鮮な玄武岩質溶岩の反射率



地球の表面と同様に、大気中のガスにも固有のスペクトル特性があります。 ただし、すべての放射線が地球の大気を通過するわけではありません。 大気を通過するスペクトル範囲は「大気窓」と呼ばれ、衛星センサーはこれらの窓で測定するように調整されています。





大気窓



最高の空間解像度の1つである衛星-WorldView-3 （空間解像度は、画像内で識別可能な最小のオブジェクトのサイズを特徴付ける値です。天底は、この点で重力の方向と一致する方向です）。

撮影機器WorldView-3の特性

スペクトル範囲	お名前	範囲、n​​m	天底の空間分解能、m
パンクロマティック（1チャンネル、スペクトルの可視部分をカバー）		450〜800	0.31
マルチスペクトル（8チャネル）	沿岸	400-450	1.24
	青	450-510
	緑色	510-580
	黄色	585-625
	赤	630-690
	レッドエッジ	705-745
	近赤外	770-895
	近赤外	860-1040
短距離マルチバンド（8チャネル）	SWIR-1	1195-1225	3.70
	SWIR-2	1550-1590
	SWIR-3	1640-1680
	SWIR-4	1710-1750
	SWIR-5	2145-2185
	SWIR-6	2185-2225
	SWIR-7	2235-2285
	SWIR-8	2295-2365

これらのチャンネルに加えて、WorldView-3には大気補正用に特別に設計された12チャンネルがあります-CAVIS（雲、エアロゾル、蒸気、氷、雪）直下での解像度30 m、波長0.4から2.2ミクロン。





WorldView-2からの例; パンクロチャンネル





異なる空間解像度のサンプル画像



他の興味深い衛星はSkySat-1とその双子のSkySat-2です。 彼らは、30フレーム/秒の周波数と1.1 mの空間解像度で、1つの領域で最大90秒続くビデオを撮影できるという点で興味深いです。





SkySat衛星からのビデオ録画



パンクロマティックチャネルの空間分解能は0.9 m、マルチスペクトルチャネル（青、緑、赤、近赤外）-2 mです。



さらにいくつかの例：

1. PlanetScope衛星コンステレーションは、赤、緑、青、近赤外線で3 mの空間分解能で調査を実行します。
   
   
2. 衛星コンステレーションRapidEyeは、赤、極端な赤、緑、青、近赤外で5 mの空間分解能で調査を実施します。
   
   
3. ロシアの一連の衛星「Resource-P」は、パンクロマティックチャネルで0.7〜1 m、マルチスペクトルチャネル（8チャネル）で3〜4 mの空間分解能でイメージングを実行します。
   
   

ハイパースペクトルセンサーは、マルチスペクトルセンサーとは対照的に、スペクトルを多くの狭い範囲（約100〜200チャネル）に分割し、10〜80 mの異なる次数の空間分解能を持ちます。



ハイパースペクトル画像は、マルチスペクトル画像ほど広範ではありません。 ハイパースペクトルセンサーが搭載されている宇宙船はほとんどありません。 その中には、NASA EO-1衛星搭載のHyperion（廃止）、欧州宇宙機関が所有するPROBA衛星搭載のCHRIS、MightySatII衛星搭載のFTHSI、米国空軍研究所、GAW（ハイパースペクトル機器）、ロシアの宇宙船リソースP」。





オンボードセンサーCASI 1500から処理されたハイパースペクトル画像。 最大228チャネル。 スペクトル範囲0.4-1 nm





EO-1スペースセンサーからの処理済み画像。 220チャンネル スペクトル範囲0.4-2.5 nm



マルチチャンネル画像の作業を簡素化するために、純粋な素材のライブラリがあります 。 純粋な材料の反射率を示しています。

画像の前処理

画像の分析を進める前に、事前処理を実行する必要があります。

1. ジオリファレンス;
   
   
2. オルソ補正;
   
   
3. 放射補正とキャリブレーション。
   
   
4. 大気補正。
   
   

衛星画像のサプライヤは、画像を前処理するために追加の測定を行い、処理済み画像と自己補正のための追加情報を含む生画像の両方を提供できます。



また、衛星画像の色補正についても言及する価値があります。自然な色（赤、緑、青）の画像を、人にとって見慣れたビューにすることができます。



おおよそのジオロケーションは、軌道上の衛星の初期位置と画像ジオメトリによって計算されます。 地理参照は、地上管理点（GCP）によって洗練されています。 これらの制御点はマップと画像で検索され、異なる座標系での座標がわかっているため、ある座標系から別の座標系への変換パラメーター（共形、アフィン、遠近法または多項式）を見つけることができます。 GCP検索は、GPS撮影を使用して実行されます [ ソース。 1秒 230、 ソース 2 ]またはGCP座標が正確にわかっている2つの画像をキーポイントで比較することによって。



オルソ画像補正は、幾何学的画像補正のプロセスであり、遠近法の歪み、回転、レンズの歪みなどに起因する歪みを除去します。 画像は計画された投影に縮小されます。つまり、地形の各ポイントが厳密に垂直方向に直下で観測される投影に縮小されます。



衛星は非常に高い高度（数百キロメートル）から撮影しているため、直下で撮影する場合、歪みは最小限に抑える必要があります。 しかし、宇宙船は直下で常に写真を撮ることはできません。さもないと、特定の地点を通過する瞬間を待つのに非常に長い時間がかかります。 この欠点を解消するために、衛星は「完成」しており、ほとんどのフレームは有望です。 撮影角度は45度に達する可能性があり、高高度ではこれが大幅な歪みにつながることに注意してください。



画像の測定および位置特性が必要な場合は、オルソ補正を実行する必要があります。 追加の操作により画質が低下します。 これは、各画像ラインのレジストレーション時にセンサーのジオメトリを再構築し、ラスター形式でレリーフを表すことにより実行されます。



衛星カメラモデルは、一般化された近似関数（有理多項式-RPC係数）の形式で表され、高度データは、地形図、立体測量、レーダーデータ、または一般的に利用可能な大まかなデジタル標高モデルからの輪郭を使用した地上測定の結果として取得できます：SRTM（解像度30-90 m）およびASTER GDEM（解像度（15-90 m）。



放射補正-使用する撮影デバイスの特性によるハードウェア放射ひずみの画像の予備準備段階での補正。



スキャナーフィルムデバイスの場合、このような欠陥は画像変調として視覚的に観察されます（垂直および水平ストライプ）。 放射補正は、不良画像ピクセルとして観察される画像欠陥も除去します。





不良ピクセルと縦縞の除去



放射補正は、次の2つの方法で実行されます。

1. 調査機器の既知のパラメーターと設定の使用（補正テーブル）;
   
   
2. 統計的に。
   
   

最初のケースでは、長時間の地上試験と飛行試験に基づいて、調査機器に必要な補正パラメータが決定されます。 統計的方法による修正は、修正される画像自体から直接欠陥とその特徴を識別することにより実行されます。



画像のラジオメトリックキャリブレーション—輝度の「生の値」を他の画像のデータと比較できる物理単位に変換します。

$$

$$B_ \ lambda = K_ \ lambda \; \ ast \; DN_ \ lambda + \; C_ \ lambda、$$

どこで $$$B_ \ lambda$ -スペクトルゾーンのエネルギー輝度 $$$\ラムダ$ ;

$$$DN_ \ lambda$ -生の輝度値;

$$$K_ \ラムダ$ -校正係数;

$$$C_ \ lambda$ キャリブレーション定数です。



衛星の電磁放射は、センサーによって検出される前に、地球の大気を2回通過します。 大気の影響には、散乱と吸収という2つの主な効果があります。 大気中のガスの粒子と分子が電磁放射と相互作用し、元の経路から偏向すると、散乱が発生します。 吸収中、放射エネルギーの一部は吸収分子の内部エネルギーに変換され、その結果、大気が加熱されます。 電磁波放射に対する散乱と吸収の影響は、スペクトルのある部分から別の部分への遷移によって変化します。





衛星センサーに入る反射日射に影響する要因



大気補正を実行するためのさまざまなアルゴリズムがあります（たとえば、 DOSメソッド-Dark Object Subtraction ）。 モデルの入力パラメーターは、太陽とセンサーの位置のジオメトリ、気体成分の大気モデル、エアロゾルモデル（タイプと濃度）、大気の光学的厚さ、表面反射係数、スペクトルチャネルです。

大気補正では、画像からヘイズを除去するためのアルゴリズムを適用することもできます- ダークチャンネル事前分布を使用した単一イメージのヘイズ除去 （ 実装 ）。





ダークチャンネル事前分布を使用した単一画像ヘイズ除去のケーススタディ

インデックス画像

マルチチャネル画像からオブジェクトを調べる場合、重要なのは絶対値ではなく、さまざまなスペクトルゾーンのオブジェクトの輝度値の特性関係であることがよくあります。 これを行うには、いわゆるインデックスイメージを構築します。 このような画像では、元の画像とは対照的に、目的のオブジェクトがよりはっきりと目立ちます。

インデックス画像の例

インデックス名	フォーミュラ	申込み
酸化鉄指数	赤/青	酸化鉄を検出するには
粘土鉱物指数	中赤外線チャネル（CIR）内の輝度値の比率。 CIK1 / CIK2、ここでCIK1は1.55から1.75ミクロンの範囲、CIK2は2.08から2.35ミクロンの範囲	粘土鉱物を検出するには
鉄鉱物指数	平均赤外線（SIK1; 1.55〜1.75μm）チャネルの輝度値と近赤外線チャネル（NIR）の輝度値の比率。 SIK1 / BIK	腺ミネラルの含有量を検出する
レッドカラーインデックス（RI）	赤（K）と緑（H）の範囲の赤色の鉱物の反射率の違いに基づきます。 RI =（K-3）/（K + 3）	土壌中の酸化鉄を検出する
正規化された微分雪指数（NDSI）	NDSIは、赤（K）と短波赤外線（KIK）の範囲の雪の反射率の差によって特徴付けられる相対値です。 NDSI =（K-KIK）/（K + KIK）	雪で覆われた領域を強調するために使用されます。 雪の場合NDSI> 0.4
水指数（WI）	WI = 0.90μm/ 0.97μm	ハイパースペクトル画像から植生の水分量を決定するために使用されます。
正規化された微分植生指数（NDVI）	植物の葉のクロロフィルは、電磁スペクトルの近赤外線（NIR）範囲の放射を反射し、赤（K）で吸収します。 これら2つのチャンネルの輝度値の比率により、野菜を他の自然物から明確に分離して分析することができます。 NDVI =（NIR-K）/（NIR + K）	植生の存在と状態を示します。 NDVI値の範囲は-1〜1です。 密な植生：0.7; 疎植生：0.5; 開いた土：0.025; 雲：0; 雪と氷：-0.05; 水：-0.25; 人工材料（コンクリート、アスファルト）：-0.5

Pythonで衛星画像を操作する

衛星画像を保存するのが一般的な形式の1つにGeoTiffがあります（これに限定されます）。 PythonでGeoTiffを使用するには、 gdalまたはrasterioライブラリを使用できます。



gdalとrasterioをインストールするには、condaを使用することをお勧めします。

conda install -c conda-forge gdal conda install -c conda-forge rasterio
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

他の衛星画像ライブラリは、pipを介して簡単にインストールできます。



gdalを介したGeoTiffの読み取り：

 from osgeo import gdal import numpy as np src_ds = gdal.Open('image.tif') img = src_ds.ReadAsArray() #height, width, band img = np.rollaxis(img, 0, 3) width = src_ds.RasterXSize height = src_ds.RasterYSize gt = src_ds.GetGeoTransform() #    minx = gt[0] miny = gt[3] + width*gt[4] + height*gt[5] maxx = gt[0] + width*gt[1] + height*gt[2] maxy = gt[3]
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

衛星画像には多くの座標系があります。 これらは、地理座標系（GCS）と平面座標系（GCS）の2つのグループに分けることができます。



GSKでは、測定単位は角度で、座標は10進度で表されます。 最も有名なHSCはWGS 84（EPSG：4326）です。



UCSでは、単位は線形であり、座標はメートル、フィート、キロメートルなどで表現できるため、線形に補間できます。 GSKからPSKに切り替えるには、地図投影法が使用されます。 最も有名な投影法の1つはメルカトル図法です。



通常、マップ（画像のマークアップ）はラスター形式ではなく、ポイント、ライン、ポリゴンの形式で保存します。 これらの幾何学的オブジェクトの頂点の地理座標は、画像のマークアップとともにファイル内に保存されます。 fionaおよびshapelyライブラリを使用して、それらを読み取り、操作できます。



ポリゴンをラスタライズするためのスクリプト：

 import numpy as np import rasterio #fiona       ,  shp  geojson import fiona import cv2 from shapely.geometry import shape from shapely.ops import cascaded_union from shapely.ops import transform from shapely.geometry import MultiPolygon def rasterize_polygons(img_mask, polygons):   if not polygons:       return img_mask     if polygons.geom_type == 'Polygon':       polygons = MultiPolygon([polygons])     int_coords = lambda x: np.array(x).round().astype(np.int32)   exteriors = [int_coords(poly.exterior.coords) for poly in polygons]   interiors = [int_coords(pi.coords) for poly in polygons                for pi in poly.interiors]     cv2.fillPoly(img_mask, exteriors, 255)   cv2.fillPoly(img_mask, interiors, 0)   return img_mask def get_polygons(shapefile):     geoms = [feature["geometry"] for feature in shapefile]   polygons = list()   for g in geoms:           s = shape(g)           if s.geom_type == 'Polygon':               if s.is_valid:                   polygons.append(s)               else:                   #                     polygons.append(s.buffer(0))           elif s.geom_type == 'MultiPolygon':               for p in s:                   if p.is_valid:                       polygons.append(p)                   else:                       #                         polygons.append(p.buffer(0))   mpolygon = cascaded_union(polygons)     return mpolygon #   geojson    src = rasterio.open('image.tif') shapefile = fiona.open('buildings.geojson', "r") #  (      ) left = src.bounds.left right = src.bounds.right bottom = src.bounds.bottom top = src.bounds.top height,width = src.shape #  mpolygon = get_polygons(shapefile) #    ,          mpolygon = transform(lambda x, y, z=None: (width * (x - left) / float(right - left),\                                                      height - height * (y - bottom) / float(top - bottom)), mpolygon) #  real_mask = np.zeros((height,width), np.uint8) real_mask = rasterize_polygons(real_mask, mpolygon)
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

画像の復号化中に、交差するポリゴンの操作が必要になる場合があります（たとえば、地図上の建物に自動的にマークを付ける、雲がある画像のそれらの場所にある建物のマークアップも自動的に削除したい）。 これにはWyler-Athertonアルゴリズムがありますが、自己交差のないポリゴンでのみ機能します。 自己交差を解消するには、すべてのエッジとポリゴンの他のエッジの交差を確認し、新しい頂点を追加する必要があります。 これらの頂点は、対応するエッジを断片に分割します。 シェイプリーライブラリには、自己交差を排除する方法があります-バッファ（0）。



GSKからUCSに転送するには、PyProjライブラリを使用できます（または、rasterioで実行できます）。

 #      epsg:4326      epsg:32637 (  ) def wgs84_to_32637(lon, lat):   from pyproj import Proj, transform     inProj = Proj(init='epsg:4326')   outProj = Proj(init='epsg:32637')     x,y = transform(inProj,outProj,lon,lat)     return x,y
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

主成分法

画像に3つ以上のスペクトルチャネルが含まれる場合、3つの主要なコンポーネントのカラー画像を作成できます。これにより、情報の顕著な損失なしにデータ量を削減できます。



このような変換は、1つまたは2つのコンポーネントで明確に現れるダイナミクスを識別するために、異なる時間に撮影された一連の画像に対しても実行され、単一の座標系に取り込まれます。



4チャンネル画像を3チャンネルに圧縮するスクリプト：

 from osgeo import gdal from sklearn.decomposition import IncrementalPCA from sklearn import preprocessing import numpy as np import matplotlib.pyplot as plt # ,    height, width, band    ,  PCA     src_ds = gdal.Open('0.tif') img = src_ds.ReadAsArray() img = np.rollaxis(img, 0, 3)[300:1000, 700:1700, :] height, width, bands = img.shape #   PCA data = img.reshape((height * width, 4)).astype(np.float) num_data, dim = data.shape pca = IncrementalPCA(n_components = 3, whiten = True) #c  (  ) x = preprocessing.scale(data) #         y = pca.fit_transform(x) y = y.reshape((height, width, 3)) #   for idx in range(0, 3):   band_max = y[:, :, idx].max()   y[:, :, idx] = np.around(y[:, :, idx] * 255.0 / band_max).astype(np.uint8)
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

PlanetScopeサテライトグループの画像（赤、緑、青、色補正なし）





PlanetScope衛星グループの画像（緑、赤、近赤外線）





主成分法を使用して撮影されたスナップショット

スペクトル混合法

スペクトル分離法は、ピクセルサイズよりはるかに小さい画像内のオブジェクトを認識するために使用されます。



この方法の本質は次のとおりです。混合スペクトルは、既知の純粋なスペクトルと比較することによって分析されます。たとえば、既に述べた純粋な材料のスペクトルライブラリからです。 この既知の（クリーンな）スペクトルと各ピクセルのスペクトル内の不純物の比率の定量的評価が行われます。 このような評価を実行した後、ピクセルの色がこのピクセルのスペクトルでどの成分が優勢であるかを意味するようにペイントされた画像を取得できます。





混合スペクトル曲線





復号化されたスナップショット

衛星画像のセグメンテーション

現時点では、バイナリイメージセグメンテーションタスクの最先端の結果は、 U-Net モデルの 変更を示しています。





U-Netモデルのアーキテクチャ（出力画像のサイズは入力画像のサイズよりも小さい。これは、ネットワークが画像の端でより悪い予測をするために行われる）



U-Netの著者は、別のモデルに基づくアーキテクチャを開発しました。これは、畳み込みワードのみを特徴とする（最大プーリングをカウントしない）完全畳み込みネットワーク（FCN）です。

U-Netは、最大プーリングがアップコンボリューションに置き換えられるレイヤーが追加されるという点でFCNと異なります。 したがって、新しいレイヤーは出力解像度を徐々に高めます。 また、エンコーダー部分からの記号は、デコーダー部分からの記号と結合されるため、追加情報によりモデルはより正確な予測を行うことができます。



エンコーダー部分からデコーダー部分への機能転送がないモデルはSegNetと呼ばれ、実際にはU-Netよりも悪い結果を示します。





最大プーリング





アップコンボリューション



U-Net、Segnet、FCNに画像のサイズに関連付けられているレイヤーがないため、異なるサイズの画像を同じネットワークの入力に送信できます（画像サイズは、最初の畳み込みレイヤーのフィルターの数の倍数でなければなりません）。



kerasでは、これは次のように実装されます。

 inputs = Input((channel_number, None, None))
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

予測と同様に、トレーニングは、GPUメモリが許可する場合、画像フラグメント（クロップ）または画像全体のいずれかで実行できます。 さらに、最初の場合：

1）バッチのサイズが大きくなり、データにノイズが多く、異種の場合、モデルの精度に大きく影響します。

2）再トレーニングのリスクが少ない フルサイズの画像でトレーニングする場合よりもはるかに多くのデータがあります。



ただし、クロップで学習する場合、エッジ効果はより顕著になります-ネットワークは、中央に近い領域よりも画像のエッジで予測の精度が低くなります（予測ポイントが境界に近いほど、次の情報が少なくなります）。 この問題は、境界上の領域をオーバーラップおよび破棄または平均化するフラグメントのマスクを予測することで解決できます。



U-Netは、バイナリセグメンテーションのタスクのためのシンプルで強力なアーキテクチャです。githubでは、DLフレームワークの複数の実装を見つけることができますが、多数のクラスをセグメント化すると、このアーキテクチャはPSP-Netなどの他のアーキテクチャに失われます。 ここでは、セマンティックイメージセグメンテーションのアーキテクチャの興味深い概要を読むことができます。

建物の高さの決定

建物の高さは、影によって決まります。 これを行うには、メートル単位のピクセルサイズ、ピクセル単位の影の長さ、および太陽（太陽）の仰角（水平線上の太陽の角度）を知る必要があります。





太陽、衛星、建物のジオメトリ



タスクの全体的な難しさは、建物の影を可能な限り正確にセグメント化し、影の長さをピクセル単位で決定することです。 画像内の雲の存在も問題に追加されます。



建物の高さを決定するためのより正確な方法があります。 たとえば、地平線上の衛星の角度を考慮することができます。







影によって建物の高さを決定するアルゴリズムの例



3 mの空間分解能を持つPlanetScope衛星グループの画像では、MAPEを使用して建物の高さを決定する際の誤差（平均絶対パーセント誤差）は約30％でした。 合計で、40の建物と1つのショットが検査されました。 しかし、サブメーター画像では、研究者はわずか4-5％のエラーを受け取りました 。

おわりに

地球のリモートセンシングは、分析のための多くの機会を提供します。 たとえば、衛星画像に基づいたOrbital Insight、Spaceknow、Remote Sensing Metrics、OmniEarth、DataKindなどの企業は、米国の生産と消費を監視し、都市化、交通、植生、経済などを分析します。 同時に、写真自体がよりアクセスしやすくなっています。 たとえば、10m以上の空間解像度を持つLandsat-8およびSentinel-2衛星からの画像はパブリックドメインにあり、常に更新されています。



ロシアでは、Sovzond、ScanEx、Racurs、Geo-Alliance、Northern Geographic Companyも衛星画像を使用して地理分析を行っており、リモートセンシング衛星オペレーター企業の公式ディストリビューターです：Russian Space Systems（ロシア）、DigitalGlobe（アメリカ）、Planet（アメリカ） ）、エアバスディフェンスアンドスペース（フランス-ドイツ）など



PS昨日、2つのタスクで構成されるオンライン衛星画像コンテストを開始しました。

1. 建物と車のセグメンテーションおよび車のカウント。
   
   
2. 建物の高さの決定。
   
   

コンピュータビジョンと地球のリモートセンシングの分野に興味のある方は、ぜひご参加ください！



私たちは、ロスコスモス、エアバスディフェンス、スペース、プラネットスコープの写真をめぐる別のコンテストを開催する予定です。