はじめに
この記事では、Tellusを使った干渉SAR解析を行います。
コードは主に、宙畑というTellusのオウンドメディアの記事、【コード付き】地盤の沈降が分かる干渉SAR解析をTellusでやってみたを参考に作成しました。
使用言語はPython 3.10.12（2023/12/7 時点のGoogle Colaboratory環境下）です。
元々上記の記事は、非常にわかりやすく、記事内のコードをそのまま使えば干渉SAR解析を誰でも無料で実行できるようになっていました。しかし、2022年度のTellusのAPI仕様の変更により、コードを一部修正する必要が出てきました。また、SAR画像として用いていたPALSAR-2（JAXAが開発したALOS-2（だいち2号）に搭載されている合成開口レーダ）のデータが、Tellus上で無料利用ができなくなってしまいました。そのため、無料で試すにはコードの変更および解析対象の変更が必要になりました。
そこで、一世代前のPALSAR-1のデータを使用し、2023年度でも無料で実施可能な干渉SAR解析のコードを提示することが、この記事の目的です。
記事の中で、間違っている点やよろしくない表現がありましたら、ご指摘をよろしくお願いいたします。
対象者
・衛星データに興味がある方
・干渉SAR解析に興味がある方
最終的にできるようになること
- 注目地点の地盤の変動調査
この記事の特徴
・データを期間と場所で絞りこむところから始めている点
・2023年度で使用可能なTellusのAPIを使用している点
・PALSAR-1のデータ（L1.1）を使用している点
・CEOSフォーマットを簡易的に解説している点
上記の４点は、参考記事 宙畑 【コード付き】地盤の沈降が分かる干渉SAR解析をTellusでやってみたとの違いでもあります。
衛星データプラットフォームTellusとは
公式ページより引用
Tellus（テルース）は、「宇宙×ITで新しい価値を創造する」というビジョンを掲げて運用する日本発の衛星データプラットフォームです。衛星データの提供をはじめとし、データを利用した新たなビジネスを創出する環境をご用意しています。
Tellusの概要
これまでの衛星データでは、各衛星ごとに異なる取得先からそれぞれ異なる特殊なデータフォーマットを扱う必要がありました。そのため、扱いづらく衛星データの活用はとても限定的でした。そこで、共通データフォーマットの整備等により、誰もが、どんな衛星データも、ワンストップで統合的に検索・閲覧したり、各種データが容易に処理できる基盤が提案されました。この基盤は、衛星データプラットフォームと呼ばれています。中でもTellusは、さくらインターネット社が経済産業省の宇宙産業政策として開発した衛星データプラットフォームです。衛星データや衛星データを扱うさまざまなツールを提供しています。
干渉SARとは
干渉SAR（InSAR）とは、衛星搭載のSAR（合成開口レーダ）が発する電波の「位相」を利用して、地表の変動を監視する技術です。同じ地点を複数回観測し、受信した電波（観測データ）を干渉させて「位相差」を計算することで、観測地点の地表の形状や変化を広範囲かつ数cm単位で調べることが可能です。
また、現地に測定装置を設置する必要がないため、干渉SARの用途として、
・地盤沈下の把握
・地震や火山活動に伴う地表の変動の把握
・山間部の災害リスクの大きい斜面の継続的監視
などがあります。
※ 厳密には、地表の形状を計測する「干渉SAR」と、地表の変化を計測する「差分干渉SAR」に分けられますが、この記事では差分干渉SAR（D-InSAR）を主に扱います。
干渉SARの処理フロー
ここからは、 宙畑 【コード付き】地盤の沈降が分かる干渉SAR解析をTellusでやってみたの干渉SARの処理フローに従って、進めていきます。
干渉SARの処理フローは以下の通りです。
- 可干渉ペアの検索
- 強度・位相画像の出力
- 画像の位置合わせ
- 干渉処理
- 軌道縞の除去
- 地形縞の除去
- 差分干渉画像の完成
干渉SARは本来、高度な専門性が必要とされています。
しかし、ここでは基本的な処理のみを行い、精緻な解析に必要な高度処理を省略しています。
先に必要なライブラリをインポートしてから、各項目についてみていきましょう。
import os
import cv2
import struct
import numpy as np
import requests
import matplotlib.pyplot as plt
import math
from io import BytesIO
from matplotlib.projections import LambertAxes
import scipy
from scipy.interpolate import interp1d
import itertools
1. 可干渉ペアの検索
まず、干渉可能な画像のペアを探します。
ここで、単に同じ地点の画像同士であれば干渉できる、というわけではないことに注意してください。PALSAR-1の場合、可干渉な画像のペアの検索では、下記の項目がすべて一致する必要があります。
- 衛星の軌道経路
- 衛星の進行方向（北行 or 南行）
- 電波照射方向（右 or 左）
- オフナディア角
- 観測範囲の中心位置
- 送受信の偏波（HH, HV, VH, VV）
- 観測モード（スポットライト、高分解能、広域観測）
最初の2つは、人工衛星が同一軌道上にあることを保証する項目です。
それ以外は、SARの観測状態・観測方法が同一であることを保証する項目です。
衛星画像の可干渉の条件を確認したところで、Tellusを通して衛星データの検索・取得を行っていきましょう。
TellusのAPIトークンの取得
まず、Tellus公式サイトからアカウント登録を行い、マイページ上でAPIトークンの発行・取得を行ってください。こちらのAPIトークンを用いて、CLI上でデータを取得していきます。TellusのAPIトークンの発行・取得画面は、次のようになります。
Tellus上の衛星データの検索・取得
希望の地点の緯度・経度情報を入力して、プロダクト（衛星データ）を検索します。この時、先ほど述べた可干渉の条件に必要な情報も出力します。
衛星データの詳細な情報は、Tellusのデータカタログから取得することができます。今回使用するPALSAR_L1.1のプロダクトIDやカタログ情報は、データセット詳細 - Tellus Travelerに記載されています。
ここで、L1.1について簡単に触れておきます。衛星から吐き出された生データは、そのまま直接利用することは難しく、意味のある数値や画像に変換する前処理が必要になります。衛星データの前処理には段階があり、おおまかにLevel-0処理→Level-1処理→高次処理となっています。SARの観測データには、電波の強度と位相情報が含まれているのですが、今回用いるL1.1は、その2つの情報が、利用可能な状態まで処理されたデータになります。衛星データの前処理についての詳細は、宙畑 【図解】衛星データの前処理とは~概要、レベル別の処理内容と解説~を参照してください。
今回使用したAPIは以下の3つです。
# dataset_idはプロダクトID、data_idは観測データID、file_idはファイルID
# PALSAR_L1.1というプロダクトから観測データ（データシーン）を検索する場合
data_scene_search_url = "https://www.tellusxdp.com/api/traveler/v1/datasets/{dataset_id}/data-search/"
# 観測データ内のファイルIDを確認する場合
data_scene_file_id_url = "https://www.tellusxdp.com/api/traveler/v1/datasets/{dataset_id}/data/{data_id}/files/"
# 所望の観測データをダウンロードする場合
data_scene_download_url = "https://www.tellusxdp.com/api/traveler/v1/datasets/{dataset_id}/data/{data_id}/files/{file_id}/download-url/"
データの取得にはプロダクトID、観測データID、ファイルIDの3つのIDが必要です。CLIでは、逐一APIを叩いてIDを確認する必要がありますが、Tellusが提供するGUIであるTellus Travelerでは、より簡単に所望のデータまでたどり着くことができます。一方、GUI操作のため大量のデータのダウンロードには向いていないことから、用途別で使い分けるのが良いと思います。TellusのAPIの仕様や使い方は、APIリファレンスやAPIリファレンスの使い方を参考にしてください。
今回は、東京ディズニーリゾート周辺の干渉SAR解析をしてみようと思います。東京ディズニーリゾートの緯度経度は、Google Mapで取得した値を用いました。また、検索期間を2006年1月1日から2011年3月18日としました。5年という期間をざっくりと決めた後、Tellus Travelerにて、存在する観測データの撮影日時を調べ、この期間としました。
一つ目のAPIを用いて以下のコードを実行すると、 PALSAR_L1.1の観測データの情報が得られます。
BASE_API_URL = "https://www.tellusxdp.com/api/traveler/v1/datasets/{}/data-search/"
dataset_id = "8836ec9a-35b5-43c3-92be-263498061e91" # PALSAR_L1.1のプロダクトID
ACCESS_TOKEN = "YOUR_TOKEN"
HEADERS = {"Authorization": "Bearer " + ACCESS_TOKEN}
url = BASE_API_URL.format(dataset_id)
# 東京ディズニーリゾートの座標
lat = 35.631 # 緯度
lon = 139.883 # 経度
# リクエストボディ
r_body= {
"intersects": {
"type": "Polygon",
# 検索範囲を領域で指定
"coordinates": [
[
[lon, lat],
[lon+0.01, lat],
[lon+0.01, lat+0.01],
[lon, lat+0.01],
[lon, lat],
]
]
},
"query": {
# 検索期間の指定
"start_datetime": {"gte": "2006-01-01T15:00:00Z"},
"end_datetime": {"lte": "2011-03-18T12:31:12Z"},
}
}
# 該当するデータの情報を取得
response = requests.post(url, json=r_body, headers=HEADERS)
res_json = response.json()
得られたres_jsonは、type、features、metaの3つのkeyを持っており、featuresが観測データごとの詳細な情報を持っています。また、観測データごとにdataset_id、geometry、id、type、propertiesの5つのkeyを持っており、propertiesが可干渉な画像に関する情報を持っています。
なお、featuresの数が、今回の検索で取得した観測データの数です。
print(res_json.keys())
print(res_json["features"][0].keys())
dict_keys(['type', 'features', 'meta'])
dict_keys(['dataset_id', 'geometry', 'id', 'type', 'properties'])
# 取得した観測データの数
print(len(res_json["features"]))
# 0番目の観測データのプロパティを表示
for k, v in res_json["features"][0]["properties"].items():
print("{:<30} + {}".format(k, v))
143
sar:frequency_band + L
processing:level + L1.1
sat:relative_orbit + 62
start_datetime + 2008-11-24T01:20:19.965000+00:00
end_datetime + 2008-11-24T01:20:29.550000+00:00
tellus:name + ALPSRP150972890
created + 2023-02-22T07:46:41.062437+00:00
tellus:can_ordered + True
tellus:sat_frame + 2890
view:off_nadir + 21.5
palsar:beam + None
sar:observation_direction + right
sar:polarizations + HH+HV+VH+VV
sar:instrument_mode + P
tellus:published_datetime + 2023-06-19T04:20:22.255696+00:00
sat:orbit_state + descending
sar:product_type + SSC
gsd + 10
updated + 2023-06-19T04:20:22.271852+00:00
次に、可干渉な画像のペアを探していきます。
今回のSAR解析では、基本的な処理のみを行うため、本来とは異なり精度よく干渉できない場合が想定されます。そこで、画像のペアを検索する段階で、画像のペア候補ができるだけ多い条件で検索し、はっきりと干渉している結果が得られるように工夫しました。
なお、候補が多くなるような条件は、事前調査によって決定しました。
※本来は日時で指定した画像ペア同士で干渉の結果を確認すべきですが、今回はそこまで精緻な干渉処理を行いませんでした。
可干渉な画像のペアに該当する観測データの取得コードは以下になります。
data_list = []
for data in res_json["features"]:
prop = data["properties"]
if (
prop["sat:relative_orbit"] == 406 # 衛星の軌道経路
and prop["sat:orbit_state"] == "ascending" # 衛星の進行方向
and prop["sar:observation_direction"] == "right" # 電波照射方向
and prop["view:off_nadir"] == 34.3 # オフナディア角
and prop["tellus:sat_frame"] == 700 # 観測範囲の中心位置
and prop["sar:polarizations"] == "HH" # 送受信の偏波
and prop["sar:instrument_mode"] == "H" # 観測モード
):
data_list.append(data["id"])
for data in data_list:
print(data)
0d3a35fb-de02-4924-9040-75f1b13995af
bc71cf24-d76d-4f02-ab4c-19fdd261ff79
94bfcb3c-5f39-424b-a1cd-00a5b4f081b0
d8dd2301-7877-44e4-8cc7-408c92696b1a
2daca5c6-b11e-4d4c-9c30-94c430f9bba0
27f04b83-eb03-4cc8-8da5-4f533f15fc90
af9f02d3-bd11-4720-84b5-bbb368bc2098
8ca7d2ce-45fd-4df9-88b5-045cfde57b69
8d51fd4f-ff6c-4c76-854c-d2152cc0a90f
276c7b18-10f8-4c25-8244-437486c5c333
e1886ee8-662d-46d0-83a9-301012774efd
c5c693a3-de28-482c-b730-0b9a30f9dd19
b663d7b5-43d6-4dc9-96a0-7e7d74eef62b
344d9e3b-973a-4517-87e6-5c9af65c451d
2e465401-8671-4065-b37b-cee7f980accc
bfb14801-bac3-4055-8060-45eb2e76bae0
可干渉な画像の候補を取得できました。
次に、各観測データのファイル情報を見てみましょう。
観測データIDを指定し、二つ目のAPIを呼び出すことで、ファイル情報を取得できます。
dataset_id = "8836ec9a-35b5-43c3-92be-263498061e91" # PALSAR_L1.1のプロダクトID
data_id = data_list[0] # 0番目の観測データID
BASE_API_URL = "https://www.tellusxdp.com/api/traveler/v1/datasets/{}/data/{}/files/"
url = BASE_API_URL.format(dataset_id, data_id)
# 0番目の観測データのファイル情報を取得
response = requests.get(url, headers=HEADERS)
res_file_json = response.json()
for res_file in res_file_json["results"]:
print(res_file)
{'size_bytes': 138298, 'mime_type': 'image/jpeg', 'name': 'EALPSRP051270700.jpg', 'id': 1, 'status': 'uploaded', 'is_downloadable': True, 'require_archived_file_download': False}
{'size_bytes': 1385423568, 'mime_type': 'binary/octet-stream', 'name': 'IMG-HH-ALPSRP051270700-H1.1__A', 'id': 2, 'status': 'uploaded', 'is_downloadable': True, 'require_archived_file_download': False}
{'size_bytes': 1262, 'mime_type': 'text/plain', 'name': 'summary.txt', 'id': 3, 'status': 'uploaded', 'is_downloadable': True, 'require_archived_file_download': False}
{'size_bytes': 1800, 'mime_type': 'binary/octet-stream', 'name': 'VOL-ALPSRP051270700-H1.1__A', 'id': 4, 'status': 'uploaded', 'is_downloadable': True, 'require_archived_file_download': False}
{'size_bytes': 960720, 'mime_type': 'binary/octet-stream', 'name': 'TRL-ALPSRP051270700-H1.1__A', 'id': 5, 'status': 'uploaded', 'is_downloadable': True, 'require_archived_file_download': False}
{'size_bytes': 12510240, 'mime_type': 'binary/octet-stream', 'name': 'LED-ALPSRP051270700-H1.1__A', 'id': 6, 'status': 'uploaded', 'is_downloadable': True, 'require_archived_file_download': False}
Tellus Travelerの取得シーン詳細ページでも、ファイル情報を確認することができます。
これらのファイル中で今回必要とするのは、HH偏波の画像（IMG-HH-ALPSRP○○-H1.1__A）とリーダファイル（LED-ALPSRP○○-H1.1__A）です。
namekeyで該当するファイルを検索し、そのファイルIDを取得します。
def get_file_info(dataset_id, data_id, pol):
BASE_API_URL = "https://www.tellusxdp.com/api/traveler/v1/datasets/{}/data/{}/files/"
url = BASE_API_URL.format(dataset_id, data_id)
response = requests.get(url, headers=HEADERS)
res_file_json = response.json()
for i in res_file_json["results"]:
if i["name"].startswith(f"IMG-{pol}"):
img = {str(i["id"]): i["name"]}
elif i["name"].startswith("LED"):
led = {str(i["id"]): i["name"]}
files = [img, led]
return files
file_id_list = []
pol = "HH" # 送受信時の偏波
for data_id in data_list:
file_id_list.append(get_file_info(dataset_id, data_id, pol))
for file_id in file_id_list:
print(file_id)
[{'2': 'IMG-HH-ALPSRP051270700-H1.1__A'}, {'6': 'LED-ALPSRP051270700-H1.1__A'}]
[{'2': 'IMG-HH-ALPSRP104950700-H1.1__A'}, {'6': 'LED-ALPSRP104950700-H1.1__A'}]
[{'1': 'IMG-HH-ALPSRP212310700-H1.1__A'}, {'5': 'LED-ALPSRP212310700-H1.1__A'}]
[{'2': 'IMG-HH-ALPSRP265990700-H1.1__A'}, {'1': 'LED-ALPSRP265990700-H1.1__A'}]
[{'2': 'IMG-HH-ALPSRP272700700-H1.1__A'}, {'3': 'LED-ALPSRP272700700-H1.1__A'}]
[{'1': 'IMG-HH-ALPSRP111660700-H1.1__A'}, {'3': 'LED-ALPSRP111660700-H1.1__A'}]
[{'1': 'IMG-HH-ALPSRP111660700-H1.1__A'}, {'3': 'LED-ALPSRP111660700-H1.1__A'}]
[{'4': 'IMG-HH-ALPSRP118370700-H1.1__A'}, {'6': 'LED-ALPSRP118370700-H1.1__A'}]
[{'4': 'IMG-HH-ALPSRP118370700-H1.1__A'}, {'6': 'LED-ALPSRP118370700-H1.1__A'}]
[{'3': 'IMG-HH-ALPSRP158630700-H1.1__A'}, {'2': 'LED-ALPSRP158630700-H1.1__A'}]
[{'5': 'IMG-HH-ALPSRP205600700-H1.1__A'}, {'1': 'LED-ALPSRP205600700-H1.1__A'}]
[{'1': 'IMG-HH-ALPSRP219020700-H1.1__A'}, {'4': 'LED-ALPSRP219020700-H1.1__A'}]
[{'1': 'IMG-HH-ALPSRP219020700-H1.1__A'}, {'4': 'LED-ALPSRP219020700-H1.1__A'}]
[{'4': 'IMG-HH-ALPSRP225730700-H1.1__A'}, {'2': 'LED-ALPSRP225730700-H1.1__A'}]
[{'4': 'IMG-HH-ALPSRP225730700-H1.1__A'}, {'2': 'LED-ALPSRP225730700-H1.1__A'}]
[{'1': 'IMG-HH-ALPSRP245860700-H1.1__A'}, {'2': 'LED-ALPSRP245860700-H1.1__A'}]
表示したファイルIDの中で、プロダクト名ALPSRP111660700、ALPSRP118370700、ALPSRP219020700、ALPSRP225730700は重複していることが分かります。これらに対応する画像は除いてしまいましょう。
これで、可干渉な画像のプロダクトID、観測データID、ファイルIDが揃いました。
それではデータをダウンロードしましょう。
以下のコードでは、３つ目のAPIを用いて可干渉な観測データとして該当した全ての観測データのファイルをダウンロードしています。ちなみに今回は、ファイル名でデータを保存しているため、重複するデータは上書きされて、一方のデータのみダウンロードされたことになっています。
def save_file(dataset_id, data_id_list, files_list):
BASE_API_URL = "https://www.tellusxdp.com/api/traveler/v1/datasets/{}/data/{}/files/{}/download-url/"
for data_id, files_data in zip(data_id_list, files_list):
for file_data in files_data:
file_id =int(list(file_data.keys())[0])
file_name = list(file_data.values())[0]
url = BASE_API_URL.format(dataset_id, data_id, file_id)
# ダウンロードURLの発行
response_post = requests.post(url, headers=HEADERS)
dl_url = response_post.json()["download_url"]
# 現在のディレクトリ下にファイルをダウンロード
with open(file_name, "wb") as f:
f.write(requests.get(dl_url).content)
save_file(dataset_id, data_list, file_id_list)
作業ディレクトリ下に、ダウンロードしたファイルが存在することを確認してください。
HH偏波の画像（IMG-HH-ALPSRP○○-H1.1__A）がSAR画像、すなわち実際に地表を観測したデータです。また、リーダファイル（LED-ALPSRP○○-H1.1__A）が観測に関する詳細なデータです。
2. 強度・位相画像の出力
次にダウンロードしたデータを画像として出力してみましょう。
SARの観測範囲は広範にわたる、かつデータサイズが大きいため、先にデータを指定の範囲で切り出してから、画像として出力していきます。なお、SAR画像は電波の強度と位相の情報を持つ複素画像です。今回は強度画像と位相画像に分解してから出力します。
SARの観測範囲が広範であることを確認するために、例として観測データ「ALPSRP051270700」のサムネイル画像を表示してみます。
観測範囲は東京都を中心として、隣接県も含むかなり広範囲であることが分かります。
今回は、赤字で囲った東京ディズニーリゾート周辺のみが解析対象範囲なので、この部分を切り出したいと思います。切り出すにあたって、緯度・経度で指定される対象領域と取得した画像のピクセル座標をマッピングする必要があります。マッピングに用いる変換式の情報など、観測データに関する詳細な情報は、リーダファイルが有しています。
このリーダファイルから所望の値を取得する際に必要となるのが、プロダクトフォーマット説明書です。今回使用するPALSAR-1のL1.1 プロダクトフォーマットは、CEOSフォーマットに準拠しています。
CEOSフォーマット
CEOSフォーマットは、各国の宇宙機関や国際機関などが国際的な協力をするために組織した地球観測衛星委員会（Committee on Earth Observation Satellites）が定めたデータフォーマットです。
バイナリデータのどの位置にどのようなデータが収められているかは、衛星やセンサごとに異なるため、既存のソフトウェアやライブラリが対応していない観測データを利用するには各データプロバイダが公開しているドキュメントを読み込む必要があります。
PALSAR-2 プロダクトフォーマット説明書（JAXA）
科学者やプロ向けのデータで、衛星データに精通していない人が使うには少々手こずるかもしれません。
出所：宙畑 【ゼロからのTellusの使い方】提供元オリジナルデータを取得する
上記がCEOSフォーマットの説明になります。初心者がCEOSフォーマットを読み解くのはハードルが高いかもしれませんが、これから丁寧に読み解いていきましょう。
PALSAR-1 L1.1 プロダクトフォーマット
プロダクトフォーマット説明書には、各ファイル内のレコード名、バイト位置、レコード内容等が記述されています。
今回参照するプロダクトフォーマット説明書は、JAXAのALOS（だいち）のデータを使うというサイトの標準成果物(レベル1)フォーマット説明書のALOS処理プロダクトフォーマット説明書 PALSARレベル1.1/1.5編 L改訂版 平成21年7月です。
p.2-1にプロダクトフォーマットの概要が記載されています。特に、2.3 ファイル構成の表2-1には、各ファイルの種別や内容等が記載されています。
出所：ALOS処理プロダクトフォーマット説明書 PALSARレベル1.1/1.5編 L改訂版 平成21年7月
p.2-3の表2-2には、処理レベルごとのリーダファイルのレコード構成が記載されています。SARリーダファイルが9種類のレコードで構成されていることが分かります。今回は、データセットサマリレコード、プラットフォーム位置データレコード、ラジオメトリックデータレコード、設備関連データレコードのデータを使用します。L1.1では、地図投影データレコードがないことに注意してください。
また、SAR画像（SARイメージファイル）は、ファイルディスクリプタレコード、シグナルデータレコード、処理済データレコードで構成されていることが分かります。シグナルデータレコードに実際の観測データが格納されています。
出所：ALOS処理プロダクトフォーマット説明書 PALSARレベル1.1/1.5編 L改訂版 平成21年7月
p.2-4の表2-3には、各レコードのレコード長やレコード数が記載されています。
表2-3は、以降の処理でもたびたび参照します。
出所：ALOS処理プロダクトフォーマット説明書 PALSARレベル1.1/1.5編 L改訂版 平成21年7月
p.3-1にはリーダファイルの各レコードの記載ページが記されています。
PALSAR-1とPALSAR-2では、レコード構成やレコード長、レコード値が一部異なっていることに注意してください。
切り出す範囲の緯度・経度をピクセル座標に変換
緯度・経度からPALSAR-1の画像のピクセル座標に変換する場合、8次多項式を用いて変換します。変換多項式の各係数は、p.3-84のリーダファイルの設備関連データレコード11に格納されています。
リーダファイルから所望のレコードを取得する際は、ファイルポインタをファイルの先頭から指定位置まで動かし、そこからレコードを読み取ります。今回は変換係数を取得したいため、リーダファイルの設備関連データレコード11の先頭まで移動します。このとき、L1.1では地図投影データレコードが存在しないことに注意してください。具体的な数値は以下のようになります。（p.2-4の表2-3参照）
fpled = open(os.path.join("LED-" + data_name + "-H1.1__A"),mode='rb')
# PALSAR-1とPALSAR-2でレコード内容やレコード長が変更している点に注意
fpled.seek(720+4096+0+4680+8192+9860+1620+0+1540000+4314000+345000+325000+325000+3072+511000+4370000+728000+15000+2064, 0)
seek関数の第一引数について、補足します。
まず、ファイルディスクリプタレコードからキャリブレーションレコードまでのバイト値が次です。
720 + 4096 + 0 + 4680 + 8192 + 9860 + 1620 + 0
L1.1のため地図投影データレコードはなし、すなわち0です。
また、キャリブレーションレコードも0でした。こちらは、実際にキャリブレーションレコードの位置までポインタを移動させて、値を確認しました。校正がある場合は、こちらに値が入るため、ここは実際に確認するしかないように思います。今回は値が0だったことから、校正がなかったと考えられます。
次に、設備関連データレコード10までのバイト値です。
1540000 + 4314000 + 345000 + 325000 + 325000 + 3072 + 511000 + 4370000 + 728000 + 15000
p.3-26, p.3-27のリーダファイルディスクリプタレコードの各データレコード長を参照しました。
最後に、設備関連データレコード11の中で、緯度・経度からピクセル座標に変換する際の多項式の係数レコードの前までのバイト値です。
2064
p.3-79からp.83までを参照しました。
ここまでで、やっと所望のレコードまでポインタを動かすことができました。
プロダクトフォーマット説明書を読み解くのは、初心者にはなかなか大変ですね。。。
次に、p.3-84を参照すると、多項式の係数のレコード長は1000バイトであることが分かります。さらに、係数はc0, c1, c2, ･･･, c24及びd0, d1, d2, ･･･, d24の計50個であることから、一つの係数につき20バイト割り当てられていることが分かります。
出所：ALOS処理プロダクトフォーマット説明書 PALSARレベル1.1/1.5編 L改訂版 平成21年7月
実際には、緯度・経度からピクセル座標（p, l）への変換係数が、pとlごとに格納されていることから、500バイトずつ読み込みました。さらに、その次2つのレコードは原点緯度と原点経度であり、変換時の基準として使用します。
緯度・経度をピクセル座標に変換する関数は次のようになります。
def get_pl(fpled,latlon):
fpled.seek(720+4096+0+4680+8192+9860+1620+0+1540000+4314000+345000+325000+325000+3072+511000+4370000+728000+15000+2064, 0)
c = np.zeros(25,dtype="float64")
d = np.zeros(25,dtype="float64")
result = fpled.read(500) # バイトNo.2065-2564の中を20バイト分ずつ、計500バイト分進む
result_str = result.decode()
arr_str = result_str.split()
coeff_list = [float(num) for num in arr_str]
for i in range(25):
c[i] = coeff_list[i]
result = fpled.read(500) # バイトNo.2065-2564の中を20バイト分ずつ、計500バイト分進む
result_str = result.decode()
arr_str = result_str.split()
coeff_list = [float(num) for num in arr_str]
for i in range(25):
d[i] = coeff_list[i]
lat0 = float(fpled.read(20)) # 原点緯度(Φo)[度] 3065-3084を読み込み、20バイト分進む
lon0 = float(fpled.read(20)) # 原点経度(Λo)[度] 3085-3104を読み込み、20バイト分進む
phi = np.zeros(2, dtype="float64")
lam = np.zeros(2, dtype="float64")
phi[0] = latlon[0] - lat0
phi[1] = latlon[1] - lat0
lam[0] = latlon[2] - lon0
lam[1] = latlon[3] - lon0
pl = np.zeros(4,dtype="float64")
for i in range(5):
for j in range(5):
id = i*5+j
pl[0] += c[id]*lam[0]**(4-j) * phi[0]**(4-i)
pl[1] += c[id]*lam[1]**(4-j) * phi[1]**(4-i)
pl[2] += d[id]*lam[0]**(4-j) * phi[0]**(4-i)
pl[3] += d[id]*lam[1]**(4-j) * phi[1]**(4-i)
return pl
次に、指定した範囲で画像を切り出し、強度画像と位相画像として出力しましょう。
SAR画像の切り出しと強度・位相画像の出力
まず、SAR画像を指定された範囲で切り出します。
実際の観測データが格納されている、SARイメージファイルを参照します。
p.3-85からイメージファイルディスクリプタレコードを順に調べていき、変換に必要なレコードのバイト位置を確認します。p.3-92の「データセット（チャンネル）あたりのライン数（境界を除く）」や「 １ラインあたりのデータグループ（ピクセル）の数」というレコードが、SAR画像の縦横ピクセル値に関係する値です。また、SAR画像のピクセル値をレコードから読み取る際には、prefixが存在することに注意する必要があります。このprefixは、PALSAR-1では412バイト、PALSAR-2では544バイトです（p.3-94）。
こちらの説明は、宙畑 【コード付き】TellusでPALSAR-2のL1.1の 画像化にチャレンジしてみたを参照してください。以下のようなわかりやすい図付きで解説しています。
出所：【コード付き】TellusでPALSAR-2のL1.1の 画像化にチャレンジしてみた
次に、シグナルデータレコードからSAR画像のデータを読み込みます。
SAR画像は、レーダから地表に照射後、跳ね返ってきた電波です。つまり、SAR画像は電波の強度と位相という２つの成分を有しています。一般的にはまとめて扱いますが、分かりやすさを優先し、今回はこの２成分を分解して使用します。
レコードから取得したSAR画像のデータは、各ピクセル座標ごとに、実数部と虚数部を持つ1次元複素データです。そのため、SAR画像の縦横幅に合わせた複素数の2次元配列として変換します。
強度成分は、校正式に従って対数変換します。校正係数（CF）は、p.3-68のラジオメトリックデータレコードに記載があります。詳細は、ALOS-2 校正・検証を参照してください。位相成分は、複素数の角度部分になります。numpyのangle関数を使って、[-π~-π]にスケーリングします。
出所：ALOS処理プロダクトフォーマット説明書 PALSARレベル1.1/1.5編 L改訂版 平成21年7月
観測データからシグナルデータレコードを読み込み、強度・位相画像に変換する関数は以下になります。
def gen_img(fpimg,off_col,col,off_row,row):
# off_col,col,off_row,rowで、対象範囲を絞り込む
cr = np.array([off_col,off_col+col,off_row,off_row+row], dtype="i4")
# データセット（チャネル）当たりのライン数（境界を除く）
fpimg.seek(236)
nline = int(fpimg.read(8))
nline = cr[3]-cr[2] # row
# 1ライン当たりのデータグループ（ピクセル）の数
fpimg.seek(248)
ncell = int(fpimg.read(8))
# PALSAR-2の場合は、Prefixが544であることに注意
prefix = 412
nrec = prefix + ncell*8
# シグナルデータレコードの前まで移動
fpimg.seek(720)
# SAR画像データを取得し、2次元配列に変換
fpimg.seek(int((nrec/4)*(cr[2])*4))
data = struct.unpack(">%s"%(int((nrec*nline)/4))+"f",fpimg.read(int(nrec*nline)))
# 1次元配列を列数「nrec/4」で分割し、それぞれ行に変換する
data = np.array(data).reshape(-1,int(nrec/4))
# 全ての行に対して、列番号が「prefix/4」以降の全ての列を選択する
data = data[:,int(prefix/4):]
# 各行の偶数列の要素を実数部とし、奇数列の要素を虚数部として、複素数の配列を作成する
data = data[:,::2] + 1j*data[:,1::2]
data = data[:,cr[0]:cr[1]]
# 位相成分と強度成分に分解
CF = -83.0
CF_offset = 32.0
sigma = 10*np.log10(abs(data)) + CF - CF_offset
phase = np.angle(data)
# 画像の輝度を、0~255の範囲で正規化 & 8ビット符号なし整数のデータ型に変換する
sigma = np.array(255*(sigma - np.amin(sigma))/(np.amax(sigma) - np.amin(sigma)), dtype="uint8")
# ヒストグラムの均一化によってコントラストを向上させる
sigma = cv2.equalizeHist(sigma)
return sigma, phase
それでは、東京ディズニーリゾート周辺の画像を切り出し、強度画像と位相画像をそれぞれ表示してみましょう。画像の切り出しにあたっては、想定通りの範囲が切り出されなかったため、ピクセル値を調整しました。
# 東京ディズニーリゾート周辺
latlon = np.array([35.644252, 35.621438, 139.897408, 139.871744], dtype="float64")
data_names = [
"ALPSRP051270700", # 観測データ0
"ALPSRP104950700", # 観測データ1
"ALPSRP212310700", # 観測データ2
"ALPSRP265990700", # 観測データ3
"ALPSRP272700700", # 観測データ4
"ALPSRP111660700", # 観測データ5
"ALPSRP118370700", # 観測データ6
"ALPSRP158630700", # 観測データ7
"ALPSRP205600700", # 観測データ8
"ALPSRP219020700", # 観測データ9
"ALPSRP225730700", # 観測データ10
"ALPSRP245860700", # 観測データ11
]
for data_name in data_names:
fpimg = open(os.path.join("IMG-HH-" + data_name + "-H1.1__A"),mode='rb')
fpled = open(os.path.join("LED-" + data_name + "-H1.1__A"),mode='rb')
pl = np.array(np.ceil(get_pl(fpled, latlon)), dtype=np.int64)
off_col = min(pl[0], pl[1])
off_row = min(pl[2], pl[3])
# ピクセル値を調整して強度・位相画像を出力
sigma, phase = gen_img(fpimg, off_col-200, 1000, off_row-300, 1000)
sigma = np.rot90(sigma, k=3)
sigma = np.rot90(sigma, k=-1)
sigma = np.fliplr(sigma) # 画像が反転しているため
phase = np.rot90(phase, k=3)
phase = np.rot90(phase, k=-1)
phase = np.fliplr(phase) # 画像が反転しているため
# 強度画像(sigma)と位相画像(phase)をそれぞれグレースケールとjetカラーマップで表示
plt.imsave('sigma_{}.jpg'.format(data_name), sigma, cmap = "gray")
plt.imsave('phase_{}.jpg'.format(data_name), phase, cmap = "jet")
np.save('sigma_{}.npy'.format(data_name), sigma)
np.save('phase_{}.npy'.format(data_name), phase)
強度画像では、なんとなく建物や幹線道路が識別できていることがわかります。しかし、位相画像は全く視認性がなく、これだけでは何を示しているのかわかりません。位相画像は、干渉させることでその価値を発揮します。なお、強度画像が歪んで見えるのは、観測時の衛星と観測地点に角度がついているためです。
では、先程ダウンロードした全てのプロダクトの強度画像と位相画像を出力してください。
その後、いよいよ干渉SAR解析を実施してきます。
3. 画像の位置合わせと干渉処理
SAR画像に対する干渉は、以下の手順で行われます。
- ２つの強度画像から画像の位置合わせを行う
- ２つの位相画像の差分をとる
1.の位置合わせは、干渉結果の出来に直結するため非常に重要です。ただし今回は、必要最低限の位置合わせのみを行います。また、2.の位相画像の差分を取るという計算が、干渉になります。
こちらの詳細は、宙畑【コード付き】地盤の沈降が分かる干渉SAR解析をTellusでやってみたの「(3)画像の位置合わせと干渉処理」を参照してください。
画像の位置合わせと干渉処理を行う関数は次のようになります。
def coregistration(A,B,C,D):
# 画像の位置合わせを行う
# A, Bが強度画像、C, Dが位相画像
py = len(A)
px = len(A[0])
A = np.float64(A)
B = np.float64(B)
# 強度画像A, Bを用いて変位量を求める
d, etc = cv2.phaseCorrelate(B,A)
dx, dy = d # 画像の変位量
print(f"画像の変位量:{d}")
# 位相画像C, Dをその変位量でカットする
if dx < 0 and dy >= 0:
dx = math.ceil(dx)-1
dy = math.ceil(dy)
rD = D[dy:py,0:px+dx]
rC = C[dy:py,0:px+dx]
elif dx < 0 and dy < 0:
dx = math.ceil(dx)-1
dy = math.ceil(dy)-1
rD = D[0:py+dy,0:px+dx]
rC = C[0:py+dy,0:px+dx]
elif dx >= 0 and dy < 0:
dx = math.ceil(dx)
dy = math.ceil(dy)-1
rD = D[0:py+dy,dx:px]
rC = C[0:py+dy,dx:px]
elif dx >= 0 and dy >= 0:
dx = math.ceil(dx)
dy = math.ceil(dy)-1
rD = D[dy:py,dx:px]
rC = C[dy:py,dx:px]
return rC, rD
def wraptopi(delta_phase):
# 値の範囲が、[-π~+π] – [-π~+π] = [-2π~+2π]となるので、それを[-π~+π]に戻す（ラップする）
delta_phase = delta_phase - np.floor(delta_phase/(2*np.pi))*2*np.pi - np.pi
return delta_phase
def get_ifgm(data_A, data_B, AB):
file_name_sigma = "sigma_{}.npy"
file_name_phase = "phase_{}.npy"
sigma1 = np.load(file_name_sigma.format(data_A))
sigma2 = np.load(file_name_sigma.format(data_B))
phase1 = np.load(file_name_phase.format(data_A))
phase2 = np.load(file_name_phase.format(data_B))
print(sigma1.shape)
print(sigma2.shape)
coreg_phase2, coreg_phase1 = coregistration(sigma1, sigma2, phase1, phase2)
# 位相画像の差分を、[-π~+π]にラップする
ifgm = wraptopi(coreg_phase2 - coreg_phase1)
np.save('ifgm{}.npy'.format(AB), ifgm)
plt.imsave('ifgm{}.jpg'.format(AB), ifgm, cmap = "jet")
return
先程得た全てのSAR画像の組み合わせ（66通り）に対して、干渉処理を実行してみましょう。
combinations = list(itertools.combinations(range(len(data_names)), 2))
for combination in combinations:
AB = str(combination[0]) + "_" + str(combination[1])
data_A = data_names[int(combination[0])]
data_B = data_names[int(combination[1])]
get_ifgm(data_A, data_B, AB)
得られた干渉画像（以下インターフェログラム）の例を示します。
縞模様の部分は干渉によって生じる干渉縞です。縞模様が判別できないところは、主に水面部分です。水面は干渉性（コヒーレンス）が低く、干渉縞が生じにくい場所になります。また、位置合わせの精度が十分でないと、縞模様は曖昧になってしまいます。
以降の処理で、インターフェログラム同士の差分を取り、ノイズ（干渉縞）を除去します。ノイズ（干渉縞）の除去方法では、3枚のSAR画像が必要です。すなわち、〇と◇というSAR画像の地表の変化を解析したい場合、△というSAR画像を基準にして、それぞれ〇-△と◇-△という組み合わせでインターフェログラムを作成すれば、(〇-△) - (◇-△) = 〇-◇ となって、所望の組み合わせの干渉SAR解析ができます。
今回12種類のSAR画像をダウンロードしたため、66通りのパターンがあります。この中で、
- 観測時期が２年間程度離れていること
- 基準となる画像と組み合わせで作成したインターフェログラムの干渉縞がはっきりと見える（干渉性が高い）こと
という条件で探したところ、基準の観測データ名ALPSRP219020700に対して、ALPSRP104950700とALPSRP212310700がよさそうでしたので採用しました。
インターフェログラムには干渉縞や地表の変化など、いくつかの情報が含まれています。
そこで、今回知りたい地形の変化のみを抽出するために、干渉縞を削除しましょう。
4. 軌道縞の除去と地形縞の除去
干渉縞には、大きく「軌道縞」と「地形縞」が含まれます。
軌道縞は、２回の観測時の衛星の位置がずれることによって生じる縞模様です。
地形縞は、地形の起伏に由来する縞模様です。
軌道縞は軌道情報から作成し、除去します。地形縞の除去には、デジタル標高モデル（DEM）を用いるのが一般的ですが、今回は簡単のため、2枚のインターフェログラムの差分を取ることで代替します。
こちらの詳細は、宙畑 「【コード付き】地盤の沈降が分かる干渉SAR解析をTellusでやってみたの「(4)軌道縞を削除し地形の変化情報を抽出する」を参照してください。
出所：【コード付き】地盤の沈降が分かる干渉SAR解析をTellusでやってみた
軌道縞の作成には、軌道情報が必要です。中でも、観測中の衛星位置が重要な情報になります。衛星位置を取得する関数は次のようになります。
基本的は、リーダファイルとSARイメージファイルから必要な情報を取得しています。
def get_sat_pos(led,img):
# シグナルデータレコードのセンサ取得ミリ秒を取得する
img.seek(720+44)
time_start = struct.unpack(">i",img.read(4))[0]/1000
# データセットサマリレコードのレーダ波長(m)を取得する
led.seek(720+500)
lam = float(led.read(16)) #m
# プラットフォーム位置データレコードのデータポイント数を取得する
led.seek(720+4096+0+140)
position_num = int(led.read(4))
# プラットフォーム位置データレコードのポイント間のインターバル時間(秒=60)を取得する
led.seek(720+4096+0+182)
time_interval = int(float(led.read(22)))
# プラットフォーム位置データレコードの第一ポイントの通算秒を取得する
led.seek(720+4096+0+160)
start_time = float(led.read(22))
# データセットサマリレコードのシーンセンタ時刻を取得する
led.seek(720+68)
center_time = led.read(32)
Hr = float(center_time[8:10])*3600
Min = float(center_time[10:12])*60
Sec = float(center_time[12:14])
msec = float(center_time[14:17])*1e-3
center_time = Hr+Min+Sec+msec
time_end = time_start + (center_time - time_start)*2
# SARイメージファイルディスクリプタのデータセット(チャネル)当たりのライン数を取得する
img.seek(236)
nline = int(img.read(8))
time_obs = np.arange(time_start, time_end, (time_end - time_start)/nline)
time_pos = np.arange(start_time, start_time+time_interval*position_num, time_interval)
pos_ary = [];
# プラットフォーム位置データレコードのデータポイント位置ベクトル(x,y,z/m)を、データポイント数分取得する。
# データポイントとは、衛星が1枚の画像を撮影している間(約10秒間)の、衛星の観測点である。
# その10秒間で衛星は100km程度移動するため、適切に補完を行い、画像の各ピクセルが取得されたときの衛星の位置を特定する。
led.seek(720+4096+0+386)
for i in range(position_num):
for j in range(3):
pos = float(led.read(22))
pos_ary.append(pos)
led.read(66)
pos_ary = np.array(pos_ary).reshape(-1,3)
fx = scipy.interpolate.interp1d(time_pos,pos_ary[:,0],kind="cubic")
fy = scipy.interpolate.interp1d(time_pos,pos_ary[:,1],kind="cubic")
fz = scipy.interpolate.interp1d(time_pos,pos_ary[:,2],kind="cubic")
X = fx(time_obs)
Y = fy(time_obs)
Z = fz(time_obs)
XYZ = np.zeros(nline*3).reshape(-1,3)
XYZ[:,0] = X
XYZ[:,1] = Y
XYZ[:,2] = Z
return XYZ
上記のコードは、宙畑 【コード付き】Tellusで干渉SARをやってみた（DEM編）のコードを参考にして、各コードに説明を付け加えました。コードの概要については、「軌道情報を使って軌道縞を除去する」の章を参照してください。
次に、軌道縞を生成する関数は以下になります。
def get_orbit_stripe(pos1,pos2,pl,led):
# データセットサマリレコードのレーダ波長(m)を取得する
led.seek(720+500)
lam = float(led.read(16))
a = []
b = []
c = []
# 設備関連データ11のピクセルとラインを緯度と経度に変換する8次多項式の係数を取得（p.3-83）
# PALSAR-2とPALDAR-1では、設備関連データの番号が異なるので注意
led.seek(720+4096+0+4680+8192+9860+1620+0+(1540000+4314000+345000+325000+325000+3072+511000+4370000+728000+15000+1024), 0)
for i in range(25):
a.append(float(led.read(20)))
for i in range(25):
b.append(float(led.read(20)))
# 原点ピクセルと原点ラインを取得
for i in range(2):
c.append(float(led.read(20)))
npix = abs(pl[0]-pl[1])
nline = abs(pl[2]-pl[3])
orb = np.zeros((nline, npix))
for i in range(npix):
if i % 100 == 0:
continue
for j in range(nline):
px = i+pl[0]-c[0]
ln = j+pl[2]-c[1]
ilat = (a[0]*ln**4 + a[1]*ln**3 + a[2]*ln**2 + a[3]*ln + a[4])*px**4 + (a[5]*ln**4 + a[6]*ln**3 + a[7]*ln**2 + a[8]*ln + a[9])*px**3 + (a[10]*ln**4 + a[11]*ln**3 + a[12]*ln**2 + a[13]*ln + a[14])*px**2 + (a[15]*ln**4 + a[16]*ln**3 + a[17]*ln**2 + a[18]*ln + a[19])*px + a[20]*ln**4 + a[21]*ln**3 + a[22]*ln**2 + a[23]*ln + a[24]
ilon = (b[0]*ln**4 + b[1]*ln**3 + b[2]*ln**2 + b[3]*ln + b[4])*px**4 + (b[5]*ln**4 + b[6]*ln**3 + b[7]*ln**2 + b[8]*ln + b[9])*px**3 + (b[10]*ln**4 + b[11]*ln**3 + b[12]*ln**2 + b[13]*ln + b[14])*px**2 + (b[15]*ln**4 + b[16]*ln**3 + b[17]*ln**2 + b[18]*ln + b[19])*px + b[20]*ln**4 + b[21]*ln**3 + b[22]*ln**2 + b[23]*ln + b[24]
ixyz = lla2ecef(ilat*np.pi/180.0,ilon*np.pi/180.0,0)
r1 = np.linalg.norm(pos1[j+pl[2],:] - ixyz);
r2 = np.linalg.norm(pos2[j+pl[2],:] - ixyz);
orb[j,i] = wraptopi(2*np.pi/lam*2*(r2-r1));
return orb
# 緯度・経度高度から地球のECEF座標系(地球を基準とした直交座標系)のXYZを出力する関数
def lla2ecef(lat, lon, alt):
# 地球の長半径
a = 6378137.0
# 地球の偏平率
f = 1 / 298.257223563
# 第一離心率の2乗(e2)
e2 = 1 - (1 - f) * (1 - f)
# 子午線曲率半径
v = a / math.sqrt(1 - e2 * math.sin(lat) * math.sin(lat))
x = (v + alt) * math.cos(lat) * math.cos(lon)
y = (v + alt) * math.cos(lat) * math.sin(lon)
z = (v * (1 - e2) + alt) * math.sin(lat)
return np.array([x,y,z])
次に、インターフェログラムから軌道縞を除去する関数は以下になります。
def del_orbit_stripe(ifgm, orbit_stripe, X, Y):
# 軌道縞の画像サイズをインターフェログラムのサイズに合わせる
orbit_stripe = orbit_stripe[0:ifgm.shape[0], 0:ifgm.shape[1]]
# 軌道縞の除去
ifgm_sub_orb = wraptopi(ifgm - orbit_stripe)
plt.figure()
plt.imsave('ifgm_sub_orb{}_{}.jpg'.format(X, Y), ifgm_sub_orb, cmap = "jet")
np.save('ifgm_sub_orb{}_{}.npy'.format(X, Y), ifgm_sub_orb)
def get_SAR_sub_orbit(A,B):
plAB = np.array(np.ceil(get_pl(fpled[A], latlon)), dtype=np.int64)
off_colAB = min(plAB[0], plAB[1])
off_rowAB = min(plAB[2], plAB[3])
x = 1000
plAB = [off_colAB, off_colAB + x, off_rowAB, off_rowAB + x]
# pos[i]は画像iの衛星の座標
orbit_stripeAB = get_orbit_stripe(pos[A], pos[B], plAB, fpled[A])
# 元画像の変換に合わせて、軌道縞も回転・反転させる
orbit_stripeAB = np.rot90(orbit_stripeAB, k=3)
orbit_stripeAB = np.rot90(orbit_stripeAB, k=-1)
orbit_stripeAB = np.fliplr(orbit_stripeAB)
np.save('orbit_stripe{}_{}.npy'.format(A,B), orbit_stripeAB)
plt.imsave('orbit_stripe{}_{}.jpg'.format(A,B), orbit_stripeAB, cmap = "jet")
ifgmAB = np.load('ifgm{}_{}.npy'.format(A,B))
del_orbit_stripe(ifgmAB, orbit_stripeAB, A, B)
では実際に、インターフェログラムから軌道縞を除去してみましょう。
select_data_names = [
data_names[1], # ALPSRP104950700
data_names[2], # ALPSRP212310700
data_names[9], # ALPSRP219020700
]
pos = []
fpled = []
for i, data_name in enumerate(select_data_names):
fpimg = (open(os.path.join("IMG-HH-" + data_name + "-H1.1__A"),mode='rb'))
fpled.append(open(os.path.join("LED-" + data_name + "-H1.1__A"),mode='rb'))
pos.append(get_sat_pos(fpled[i], fpimg))
combinations = list(itertools.combinations(range(len(data_names)), 2))
for combination in combinations:
get_SAR_sub_orbit(int(combination[0]), int(combination[1]))
観測データ1と9の除去後のインターフェログラムと軌道縞
観測データ2と9の除去後のインターフェログラムと軌道縞
縞模様がきれいになくなり、いい感じのインターフェログラムが得られました。２つのインターフェログラムで全体の色調が異なっているのは、ピクセルの値が異なっているためです。ピクセルの値は、[-π~π]の範囲の周期的な値なので、一律な値を引くことで、色調は調整できます。
次に、地形縞を除去して、最終的な地表の変化のみを抽出する関数は以下になります。
def get_insar(pairs):
pair1 = pairs[0]
pair2 = pairs[1]
ifgm_pair1 = np.load('ifgm_sub_orb{}_{}.npy'.format(pair1[0], pair1[1]))
ifgm_pair2 = np.load('ifgm_sub_orb{}_{}.npy'.format(pair2[0], pair2[1]))
# 地形縞の除去
# 色彩調整用の + np.pi*2/4.0 および-ifgm_trueは画像ペアによって適宜修正
ifgm_true = wraptopi(ifgm_pair1[0:999, 0:999] - ifgm_pair2[0:999, 0:999] + np.pi*2/4.0)
ifgm_true = -ifgm_true
plt.figure()
plt.imsave('D-InSAR_{}_{}_{}_{}.jpg'.format(pair1[0], pair1[1], pair2[0], pair2[1]), ifgm_true, cmap = "jet")
plt.imshow(ifgm_true, cmap = "jet")
plt.colorbar()
では最後に、地形縞を除去しましょう。
pairs = [[1,9], [2,9]] # +2/4 -> 全体を-
get_insar(pairs)
遂に差分干渉画像が完成しました！
なお、元画像は観測角度の影響を受けて歪んでいるため、横方向に少し引き延ばした画像を表示しています。
強度画像と比べると、陸地と水辺の境界が適切に表現されているのが分かります。また、ディズニーシー周辺が濃い青色になっており、数cm程度地盤が沈下しているように見受けられます。ちなみに、観測データ1と2の観測年月日はそれぞれ2008年1月13日と2010年1月18日です。観測データ1と、そこから約2年後の結果を比較していると言えます。
さらに、別の観測データのペアの解析結果や他の機関が行った解析結果と比較してみましょう。
結果の比較
国土地理院の解析結果との比較
まずは、国土地理院の解析結果と比較してみましょう。
国土地理院 平成 23 年 高精度地盤変動測量 高精度地盤変動測量（干渉 SAR）
国土地理院の解析の比較期間は、2008年4月14日から2011年1月21日の約3年となっています。観測データ1と、そこから約3年後の結果を比較していると言えます。こちらの結果からでも、東京ディズニーリゾートの一部では、数cmの地盤沈下が起きているという結果が読み取れます。比較してみると、今回得られた結果は、国土地理院の解析結果にかなり近い結果が得られたと言えるのではないでしょうか。
上手く干渉するような画像のペアをわざと選んではいるものの、簡易的な処理だけでも良い結果が得られることを示せました。
比較期間が異なる場合との比較
次に、比較期間が異なる干渉SARの解析結果と比較してみましょう。
左の2010年3月5日から2010年4月20日の約1か月間の干渉SARでは、地表にほぼ変動がないことが分かります。このことから、観測データ1と2との解析結果における一部の地盤沈下が、特定の期間に生じた現象であることが確かめられました。
まとめ
少々長くなってしまいましたが、無事PALSAR-1の画像を干渉させることができました。
この記事を通して、干渉SAR解析に興味を持たれた方が、PALSAR-2など最新のSARデータを用いて解析していくきっかけになればと思います。
参考文献
- 宙畑 【コード付き】地盤の沈降が分かる干渉SAR解析をTellusでやってみた
- 宙畑 【コード付き】TellusでPALSAR-2のL1.1の 画像化にチャレンジしてみた
- 宙畑 【コード付き】Tellusで干渉SARをやってみた（DEM編）
- 宙畑 【図解】衛星データの前処理とは~概要、レベル別の処理内容と解説~
- 宙畑 【ゼロからのTellusの使い方】提供元オリジナルデータを取得する
- 国土地理院 SAR干渉画像の作成手順と見方
- 国土地理院 平成 23 年 高精度地盤変動測量 高精度地盤変動測量（干渉 SAR）
- ALOS処理プロダクトフォーマット説明書 PALSARレベル1.1/1.5編 L改訂版 平成21年7月
- TellusのAPIからPALSAR2の生データを取得する