高解像度光学衛星画像を買って評価！～購入手順から質の評価方法まで一挙に解説～
光学衛星のインターネット上での購入手順や購入後の評価方法について、詳しくご紹介します。
1.はじめに
最近、「衛星画像」という言葉を耳にすることがずいぶん増えてきました。
政府が衛星データの利活用を積極的に推進していることもあり、テレビや新聞などのメディアでも衛星画像を使った新たな取り組みや活用事例が以前にも増して取り上げられるようになったからでしょう。
以前は、衛星画像といえば海外のプロバイダーが中心でしたが、日本国内でも衛星データを提供する企業が登場したことで、私たちにとって衛星画像がぐっと身近な存在になったように感じます。
こうした流れを背景に、衛星画像、とりわけ光学衛星画像について基礎的な知識を持つ人が徐々に増えてきました。しかし同時に、海外のプロバイダーを中心に様々な光学衛星が次々と打ち上げられているため、「どの衛星の画像を選べばいいのかよく分からない」と迷う方も多いのではないでしょうか。
例えば、「空間分解能30cm級」の衛星だけを見ても、Vantor社のWorldView-3やWorldView Legion、Planet社のPelican、Airbus社のPléiades Neoなど、複数の候補があります。一言で「光学画像」といっても、衛星毎にセンサー特性や撮影条件が異なるため、カタログスペックが同じでも、実際に得られる画像の品質や見え方には差があったりします。
こうした「スペックは似ているように見えるけれど、実際には何が違うの？」という疑問を解消するために、本記事では光学衛星画像の質の評価方法を解説します。
記事の前半では、光学衛星画像の仕組みを簡単に振り返りながら、光学画像の品質を評価するための方法を紹介します。続いて後半では、海外プラットフォームを通じて実際に高解像度光学衛星画像を購入し、前半で紹介した評価方法に沿って評価を行い、どんな結果が得られるのかをお見せします。
光学衛星画像をより深く理解して、実務や研究でしっかり活用していきたい――皆さまのお役に立てば幸いです。
2. 光学衛星の仕組みをコンパクトに復習
光学衛星とは、太陽光を光源として地表に反射した光（波長）を捉えるセンサを搭載した衛星のことです。センサが捉える波長には、私たち人間の目に見える可視光だけでなく、目に見えない赤外線も含まれます。
波長について復習したい方は下記をご覧ください。
光の波長って何？ なぜ人工衛星は人間の目に見えないものが見えるのか
光学衛星から得られる画像は、主に可視光で表示されるため、私たちが日常的に見ている「写真」とよく似ています。色や形が直感的に分かるため、建物や車両を見分けたり、都市や農地の様子を把握したりすることができ、幅広い場面で活用されています。
さらに、近赤外線など可視光以外の波長にも対応した撮影（マルチスペクトルと呼ばれます）が可能な衛星を使うと、目には見えない植物の活性や地表の微妙な変化を捉えることができます。これは、物体ごとに波長を反射する特性（反射特性）が異なることを上手く利用した技術です。
「波長の反射特性」についてもう少し詳しく知りたい方は、下記の記事も参考にしてください。
衛星データのキホン~分かること、種類、頻度、解像度、活用事例~
このように、光学衛星は搭載しているセンサを通じて観測を行うため、「センサがどのような波長を捉えるのか」「センサやレンズがどのように設計されているのか」によって画像の質が大きく変わります。例えば、センサの大きさやレンズの性能、焦点距離などの光学系の特性は、最終的な分解能やノイズの少なさに直結します。
そのため、空間分解能が同じ「30cm級」と表記されている衛星でも、搭載されているセンサの違いによって実際に撮れる画像の質には差が生じます。衛星データを選ぶ際は、どの衛星がどの程度の品質で画像を提供できるのかをしっかり把握しておくことが重要です。
3. 光学衛星画像の「質」を評価する4つの指標
では、「画像の質をしっかり見極めるには、実際に何をチェックすればよいのか？」という疑問に答えるのが、これからのパートです。ここでは、光学衛星画像の質を左右する以下の4つの評価指標を紹介します。
- 空間分解能（Spatial Resolution）
- 信号対雑音比（Signal-to-Noise Ratio）
- 放射分解能（Radiometric Resolution）
- 幾何精度（Geolocation Accuracy）
これらの指標がそれぞれ何を意味し、衛星データの活用時にどのような影響を及ぼすのかを、順を追って解説していきます。
3.1 空間分解能（Spatial Resolution）
まず初めに、空間分解能（Spatial Resolution）を紹介します。
空間分解能とは、衛星画像が地上の空間的情報をどの程度まで細かく認識できるかを示す指標です。
リモートセンシング分野では、一般的に「空間分解能」と言うと「地上サンプリング距離（Ground Sampling Distance: GSD）」という概念を指す傾向にあります。一方、より実際の識別能力に着目した「地上識別距離（Ground Resolved Distance: GRD）」という概念も存在します。
文脈によって、「空間分解能」がGSDを意味するのか、またはGRDを意味するのか異なる場合があるため、単に「空間分解能」と聞いた時には注意が必要です。
ここでは、空間分解能に対する理解を深めるため、GSDとGRDという2つの概念について解説します。
3.1.1 地上サンプリング距離（GSD）
地上サンプリング距離（GSD）とは、衛星が地上を撮影する際のサンプリングの間隔を地上距離で表したものです。例えば、「GSD=1m」であれば、地上1m四方が1ピクセルとして記録されることを意味します（下図）。GSDの値が大きいほど、地上の空間的情報をより粗く観測することになります。
ただし、GSDは不変というわけではありません。センサの画素サイズや視野角に加え、衛星の高度や撮影時の角度（直下視か斜め視か）など、複数の要因から影響を受けます。そのため、同一センサでも、撮影の角度や軌道が異なれば、画像ごとにGSDは変化します。このような背景から、多くの衛星画像のカタログでは「直下視でGSD＝30cm」などと、標準条件下での代表値が記載されているのが一般的です。
なお、GSDは「解像度（Resolution）」と呼ばれることもありますが、本記事では両者を厳密に区別します。
GSDはセンサ設計に基づいて定まる物理的かつ理論的な値であるのに対し、解像度は画像処理の結果として得られる表示上のピクセル間隔です。後者はユーザーが任意に設定できるため、例えば30cmGSDの画像を15cmや45cmの解像度で出力することも可能です。しかし、表示上のピクセル数を増やしても、実際に取得された空間的情報が増えるわけではありません。
また、通常、GSDと解像度は一致しており、その点が両者の概念を混同されやすい原因となります。ただ、上記で説明したとおり、GSDと解像度は密接に関連しつつも異なる概念となります。
詳しくは以下の宙畑の記事でも解説しています。
解像度と分解能の違いとは？ 地球観測含む様々な業界の言葉の定義と違いを知る意義
3.1.2 地上識別距離（GRD）
地上識別距離（GRD）は、センサが地上の空間的情報を「実際にどこまで細かく見分けられるか」を示す指標です。例えば、「GRD=1m」であれば、大きさが1mほどあれば識別できることを意味します（下図）。GRDの値が大きいほど、地物を細かく識別することが難しくなります。
GRDは、GSDよりも多くの要因に左右されます。具体的には、望遠鏡の口径や焦点距離といった光学的設計要素に加え、撮影中の衛星の微小な揺れ、大気のかすみ、焦点のズレ、更には画像処理の手法方法などが複合的に影響します。
このため、GRDは実際の画像を用いて評価される実測値であり、たとえセンサ側でGSDが細かく設定されていたとしても、画像上で得られる「見える細かさ」（GRD）はGSDよりも劣る場合が多々あります。衛星画像のカタログに「直下視でGSD＝30cm」と記載されていても、GRDは30cmより大きい値になっている場合もあるので注意が必要です。
GSDとGRDの違いを直感的に理解するため、自動車を撮影した衛星画像を模した下図を用意しました。
いずれも同一GSD及び同一解像度を前提としていますが、左から右にかけて画像のぼけ具合を段階的に強めています。ご覧のとおり、ぼけが強くなるにつれて自動車の形状が次第に不明瞭になっていきます。実際、同じGSDの衛星画像でも、ここまで極端ではありませんが、物体の見え方は異なります。GRDは、こうした「ぼけの程度」を数値として反映するための指標です。
GRDの測定には、「ISO 12233」などの国際規格や、NASAの技術文書に基づく手法が用いられます。
GRDは一般的に公開されていませんが、公開している衛星プロバイダーも存在し、例えば、Satellogic社は自社衛星のGRDの測定結果をこちらで公表しています。
GRDは、画像の鮮明さを定量的に評価するための重要な指標であり、衛星画像から物体を検出するようなタスクを中心に非常に有用な情報となります。
3.1.3 エッジのぼけからGRDを測定する方法
GRDの概念を理解したところで、その代表的な測定方法を紹介します。GRDは一般的に、画像内に含まれるエッジ（境界線）のぼけ具合を分析することで求められます。
まず、「ぼけ」とは何か、またそれが画像にどのように影響を与えるかを説明します。
下図に示すように、画像の「ぼけ」とは、ある1点を撮影した際に、その像が一点に集中せず、広がりをもって写る現象を指します。画像中のあらゆる像は、こうした点像の集まりで構成されているため、ぼけが生じると、線やエッジの境界もにじんで見えるようになります。
つまり、ぼけが大きいほど、線やエッジの輪郭は曖昧になり、拡がりも大きくなります。この関係を逆に利用すれば、線やエッジの像（拡がり具合）を調べることで、画像に含まれる「ぼけ」の特性を定量的に推定することが可能です。今回紹介するGRDの測定方法は、まさにこの考え方に基づいています。
ぼけを測定する手法として、一般的なカメラの画質評価では、線やエッジを含む解像度チャートと呼ばれる専用パターンを撮影し、ぼけの特性を評価します。解像度チャートには、様々な線やエッジが並んだパターンが含まれており、ぼけの方向や場所の依存性も測定できるよう設計されています。
一方、光学衛星の画質評価では、「チャートサイト」と呼ばれる、解像度チャートの地上版とも言える場所を撮影します。このチャートサイトには、1辺が最大で数十メートルにも及ぶ白黒のパターンが整然と設置されています。チャートサイトには、設置可能なパターンの種類には制限があるものの、エッジの鮮明さを高精度に測定できる最適な場所として広く利用されています。
それでは、具体的にGRDをどのように測定するのか、その手順を見ていきましょう。下図に、全体の流れを示します。
3.1.3.1 エッジ拡がり関数（Edge Spread Function: ESF）
エッジ拡がり関数（Edge Spread Function: ESF）とは、実画像におけるエッジ部分の明るさの変化（信号強度）を表す関数です。
この概念を直感的に理解するために、下図をご覧ください。
衛星がチャートサイトに設置された黒と白のパターンのエッジを撮影した場合、黒から白へと変化する方向に沿って画素ごとの明るさ（信号強度）をプロットすると、図中央下の青線のような明るさの変化プロファイルが得られます。
ぼけのない理想的な画像では、このプロファイルは90°の垂直線になります（図の中央下）。しかし、ぼけを含む実画像では、黒から白への変化はなだらかな曲線となります（図の右側下）。
この明るさの変化パターンのことをESFと呼びます。
3.1.3.2 線拡がり関数（Line Spread Function: LSF）
ESFは、「面」のエッジに対する明るさの変化を表す関数でした。これに対して、「線」上での明るさの変化の速さ（変化率）に着目したものが、線拡がり関数（Line Spread Fucntion: LSF）です。
変化率を計算するということで、ESFの微分によりLSFが得られます。これは、黒から白へと変化するエッジにおいて、どの位置でどれだけ急激に明るさが変化しているか、つまり線上でどの程度ぼけているかを定量的に評価するための指標です。
LSFを用いることで、エッジの鮮明さや画像の解像性能を、より詳細に分析することが可能になります。
3.1.3.3 放射伝達関数（Modulation Transfer Function: MTF）
LSFは、エッジ上における明るさの変化の急峻さを示すことで、画像の鮮明さを評価する指標でした。
しかし、画像の評価においては、エッジ周辺だけでなく画像全体に含まれる模様や構造も定量的に把握する必要があります。こうした模様や構造は、数学的には様々な細かさのパターンが重なり合って出来ていると考えられます。例えば、太い白黒の縞模様や、非常に細かいノイズのような模様も、すべて何かしらの空間的パターンの重ね合わせとして表現できます。
これらのパターンの細かさは、空間周波数という尺度で表されます。そして、画像に含まれる、これらの周波数成分を取り出す変換が、フーリエ変換です。 画像をフーリエ変換（=周波数領域への変換）すると、どのくらい細かいパターン（=空間周波数）が、どのくらいの強さで含まれているかが分かります。
—-
画像と縞模様とフーリエ変換の関係性について詳しく知りたい方は、こちらの動画が分かりやすいです。
—-
この考え方をLSFに適用し、LSFをフーリエ変換することで得られるのが、放射伝達関数（Modulation Transfer Function: MTF）となります。
もう少し踏み込んでお話すると、LSFはセンサが空間的に信号をどのように「ぼかす」かを表しています。これをフーリエ変換することで、模様の細かさ（空間周波数）ごとに、信号がどれだけ減衰しているか（コントラストがどの程度維持されるか）を評価することが可能です。このとき得られるのがMTFで、MTFは「細かさ（空間周波数）ごとの伝達効率」を表します。
MTFは、通常「MTFカーブ」として次のように表現されます（上図）：
・横軸：空間周波数 f（模様の細かさ）
単位：[cycles/m]（1mあたりに白黒の繰り返しが何回あるか）
・縦軸：コントラストの再現率（0〜1）
一般的な傾向として：
・模様が太い（低周波）→ MTF ≒ 1（コントラストが保たれる）
・模様が細かい（高周波）→ MTF が徐々に低下（再現が難しくなる）
となります。
ただし、MTFカーブは連続的な関数であり、そのままではセンサの空間分解能を単一の数値で定量的に示すことができません。つまり、MTFカーブのどの点を基準としてセンサの性能を評価するか、あらかじめ基準を設定する必要があります。
その代表的な基準がMTF = 0.5です。これは、「センサが模様の明暗差（コントラスト）を元の50%まで再現できる限界」を意味し、この時の衛星画像の視点における「センサがかろうじて構造として識別できる最小の線間隔（地上距離）」がGRDに該当します。ここでようやく、エッジのぼけの話からGRDへと繋がりました。
では、MTF = 0.5 に相当する横軸（空間周波数)の模様とは、地上ではどのような構造になるでしょうか。
ここで、下図のような等間隔の白黒のストライプパターンを考えます。
このストライプには、1周期の中に明暗の境界（エッジ）が2つ含まれています。隣り合うエッジ間の地上距離を d とすると、 1周期の長さは 2d になります。この時、対応する空間周波数fは次のように表されます：
この空間周波数において、センサのMTFがちょうど0.5になるような d が、GRDとなります。
式で表すと、以下のようになります：
この式でGRDが意味するところは：
センサがコントラストを半分まで保って再現できる、隣り合う線の地上間隔
となります。つまり、繰り返しになりますが、GRDとは、「センサがかろうじて構造として識別可能な最小の線間隔」を、地上距離として定義した指標です。
なお、実際の衛星画像では、センサの走査方向によってぼけの特性が異なる場合があります。
特にプッシュブルーム方式やTDI（Time Delay Integration）型センサでは、アロングトラック方向（衛星進行方向）とクロストラック方向（直交方向）でMTFやGRDの値が異なることが知られています。
また、エリアセンサであっても、光学系の収差や設計上の非対称性によって方向依存性が生じる場合があります。
そのため、GRDを評価する際には、可能であれば複数方向での測定を行うことが望まれます。
加えて、実際の評価では、十分な数の有効なエッジを統計的に確保する必要があり、条件を満たしたデータを集めること自体が容易ではありません。
本章の後半で紹介する実測例では、簡潔化のため一方向の測定結果を用いていますが、現実にはこうした多方向性やサンプル数確保の難しさを伴うことを理解しておくと良いでしょう。
補足：GRDの定義について
今回、MTFが0.5になる空間周波数を基準にしたGRDの定義を採用しました。この定義はこの文献を参考基づいており、一般的な解像度指標としても広く採用されている「MTF50」という基準に該当します。
https://ntrs.nasa.gov/api/citations/20240003473/downloads/ASemple-JACIE-2024.pdf
一方で、Satellogic社の公開情報では、GRDを、LSFの半値幅 (ピークから半分になる時の幅: Full Width Half Max: FWHMなどと呼ばれます)として定義しています。
https://satellogic.com/2023/03/20/optimal-spatial-resolution-imagery/
MTF50とFWHMは、どちらも線の「ぼけ」を定量化する指標であり、いずれも線の拡がりを表していますが、数値として得られる値は異なります。そのため、異なる文献やプロバイダの数値を比較する際には、どちらの定義に基づいているかを確認することが重要です。
なお、GRDとは別に、FWHMの値だけを掲載している資料もあります。実務上は、一度LSFを求めておけば、FWHWもMTFも求めることが可能なので、両方の値を算出しておき、比較対象の仕様に応じて使い分けるのが望ましいでしょう。
3.2 信号対雑音比（Signal-to-Noise Ratio）
信号対雑音比（Signal to Noise Ratio: SNR）とは、衛星画像において各ピクセルに記録された値が、どれだけ信頼できる測定値として扱えるかを判断するための指標の一つです。
SNR の「信号（Signal）」は、センサー設計の段階では検出器が出力するデジタル値（DN）を指すこともあります。しかし本記事では、ユーザーが実際に画像解析を行う場面を想定し、キャリブレーション後の物理量（放射輝度や反射率）を信号として説明します。SNR とは「本来得たい信号の強さ」と「それに含まれるノイズ（ばらつき）」との比率を表し、SNR が高いほどノイズの影響は小さく、ピクセル値の信頼性が高いと考えられます。
このSNRの影響は、統計的な解析や物理量の抽出のみならず、AIを用いた画像解析にも及びます。例えば、物体検出やセグメンテーションを行うAIモデルでは、地物の輪郭や形状、反射特性などのパターンを学習して認識を行います。ところが、画像にノイズが多いと輪郭が不明瞭になったり、画素値の一貫性が失われたりして、学習そのものが不安定になり学習精度が著しく低下します。その結果、誤検出や見逃しといったエラーの増加を招く可能性が高くなります。
さらに、SNRが低い場合、時間差で撮影された画像の比較による変化検出にも影響が出ます。例えば、洪水による水域の拡大や建物の損壊といった実際の変化が、ノイズによるランダムなばらつきと区別しにくくなり、判断を誤るリスクが高まります。
このように、SNRは画像を「見る」「計測する」だけでなく、「AI に理解させる」うえでも極めて重要な基礎指標であり、衛星画像の品質を評価する際に必ず考慮すべき要素となります。
SNRの評価方法は、次のようになります。
下図に、大きさμの信号成分と、ノイズ成分を用意しました。実際の撮影画像は、信号成分とノイズ成分の合成画像となります。ノイズに偏りがない場合、撮影画像は信号成分の大きさμの値を中心にばらついたものとなるでしょう。このばらつきの大きさがノイズの大きさと考えることができます。これを数式的には以下で定義される標準偏差σによって評価します。
vi, ：各点での画素値
n ：全画素数
このようなノイズの大きさの評価にあたって、本来は光の入らない状態で撮影されたダーク画像を使用してセンサーによるノイズを純粋に評価するのが理想です。これは信号を含まない「ノイズのみ」の状態を測るため、最も厳密な方法と言えます。
さらに、ノイズの特性は入射輝度（明るさ）によっても変化します。暗い領域（穏やかな海面など）では信号自体が弱いため、ノイズの影響が相対的に大きくなります。一方、明るい領域（砂漠や雪原など）では、センサが受け取れる明るさの上限に達してしまったり（飽和）、明るさを数値に変換する際のきざみ誤差（量子化誤差）が支配的になったりすることがあります。そのため、SNRを評価する際には、複数の輝度条件においてノイズ特性を測定し、センサがどの明るさ範囲で安定して動作しているかを確認することが重要です。
3.3 放射分解能（Radiometric Resolution）
放射分解能（Radiometric Resolution）は、衛星画像の1ピクセルが持つ明るさ（信号の強さ）を、どれだけ細かく区別できるかを示す指標です。
例えば、8ビットなら256段階、10ビットなら1024段階の階調で、地表からの反射や放射の強さを表現できます。段階の数が多いほど、微妙な明るさの違いを捉えやすくなり、例えば作物の反射率や植生の健康状態など、わずかな変化・違いの解析精度が向上します。
この放射分解能は、どの範囲の信号（「ダイナミックレンジ」と言います）を、何段階に分けて記録できるかによって決まります。たとえ同じビット数であっても、ダイナミックレンジが広ければ 1 段階あたりの幅は大きくなり、微小な変化は表現しづらくなります。逆に範囲が狭すぎると、明るい部分の信号が上限を超えて「飽和」し、詳細な情報が失われるおそれがあります。
こうしたダイナミックレンジと放射分解能のバランスは、センサの運用方針や設計時にあらかじめ決められていることが多く、観測のたびに自由に変更できるものではありません。ただし、機種によっては撮影条件や対象に応じて自動補正・調整が行われる場合もあります。例えば、雲や雪面など、高反射の対象を撮影する際には、飽和を防ぐ処理が自動的に働くことがあります。
さらに、放射分解能の評価はデータの形式にも影響を受けます。衛星画像が輝度値や反射率などの物理量に変換された形で提供されている場合、階調幅や信頼性をより直接的に評価できます。
一方、近年の商用衛星画像では、「デジタルナンバー（Digital Number）」のみで提供されるケースも多く、その場合は実際の明るさを直接知ることはできません。この場合は、画像全体のダイナミックレンジやヒストグラムの広がり、画素値の分布を基に、階調の細かさや信頼性を画像内で相対的に判断することになります。
放射分解能は、画像の見た目のなめらかさだけでなく、ピクセルの値をどの程度信頼して数値解析や機械学習に活用できるかを左右する、重要な指標の一つです。
補足：実効放射分解能について
ここまで説明してきた放射分解能は、信号をどの程度の階調でサンプリングし記録できるかという、いわば量子化分解能の観点に基づくものでした。
一方、実際の観測では、センサのノイズ特性（SNR）によって識別できる階調数が制約されることがあります。
このように、ノイズの影響を含めて実際に区別可能な階調の細かさを示す概念を、実効放射分解能と呼びます。
理論上のビット深度よりも、実際の有効階調数は低くなる場合が多く、センサのSNRや光学系の安定性、入射輝度条件などによって変動します。
詳細なノイズ特性やその明るさ依存性については、3.2節のSNRの項とあわせて理解するとよいでしょう。
3.4 幾何精度（Geolocation / Geometric Accuracy）
これまで、空間分解能や信号対雑音比、放射分解能といった「画質」に関する指標をご紹介しました。しかし、実際に解析を行う際には「幾何精度」も非常に重要となります。
例えば、画像上で検出した地物を地図上に正しく載せたり、他のデータと組み合わせて比較や解析したりする場合、画像の位置情報がずれていると、得られる結果に誤差が生じてしまいます。
そこで、ここからは衛星画像が「どの程度正しい位置を示しているのか」を表す幾何精度について、より詳しく解説します。
3.4.1 幾何精度とは何か
幾何精度とは、画像上に写っている対象物が実際の地表とどのくらい正確に一致しているかを示す指標です。都市計画やGISデータの整備、防災対策など、位置情報が重要な業務では、この幾何精度が低いと数メートルから数十メートル単位のずれが発生し、実務に大きな影響を及ぼす可能性があります。
3.4.2 絶対精度とバンド間精度
幾何精度を理解するうえで押さえておきたいのが、「絶対精度」と「バンド間精度」という2つの概念です。
絶対精度は、GPSや地上基準点（Ground Control Point: GCP）などと比較して画像上の座標誤差を数値化したものです。例えば、「平均5m誤差」などの具体的な数値で示され、衛星画像を位置情報に紐づけて利用するうえでの基本的な評価項目です。位置情報の正確さが求められる解析や利用時には、この絶対精度をしっかり確認しておく必要があります。
Vandor社が公開しているWorldView-3 Data Sheetには、「Predicated < 3.5m CE90 without ground control」と明記されています。これは、検証に使用した地点の90%が3.5m以内の範囲に収まることを意味します。
次に、マルチスペクトル画像を扱う場合には、バンド間精度も見逃せません。バンド間精度とは、可視光や近赤外など複数の波長バンドが正しく重なり合っているかを評価する指標です。例えば、植生指数（NDVI）の解析やバンドを組み合わせて行う画像解析では、バンド間の位置ずれが大きいと結果に大きな誤差が生じます。衛星画像を用いて高い分析精度を目指す場合は、バンド間精度が良好な衛星画像を選ぶ必要があります。
米国地質調査所（the United States Geological Survey: USGS）傘下の地球資源観測科学センターによる調査によれば、WorldView-3のバンド間精度は0.062m以内と、1ピクセル内に収まっています。
幾何精度は衛星やプロバイダーによって異なり、高精度補正を売りにする衛星ほど、絶対精度・バンド間精度ともに優れている傾向にあります。
また、衛星プロバイダーが提供する以上に厳密な幾何精度が必要な場合には、ユーザー自身が追加のGCPを使って独自に再補正を行うケースもあります。ただし、GCPを使っての独自補正を行ったとしても、多少の誤差が残り、そうした誤差のことを上記「絶対精度」と対になるものとして、「相対精度」と言います。
幾何精度は、画像が持つ位置情報の正確さを評価するうえで欠かせない指標です。いくら高解像度な画像であっても位置情報が不正確であれば精密な解析は困難になります。したがって、衛星画像を利用する際には事前に幾何精度を必ず確認し、用途に合った精度を満たしているかをチェックしましょう。
以上、光学衛星の画質評価に用いられる4つの指標について紹介しました。各衛星プロバイダーはこれらの指標に基づいて性能評価を実施し、その結果を反映させながら、観測された生データを標準プロダクトなどの利用可能な形式へ変換する際に用いるセンサーモデルの最適化および更新を行っています。
4. 実光学衛星画像での評価
それでは、ご紹介した4つの指標に基づき、実際の光学衛星画像の評価に進んでいきましょう。
まずは①今回評価する光学衛星画像を紹介し、続いて②今回利用する海外プラットフォームの紹介し、その後③そのプラットフォームでの購入手順を確認します。最後に、④実際にどのように画像の質をチェックするのかを具体的に解説します。
4.1 評価する光学衛星画像の紹介
購入する光学衛星は、中国の光学衛星「SuperView Neo-1」となります。
SuperView Neo-1の情報は以下となります。
| 別名 | Siwei Gaojing |
|---|---|
| 打ち上げ日 | 1&2号機：2022年4月29日 3&4号機：2025年2月27日 |
| 運用状態 | 1&2号機：運用中（2025年3月時点） 3&4号機：運用中（2025年3月時点） |
| 実運用機数／ 将来運用機数 |
4/16 |
| 空間分解能（GSD） | 0.3m |
| 観測幅 | 12km |
| 製造機関 | the China Academy of Space Technology (CAST) |
| 運用機関 | China Siwei Survey and Mapping Technology Co. Ltd. |
| 軌道 | 太陽同期軌道 |
| 高度 | 約530km |
| 重量 | 540kg |
SuperView Neo-1は、SuperView Neo-2及びSuperView Neo-3とコンステレーションを組んで運用されています。このうち、SuperView Neo-2はSAR衛星で、0.5m空間分解能を持っています。一方のSuperView Neo-3は光学衛星で、0.7m空間分解能に加えて、100kmもの広い範囲を一度に撮影できるのが特徴です。
最終的には、SuperView-1を16機、SuperView Neo-2を8機、SuperView Neo-3を4機の、合わせて28機体制のコンステレーションを目指しているとのことです。
4.2 利用する海外プラットフォームの紹介
衛星画像を入手する際、衛星画像の提供元（プロバイダー）と直接ごとやり取りする方法がありますが、複数のプロバイダーの画像をまとめて検索・購入できるプラットフォームを利用するのも便利です。プラットフォームでは、複数のプロバイダの衛星画像や過去の画像（アーカイブ）を簡単に検索・確認できます。
利用する海外プラットフォームは、米国のプラットフォーム「SkyFi」となります。
SkyFiは、ウェブで検索すると「SkyFi: Earth Observation, Made Easy」というフレーズが目に入るように、衛星データのアクセス・利用を手軽にすることを目指すプラットフォームです。パソコンだけでなくスマートフォン向けのアプリも提供しているため、誰でもどこでも簡単に衛星画像を取得し、活用できます。
さらに、SkyFiのエンドユーザーライセンス（End User License Agreement）はユーザーフレンドリーな設計となっており、プロバイダーとSkyFiのクレジットを表示さえすれば、ウェブやソーシャルメディアへの掲載はもちろん、画像解析の成果物を販売することも可能となっています。
2025年3月時点で連携している衛星プロバイダー（拠点）は、以下となっています。
| 光学 | ・Planet（米国） ・Satellogic（ウルグアイ） ・IMPRO（トルコ） ・GEOSAT（スペイン） ・Siwei（中国） |
|---|---|
| SAR | ・Umbra（米国） |
| その他 | ・Vexcel（米国） ・Urban Sky（米国） ・Near Space Labs（米国） |
4.3 今回の光学衛星画像の購入手順
ここからは、実際に海外プラットフォームの「SkyFi」を使って光学衛星画像を購入する手順を見ていきましょう。
まず、SkyFiのトップページにアクセスし、画面左に配置されているメニューリストの下部にある人型ボタン（赤枠）からアカウント登録を行います。
アカウント登録が済んだら、トップページ中央の「Explore existing images」から「COMMERCIAL」をクリックし、アーカイブ画像の検索画面に移動します。
今回は、インドのテランガーナ州の州都ハイデラバード近郊のShadnagarを購入エリアに指定します。このエリアは、衛星画像の校正地の一つとして有名な場所です。
画面の左上にある検索欄に、経度緯度（17.034288914469386, 78.1830843717838）を入力すると、地図がその場所に移動します。
その後、検索欄の下にあるフィルターボタンをクリックし、衛星の種類、Resolution（GSD）、バンド数、プロバイダーを順に設定し、「APPLY FILTERS」を押します。今回は、下図の設定にしています。
条件に合致したアーカイブ画像が画面左側に一覧で表示されます。ここで注意するのが購入面積です。今回購入する画像の最小購入面積は5km²と決まっているので、初期の1km²のままだと購入ができません。そのため、購入面積が5km²以上になるように購入エリアを調整します。範囲を広げたら、表示されているアーカイブ画像から希望のものを選びましょう。
希望する画像をクリックすると、観測日時、雲の割合、バンド数、オフナディア角度、提供形式、購入価格などの情報が確認できます。
購入価格を確認すると、112.50ドル、つまり約17,000円（1ドル150円換算）となっており、非常にリーズナブルな価格になっています。
宙畑メモ
一般的に、WorldViewなどの高解像度衛星のアーカイブ画像を購入する場合、最小購入面積が25km²と設定されているなどし、価格は今回の約10倍になることも珍しくありません。価格面だけを考えると、今回の価格設定はかなりお得だと言えるでしょう。
内容に問題がなければ「ADD TO CART」をクリックしてカートに入れ、「Checkout」ボタンから決済画面に進み、クレジットカードの情報を入力すれば購入は完了となります。
データの提供準備が整い次第、メールで通知が届きますので、それまでしばらく待ちましょう。なお、細かな情報にはなりますが、プラットフォーム上でオーダーの進行状況を随時確認することも可能です。
今回のアーカイブ画像は48時間以内には手元に届く予定となっており、今回は2時間6分で届きました。
では最後に、受け取ったデータの内容を確認してみましょう。
下図に示されているとおり、視覚化に適したPNG画像と、TIFFファイルやメタデータが格納されたZipファイルの2種類が提供されます。メタデータには、観測日や衛星天頂角、処理レベルなどの情報が含まれています。オルソ補正は既にされているプロダクトのようです。
Zipファイル内のTIFFファイルをQGISで表示してみたところ、正しく表示されることが確認できました。
これでSkyFiを通じた高解像度光学画像の購入手続きは以上となります。オンラインショッピングのように衛星画像を手軽に購入することができました。
ここからは、購入した画像の質を評価していきます。
4.4 空間分解能
まずは、空間分解能の一つであるGRDを測定します。GSDについては、既にプロバイダーより「0.3m」との情報が提供されているため、改めて測定する必要はありません。
GRDを求めるには、まずエッジ拡がり関数（ESF）の算出から取り掛かります。その前に、購入した画像における1ピクセルあたりの地上距離を確認しておきましょう。
位置情報が付与されたGeoTIFF形式の画像であれば、以下のコードで取得できます。
＜入力＞
import rasterio
tif_path = "img_SuperviewNeo1.tif"
with rasterio.open(tif_path) as src:
transform = src.transform
resolution_x = abs(transform.a) # pixel width (m/pixel)
resolution_y = abs(transform.e) # pixel height (m/pixel)
crs = src.crs
print(f"Resolution (X): {resolution_x:.4f} m/pixel")
print(f"Resolution (Y): {resolution_y:.4f} m/pixel")
print(f"CRS: {crs}")
＜出力＞
Resolution (X): 0.3000 m/pixel
Resolution (Y): 0.3000 m/pixel
CRS: EPSG:32644
出力から、解像度は0.3mであることがわかります。この値は、後ほど、ピクセル値から距離に換算する際に使用します。
次に、チャートサイト付近の画像を確認しておきます。
import tifffile
import numpy as np
import matplotlib.pyplot as plt
# tifffileを開く
# memmapでメモリ節約
with tifffile.TiffFile(tif_path) as tif:
mmap = tif.asarray(out="memmap")
chart_h = 3000
chart_hsize = 800
chart_w = 3250
crop_wsize = 800
chart_img = np.asarray(
mmap[chart_h : chart_h + chart_hsize, chart_w : chart_w + crop_wsize, :4]
)
chart_img = chart_img.astype(np.float32) / 1500
chart_img = chart_img.clip(0,1)
plt.imshow(chart_img[:, :, :3], vmin=0, vmax=1)
plt.axis("off")チャートサイトが綺麗に撮れていることが確認できました。
余談ですが、TIFF ファイルから特定領域を取り出す際、ファイルサイズの大きな GeoTIFF を処理する場合には、メモリマップ（memory map）を使って読み込むことで物理メモリの消費を抑えることができます。特に大容量の GeoTIFF を扱う際には、この方法が有効です。
では、チャートのエッジ部分を抽出します。画像中央に見える白黒の境界線を対象にします。ただし、今後の説明の都合上、上下を反転している点には注意してください。
crop_h = 3220
crop_hsize = 180
crop_w = 3630
crop_wsize = 70
img_org = np.asarray(
mmap[crop_h : crop_h + crop_hsize, crop_w : crop_w + crop_wsize, :4]
)
crop = img_org.astype(np.float32) / (1500) * 255
crop = crop[::-1, :, :]
plt.imshow(crop[:, :, :3].astype(np.uint8), vmin=0, vmax=255)
ESFを算出するために、エッジ画像からエッジのピクセル値の変化（エッジプロファイル）をサンプリングしていきます。
ここで、1つテクニックがあります。エッジをサンプリングする際、ある1つの列に沿って黒側から白側への値を取得してもピクセル値の変化を確認することはできますが、ピクセルの間隔以上に細かくサンプリングすることはできません。サンプリングの粗さは後のGRDの測定精度に直結します。そのため、今回は、エッジが画像上でわずかに傾いて写っていることを利用し、複数行にまたがるピクセル値からサブピクセル単位でサンプリングする手法を用います。詳しくは、こちらをご覧ください。
この手法では、下図のように、まず黒と白の境界となるエッジの線を直線としてモデル化します。続いて、エッジ付近の各ピクセルについて、画像座標からエッジ直線までの距離dを計算します。このdと対応するピクセル値のペアは、画像内のすべてのピクセルから得られるので、異なるdにおけるピクセル値を高密度にサンプリングできます。エッジが傾いていることで、dに微妙な違いが生まれ、結果としてサブピクセル単位の滑らかなプロファイルが得られる、という仕組みです。
上記をコードに落とし込むと以下になります。
from matplotlib.colors import LogNorm
def line_coefficients(x1, y1, x2, y2):
# エッジ上の2点の座標から直線の係数を計算
a = y2 - y1
b = x1 - x2
c = x2 * y1 - x1 * y2
return a, b, c
def signed_distance_and_value(img, x1, y1, x2, y2, starty=0, endy=None):
H, W = img.shape[:2]
if endy is None:
endy = H # デフォルトは画像全体
# 行方向の範囲でマスク作成
yy, xx = np.meshgrid(
np.arange(starty, endy), np.arange(W), indexing="ij"
) # shape: (endy-starty, W)
a, b, c = line_coefficients(x1, y1, x2, y2)
signed_dist = (a * xx + b * yy + c) / np.sqrt(a**2 + b**2)
values = img[starty:endy, :].astype(np.float32)
print(signed_dist.shape)
print(values.shape)
return np.stack([signed_dist.ravel(), values.ravel()], axis=1)
# ピクセル間隔 (cm/px)
PIXEL_SPACING_CM = 30
# チャンネルの定義
CH = 0
COLOR = ["R", "G", "B", "IR"][CH]
crop_ch = crop[:, :, CH]
# エッジ上の2点の座標
x1, y1 = 27, 1
x2, y2 = 49, 178.5
# エッジ上の2点の座標から距離と画素値の取得
dvals = signed_distance_and_value(crop_ch, x1, y1, x2, y2) # , starty=0, endy=100)
# 距離と画素値の2次元ヒストグラムを作成
plt.figure(figsize=(8, 4))
plt.hist2d(dvals[:, 0], dvals[:, 1], bins=200, cmap="viridis", norm=LogNorm())
plt.colorbar(label="Log-scaled Counts")
plt.xlabel("Signed Distance (pixels)")
plt.ylabel("Pixel Intensity")
plt.title("Raw ESF Samples")
plt.grid(True, linestyle=":", alpha=0.3)
plt.tight_layout()
plt.show()
上図に示すように、エッジからの距離とピクセル値の関係性を示すヒストグラムが作成できました。ここでは、各点の色がカウント（出現頻度）を表しており、青色から黄色へと変化するほど、その値のカウントが多いことになります。
次に、得られたサンプル群を距離ごとに平均化し、エッジプロファイルを求めます。このプロファイルがESFとなります。
from numpy.polynomial import Polynomial
def bin_esf(dvals, bins=500, range_x=(-10, 10)):
"""
距離-画素値の散布図（dvals）からESFをビン平均で求める
"""
x = dvals[:, 0] # 距離
y = dvals[:, 1] # 画素値
bin_edges = np.linspace(*range_x, bins + 1)
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
digitized = np.digitize(x, bin_edges) - 1
esf_vals = np.full_like(bin_centers, np.nan, dtype=np.float32)
for i in range(bins):
vals = y[digitized == i]
if len(vals) > 0:
esf_vals[i] = np.mean(vals)
return bin_centers, esf_vals
range_min, range_max = -20, 20
xx, esf_val = bin_esf(dvals, bins=400, range_x=(range_min, range_max))
plt.figure(figsize=(8, 4))
plt.plot(xx, esf_val, ".", color="black", label="Binned ESF", markersize=3)
plt.xlabel("Signed Distance (pixels)")
plt.ylabel("Pixel Intensity")
plt.title("Edge Spread Function (ESF)")
plt.grid(True)
上図に示すとおり、無事にESFができました。
次は、LSFを求めます。LSFは、ESFを微分することで求めることができます。
from scipy.optimize import curve_fit
from scipy.special import j1
from scipy.interpolate import interp1d
def bessel_lsf(x, a, b, x0):
"""
a * J1(b*(x - x0)) / (x - x0)
- a: 振幅
- b: スケール（波長的スケール）
- x0: 中心位置（ピークの位置）
"""
x_shifted = x - x0
x_shifted = np.where(x_shifted == 0, 1e-6, x_shifted) # avoid divide by zero
return a * (j1(b * x_shifted) / x_shifted)
def fit_lsf_with_bessel(xx, lsf_vals, fit_range=(-10, 10), p0=None):
# フィット範囲をマスク
mask = (xx >= fit_range[0]) & (xx <= fit_range[1])
x_fit = xx[mask]
y_fit = lsf_vals[mask]
# 初期値を設定
if p0 is None:
a0 = np.max(y_fit) - np.min(y_fit)
b0 = 1.0 / (fit_range[1] - fit_range[0])
x0 = x_fit[np.argmax(y_fit)]
p0 = [a0, b0, x0]
popt, _ = curve_fit(bessel_lsf, x_fit, y_fit, p0=p0)
return popt
def compute_lsf_from_esf(x, esf_vals):
"""
ESFの値をx方向で数値微分してLSFを得る。
x: 距離軸（等間隔であること）
esf_vals: ESFの値 (NaNを含む可能性あり)
"""
# NaNを線形補間で埋める
mask = ~np.isnan(esf_vals)
x_valid = x[mask]
esf_valid = esf_vals[mask]
# 再サンプリングして等間隔に整える（念のため）
interp_func = interp1d(x_valid, esf_valid, kind="linear", fill_value="extrapolate")
x_uniform = np.linspace(x_valid.min(), x_valid.max(), len(x)) # 元と同じ長さ
esf_uniform = interp_func(x_uniform)
# 微分（中心差分）
dx = np.mean(np.diff(x_uniform))
lsf_vals = np.gradient(esf_uniform, dx)
return x_uniform, lsf_vals
def compute_fwhm(x, y):
"""
任意の関数 y(x) に対して FWHM（半値全幅）を計算
- x: 等間隔であることが望ましい
- y: 関数値（最大値に対して正規化される）
"""
y = y / np.max(y) # 最大値で正規化
half_max = 0.5
# 正→負になる区間を探す
indices = np.where(np.diff(np.sign(y - half_max)))[0]
if len(indices) < 2:
return None # FWHMを定義できない
x0 = x[indices[0]]
x1 = x[indices[0] + 1]
x2 = x[indices[1]]
x3 = x[indices[1] + 1]
# 線形補間で厳密な交点を求める
fwhm_left = x0 + (x1 - x0) * (half_max - y[indices[0]]) / (
y[indices[0] + 1] - y[indices[0]]
)
fwhm_right = x2 + (x3 - x2) * (half_max - y[indices[1]]) / (
y[indices[1] + 1] - y[indices[1]]
)
return abs(fwhm_right - fwhm_left)
xx_lsf_val, lsf_val = compute_lsf_from_esf(xx, esf_val)
# フィット
bessel_params = fit_lsf_with_bessel(xx_lsf_val, lsf_val, fit_range=(-10, 10))
lsf_fit_bessel = bessel_lsf(xx_lsf_val, *bessel_params)
# FWHM
fwhm_raw = compute_fwhm(xx_lsf_val, lsf_val)
fwhm_fit = compute_fwhm(xx_lsf_val, lsf_fit_bessel)
fwhm_raw_cm = fwhm_raw * PIXEL_SPACING_CM
fwhm_fit_cm = fwhm_fit * PIXEL_SPACING_CM
# プロット
plt.figure(figsize=(8, 4))
title_txt = f"LSF (ch:{COLOR}, {PIXEL_SPACING_CM} cm/px)\n"
title_txt += f"FWHM-raw: {fwhm_raw:.3f} px ({fwhm_raw_cm:.1f} cm) | FWHM-fit: {fwhm_fit:.3f} px ({fwhm_fit_cm:.1f} cm)"
plt.title(title_txt)
plt.plot(xx_lsf_val, lsf_fit_bessel, label="Fit", linestyle="--", color="red")
plt.plot(xx_lsf_val, lsf_val, label="Raw", alpha=0.5, color="black")
plt.xlabel("Distance (pixels)")
plt.ylabel("LSF")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.xlim(-10, 10) # 拡大表示
plt.ylim(-25, 60) # 拡大表示
plt.show()
上図が、LSFの結果となります。灰色の線はESFをそのまま微分して得られた生データ（Raw）であり、それに対して滑らかにフィッティングを行った結果が赤線（Fit）になります。
なお、今回のフィッティングには、ベッセル関数（以下の式）を用いています。
この関数は、レンズなどの光学的な特性からくる理想的な光のぼけ方 (回折限界) の近似式として用いられるものです。この関数はこうした根拠のある式ではありますが、あくまで理想的な条件を仮定しているので、実際の測定では必ずしもこの関数をLSFのフィットに用いなければいけないというわけではありません。この関数以外に、ガウスの正規分布や、高次の多項式をなどが使われる場合があります。また、LSFではなく、微分する前のESFの段階でフィッティングしているものもあります。
今回は、LSF にした際に見られるピークの裾のスパイク成分も考慮し、それに適した特性を持つベッセル関数をフィッティング関数として採用しました。
さらに、ここで LSF の半値幅（FWHM） を計算しています。今回の結果では、FWHMは約 2.205 ピクセルとなりました。これを地上距離に換算すると、2.205 × 0.3m = 約0.66m（66cm） となります。
次に、LSFをフーリエ変換してMTFを求めます。
from numpy.fft import fft, fftfreq
def compute_mtf_and_mtf50(x, lsf, pad_factor=2):
dx = np.mean(np.diff(x))
N = len(lsf)
N_fft = 2 ** int(np.ceil(np.log2(N)) + pad_factor) # パディングで高分解に
mtf = np.abs(fft(lsf, n=N_fft))[: N_fft // 2]
mtf /= mtf[0]
freqs = fftfreq(N_fft, d=dx)[: N_fft // 2]
idx = np.where(np.diff(np.sign(mtf - 0.5)))[0]
if len(idx) == 0:
mtf50 = None
else:
i = idx[0]
x0, x1 = freqs[i], freqs[i + 1]
y0, y1 = mtf[i], mtf[i + 1]
mtf50 = x0 + (0.5 - y0) * (x1 - x0) / (y1 - y0)
return freqs, mtf, mtf50
# フィットしたLSFをより細かくサンプリング
x_fine = np.arange(-20, 20, 0.001)
lsf_fit_fine = bessel_lsf(x_fine, *bessel_params)
# x ∈ [-10, 10] に切り取り
mask = (x_fine >= -10) & (x_fine <= 10)
x_crop = x_fine[mask]
lsf_crop = lsf_fit_fine[mask]
# MTF計算
freqs, mtf_vals, mtf50 = compute_mtf_and_mtf50(x_crop, lsf_crop, pad_factor=0)
# プロット
plt.figure(figsize=(8, 4))
plt.plot(freqs, mtf_vals, label="MTF")
plt.axhline(0.5, color="gray", linestyle="--", alpha=0.5)
if mtf50 is not None:
plt.axvline(
mtf50, color="red", linestyle="--", label=f"MTF50 = {mtf50:.3f} cyc/pix"
)
plt.xlabel("Spatial Frequency (cycles/pixel)")
plt.ylabel("MTF")
grd = 1 / (2 * mtf50) * PIXEL_SPACING_CM
title_txt = f"MTF from LSF (ch:{COLOR}, {PIXEL_SPACING_CM} cm/px)" if mtf50 else "MTF from LSF"
title_txt += f"\nGRD: {grd:.1f} cm"
plt.title(title_txt)
plt.grid(True)
plt.legend()
plt.xlim(0, 1.1)
plt.ylim(0, 1.1)
plt.tight_layout()
plt.show()
上図のようなMTF曲線を得られました。この曲線からMTFが0.5になる時、x軸の線は0.277 [cycles/pixel] となることが分かります。
この場合、GRDを算出する式は以下でしたので、
したがって、
となります。この結果から、購入画像では
「54.2 cmの間隔の線であればコントラストを半分まで保って再現できる
(それより短い間隔ではコントラストが弱くなり見分けが付きにくくなる)」
ということになります。これはGSDの30cmよりも大きな値であることが分かります。あくまで一般論になりますが、衛星の光学系では、光学的なぼけを十分にサンプリングできるようなGSDに設計される場合が多いです。その結果、GRDはGSDよりも大きくなる傾向にあります。したがって、今回の測定値もおよそ妥当だと考えられます。
また、LSFのFWHMについても、今回測定された値は約2.2ピクセルでしたが、こちらの文献に記載されている様々な衛星のFWHMの値と比較すると、Super Doveで3.3ピクセル、 Black Skyで 2.7ピクセル、WorldView-2 で1.5ピクセルとなっており、おおよそ1ピクセル～3 ピクセルに収まることが分かります。これらの衛星はGSDが異なるため、単純に画質の優劣を比較することはできませんが、GSDに対するFWHMの比率を参考にすることで画質の相対的な傾向はつかむことが出来ると思います。
以上までは購入画像の赤バンドでの結果となりますが、他のバンド(青、緑、近赤外)での結果も合わせて以下にまとめます。
| 青 | 緑 | 赤 | 近赤外 | |
|---|---|---|---|---|
| FWHM [cm] | 65.6 | 65.8 | 66.2 | 65.1 |
| GRD [cm] | 53.8 | 54.0 | 54.2 | 53.5 |
バンドごとの差は 0.3% 以内に収まり、いずれのバンドでもほぼ同等の解像度が得られていることが確認できました。
一般に、光は波長が長くなるほど回折しやすくなるため、分解能は低下し（＝ぼけが大きくなる）、波長によって解像度に差が生じることが考えられます。
しかし今回、そのような傾向は見られませんでした。この理由として、測定精度や手法による不定性も考えられますが、もし購入画像が高解像度の単一バンド（パンクロ）画像を基に各バンドをシャープ化するパンシャープン処理などの画像処理が適用していた場合、それによってバンド間の解像度が均質化された可能性もあります。
以上が、GRDの基本的な測定方法の紹介となります。
この評価の発展として、白黒のチャートが左右ではなく上下に並んだ部分から水平方向のエッジも同様に測定してエッジの方向の依存性を評価する事もできます。さらに、各バンドでのエッジの中心の位置を評価し、バンドごとの位置ずれを評価することもできます。
4.5 信号対雑音比（SNR）
では次に、SNRの観点から画像の質評価を行います。
SNRは、「本来得たい信号の強さ」と、それに含まれる「ノイズ（ばらつき）」との比率を表していました。高いSNRは、画像の画素値が物理的な意味（例：反射率や輝度）として信頼できることを示します。
SNRの評価にあたっては、ダーク画像を用いることが理想でした。しかし、衛星画像の多くではダーク画像は非公開であることが多く、現実には利用できません。
そこで今回は、実務の中で利用可能なものとして、以下の代替手法を採用します：
・画像中の均一な領域を抽出
見た目に均質な明るさを持つ領域を複数ピックアップし、そこを局所的に「信号 + ノイズ」とみなします。領域が十分に大きく均質であれば、シグナル成分の大きさは領域の全画素の平均値とみなせます。そして、SNRをシグナル成分の大きさ（μ）とノイズ成分（σ）の大きさの比
として表します。
・バンド別の処理
各領域について、R（Red）、G（Green）、B（Blue）、NIR（近赤外）のバンドごとに平均値と標準偏差（ノイズ量）を算出します。
・画素値ごとのSNR傾向の評価
ノイズの性質は信号強度によって異なるため、明るさの異なるいくつかの領域を選定します。これにより、SNRの変化傾向を確認できます。
以上の手順を基に、衛星画像からSNRを求めていきましょう。
import tifffile
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
tif_path = "img_SuperviewNeo1.tif"
with tifffile.TiffFile(tif_path) as tif:
mmap = tif.asarray(out="memmap")
# 表示用の画像のスケールをおおまかに決める
# 画像を間引いて各チャンネルの0をのぞいた最小値と最大値を取得
img_subsampled = mmap[::50, ::50, :4]
img_mins = img_subsampled.min(axis=(0, 1))
img_maxs = img_subsampled.max(axis=(0, 1))
min_max_list = []
for ch in range(4):
nonzero = img_subsampled[:, :, ch] > 0
vmin = img_subsampled[nonzero, ch].min()
vmax = img_subsampled[nonzero, ch].max()
min_max_list.append((int(vmin), int(vmax)))
# 評価に使う領域のh, w, sizeを指定
h_w_s_list = [
(3662, 4079, 30),
(3333, 4189, 30),
(3406, 4002, 30),
(3246, 3522, 20),
(3249, 3524, 20),
(3300, 3483, 20),
]
results = []
for crop_h, crop_w, crop_size in h_w_s_list:
# 評価に使う領域のcrop
hstart = crop_h
hend = crop_h + crop_size
wstart = crop_w
wend = crop_w + crop_size
img = np.asarray(mmap[hstart : hend, wstart : wend, :4])
# 表示用のfigの作成
fig = plt.figure(figsize=(10, 5))
ax_1 = fig.add_subplot(2, 5, 1)
ax_2 = fig.add_subplot(2, 5, 2)
ax_3 = fig.add_subplot(2, 5, 3)
ax_4 = fig.add_subplot(2, 5, 4)
ax_5 = fig.add_subplot(2, 5, 5)
axes = [ax_1, ax_2, ax_3, ax_4, ax_5]
cols = ["red", "green", "blue", "brown"]
channels = ["R", "G", "B", "NIR"]
ax_5 = fig.add_subplot(2, 5, (6, 8))
ax_6 = fig.add_subplot(2, 5, (9, 11))
margin = 80
img_roi = np.asarray(mmap[hstart - margin : hend + margin,
wstart - margin : wend + margin ,
:4])
# RGBの表示用の画像の作成
img_plot = np.zeros_like(img_roi).astype(np.float32)
for ch in range(4):
img_plot_ch = img_roi[:, :, ch].copy()
ch_min, ch_max = min_max_list[ch]
img_plot_ch = np.clip(img_plot_ch, ch_min, ch_max)
img_plot[:, :, ch] = (img_plot_ch - ch_min)/(ch_max - ch_min)
axes[0].imshow(img_plot[:, :, :3])
axes[0].add_patch(
Rectangle((margin, margin), crop_size, crop_size, edgecolor="red", fill=False, lw=1)
)
axes[0].axis("off")
axes[0].set_title("RGB")
res_cols = []
for i in range(4):
# チャンネルごとにsnrを計算
img_col = img[:, :, i]
mean = img_col.mean()
std = img_col.std()
snr = mean / std
# mean, snrの誤差を計算
N = img_col.size
mean_error = std * np.sqrt(1/N)
snr_error = snr * np.sqrt(3/(2*N))
print(f"{channels[i]} mu: {mean:.1f}, sigma: {std:.1f}")
axes[i+1].set_title(f"{channels[i]}")
# チャンネルごとの画像の表示
img_ch_plot = img[:, :, i]
ch_min, ch_max = min_max_list[i]
img_ch_plot = np.clip(img_ch_plot, ch_min, ch_max)
img_ch_plot = (img_ch_plot - ch_min)/(ch_max - ch_min)
axes[i+1].imshow(img_ch_plot, cmap="gray", vmin=0, vmax=1)
axes[i+1].axis("off")
kwargs = dict(bins=20,
alpha=0.5,
histtype="stepfilled",
edgecolor=cols[i],
color=cols[i],
label=channels[i])
# チャンネルごとのヒストグラムの表示
ax_5.hist(img_col.ravel(), **kwargs)
ax_5.set_xlabel("DN")
ax_5.set_ylabel("count")
res_cols.append({"mu": mean, "sigma": std, "snr": snr, "mu_error": mean_error, "snr_error": snr_error})
ax_5.legend(loc="center left", bbox_to_anchor=(1, 0.75))
results.append(res_cols)
# 結果の表示
text = ""
for c, ch in enumerate(channels):
text += f"{ch} $\mu$: {res_cols[c]['mu']:.1f}, "
text += f"$\sigma$: {res_cols[c]['sigma']:.1f}, "
text += f"SNR: {res_cols[c]['snr']:.1f} "
text += f"$\pm$ {res_cols[c]['snr_error']:.1f}"
text += "\n"
text += "\n"
ax_6.text(0., 0.4, text, va="top")
ax_6.axis("off")
plt.show()この処理で、以下のような結果が各領域ごとに得られます。図の上段 “RGB”の赤枠で抽出した領域に対して、Red、 Green、 Blue、NIR (近赤外)のそれぞれのバンドで輝度値のヒストグラムを作成しました (図下段)。それぞれのバンドに対して、平均値μ、標準偏差σ、SNR (μ/σ) を記載しています。
各領域での画素値に対するSNRの値をプロットしてみましょう。
channels = ["R", "G", "B", "NIR"]
colors = ["red", "green", "blue", "brown"]
fig, ax = plt.subplots(figsize=(8, 6))
for ch_idx, (ch, color) in enumerate(zip(channels, colors)):
mu_list = []
snr_list = []
mu_err_list = []
snr_err_list = []
for patch in results:
mu_list.append(patch[ch_idx]["mu"])
snr_list.append(patch[ch_idx]["snr"])
mu_err_list.append(patch[ch_idx]["mu_error"])
snr_err_list.append(patch[ch_idx]["snr_error"])
ax.errorbar(
mu_list,
snr_list,
xerr=mu_err_list,
yerr=snr_err_list,
fmt="o",
capsize=3,
label=ch,
color=color,
alpha=0.8,
)
ax.set_xlabel(r"$\mu$ (Mean pixel value)", fontsize=13)
ax.set_ylabel("SNR", fontsize=13)
ax.set_title("SNR vs $\mu$ by Channel", fontsize=14)
ax.grid(True)
ax.legend(title="Channel")
plt.tight_layout()
plt.show()結果が以下のプロットになります。
ここでは、なるべく均質でノイズ評価に適した領域を複数選び、SNRを各チャンネルごとに調べました。ただし、画素値が1000を超えるような明るい場所では、そうした均一な領域を見つけるのが難しく、分析は比較的暗めの領域（画素値500～900程度）に限定されています。プロットは横軸にその領域の平均値μ、縦軸にSNRをとり、明るさに対してSNRがどう変わるかを見ています。
なお、プロットのSNRの誤差は、抽出した領域の画素数 n によるμとσの測定のばらつきを考慮して
としています。
上図のプロットを見ると、すべてのチャンネルにおいて、画素値（明るさ）が大きくなるにつれてSNRも高くなる傾向が確認できました。特に青（B）チャンネルではその関係が比較的はっきりしており、明るくなるほどノイズの影響が相対的に小さくなるという、Poissonノイズ型のセンサ特性に沿った挙動が見られます。一方で、一部のデータ点ではこうした傾向から外れるケースもあり、地物の微妙な構造や領域の不均質さなどが影響している可能性があります。
まとめると、画像におけるSNR全体としては、得られたSNR値はおおよそ120〜210の間に収まっており、信号に対するノイズの影響は比較的抑えられていると考えられます。SNRの逆数をとって、信号に対するノイズの割合に換算すると概ね0.5～0.8%程度であり、この範囲であれば、定量的な画像解析（例えば、反射率の比較や指数計算、変化検出）にも利用できる水準と言えます。ただし、ノイズの許容範囲は解析の目的や精度要件によって異なるため、ケースバイケースでの判断が必要です。
4.6 放射分解能
放射分解能（Radiometric Resolution）は、衛星画像の1ピクセルが持つ明るさ（信号の強さ）を、どれだけ細かく区別できるかを示す指標でした。
今回の画像データではデジタルナンバー（DN）のみが提供されており、物理量としての反射率や放射輝度への変換情報は含まれていません。
そのため、放射分解能を厳密に定量化することはできませんが、実際の画素値の分布から定性的に評価することは可能です。
例えば、画素値が11ビット=2048段階に対応している画像において特定のバンドの画素値がごく一部の範囲に集中していると、限られた階調しか使われていないことになります。つまり、わずかな明暗差の表現は制限されている可能性があります。逆に、画素値が広い範囲にわたって分布していれば、それだけセンサがシーン内の明暗を細かく記録しており、階調を十分に活かせている状態だと考えられます。
購入した画像の画素値の範囲を調べていきましょう。
まず、階調については、画像とともに提供されるメタデータよりビット深度が11bitと記載されているので、2048階調となります。
次に、画像の各バンドのヒストグラムと最大、最小値を見ていきましょう。
import tifffile
import numpy as np
import matplotlib.pyplot as plt
tif_path = "SkyFi_2513LIIU-1_2024-06-17_0548Z_DAY_SUPER-HIGH_Telangana-India/SkyFi_2513LIIU-1_2024-06-17_0548Z_DAY_SUPER-HIGH_Telangana-India.tif"
with tifffile.TiffFile(tif_path) as tif:
mmap = tif.asarray(out="memmap")
# 再描画（前提：mmapが [H, W, 4] で定義済み）
band_names = ["Red", "Green", "Blue", "NIR"]
colors = ["red", "green", "blue", "brown"]
channel_stats = []
plt.figure(figsize=(10, 6))
# ヒストグラム描画と統計計算
for i in range(4):
band = mmap[:, :, i]
band_nonzero = band[band > 0]
min_val = int(band_nonzero.min())
max_val = int(band_nonzero.max())
channel_stats.append((min_val, max_val))
plt.hist(
band_nonzero.flatten(),
bins=256,
alpha=0.5,
label=band_names[i],
color=colors[i],
range=(0, 2**11-1),
)
plt.xlabel("DN Value")
plt.ylabel("Count")
plt.legend(fontsize=12)
plt.grid(True, linestyle="--", linewidth=0.5)
# 軸の枠線
for spine in plt.gca().spines.values():
spine.set_visible(True)
# プロット右側に統計情報を表示
xmax = plt.xlim()[1]
ymax = plt.ylim()[1]*0.7
line_spacing = ymax * 0.07
for i, (name, (min_val, max_val)) in enumerate(zip(band_names, channel_stats)):
plt.text(
xmax * 0.97,
ymax - i * line_spacing,
f"{name}: Min, Max = {min_val}, {max_val}",
verticalalignment="top",
horizontalalignment="right",
fontsize=14,
bbox=dict(facecolor="white", edgecolor="none", pad=1.5),
)
plt.tight_layout()
plt.show()
上図より、全バンド含めてDN値はおおよそ368〜1943の範囲に分布しており、最大値から最小値が、およそ1600階調であらわされている事が分かります。使用されているビット深度が 11ビット（= 2048階調）であることを考慮すると、1600階調は約8割に当たります。
このことは、画像上で明部から暗部まで比較的広いダイナミックレンジが確保されており、自然地物や植生などの微妙な輝度差を視覚的に、ある程度滑らかに再現できていることを意味します。
ただし、観測条件や対象によっては、一部の領域でDN値が飽和する可能性もあります。そのため、放射分解能の絶対的な評価には注意が必要です。今回のように、画像のヒストグラムや表示時の調整値を参考にして、データの扱いやすさや解析の信頼性を見積もることに留めておくことが重要です。
4.7 幾何精度（Geolocation accuracy）
では最後に幾何精度を確認します。
本来であれば、GCPを利用して厳密に幾何精度を検証するのが理想的です。しかし今回の観測範囲に対応するGCPを持ち合わせていないため、代替手段としてGoogle Mapsとの比較を行います。 ただし、Google Maps自体にも若干の誤差が含まれている可能性があるため、あくまで簡易的な形での検証と捉えてください。また、筆者の目視・手作業での検証となりますので、その点にもご留意ください。
では早速、QGIS上で20か所のサンプルポイントを配置し、それぞれの位置がGoogle Maps上の地物とどれほど一致しているかを目視で確認します。
結果として、20か所すべてで明確な位置ずれが確認できました。加えて、各ポイントのずれ幅を計測すると、平均で約6.2mの誤差があることが分かりました。この程度の誤差であれば、1枚の画像から物体を検出して地図に載せる分には大きな問題にはならないかもしれません。一方、複数の画像を組み合わせて解析する際には幾何補正をしないまま進めるのは注意が必要と考えられます。
Credit : China Siwei Survey and Mapping Technology Co. Ltd. © SkyFi © Maps Data: Google
5. まとめ
本記事では、光学衛星画像の仕組みを振り返りながら、画像の質を評価する方法をご紹介しました。また、海外プラットフォーム「SkyFi」を通じて、中国の光学衛星「SuperView Neo-1」の画像を実際に購入し、評価を試みました。
今回はアーカイブ画像の購入手順を紹介しましたが、新規画像の購入方法についてはSkyFiがTutorialを公開しています。英語に抵抗感がある方は日本語字幕で閲覧してみてください。
How to Task a Satellite in 90 Seconds with SkyFi
あわせて、本記事の中盤でも述べたとおり、こうした商用衛星画像の品質評価は、米国地質調査所（USGS）の地球資源観測科学センター（特にJoint Agency Commercial Imagery Evaluation）が実施しています。近年は、様々な企業や組織から多数の地球観測衛星が打ち上げられており、それぞれの特徴や特性を理解したうえで衛星画像を適切に選定・活用することがますます重要になっています。Joint Agency Commercial Imagery Evaluationは、こうした目的に対応するために設立され、独立的な立場から画像の評価を行う組織です。
今回選定したインド・テランガーナ州の州都ハイデラバード近郊（Shadnagar）の校正サイトは、USGSが公表している校正サイトカタログの中からピックアップしたものになります。校正サイトには多様な種類があり、目的に応じて使い分けができるようになっています。ご興味のある方は「どのような校正サイトが、どのような目的で利用されているか」をぜひチェックしてみてください。
衛星画像の利用が広がり、様々な衛星プロバイダーが登場している今、衛星画像の質とは何か？をしっかりと理解し、どのデータを選ぶ・活用するかを考えることは一層重要になってきています。本記事が、皆さまの衛星データ活用に少しでもお役に立てれば幸いです。
