480 likes | 1.06k Vues
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:
E N D
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 xzvfsuppl.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データの解凍とディレクトリの作成link_raw_pixelデータの解凍とディレクトリの作成 • (PIXELアーカイブのファイル名フォーマットを読み込むように、オリジナルのlink_rawというシェルスクリプトを改変してある。) • このコマンドにより、解析ディレクトリの下にフレーム番号(7020)が作成され、 その下にyymmddという名前のディレクトリが作成される。その中に、RAWデータ のリンクが置かれる。 使用法: link_raw_pixeldata_diranalysis_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の作成・編集最終的な解析領域を設定master_crop.inの作成・編集最終的な解析領域を設定 • cp $MY_SCR/master_crop.in . • vi master_crop.in (もちろんviでなくgedit等でも可) first_l1501 last_l4000 first_p1501 last_p2750 • first_lとlast_lは、アジマス方向のライン数。切り出された後のライン数を基準に設定。 • (元々切り出すと決めた領域は1000,4000,1401,2750に相当するが、解析領域を絞るため若干変えた)
step_master_setupマスターの準備 • パラメタファイルの作成などとともに、20080718_crop.ras ができる。この画像が、最終的な解析領域。解析したい領域が適切に切り出されているかチェックし、必要に応じmaster_crop.inを編集してやり直す。 • もし、20080718_crop.rasが真黒な場合、切り出しに失敗している。master_crop.inで設定した領域がはじめに切り出した領域外になっているとこうなる。 cd ~/tambora/7020/SLC/20080718 step_master_setup
make_slc_alosすべてのスレーブ画像のSLCを作成 • しばらく時間がかかる。 • 各データのディレクトリ内にras形式の強度画像が作成されているので、 問題なくできているかを確認する。 • Toyviewerを使うなら、rasファイルはToyViewerで開くことを定義しておき、open */*.rasとやると簡単。 cd ~/tambora/7020/SLC make_slcs_alos
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を使って、以下のようにデータを切り出しする(2バイトベタ書きフォーマットのものを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'); whosb 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 //rowscols 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_timingTiming(アジマス方向)Errorの推定step_master_timingTiming(アジマス方向)Errorの推定 • DEMと位置合わせをすることにより、軌道データから推定されるアジマス方向の位置のずれを推定。 cd ~/tambora/7020/INSAR_20080718 step_master_timing
(続き)以下のような出力が得られれば正常 Number posLposPoffsetLoffsetP 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_orbitsINSAR_20080718のディレクトリ内に各スレーブ用のディレクトリが作成されるmake_coarse大雑把な位置合わせmake_orbitsINSAR_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_prepPS候補点の選定 cd ~/tambora/7020/INSAR_20080718/ mt_prep 0.4
(mt_prepのオプション) (意味) 0.4: Amplitude dispersionのしきい値(これ以下を用いる) 2: レンジ方向のパッチ数(デフォルト1) 2: アジマス方向のパッチ数(デフォルト1) 50: レンジ方向のパッチ間のオーバーラップピクセル数(デフォルト50) 200: アジマス方向のパッチ間のオーバーラップピクセル数(デフォルト200) StaMPSの処理はパッチごとに行われ、あとでマージされる。解析エリアに 対してメモリに余裕があれば、パッチに分ける必要はない。 mt_prep 0.4 2 2 50 200
以下の処理は、matlab上で実行する。matlabを起動して、処理ディレクトリに移動する。その後、処理を実行するという流れが基本。以下の処理は、matlab上で実行する。matlabを起動して、処理ディレクトリに移動する。その後、処理を実行するという流れが基本。
(getparm)処理パラメタの確認 • ~/tambora/7020/INSAR_20080718/ で、 >> getparm と打つと、処理パラメタが出力される。 • コマンドsetparmでパラメタを変更することができる。通常、変更する必要はない。パラメタの説明は、マニュアルには書いてない。Hooper et al. (2007)のアルゴリズムの説明を読めば、だいたいは直感的に理解できる。
stampsPS処理の実行 • 何も指定しなければ、全ステップの処理を自動的にしてくれる。 • stampsの処理は8つのステップに分かれており、オプションを指定することにより各々のステップのみを実行することも可能。詳しくは、help stamps 参照。 >> cd ~/tambora/7020/INSAR_20080718 >> stamps
ps_plot結果のプロット・エクスポート • このmatlab関数には、非常の多くの機能がある。詳しくは、 help ps_plotや、StaMPSマニュアルの Chapter 9 参照。
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の2つ目の変数に-1を入れると、プロットされるはずのデータがmatファイルにエクスポートされる。matファイルの名前は、ps_plot_処理方法指定.mat となる。 • この例では、一度エクスポートしたものをロードし、ロードされたph_dispという変数をさらにkmlファイルにエクスポートしている。
LOSベクトルの計算 • heading角と入射角が必要。 • heading角は、stamps処理の最中に、heading.1.in という名前のファイルに書き出される。ROI_PACのパラメタファイル(*.rsc)にも書いてある。 • 入射角は、干渉画像が置かれているスレーブのディレクトリ(INSAR_20080718/yyyymmdd、どれでもよい)にあるdorisのファイルcoreg.out, dem.out, geocode.outに書かれている。(grepincidence 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 scattererInSAR for crustal deformation analysis, with application to VolcanAlcedo, 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.