Transcript 講習スライド（事後補足・修正付）

StaMPS講習会
解析処理 Step by Step
福島 洋・小澤 拓
2011/12/20-21
事後微変更済み by 福島
（最終更新日：2011/12/31）
解析データ：タンボラ火山
• Sensor: ALOS/PALSAR
• Path: 417, Frame:7020 (Ascending, HH,
Offnadir:34.3deg.)
• Date:
– 20070113, 20070228,
20070716, 20070831, 20071016, 20080302,
20080417, 20080602, 20080718, 20080902,
20081018, 20090118, 20090721, 20090905,
20091021, 20100308
ソフトインストール状況の簡易チェック
•
•
•
•
•
•
roi
make_raw_alos.pl
doris
mt_prep
matlab
(in matlab) which stamps
• あとは、途中で問題が発生した時点で解決しま
しょう。
以下のようにパスを確認し、パスが
通っていなければ修正してください。
• echo $FFTW_LIB
/Users/yo/FUKU/programs/ROI_PAC_3_0_1/ROI_PAC/NetInst
/fftw-111104-1407/lib
• echo $MY_BIN
…./ALOS_preproc/bin/i386
• echo $STAMPS
/Users/yo/FUKU/programs/StaMPS_v3.2.1
• echo $MY_SCR
$STAMPS/ROI_PAC_SCR
• echo $DORIS_SCR
$STAMPS/DORIS_SCR
• echo $MATLABPATH
$STAMPS/matlab
データ準備
• データ（path-frame-yyyymmdd.tar.gz）を
~/tambora/RAW 内にまとめて置いてください。
suppl.tar.gz（赤は講習会後変更部分）
•
まず、解凍＆展開してください（tar xzvf suppl.tar.gz）。ファイルを次に指定する場
所に移動させてください。ここに含まれているM関数ファイルを使えば、stampsは
Matlab本体のみ（ツールボックス不要）で処理することができます。
•
•
•
•
link_raw_pixel
$STAMPS/bin
prep_alos_raw_pixel
$STAMPS/bin
link_alos_frame_raw_yf
$STAMPS/bin
ps_est_gamma_quick.m
$STAMPS/matlab
（このファイルは、同名のファイルの247行目あたりを修正したものです）
gausswin.m
$STAMPS/matlab
dem_le.r4
~/tambora/DEM
•
•
以下は、自分でDEMを準備するときのために配布します。DEM作成用関数
prep_demcgiara を使うときは、すべてのM関数ファイルをmatlabのパスが通っている
場所に置いてください。
• prep_demcgiara.m, loadreal.m, savereal.m, cartprod.m, ind2subVect.m
ロケールの変更
• 不具合が生じないように、ターミナルのlocale
を英語に変えておく。（日本語環境では、処理
中に用いられるdateコマンドの挙動が異なる
ため。）
• 解析を行うターミナル上で以下のいずれかを
実行。
export LANG=en_US.UTF-8（bashの場合）
setenv LANG en_US.UTF-8（csh, tcshの場合）
link_raw_pixel
データの解凍とディレクトリの作成
• （PIXELアーカイブのファイル名フォーマットを読み込
むように、オリジナルのlink_rawというシェルスクリプ
トを改変してある。）
• このコマンドにより、解析ディレクトリの下にフレーム
番号（7020）が作成され、 その下にyymmddという名
前のディレクトリが作成される。その中に、RAWデー
タ のリンクが置かれる。
使用法：
link_raw_pixel data_dir analysis_dir
使用例：
link_raw_pixel /Users/yo/tambora/RAW /Users/yo/tambora
（参考）基線長
撮像日 垂直基線長（km）
20070113 0
20070228 1.585
20070716 1.5035
20070831 1.6267
20071016 1.3662
20080302 1.129
20080417 1.1105
20080602 0.7809
20080718 1.5286
20080902 1.815
20081018 2.3616
20090118 2.1158
20090721 1.6222
20090905 1.7799
20091021 1.6994
20100308 1.3613
step_slc_alos
マスター画像のSLC作成
• 適当なマスターデータを決める。今回は、お
およそBperpおよび観測日が中間値となる
データとして、20080718とする。
• そこのディレクトリに入り、step_slc_alosを実
行。
cd ~/tambora/7020/SLC/20080718
step_slc_alos
注）ALOS_preprocの`setenv MY_BIN (hoge)` (tcsh) 、`export MY_BIN=(hoge)`
(bash) を忘れずに！
SLCの強度画
像
（ルック数
5 x 15）
レーダ座標系なので、
上下逆であることに
注意。
強度画像image.slc_5l.rasを開く
• GIMP (Linux/Mac) 、ToyViewer (Mac) など
• ここでやること：解析範囲（最終的な干渉画像
の範囲）を決定する
• そのために、まずras画像中で切り出したい領
域の画素番号をしらべる（今回は、X:281-550,
Y:1201-1400あたりを切り出すことにする）
• シングル・ルックの場合にどうなるかを計算。
X, Yそれぞれ5と15をかけた値をメモる。
../roi.proc の編集
切り出し領域を設定
• roi.procの中身
use1dopp=1
mean_pixel_rng = 2075 ←レンジ方向の解析領域中心
ymin = 17000 ←アジマス方向の解析領域下限
ymax = 22000 ←アジマス方向の解析領域上限
• アジマス方向のライン数は、1000以上余裕を与えて
設定する。つまり、最終的な解析領域を基準に、ymin
は1000以上小さい値、ymaxは1000以上大きい値に
設定する。
• ファイルの編集が終わったら、再度step_slc_alosを実
行。切り出された領域でSLCが作成される。
image.slc_5lk.ras
アジマス方向に切り出された画像ができる。
master_crop.inの作成・編集
最終的な解析領域を設定
• cp $MY_SCR/master_crop.in .
• vi master_crop.in
（もちろんviでなくgedit等でも可）
first_l 1501
last_l 4000
first_p 1501
last_p 2750
• first_lとlast_lは、アジマス方向のライン数。切り出され
た後のライン数を基準に設定。
• （元々切り出すと決めた領域は1000,4000,1401,2750
に相当するが、解析領域を絞るため若干変えた）
step_master_setup
マスターの準備
cd ~/tambora/7020/SLC/20080718
step_master_setup
• パラメタファイルの作成などとともに、
20080718_crop.ras ができる。この画像が、最終的な
解析領域。解析したい領域が適切に切り出されている
かチェックし、必要に応じmaster_crop.inを編集してや
り直す。
• もし、20080718_crop.rasが真黒な場合、切り出しに失
敗している。master_crop.inで設定した領域がはじめ
に切り出した領域外になっているとこうなる。
make_slc_alos
すべてのスレーブ画像のSLCを作成
cd ~/tambora/7020/SLC
make_slcs_alos
• しばらく時間がかかる。
• 各データのディレクトリ内にras形式の強度画像
が作成されているので、 問題なくできているかを
確認する。
• Toyviewerを使うなら、rasファイルはToyViewerで
開くことを定義しておき、open */*.rasとやると簡
単。
Tips
• StaMPSのシェルスクリプトには、step… という
名前と make… という名前のものがある。
• step… のスクリプトは、あるひとつの画像に対
する処理をおこなう。
• make… は、すべての画像に対する処理をお
こなう。step… を順番に実行している。
DEMの準備
• この講習会では、~/tambora/DEM に置いた
dem_le.r4 を使う。
• このDEMは、hole-filled SRTM DEM を切り出し、
オーバーサンプリングしたものです。StaMPS
では、ジオコードした干渉画像にPS処理を適
用するため、ピクセル間隔の細かいDEMを用
意する必要がある。
• hole-filled SRTM: http://srtm.csi.cgiar.org/
参考：DEMの準備（小澤メソッド）
• Photoshopなどでtifファイル（srtm_60_14.tif）を読み出し、汎
用フォーマット（RAW（バイト順序：IBM PC）で書き出すことに
より、2バイトベタ書きの形式に変換する。そこで、GMTを使っ
て、以下のようにデータを切り出しする（２バイトベタ書き
フォーマットのものをsrtm_60_14.raw、切り出すDEMファイル
の名前をdem_le.r4とする）。
xyz2grd srtm_60_14.raw -Gtemp1.grd -R115/120/-10/-5 -I3c -Zh -fg -V
grdmath -V temp1.grd -10 GE temp1.grd x = temp2.grd
grdsample temp2.grd -Gtemp3.grd -R117.7/118.3/-8.5/-8-7 -I0.75c -V
grd2xyz temp3.grd -Zf > dem_le.r4
参考：DEMの準備（福島メソッド）
• まず、suppl.tar.gz に含まれている loadreal.m、
prep_demcgiara.m を各自のmatlabパスに保存。
• 次に、下のようにmatlab上で実行。wgetかcurlがインストール
＆セッティングされていなければなりません（which wgetや
which curlとやってみてください）
lim = [117.7,118.3,-8.5,-8.0];
[width,height,ullat,ullon,dlat,dlon] = prep_demcgiara(lim,'online','dem.r4');
a = loadreal('dem.r4','float32',width,'b');
a = interp2(a,2);
savereal('dem_le.r4','float32',a,'l');
b = loadreal('dem_le.r4','float32',(width-1)*4+1,'l');
whos b
format long
dy = dlat./4
dx = dlon./4
ymax = ullat
xmin = ullon
← などとやると、必要なパラメタが得られます
timing.dorisin の編集
DEMパラメタの設定
• ~/tambora/7020/INSAR_20080718内に
timing.dorisinというファイルがある。これはDORIS
（M_SIMAMP, M_TIMING）用のパラメータファイル。
この中のDEMに関するパラメータを書き換える。
SAM_IN_FORMAT real4 // Type of format
SAM_IN_DEM /Users/yo/tambora/7020/DEM/dem_le.r4 // Path to
DEM file
SAM_IN_SIZE 2397 2877
//rows cols
SAM_IN_DELTA 0.00020833333 //in degrees
SAM_IN_UL
-8.000416667 117.7004167 //lat and lon of upperleft
SAM_IN_NODATA 0 // No data value
SAM_OUT_FILE master_sam.raw // synthetic amplitude
SAM_OUT_DEM dem_sam.raw // cropped DEM
step_master_timing
Timing（アジマス方向）Errorの推定
• DEMと位置合わせをすることにより、軌道
データから推定されるアジマス方向の位置の
ずれを推定。
cd ~/tambora/7020/INSAR_20080718
step_master_timing
（続き）以下のような出力が得られれば正常
Number posL posP offsetL offsetP correlation
0 2318 1782 -8 -6 0.467845
1 2318 1971 -8 -6 0.450781
2 2318 2160 -8 -6 0.436478
.
.
.
27 2681 2089 -8 -6 0.450372
28 2681 2278 -8 -6 0.426846
29 2681 2467 -8 -6 0.418827
Estimated total offset (l,p): -8, -6
make_orbits
INSAR_20080718のディレクトリ内に各ス
レーブ用のディレクトリが作成される
make_coarse
大雑把な位置合わせ
cd ~/tambora/7020/INSAR_20080718
make_orbits
make_coarse
make_coreg
精密位置合わせ
• このステップは時間がかかるため、時間節約のため
にパッチ数を調整してから実行する。
• パッチ数の変更：$STAMPS/matlab/coreg_pos.mのファ
イルを編集。nWinxとnWinyをそれぞれ20と50に変更
する。
• 注）MATLABPATHの設定`setenv MATLABPATH
/home/hoge/StaMPS_v3.2.1/matlab`を忘れずに！
• このあとで以下を実行。
cd ~/tambora/7020/INSAR_20080718
make_coreg
make_dems
地形縞のシミュレート画像を作成
• マニュアルによると、地形の起伏が大きくなけ
れば、実行しなくてよいとのこと。
• refdem_11.raw などのファイルが各スレーブ
に対して作られる。
cd ~/tambora/7020/INSAR_20080718
make_dems
make_resample
スレーブSLCのリサンプリング
• マスター画像と同じグリッドのスレーブ画像が
作成される
• このあと、リサンプルされた画像のファイルサ
イズがすべて同じになっているかを確認する。
cd ~/tambora/7020/INSAR_20080718
make_resample
ls –l */slave_res.slc
make_ifgs
干渉画像の作成
• 干渉画像 cint.minrefdem_5l.ras が各スレー
ブディレクトリ内にできる。
• できた干渉画像は、gimpやToyViewerなどで
確認する。（open */cint.minrefdem_5l.ras）
• 全く縞が確認できないようなものは、失敗して
いる可能性が高い。本演習の例題の場合、
処理がうまくいけば干渉性は高くなるはず。
step_geo
ジオコーディング
• これは、後のPS処理のために必要。 スレーブ
画像のディレクトリのどれか一つだけに移動
し、コマンドを実行する。
• 実行後、InSAR_20080718にlat.raw、lon.raw
が作られる。
cd ~/tambora/7020/INSAR_20080718/20070113
step_geo
ここからPS処理
mt_prep
PS候補点の選定
cd ~/tambora/7020/INSAR_20080718/
mt_prep 0.4
（mt_prepのオプション）
mt_prep 0.4 2 2 50 200
（意味）
0.4: Amplitude dispersionのしきい値（これ以下を用いる）
2: レンジ方向のパッチ数（デフォルト1）
2: アジマス方向のパッチ数（デフォルト1）
50: レンジ方向のパッチ間のオーバーラップピクセル数（デフォルト50）
200: アジマス方向のパッチ間のオーバーラップピクセル数（デフォルト
200）
StaMPSの処理はパッチごとに行われ、あとでマージされる。解析エリアに
対してメモリに余裕があれば、パッチに分ける必要はない。
以下の処理は、matlab上で実行す
る。matlabを起動して、処理ディレク
トリに移動する。その後、処理を実
行するという流れが基本。
(getparm)
処理パラメタの確認
• ~/tambora/7020/INSAR_20080718/ で、
>> getparm
と打つと、処理パラメタが出力される。
• コマンド setparm でパラメタを変更することがで
きる。通常、変更する必要はない。パラメタの説
明は、マニュアルには書いてない。Hooper et al.
(2007)のアルゴリズムの説明を読めば、だいた
いは直感的に理解できる。
stamps
PS処理の実行
• 何も指定しなければ、全ステップの処理を自
動的にしてくれる。
• stampsの処理は８つのステップに分かれてお
り、オプションを指定することにより各々のス
テップのみを実行することも可能。詳しくは、
help stamps 参照。
>> cd ~/tambora/7020/INSAR_20080718
>> stamps
ps_plot
結果のプロット・エクスポート
• このmatlab関数には、非常の多くの機能があ
る。詳しくは、 help ps_plot や、StaMPSマニュ
アルの Chapter 9 参照。
ps_plot(‘u-dm’,0,0,1) の結果
ps_plot 使用例：kmlファイル出力
>> ps_plot('v-d',-1)
>> load ps_plot_v-d
>> ps_gescatter('tambora-psinsar.kml',ph_disp,10,0.4)
Tips:
• ps_plotの２つ目の変数に-1を入れると、プロットされ
るはずのデータがmatファイルにエクスポートされる。
matファイルの名前は、ps_plot_処理方法指定.mat
となる。
• この例では、一度エクスポートしたものをロードし、
ロードされたph_dispという変数をさらにkmlファイル
にエクスポートしている。
Google Earthで見た結果（平均速度）
LOSベクトルの計算
• heading角と入射角が必要。
• heading角は、stamps処理の最中に、
heading.1.in という名前のファイルに書き出され
る。ROI_PACのパラメタファイル（*.rsc）にも書いて
ある。
• 入射角は、干渉画像が置かれているスレーブの
ディレクトリ（INSAR_20080718/yyyymmdd、どれ
でもよい）にあるdorisのファイルcoreg.out,
dem.out, geocode.outに書かれている。
（grep incidence coreg.out などとやると簡単）
その他諸々
• ps_plotには、時系列をプロットするオプション
がある。少し試してみたところ、プロット領域を
大きくしたときは異常終了してしまった。
• データのエクスポートは、ps_plotだけでなく、
ps_outputでもできる。
• ref_lon, ref_lat（またはref_centre_lonなど）で
不動域を指定することができます（ちゃんとで
きるかは未確認です）
参考文献など
•
•
•
•
•
•
•
StaMPSマニュアル（StaMPSパッケージに含まれる。下のウェブサイトからも入手可能）
StaMPSウェブサイト
http://radar.tudelft.nl/~ahooper/stamps/index.html
ユーザーグループ（過去ログの検索可）
http://groups.google.com/group/mainsar
StaMPSアルゴリズム
Hooper, A., P. Segall, and H. Zebker (2007): Persistent scatterer InSAR for crustal
deformation analysis, with application to Volcan Alcedo, Galpagos, J. Geophys. Res., 112,
B07407, doi:10.1029/2006JB004763.
StaMPSを用いた解析処理の流れの説明
福島 洋(2011): StaMPS パッケージを用いた PS 干渉 SAR 解析, 測地学会誌, 57, 41-48.
PSとSB解析の融合アルゴリズム
Hooper, A. (2008): A multi-temporal InSAR method incorporating both persistent
scatterer and small baseline approaches. Geophys. Res. Lett., 35:L16302.
アンラッピングアルゴリズム
Hooper, A. (2009): A statistical-cost approach to unwrapping the phase of InSAR time
series, Proceedings FRINGE Workshop, Frascati.