🛰️

SAR衛星画像のRAWデータ再生処理を実装してみる (ALOS1)

2025/02/01に公開

https://github.com/hiydee/sar_raw_process

概要

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)より

E_t(\tau) = E'_0\exp(i\pi f_c \tau + i \alpha t^2)

f_c はオフセット周波数で、ここでは 0としています。
\alphaについてですが、\frac{\alpha}{\pi} がチャープ率に相当します。 処理上ではチャープ率を変数 RANGE_CHIRP_RATE (\alpha') としていて、結局、

E_t(\tau) = E'_0\exp(i \pi \alpha' t^2)

となります。

上の式で計算したチャープ信号をフーリエ変換して参照信号としています。
(大内本では定常位相近似でフーリエ積分を計算していますが、ここでは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)にアジマス参照信号の以下のような式があります。

E_r(t) = rect(t/T_A) \exp(i2\pi f_{DC}t -i\beta t^2)

rect(t/T_A)は矩形関数で、実装上はpaddingで調整します。f_{DC}はドップラー中心周波数ですが、今は0としています。(厳密には計算が必要だと思います。)
上の式は以下のような形となります。

E_{doppler} = e^{-i \beta t^2}

\betaはドップラー定数で

\beta = -2 \pi \frac{v^2}{\lambda R}

となります。

実効衛星速度の導入

衛星速度の情報はRAWデータセットからECR座標系での衛星速度(V_{x}, V_{y}, V_{z})から得られますが、このノルム v_s = \sqrt{{V_x}^2 + {V_y}^2 + {V_z}^2} をそのまま使っても結像はうまくいきません。
以下の文献でも紹介されているような衛星の実効速度 v_r を導入することで、アジマス方向の結像がより鮮明になります。
論文:SAR image formation enhancement using effective velocity estimation method

概要を簡単に説明します。
さきほどのドップラー定数 \beta の式は、衛星から散乱体の距離(スラントレンジ) r の時間変化から来ています。この r の計算の際に、衛星の軌道上の速度に対応するレーダーの地表面上での照射位置の変化を考慮に入れます。

上図で、R_0 は衛星が散乱体の真横に来た時のレンジ距離になります。図中央の考え方は最もシンプルで、スラントレンジを三平方の定理から求めています。一方、図の右側の表式では、地表の球面に沿った形で衛星が軌道上を動いていると仮定しています。この場合、レーダーが照射される位置も衛星の速度に応じて変化します。この変化 v_g を考慮すると、近似的に、\sqrt{{R_0}^2 + (v_st) (v_gt)} と計算できます。これを中央の表式 \sqrt{{R_0}^2 + (v_st)^2} と比べると、

v_r = \sqrt{v_s v_g}

が衛星の実効的な速度であると考えることができます。

では v_g はどのように計算すれば良いでしょうか。以下のように考えます。

図では、衛星が地表面から高さ h の位置にあり、R_eは地球の半径になります。 ここで衛星は自身から R_0 の距離にある地表面にレーダーを照射しています。衛星はこの図を奥に突き抜ける方向に v_s で移動しています。 地球中心の座標系で衛星の高さ方向に軸を取ると、v_g は高さ成分に応じて v_s からの寄与を受けると考えます。高さR_h で衛星速度が v_s で動いており、高さ成分が R_e\cos{\alpha} で 地表点が v_g で動いているという関係を相似形として、

v_g = \frac{R_e}{R_h} v_s \cos{\alpha}

が得られます。
\cos{\alpha} は、図から R_e, R_h, R_0 を用いて余弦定理で求めることができます。

今回、2次元平面を使ってこのような説明になっていますが、わかりにくかった場合は3次元の座標系で位置ベクトルをとって考えてもよいかと思います。

以上から、v_r = \sqrt{v_s v_g} を計算すると、

v_r = \frac{v_s}{R_h} \sqrt{\frac{ {R_e}^2 + {R_h}^2 - {R_0}^2 }{2} }

となります。

実装

実装では、v_{s} は各アジマスラインの時刻 (アジマス時刻) に対応する衛星の速度データを使うことになります。はじめに satpos_df として用意しておいたものです。
このデータは各時刻カウントとそのときの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