Transcript StaMPS*** **** Step by Step

StaMPS講習会
解析処理 Step by Step
（講習会後改訂版）
福島 洋・小澤 拓
2012/9/11-12
（最終更新日：2012/9/14）
1
StaMPS解析のコツ
• ひとつひとつのステップを確実に！間違えた処理をし
た場合、その正しくステップをやり直しても失敗するこ
とがあるため、最悪、最初からやり直ししなければなら
なくなります。
• このスライド中、「正常な処理が行われればこうなる」
と書いてあるところは特に重要ポイントです。そのス
テップの結果を注意して確認してください。あと、
Warning表示は無視してよい場合が多いですが、
ERROR表示が出たら何かがおかしい可能性が大きい
です。
• このスライドの通りにやったのにうまくいかない場合、
まずマニュアル（StaMPS_Manual_v3.2.pdf）を参照して
ください。うまくいかなかった場合にどうするか、再処
理のやり方（特に5.2と5.3）などが書いてあります。
2
解析データ：九重火山
• Sensor: ALOS/PALSAR
• Path: 73, Frame:2950 (Descending, HH,
Offnadir:34.3deg.)
• Dates:
– 20070107, 20071125, 20080527, 20090530,
20100417, 20110420 の６枚
3
九重
阿蘇
4
ソフトインストール状況の簡易チェック
• 以下のコマンドを打ち込み、しかるべき反応が返ってく
るかチェック
–
–
–
–
–
–
–
–
–
roi
講習会中、追加インストールが必要に
ALOS_fbd2fbs
なったパッケージ：gawk, csh, tcsh
make_raw_alos.pl
doris
mt_prep
matlab
(in matlab) which stamps
(in matlab) !triangle ← stampsの正常実行のため必要
(in matlab) !snaphu ← stampsの正常実行のため必要
5
以下のようにパスを確認し、パスが
通っていなければ修正してください。
• 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
6
データ・スクリプトの準備
•
•
•
データ（path-frame-yyyymmdd.tar.gz）を ~/kuju/RAW 内にまとめて置いてください。
4931_ds.dehm を ~/kuju/DEM に置いてください。
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
– make_cohs
$STAMPS/bin
– step_coh_sb
$STAMPS/bin
– ps_est_gamma_quick.m
$STAMPS/matlab
（このファイルは、同名のファイルの247行目あたりを修正したものです）
– mt_ml_select.m
$STAMPS/matlab
– repair_lonlat_5lraw.m
$STAMPS/matlab
– loadreal.m, savereal.m
$STAMPS/matlab
– gausswin.m
$STAMPS/matlab
7
ロケールの変更
• 不具合が生じないように、ターミナルのlocale
を英語に変えておく。（日本語環境では、処理
中に用いられるdateコマンドの挙動が異なる
ため。）
• 解析を行うターミナル上で以下のいずれかを
実行。
export LANG=en_US.UTF-8（bashの場合）
setenv LANG en_US.UTF-8（csh, tcshの場合）
これを確実にやらないとハマります！
特に、複数のターミナルを開く場合に注意。
8
link_raw_pixel
データの解凍とディレクトリの作成
• （PIXELアーカイブのファイル名フォーマットを読み込
むように、オリジナルのlink_rawというシェルスクリプ
トを改変してある。）
• このコマンドにより、解析ディレクトリの下にフレーム
番号（2950）が作成され、 その下にyymmddという名
前のディレクトリが作成される。その中に、RAWデー
タ のリンクが置かれる。
使用法：
link_raw_pixel data_dir analysis_dir
使用例：
link_raw_pixel /Users/yo/kuju/RAW /Users/yo/kuju
（注意：絶対パスで与えること）
9
（参考）基線長
撮像日 垂直基線長（km）
2007/1/7
0
2007/11/25 -0.1407
2008/1/10 -0.3911
2008/2/25 -0.9280
2008/4/11 -1.2814
2008/5/27 -1.6176
2008/7/12
0.9758
2008/10/12 3.0673
2008/11/27 2.8823
2009/1/12
2.6156
2009/4/14
1.7866
2009/5/30
1.1259
2009/8/30
1.8113
2009/10/15 1.1608
2009/11/30 0.9332
2010/1/15
0.4644
2010/3/2
0.0161
2010/4/17 -0.4705
2010/12/3 -1.3912
2011/1/18 -1.6955
2011/3/5
-2.1100
2011/4/20 -2.6241
◯
◯
◯
◯
◯
◯
10
step_slc_alos
マスター画像のSLC作成
• 適当なマスターデータを決める。マニュアルには、おおよそBperpおよび
観測日が中間値となるデータをマスターとするように書いてある。今回
は、20090530とする。注：マスターを変えると、（多分画像枚数が少なけ
れば少ないほど）結果が変わります。
• 20090530ディレクトリに入り、step_slc_alosを実行。
cd ~/kuju/2950/SLC/20090530
step_slc_alos
注）ALOS_preprocの`setenv MY_BIN (hoge)` (tcsh) 、`export MY_BIN=(hoge)`
(bash) を忘れずに！
以下のようなメッセージが最初に表示されるが気にしない。
ls: No match.
ls: IMG-HH-ALP*1.0__A: No such file or directory
11
Mac固有の問題
• step_slc_alosがうまくいかない（20070107.raw
not found のようなエラーメッセージが出る）ケー
スが講習会で複数あった。
• ROI_PACのインストールがうまくいっていないた
めと思われる。
• 現在のところ（2012/9/11）、問題の解決方法は
見つかっていない。
• http://www.roipac.org/Installation の下のほうに、
Mac上でのコンパイルに関する情報が載ってい
る。ここを参考に（できれば複数の）コンパイラを
試し、まずはROI_PACで配布されているテスト
データでちゃんと解析できるようにする。
12
Linux固有の問題
• Ubuntu Linux 11を32 bitマシンの場合で、
step_slc_alosにすごく時間がかかるケースが
複数確認された。
• Ubuntu 10 & 32 bitのレンタルPCでは問題な
かった。また、Ubuntu 11 & 64 bitのマシンで
も問題なかった。
• 原因の特定のためには、問題のあったPCに
Ubuntu 10を入れて試してみる等が必要。
13
SLCの強度画
像
（ルック数
5 x 15）
レーダ座標系なので、
左右逆であることに
注意。（北行軌道の
データだと、上下逆
になります。）
14
強度画像image.slc_5l.rasを開く
• GIMP (Linux/Mac) 、ToyViewer (Mac) など
• ここでやること：解析範囲（最終的な干渉画像
の範囲）を決定する
• そのために、まずras画像中で切り出したい領
域の画素番号をしらべる（今回は、X:230-700,
Y:1270-1600あたりを切り出すことにする）
• シングル・ルックの場合にどうなるかを計算。
X, Yそれぞれ5と15をかけた値をメモる。
（x=1150-3500, y=18000-24000)
15
../roi.proc の編集
切り出し領域を設定
• roi.procの中身
use1dopp=1
mean_pixel_rng = 2325 ←レンジ方向の解析領域中心
ymin = 17000 ←アジマス方向の解析領域下限
ymax = 25000 ←アジマス方向の解析領域上限
• アジマス方向のライン数は、1000以上余裕を与えて
設定する。つまり、最終的な解析領域を基準に、ymin
は1000以上小さい値、ymaxは1000以上大きい値に
設定する。
• ファイルの編集が終わったら、再度step_slc_alosを実
行。切り出された領域でSLCが作成される。
16
image.slc_5l.ras
アジマス方向に切り出された画像ができる。
（合成開口長より短くはならなさそう）
17
master_crop.inの作成・編集
最終的な解析領域を設定
• cp $MY_SCR/master_crop.in .
• vi master_crop.in
（もちろんviでなくgedit等でも可）
first_l 2001
last_l 6500
first_p 1151
last_p 3500
• first_lとlast_lは、アジマス方向のライン数。
• 切り出された後のライン数を基準に設定。
18
step_master_setup
マスターの準備
cd ~/kuju/2950/SLC/20090530
step_master_setup
• パラメタファイルの作成などとともに、20090530_crop.ras
ができる。この画像が、最終的な解析領域。解析したい領
域が適切に切り出されているかチェックし、必要に応じ
master_crop.inを編集してやり直す。
• もし、20090530_crop.rasが真黒な場合、切り出しに失敗し
ている。master_crop.inで設定した領域がはじめに切り出
した領域外になっているとこうなるので、master_crop.inを
修正してやり直すこと。
19
make_slc_alos
すべてのスレーブ画像のSLCを作成
cd ~/kuju/2950/SLC
make_slcs_alos
• しばらく時間がかかる。（今回のデータセット（スレーブ
５枚）で10〜30分）
• 各データのディレクトリ内にras形式の強度画像が作
成されているので、 問題なくできているかを確認する。
• Toyviewerを使うなら、rasファイルはToyViewerで開くこ
とを定義しておき、open */*.rasとやると簡単。
20
Tips
• StaMPSのシェルスクリプトには、step… という
名前と make… という名前のものがある。
• step… のスクリプトは、あるひとつの画像に対
する処理をおこなう。
• make… は、すべての画像に対する処理をお
こなう。step… を順番に実行している。
21
DEM（DEHM）の準備
• 今回は、国土地理院基盤地図情報10m（0.4秒）メッ
シュ標高とEGM96から作成したDEHM(Digital
Ellipsoildal Height Model)を解析用に整形したものを
用いる。
• DEHMについて：pixelメンバー（共有データアクセス可
能な方）は、pixelサーバーからダウンロード可能。小
澤さんが数年前に準備（[pixel: 986] のメールでアナウ
ンス）。http://pixel.eri.u-tokyo.ac.jp/data/DEHM-10m/
• 今回使うDEHM（DEM/4931_ds.dehm）は、データを半
分に間引き（ピクセル間隔を倍にし）、範囲を狭く切り
取ったもの。
22
用いるDEHM
23
timing.dorisin の編集
DEMパラメタの設定
• ~/kuju/2950/INSAR_20090530内にtiming.dorisinと
いうファイルがある。これはDORIS（M_SIMAMP,
M_TIMING）用のパラメタファイル。この中のDEMに
関するパラメタを書き換える。
SAM_IN_FORMAT real4
SAM_IN_DEM /Users/yo/kuju/DEM/4931_ds.dehm
SAM_IN_SIZE 2000 2500
// rows cols
SAM_IN_DELTA 0.0002222222 // in degrees
SAM_IN_UL
33.3332778 131.0000556 // lat and lon of upper left
SAM_IN_NODATA 0
SAM_OUT_FILE master_sam.raw // synthetic amplitude
SAM_OUT_DEM dem_sam.raw // cropped DEM for
24
step_master_timing
Timing（アジマス方向）Errorの推定
• DEMと位置合わせをすることにより、軌道
データから推定されるアジマス方向の位置の
ずれを推定。
cd ~/kuju/2950/INSAR_20090530
step_master_timing
25
step_master_timing中、下のようなWarningが出るとうまくいかな
い場合がある。もし、次のページのような結果が出なければ、
失敗していることになるので、修正してやり直す。
こうなる原因は二つ（かそれ以上）ある。１：「時間」がおかしい
（ロケールの変更をしていないとこうなる）、２：DEMの範囲が狭
い。これは、DEMの範囲を拡げれば解決する。
•
•
•
•
•
•
WARNING : DEM does not cover entire interferogram.
WARNING : input DEM should be extended to the North.
WARNING : indexphiNDEM: 2085
WARNING : DEM does not cover entire interferogram.
WARNING : input DEM should be extended to the South.
WARNING : indexlambda0DEM: -126
DEMを変えてやり直す場合、step_master_timingからで多分大丈夫。ロケール
の問題があるかどうかは、master.res（dorisでの処理に使うパラメタファイル）
の中の「First_pixel_azimuth_time (UTC):」をチェックすればよい。正常であれ
ば、07-Jan-2007 01:45:44.76455881 というようなフォーマットになっている。もし
このようになっていなければ、make_slcs_alosからやり直すことをお勧めする。26
step_master_timing最後に出るメッセージ
Number
0 3044
1 3044
2 3044
posL posP offsetL offsetP correlation
1682 24 0 0.427386
2037 22 0 0.321141
2391 20 0 0.294289
.
.
.
27 5455 2258 22 0 0.401071
28 5455 2612 22 0 0.370528
29 5455 2967 20 0 0.379082
Estimated total offset (l,p): 20, 0
Coherence NaN values are disregarded in the analysis.
*********************************************************
**********
correlationが0.1以下だったり、Estimated total offset (l,p)が数100以上？
だったりした場合、多分うまくいっていない。前頁に従い問題を解決する。
27
make_orbits
INSAR_20090530のディレクトリ内に各ス
レーブ用のディレクトリ＆軌道データを含
むパラメタファイルが作成される
make_coarse
大雑把な位置合わせ
cd ~/kuju/2950/INSAR_20090530
make_orbits
make_coarse
28
make_coreg
精密位置合わせ
• 通常、このステップは時間がかかる。
• 時間節約のためにパッチ数を調整してから実行してもよい。
• パッチ数の変更：$STAMPS/matlab/coreg_pos.mのファイルを編
集。nWinxとnWinyをそれぞれ20と50等に変更する。
• coregというディレクトリの下に位置合わせ結果ファイルができ
る。やり直す場合、CPM_Data.x.x を消すかcoreg自体を消す。
cd ~/kuju/2950/INSAR_20090530
make_coreg
MacでR2011?以降のmatlabを使用している人は、以下の様なメッセージが出るかもしれ
ないが、無視してよい。
2012-08-27 00:03:55.953 MATLAB[16905:8103] This process is attempting to exclude an
item from Time Machine by path without administrator privileges. This is not supported.
29
Tips
• make_coregでエラーになるときは、matlabの
パスが通っていない疑いがある。
• コマンドラインからmatlab と打ち、matlabが起
動される状態になっているか確認し、そうなっ
ていなければパスを通してください。
30
make_dems
地形縞のシミュレート画像を作成
• refdem_11.raw などのファイルが各スレーブ
に対して作られる。
• この処理も結構時間がかかる。
• マニュアルによると、make_coregと並行して
走らせて問題ない。（が、人為的ミスを防ぐた
めにも順番に走らせることをお勧めする）
cd ~/kuju/2950/INSAR_20090530
make_dems
31
make_resample
スレーブSLCのリサンプリング
• マスター画像と同じグリッドのスレーブ画像が
作成される。所要時間５分程度。
• このあと、リサンプルされた画像のファイルサ
イズがすべて同じになっているかを確認す
る。
cd ~/kuju/2950/INSAR_20090530
make_resample
ls -l */slave_res.slc
32
make_ifgs
干渉画像の作成
• 干渉画像 cint.minrefdem_5l.ras が各スレー
ブディレクトリ内にできる。すぐ終わる。
• できた干渉画像は、Matlabでplot_all_ifgsを使
うか、gimpやToyViewerなどで確認する。
• 全く縞が確認できないようなものは、失敗して
いる可能性が高い。本演習の例題の場合、
処理がうまくいけば干渉性は高くなるはず。
33
plot_all_ifgs(0,3)
07 Jan 2007
25 Nov 2007
27 May 2008
30 May 2009
17 Apr 2010
20 Apr 2011
MASTER
34
step_geo
ジオコーディング
• これは、後のPS処理のために必要。 スレーブ
画像のディレクトリのどれか一つだけに移動
し、コマンドを実行する。
• 実行後、InSAR_20090530にlat.raw、lon.raw
が作られる。
cd ~/kuju/2950/INSAR_20090530/20070107
step_geo
35
ここからPS処理
mt_prep
PS候補点の選定
cd ~/kuju/2950/INSAR_20090530/
mt_prep 0.4
36
（mt_prepのオプション）
mt_prep 0.4 2 2 50 200
（意味）
0.4: Amplitude dispersionのしきい値（これ以下を用いる）
2: レンジ方向のパッチ数（デフォルト1）
2: アジマス方向のパッチ数（デフォルト1）
50: レンジ方向のパッチ間のオーバーラップピクセル数（デフォルト50）
200: アジマス方向のパッチ間のオーバーラップピクセル数（デフォルト
200）
StaMPSの処理はパッチごとに行われ、あとでマージされる。解析エリアに
対してメモリに余裕があれば、パッチに分ける必要はない。
37
ディスクスペース
• 今回の講習会の解析で、最終的に12.5GB程
度のディスク容量を占有します。より多くの画
像を用い、フルシーンの解析をやると、サイズ
がずっと大きくなります。
• stampsのステップ１が終わったら、以下のスク
リプトを使って不要となったファイルを消去で
きます。これらのファイルは簡単に作り直せま
す。
38
ディスクスペース節約のためのスクリプト
（すべて INSAR_xxxxx で実行）
• make_clean_ifgs
– 干渉画像の消去。復活させたければ、make_ifgs
を再実行。
• make_clean_resample
– リサンプルされたスレーブSLCの消去。復活させ
たければ、make_resampleを再実行。
• make_clean_raw
– オリジナルのSLC（ROI_PACにより作成）の消去。
復活させるには、make_slcs_xxxxを再実行。
これらは、最終結果が得られ、それ以上立ち戻る必要がなくなってから実行するこ
とをお勧めします。ディスクに余裕があれば、もちろん実行の必要はありません。 39
以下の処理は、matlab上で実行す
る。matlabを起動して、処理ディレク
トリに移動する。その後、処理を実
行する。
40
(getparm)
処理パラメタの確認
• ~/kuju/2950/INSAR_20090530/ で、
>> getparm
と打つと、処理パラメタが出力される。
• コマンド setparm でパラメタを変更することが
できる。通常、変更する必要はない。パラメタ
の説明はマニュアルにあるが、すべてではな
いかもしれない。
41
R2012aだとエラーになる部分あり。
その回避方法：
• matlabのコマンドラインで、”edit uw_interp”と
打ち、uw_interp.mをエディターで開く
• ５５行目を以下のように変更。
（変更前）
Z=dsearch(x,y,ele(:,2:4),X,Y); %index from grid to pixel node
（変更後）
Z=dsearchn([x,y],ele(:,2:4),[X(:),Y(:)]); %index from grid to pixel node
Z = reshape(Z,nrow,ncol);
42
新しいMATLABのバージョンだと、以下のエラーが出る（少し古いバージョンだと
Warning）。
Error using griddata (line 62)
GRIDDATA no longer requires Qhull-specific options.
Please remove these options when calling GRIDDATA.
Error in ps_load_initial (line 258)
la0_ps=griddata(gridX,gridY,la0,ij(:,3),ij(:,2),'linear',{'QJ'});
Error in stamps (line 116)
ps_load_initial
このエラーを回避するために、ps_load_initial.m でgriddataを呼んでいる行で、オプ
ションの ,{‘QJ’} を消す。
（修正前）griddata(gridX,gridY,la0,ij(:,3),ij(:,2),'linear',{'QJ'});
（修正後）griddata(gridX,gridY,la0,ij(:,3),ij(:,2),'linear');
43
stamps
PS処理の実行
• 何も指定しなければ、全ステップの処理を自
動的に実行してくれる。
• stampsの処理は８つのステップに分かれてお
り、処理ステップを指定可能。例：stamps(1,3)
だと、ステップ１〜３を実行。stamps(4,8)だと、
ステップ４〜８を実行。
>> cd ~/kuju/2950/INSAR_20090530
>> stamps
44
Tips
• ファイル修正や設定のし忘れ等により、途中でエ
ラーが発生しやり直す場合、pwdでディレクトリを
確認すること。
• INSAR_20090530の下のPATCH_xにいる可能性
があります。そこのディレクトリでstampsを実行す
ると、そのパッチでしか処理がされません。
• INSAR_20090530にいることを確認してstampsを
走らせてください。
• stampsは、一度失敗しても、やり直しがききます。
45
ps_plot
結果のプロット・エクスポート
• このmatlab関数には、非常の多くの機能があ
る。詳しくは、 help ps_plot や、StaMPSマニュ
アルの Chapter 9 参照。
46
ps_plot(‘u-dm’,0,[-4 4],1) の結果
このあたりアンラップエラー？
（時間的にも空間的にも不連続の場合疑わしい）
stamps ステップ6以降をもう一度やると改善する
可能性がある（詳しくはマニュアル参照）
47
結果を少し細かく見てみましょう。
• ps_plot(‘w’)
– もしノイズ大きすぎなら、weed_standard_devを小さくしてstep 4以降を再実
行。
– もしピクセル少なすぎなら、weed_standard_devを大きくしてstep 4以降を再実
行。
• ps_plot(‘u’)
– もし時間的かつ空間的に不連続な場合、アンラップエラーが考えられる。
steps 6から8を再実行すれば改善するかもしれない。[-2*pi 2*pi]で表示させ
るとわかりやすい。
• ps_plot(‘d’)
– Spatially Correlated Look Angle Error (SCLA)
• ps_plot(‘m’)
– Master Atmosphere and Orbital (AOE) Error
• ps_plot_(‘u-dms’)
– DEM（d）、マスター画像の誤差（m）とともに、各スレーブ画像の誤差（s）も除
去。時間軸上のローパスフィルタを使うので、非定常な変動があれば前後に
染み出すことになるので注意。
• Step 7の前に setparm(‘scla_deramp’,’y’) とすると、主に軌道残存縞である
位相トレンドを除去して解析する。推定された位相トレンドは、ps_plot(‘o’)
48
で図示化される。
結果のrefinement
• ps_plot(‘u’) とps_plot(‘u-dm’) を比べてみる。もし、後者が前者
よりなめらかな時系列になっていなければ、どこかにアンラップ
エラーがある可能性が高い（通常、Bperpが大きいデータ）。
• このような場合、とりあえず、エラーがありそうな干渉画像を除
いて（scla_drop_indexに設定）step 7を再実行。
• ps_plot(‘u-dm’)がなめらかになったら、step 6を再実行。
• 結果がよくなっていたら、先ほど除いた干渉画像を元に戻して
step 7を再実行。
• アンラップが正常にできるか、改善がもう見られないところまで
上の手順を繰り返す。
• さらにアンラップの精度を上げるための方策は、マニュアルp20
を参照。
49
Small Baseline Processing
短基線長（SB）解析
• PS解析とは別の時系列解析を行う手法。
• 「StaMPS」は、正確にはPS解析の部分を指す。
• StaMPS/MTI（Multi-temporal InSAR）は、StaMPSの拡張版で、SB
解析や、PS解析とSB解析を合わせた解析を行うことができる。
• StaMPS/MTIのSB解析は、single-look干渉画像で処理を行うの
で、基本的にPS解析と類似した結果になる。
• SB解析のメリット：アンラッピングの信頼度評価ができる。PS解
析と違うアルゴリズムで解くので、結果の比較が可能。PS+SB解
析することにより、信頼性が向上、ピクセル数増加。
50
Tips
• SB解析は、INSAR_20090530のディレクトリの
下で実行しますが、いろいろ試行錯誤をして
いるうちに、元のPS解析結果の重要なファイ
ルが消されたりすることがあります。
• PS解析まで終わったら、結果のディレクトリを
保存しておくことをオススメします。
51
（PS解析をせずSB解析をする場合）
• MATLABで基線長に関する情報をロードする
ために以下を実行。
mt_extract_info
matlab
>>ps_load_info
52
sb_find
短基線長ペアの選定
• コヒーレンス、時間基線長、垂直基線長を基
準に短基線長ペアを選定する。
(in matlab)
>>cd ~/kuju/2950/INSAR_20090530
>> sb_find(0.5, 5000, 4000)
（0.5: コヒーレンスのしきい値、5000: 臨界時間基線長、
4000: 臨界垂直基線長）
コヒーレンス ＝ (1 – Bt/Btc)*(1-Bp/Bpc) というモデルから計算
したコヒーレンスがしきい値を超えるものを選出する。
表示される図を見てみて、線でつながっていない点があった
り複数のクラスタに分かれたりする場合は、パラメタを調整し
てつながるようにする。
53
sb_find（続き）
• sb_findを実行すると、small_baselines.listとい
うファイルができる。その中にペアが書かれ
ている。
• 自分でこのファイルにペアを書き込んでもよ
い。
• plot_sb_baselinesを実行すると、
small_baselines.listに含まれるペアが図示化
される。
54
55
make_small_baselines
干渉処理
• PS干渉画像との違いは、common bandpath
filterをかけるというところ。range and/or azimuth
フィルタをかけないようにするというオプションも
ある（マニュアル参照）。
• この処理により、INSAR_20090530の下に
SMALL_BASELINESというディレクトリができ、そこ
に干渉画像が格納される。
cd ~/kuju/2950/INSAR_20090530
make_small_baselines
56
in SMALL_BASELINES
plot_all_ifgs(0,4,10)
57
mt_prep
候補点の選定
cd ~/kuju/2950/INSAR_20090530/SMALL_BASELINES
mt_prep 0.6
58
mt_prepのオプション:SB解析編
(PSの場合とほとんど意味は同じ）
mt_prep 0.6 2 2 50 200
（意味）
0.6: Amplitude difference dispersionのしきい値（これ以下を用いる）
2: レンジ方向のパッチ数（デフォルト1）
2: アジマス方向のパッチ数（デフォルト1）
50: レンジ方向のパッチ間のオーバーラップピクセル数（デフォルト50）
200: アジマス方向のパッチ間のオーバーラップピクセル数（デフォルト200）
StaMPSの処理はパッチごとに行われ、あとでマージされる。解析エリアに
対してメモリに余裕があれば、パッチに分ける必要はない。
59
stamps
SB処理の実行
• PS処理と同様、何も指定しなければ、全ステップの
処理を自動的にしてくれるし、処理ステップの範囲を
指定可能。
• マニュアルによると、PS処理よりメモリをだいぶ食う
とのこと。パッチに分けて解析する必要性が生じる
かもしれない。
>> cd ~/kuju/2950/INSAR_20090530/SMALL_BASELINES
>> stamps
60
新しいMATLABのバージョンだと、以下のエラーが出る（少し古いバージョンだと
Warning）。
Error using griddata (line 62)
GRIDDATA no longer requires Qhull-specific options.
Please remove these options when calling GRIDDATA.
Error in sb_load_initial (line 227)
la0_ps=griddata(gridX,gridY,la0,ij(:,3),ij(:,2),'linear',{'QJ'});
Error in stamps (line 116)
sb_load_initial
このエラーを回避するために、sb_load_initial 116行目、117行目、254行目
の ,{‘QJ’} を消す。
la0_ps=griddata(gridX,gridY,la0,ij(:,3),ij(:,2),'linear');
la1000_ps=griddata(gridX,gridY,la1000,ij(:,3),ij(:,2),'linear');
bp0_ps=griddata(gridX,gridY,bp0,ij(:,3),ij(:,2),'linear');
61
ps_plot(‘u-dm’,0,[-4 4],1) の結果
(default parameters)
62
ps_plot(‘usb’,0,[-4 4],1)
（干渉画像）
63
ps_plot(‘rsb’,0,0,1)
（アンラップ誤差）
+-p以内の周囲と無相関な誤差は気にしな
くてよいとのこと。このrsbプロットに空間的
相関がある場合、悪さをしていると思われ
る干渉画像を除いて再解析してみる
（Manual p22）。
64
ピクセル数を増やすため、やり直し。
• SBの場合、’density_rand’がデフォルトで2に設定され
ている（PSは20）。
• setparm(‘density_rand’,20) とすると、PSと同じパラメタ
値でやることになる。
• このパラメタは、km^2あたりに許容されるノイズ大の
ピクセル数。
• このパラメタが関係するのは、Step 3なので、その先
からやり直せばよい（stamps(3,8))。
• なお、先に”setparm”をしてから”stamps”をすると、自
動的にdensity_randが2に書き換えられてしまった。
stamps(1,2)→setparm(‘density_rand’,20)→stamps(3,8)
とやれば多分OK。
65
ps_plot(‘u-dm’,0,[-4 4],1) の結果
(density_rand = 20)
66
参考 ps_plot(‘u-dm’,0,[-4 4],1) の結果
(default parameters)
67
PSとSB解析結果の融合
>> cd ~/kuju/2950/INSAR_20090530
>> ps_sb_merge
これにより、INSAR_20090530の下に MERGED というディレクトリができ、
そこにマージしたデータが置かれる。これは、stamps処理のステップ５ま
で（アンラッピングの手前まで）の処理結果を融合したものである。よっ
て、最終結果を得るためには残りのステップ６〜８を実行する。
>> cd MERGED
>> stamps(6,8)
これらの処理は、すぐ終わります。
68
ps_plot(‘u-dm’,0,[-4 4],1) in MERGED
69
ps_plot(‘u-dms’,0,[-4 4],1) in MERGED
（スレーブ画像のノイズも軽減）
70
マルチルック処理をしたSB解析
• この機能は、公開パッケージには含まれてい
ませんが、開発者（Andy Hooper）に配布の許
可をもらいました。
• 通常のSB解析と同様に、マルチルック処理を
した干渉画像を作成し、それを用いて最小二
乗法で解きます。
71
ルック処理版SB解析
•
•
•
•
•
cd ~/kuju/2950/INSAR_20090530
mv SMALL_BASELINES SMALL_BASELINES_slk （Single lookのSB解析の結果を上書きされな
いよう別名にしておくだけ）
make_small_baselines
cd ~/kuju/2950/INSAR_20090530/SMALL_BASELINES
ここで make_ifgs.listを作成
/Users/yo/kuju/2950/INSAR_20090530/SMALL_BASELINES/20070107_20071125
/Users/yo/kuju/2950/INSAR_20090530/SMALL_BASELINES/20070107_20080527
…（省略）…
/Users/yo/kuju/2950/INSAR_20090530/SMALL_BASELINES/20080527_20110420
/Users/yo/kuju/2950/INSAR_20090530/SMALL_BASELINES/20090530_20100417
• make_cohs （！今講習会で配布するバージョン）
• mt_extract_info
(in matlab)
• >> cd ~/kuju/2950/INSAR_20090530
• >> repair_lonlat_5lraw （今講習会で配布。バグに対応するため福島作）
• >> cd SMALL_BASELINES
• >> mt_ml_select （今講習会で配布）
• >> stamps(6,8)
72
ルック数の指定
• ルック数は、デフォルトで5 x 15。
• これを変えたい場合、mt_ml_selectを実行す
る前に looks.txt と ar.txt を編集する。
• looks.txtには、レンジ方向のルック数を書く。
ar.txtは、アスペクト比が書いてあるが、デフォ
ルトでは3（つまり、アジマス方向のルック数＝
レンジ方向のルック数（5）×アスペクト比（3）
で15）。
73
ps_plot(‘u-dm’,0,[-4 4],1)
SMALL_BASELINES (マルチルック版）
74
ps_plot(‘u-dm’,0,[-4 4],1) in MERGED
75
ps_plot(‘u-dm’,0,[-4 4],1) の結果
(SB only, density_rand = 20)
76
ps_plot(‘u-dm’,0,[-4 4],1) の結果
(PS only, density_rand = 20)
77
ps_plot 使用例：kmlファイル出力
>> ps_plot('v-d',-1,[-20 20])
>> load ps_plot_v-d
>> ps_gescatter('kuju-psinsar.kml',ph_disp,1,1)
Tips:
• ps_plotの２つ目の変数に-1を入れると、プロットされ
るはずのデータがmatファイルにエクスポートされ
る。matファイルの名前は、ps_plot_処理方法指
定.mat となる。
• この例では、一度エクスポートしたものをロードし、
ロードされたph_dispという変数をさらにkmlファイル
にエクスポートしている。
78
Google Earthで見た結果（MERGEDの
平均速度） [-20 20] mm/year
79
80
結果をファイルに落とす
>> save filename.dat ph_disp –ascii
などとすれば、ファイルに結果を落とすことがで
きる。→GMT等でプロット可能。
ピクセル位置はps2.mat内に書かれているの
で、>> load ps2
>> save lonlat.dat lonlat –ascii
などとやる。
81
時系列のプロット
（PS解析の結果以外（SBやMERGED）
ではエラーが出てできなかった）
>> ps_plot('v-d', 'ts')
82
時系列のプロット （PS解析の結果以外（SBや
MERGED）ではエラーが出てできなかった）
>> ps_plot('v-d', 'ts’)
これが時系列の図ですが、図の意
味がわからない部分あり。わかった
ら教えてください。
83
LOSベクトルの計算
• heading角と入射角が必要。
• heading角は、stamps処理の最中に、
heading.1.in という名前のファイルに書き出され
る。ROI_PACのパラメタファイル（*.rsc）にも書いて
ある。
• 入射角は、干渉画像が置かれているスレーブの
ディレクトリ（INSAR_20090530/yyyymmdd、どれ
でもよい）にあるdorisのファイルcoreg.out,
dem.out, geocode.outに書かれている。
（grep incidence coreg.out などとやると簡単）
84
その他諸々
• データのエクスポートは、ps_plotだけでなく、
ps_outputでもできる。
• ref_lon, ref_lat（またはref_centre_lonなど）で
不動域を指定することができます（ちゃんとで
きるかは未確認です）
85
参考文献など
•
•
•
•
•
•
•
•
StaMPSマニュアル（StaMPSパッケージに含まれる。下のウェブサイトからも入手可能）
StaMPSウェブサイト
http://radar.tudelft.nl/~ahooper/stamps/index.html
ユーザーグループ（過去ログの検索可）
http://groups.google.com/group/mainsar
時系列解析方法レビュー
Hooper, A., D.Bekaert, K. Spaans, and M. Arıkan (2012): Recent advances in SAR interferometry time
series analysis for measuring crustal deformation, Tectonophysics, Volumes 514–517, 1-13, doi:
10.1016/j.tecto.2011.10.013.
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.
http://earth.eo.esa.int/workshops/fringe09/proceedings/papers/p1_26hoop.pdf
86