SAR衛星画像のRAWデータ再生処理を実装してみる (ALOS1)
概要
SARの理解のため、以下のような再生処理のざっくりした実装をしました。
- ALOS/PALSARのL1.0の画像再生処理
- レンジドップラーの周波数領域での処理
- レンジマイグレーションは未実装
- 数式は以下の本を参考にしています。
『リモートセンシングのための合成開口レーダの基礎』
大内和夫
東京電機大学出版局 第2版 (2009/6/20)
(本記事で"大内本"と言っているのはこの書籍になります)
データの取得
ここで使うデータはAlaska Satellite Fasility (ASF)のサイトから無料でダウンロードしたものです。(要登録)
ASFのサイトで、再生処理前のLEVEL1.0のプロダクトを含むシーンを検索し、ダウンロードします。
LEVEL1.1や2.2は今回は使いませんが、結果の比較のためダウンロードしておくと良いと思います。
ここでは、 scene_name: 273720840 で説明していきます。(ASFリンク)
(Dataset: ©JAXA/METI ALOS PALSAR L1.0 2011. Accessed through ASF DAAC 6 Jan 2025.)
ダウンロードしたzipファイルを展開すると、以下のようなファイル群が得られます。
- ALPSRP273720840.l0.workreport
- IMG-HH-ALPSRP273720840-H1.0__A
- LED-ALPSRP273720840-H1.0__A
- TRL-ALPSRP273720840-H1.0__A
- VOL-ALPSRP273720840-H1.0__A
- workreport
これらはCEOS スーパーストラクチャというフォーマットに準拠した複数のファイルから構成されており、"IMG-HH..."がSARイメージのRAWデータになります。
再生処理に必要なパラメータの設定
リポジトリでは、プロダクトフォーマットを参考にして、データファイル群から画像再生に必要な情報を抽出して定義しています。抽出の処理については、別の場所で改めて書きたいと思います。(→ 書きました)
ベタ書きで書くと以下のようなパラメータの定義になります。
# imageファイルヘッダ領域長
N_HEADER = 720
# レンジラインのレコード長
RECORD_LENGTH = 21100
# レンジラインの実際のデータ数
RANGES_FULL = 20608
# レンジラインのprefixデータ領域の長さ
RANGE_PREFIX = 412
# アジマス方向のフルのデータ数
AZIMUTH_FULL = 35345
# 送信パルス幅
PULSE_DUR = 2.7000000000000002e-05
# サンプリング周波数
SAMPLE_FREQ = 32000000.0
# チャープ帯域
CHIRP_BANDWIDTH = 28000000.0
# チャープ率 (大内本でいうところの alpha/piで B/tau0 で計算)
RANGE_CHIRP_RATE = 1037037037037.037
# センサー取得開始時間
DATA_START_SEC = 46611.693
# Pulse Repetition Frequency
PRF = 2155.172
# 地球半径 (GRS80の長軸と短軸の平均としている)
R_EARTH = 6367444.6570500005
# 光速
VC = 299792458
# (アジマスライン0番目の)ニアレンジ長
NEAR_RANGE_0 = 848665
# レーダ波長
LAMBDA_RADAR = 0.2360571
# アンテナ直径
D_ANTN = 8.9
# オフナディア角
OFF_NADIR_ANGLE = 34.3
# 軌道方向
LINE_TIME_DIRECTION = 'ASCEND'
衛星の速度情報も再生処理に使いますが、各時刻における衛星位置・速度のデータをdataframeとして持っておきます。今回のシーンでの抽出後のcsvファイルはこちら。
satpos_df = pd.read_scv(f"satpos.csv")
としておきます。
信号データの読み込み
今回使うモジュール
import os
import numpy as np
import matplotlib.pyplot as plt
import cv2
import pandas as pd
from scipy.interpolate import interp1d
from pymap3d.ecef import ecef2geodetic
データの読み込み
# 今回の処理用のクロップ領域のパラメータを定義
RANGE_OFFSET = 0
N_RANGES = 7000
AZM_OFFSET = 13000
N_AZMS = 10000
image_filepath = "IMG-HH-ALPSRP273720840-H1.0__A"
f = open(image_filepath, "rb")
rb = f.read()
r = np.frombuffer(rb, dtype=np.ubyte)
# ヘッダー領域を省く
# h,w がそれぞれ アジマス方向、レンジ方向のデータ数の行列にする
a_r_mtx = r[N_HEADER:].reshape((AZIMUTH_FULL, RECORD_LENGTH))
# レンジ方向のデータのみを抽出
a_r_mtx = a_r_mtx[:, RANGE_PREFIX : RANGE_PREFIX + RANGES_FULL]
# レンジ行は実部0, 虚部0, 実部1, 虚部1, ... の順番で入っているので
# 実部と虚部を足し合わせて複素数にする
a_r_mtx = a_r_mtx[:, ::2] + a_r_mtx[:, 1::2] * 1j
# crop
a_r_mtx = a_r_mtx[AZM_OFFSET:AZM_OFFSET+N_AZMS, RANGE_OFFSET:RANGE_OFFSET+N_RANGES]
rawデータを強度として見てみます。
# 適当にスケールしている
a_r_mtx_save = np.abs(a_r_mtx) / np.abs(a_r_mtx).max() * 255
a_r_mtx_save[a_r_mtx_save > 255] = 255
plt.imshow(a_r_mtx_save.astype(np.uint8), cmap="gray", vmin=100, vmax=200)
cv2.imwrite(f"rawimage.png", a_r_mtx_save.astype(np.uint8))
アジマス0番目のレンジ行のデータ(実部と虚部)を見てみると、
fig, axes = plt.subplots(2, 1, figsize=(9, 9))
axes[0].plot(a_r_mtx[0,:].real)
axes[0].set_xlabel("Range Count")
axes[0].set_ylabel("Real")
axes[1].plot(a_r_mtx[0,:].imag)
axes[1].set_xlabel("Range Count")
axes[1].set_ylabel("Imag")
こんな感じになります。
レンジ圧縮
レンジのチャープ信号の取得
レンジのチャープ信号の取得は
- 式から直接計算する
- メタデータのチャープレプリカを使う
などの異なる方法がありますが、ここでは式から直接計算してみます。
チャープ信号は大内本(第二刷) P138の式(5.24)より
となります。
上の式で計算したチャープ信号をフーリエ変換して参照信号としています。
(大内本では定常位相近似でフーリエ積分を計算していますが、ここではnumpy.fft をするだけにしています。)
times = np.arange(0, PULSE_DUR, 1/SAMPLE_FREQ)
times = times-times.mean() # 時間の中心を0に
phase = np.pi*RANGE_CHIRP_RATE*times**2
chirp = np.exp(-1j*phase)
# edge_zero_pad
n_zero_pad = a_r_mtx.shape[1] - len(times)
if n_zero_pad %2 ==0:
chirp_padded = np.pad(chirp, (int(n_zero_pad/2), int(n_zero_pad/2)))
else:
chirp_padded = np.pad(chirp, (int(n_zero_pad/2)+1 , int(n_zero_pad/2)) )
# チャープ信号をずらす
# 信号成分を両端にしないと生成後の画像がずれる
chirp_padded = np.roll(chirp_padded, int(len(chirp_padded)/2 ))
fig, axes = plt.subplots(2, 1, figsize=(8, 4))
axes[0].set_title("Range Chirp Signal Time")
axes[0].plot(chirp_padded.real)
axes[1].plot(chirp_padded.imag)
axes[1].set_xlabel("Time Point")
# FFT
fchirp_padded = np.fft.fft(chirp_padded)
fchirp_padded_shift = np.fft.fftshift(fchirp_padded)
fig, axes = plt.subplots(2, 1, figsize=(8, 4))
axes[0].set_title("Range Chirp Signal Freq")
axes[0].plot(fchirp_padded_shift.real)
axes[1].plot(fchirp_padded_shift.imag)
axes[1].set_xlabel("Freq Point")
チャープ信号をずらしている部分がありますが、これは生成後の画像がずれてしまうため後から調節してこうなりました。信号の位相がどうなるかもう少し理解する必要がありますが、とりあえずこれでいきます。
画像はチャープ信号の時間ドメイン(上)と周波数ドメイン(下)での表示です。
横軸の単位は配列の何番目のポイントかというもので、時刻 [s] や周波数 [Hz] などにはなっていません。
レンジの相関処理
# レンジラインでのFFT
a_fr_mtx = np.fft.fft(mtx_sub, axis=1)
a_fr_mtx_shift = np.fft.fftshift(a_fr_mtx, axes=(1,))
# 相関処理の実行
comp = fchirp_padded_shift[np.newaxis,:] * np.conj(a_fr_mtx_shift)
# レンジ成分を逆フーリエ変換して時間領域に戻す
a_r_rangeimage = np.fft.ifft(comp,axis=1)
# 相関処理後、レンジ方向の画像の反転
# 生成画像の方向調整のために追加
a_r_rangeimage = np.fliplr(a_r_rangeimage)
できた画像を見てみます。
# 適当にスケールしている
a_r_rangeimage_save = np.abs(a_r_rangeimage) / np.abs(a_r_rangeimage).max()*255
a_r_rangeimage_save[a_r_rangeimage_save > 255] = 255
plt.imshow(a_r_rangeimage_save.astype(np.uint8), cmap="gray", vmin=40, vmax=80)
cv2.imwrite(f"{out_dir}/rangeimage.png", a_r_rangeimage_save.astype(np.uint8))
アジマス圧縮
概要
大内本の式(6.20)にアジマス参照信号の以下のような式があります。
上の式は以下のような形となります。
となります。
実効衛星速度の導入
衛星速度の情報はRAWデータセットからECR座標系での衛星速度(
以下の文献でも紹介されているような衛星の実効速度
論文:SAR image formation enhancement using effective velocity estimation method
概要を簡単に説明します。
さきほどのドップラー定数
上図で、
が衛星の実効的な速度であると考えることができます。
では
図では、衛星が地表面から高さ h の位置にあり、
が得られます。
今回、2次元平面を使ってこのような説明になっていますが、わかりにくかった場合は3次元の座標系で位置ベクトルをとって考えてもよいかと思います。
以上から、
となります。
実装
実装では、
このデータは各時刻カウントとそのときのECR座標系での衛星位置(X, Y, Z)と速度(dx, dy, dz)が記録されています。ただし、この時刻カウントは、アジマス時刻に対応する時刻ではありません。速度データの時間変化を補間して、アジマス時刻での速度を内挿値として計算する必要があります。
def get_psat_vsat_from_df(df, azm_offset, azm_length, vnorm=True):
rel_time = df["Time s"]
data_end_time = DATA_START_SEC + AZIMUTH_FULL/PRF
azm_times = np.arange( DATA_START_SEC, data_end_time , 1/PRF )
# sat postion
finterp_x = interp1d(rel_time, df["x"], kind='cubic')
finterp_y = interp1d(rel_time, df["y"], kind='cubic')
finterp_z = interp1d(rel_time, df["z"], kind='cubic')
sat_X = finterp_x( azm_times )
sat_Y = finterp_y( azm_times )
sat_Z = finterp_z( azm_times )
sat_lla = ecef2geodetic(sat_X, sat_Y, sat_Z)
# sat velocity
finterp_dx = interp1d(rel_time, df["dx"], kind='cubic')
finterp_dy = interp1d(rel_time, df["dy"], kind='cubic')
finterp_dz = interp1d(rel_time, df["dz"], kind='cubic')
sat_Vx = finterp_dx(azm_times)
sat_Vy = finterp_dy(azm_times)
sat_Vz = finterp_dz(azm_times)
sat_Vnorm = np.sqrt(sat_Vx**2 + sat_Vy**2 + sat_Vz**2)
sat_lla = np.array(sat_lla)[:,azm_offset: azm_offset + azm_length]
if vnorm is True:
sat_Vnorm = sat_Vnorm[azm_offset: azm_offset + azm_length]
return sat_lla, sat_Vnorm
else:
sat_V = np.array([sat_Vx, sat_Vy, sat_Vz])
sat_V = sat_V[:,azm_offset: azm_offset + azm_length]
return sat_lla, sat_V
lla, vsat = get_psat_vsat_from_df(satpos_df, AZM_OFFSET, N_AZMS, vnorm=True)
速度の取得が準備できたら、参照信号を計算していきます。
# 衛星の高度
satH = lla[2,:]
# アジマス周波数
azm_freqs = np.arange(0, PRF, PRF/ N_AZMS) - PRF/2
# アジマス時間
t_az = np.arange((N_AZMS))/PRF - N_AZMS/PRF/2
# レンジのサンプリング間隔
dr = VC / (2*SAMPLE_FREQ)
# レンジ距離
range_dis = NEAR_RANGE_0 + np.arange(N_RANGES)*dr
# 衛星の実効速度の計算
Rh = satH + R_EARTH
v_eff = (
np.sqrt((R_EARTH**2 + Rh[:, np.newaxis] ** 2 - range_dis**2) / 2)
* vsat[:, np.newaxis] / Rh[:, np.newaxis]
)
# ドップラースペクトラムの計算
beta = -2 * np.pi * (v_eff) ** 2 / LAMBDA_RADAR / range_dis
dplr_spec_phase = beta * t_az[:,np.newaxis]**2
dplr_spec = np.exp(-1j*dplr_spec_phase)
# ビーム幅 (合成開口幅)の計算
# 速度は平均値を使っている
fig, axes = plt.subplots(2,1, figsize=(12,4))
beam_width = NEAR_RANGE_0 * LAMBDA_RADAR/ D_ANTN
beam_time = beam_width/v_eff.mean()
beam_points = int(beam_time* PRF)
not_beam_points = N_AZMS - beam_points
# 合成開口の領域以外を0にする
dplr_spec_synth = dplr_spec.copy()
dplr_spec_synth[ :not_beam_points//2 ,:] = 0
dplr_spec_synth[ -not_beam_points//2: ,:] = 0
# 信号の位相の調整
dplr_spec_synth = np.roll(dplr_spec_synth, len(dplr_spec_synth)//2 ,axis=0)
# ドップラースペクトルのプロット
axes[0].plot(dplr_spec_synth[:, 0].real)
axes[0].set_ylabel("Real")
axes[1].plot(dplr_spec_synth[:,0].imag)
axes[1].set_xlabel("Point")
axes[1].set_ylabel("Imag")
下のような参照信号が得られます。
相関処理を実施していきます。
# 画像のアジマス成分を平均値で引く
azm_means = np.mean(a_r_rangeimage, axis=0)
a_r_rangeimage = a_r_rangeimage - azm_means[np.newaxis, :]
# 画像のアジマス両端は合成開口に入らないので0にする
a_r_rangeimage[:n_synth//2, :] = 0
a_r_rangeimage[-n_synth//2:, :] = 0
# アジマス方向のfft
fazm_image = np.fft.fft(a_r_rangeimage, axis=0)
fazm_image = np.fft.fftshift(fazm_image, axes=(0,))
# レンジマイグレーション(未実装)
# fazm_image = sardata.range_migraion(rangeimage_fa_r, beta)
f_dplr_spec_synth = np.fft.fft(dplr_spec_synth, axis=0)
f_dplr_spec_synth = np.fft.fftshift(f_dplr_spec_synth, axes=(0,))
comp = f_dplr_spec_synth * np.conj(fazm_image)
image = np.fft.ifft(comp, axis=0)
アジマス成分のオフセットの減算について
上の処理で、アジマス成分のオフセットを平均値で引いています。
このオフセット除去は、この後の処理で合成開口長に合わせてアジマスラインの両端の一部を0とする際に、オフセットがあるままだとフーリエ変換後に余分な周波数成分が入るの避けるために入れています。
信号0にするのは信号の折り返しを防ぐ意味があります。(生成してから折り返しを含む部分をカットしてもよいです。)
また、オフセットの減算をせずに、両端を0ではなくアジマスラインの平均値で置き換えるても結果は良くなりませんでした。むしろ周期的な模様が入って悪化しているようでした。
ただし、信号の絶対値を保持し対場合はこの方法の改善が必要になります。
結果を見てみます。
# 適当に規格化している
image_save = np.abs(image) / np.abs(image).max() * 15 * 255
print(image_save.max())
image_save[image_save > 255] = 255
plt.imshow(image_save.astype(np.uint8), cmap="gray")
cv2.imwrite(f"{out_dir}/image.png", image_save.astype(np.uint8))
一応は地形が見える程度の結像はできました。
レンジマイグレーション処理は未実装なので、細かいところを見ると像が若干カーブがかっているところも見えます。
画像が縦長になっているのは、アジマス方向、レンジ方向で1ピクセルに相当の距離が異なるためです。 簡単にですが、調整して見ます。
# シーン中心付近のレンジ距離、グランドレンジ距離の変化量を取得
center_ranges = range_dis[N_RANGES // 2 : N_RANGES // 2 + 2]
center_ranges_ground = center_ranges * np.sin(np.deg2rad(OFF_NADIR_ANGLE))
drange = center_ranges_ground[1] - center_ranges_ground[0]
dazm = 1/PRF*vsat.mean()
print(f"{drange=}, {dazm=}")
image_save_corr = cv2.resize(image_save, None, fx=1, fy=drange/dazm)
plt.imshow(image_save_corr.astype(np.uint8), cmap="gray")
cv2.imwrite(f"{out_dir}/image_corr.png", image_save_corr.astype(np.uint8))
シーン中心で取得した位置変化量を一律に適用しており、まだそこまで正確ではありません。こちらでL1.1の画像と比べても違いがありますが、ひとまずここまでとします。
まとめと課題
ALOSのPALSARのL1.0データからひとまず再生するというところまでをやりました。
以下のような気になる点はいくつかあり。もう少し調べていきます。
- 本来はアジマス周波数帯域の中心値を考慮すべきだが、今回は簡略化して 0 にしている。
- 信号のメタデータはレンジラインごと定義されているものもある。今回は一行で取得したパラメータを一律で使っていたので、ラインごとに定義すると良い部分があるかもしれない
- オフセットを引かずにアジマスの相関処理をする方法
- 画像の位相合わせの理解
- レンジマイグレーション処理
Discussion