為解決 SAR 影像應用於淹水實際操作之複雜性，目前建構一套由以腳本操 作且相對自動化的緊急災害偵測標準流程，需要輸入事件前與事件後影像的檔 名，處理的結果可產生一張由藍色系突顯可能淹水變異區域的變遷偵測影像，做 為緊急應變的依據。
整個處理程式包括以下檔案：
⚫
runGPT_General_1.BAT - 使用 Sentinel-1 以外 SAR 影像的串聯程序之主程 式
⚫
runGPT_S1_1.BAT -使用 Sentinel-1 單張 SAR 影像的串聯程序的主程式
⚫
runGPT_S1_2.BAT -使用 Sentinel-1 兩張 SAR 影象的串聯程序的主程式
⚫
orb.cal.xml - Sentinel-1 專用，給予精密軌道與輻射校正的 xml 參數檔
⚫
Coregistration.xml -事件前後影像套疊 xml 參數檔
⚫
runGRASSGIS.BAT - GRASS-GIS 工作的啟動程式
⚫
GRASS_JOB.BAT - GRASS-GIS 計 算 height above the nearest drainage
（HAND）的內容
runGPT*.BAT 內需要設定一些參數，包括影像檔案位置（master* and slave*） 、
裁切範圍（subregion） 、前處理 xml 參數檔與 GRASS-GIS 程式所在目錄（graph） 、
輸出資料夾（project）與影像種類（sar_type）的選擇。以下是每個時間段兩張
Sentinel-1 影像為輸入、興趣範圍（W/E/S/N）為 138.576/140.994/35.62/37.368 與
輸出產品名稱為 output 的範例：
設定好之後，打開 cmd，執行
>runGPT_S1_2.BAT
runGPT*.BAT 將執行以下步驟來完成淹水變異區域的標示
1.
當興趣區間跨兩張共軌道影像的區域，需先將共軌道影像合併。
2.
裁切至興趣區間來降低所需處理的資料量，並減少處理時間。
3.
給予精密軌道資訊與輻射校正
4.
斑駁濾波降低雜訊
set master1='..\IMAGE\S1A_IW_GRDH_1SDV_20191006T204300_20191006T204325_029343_035 5EA_E1F1.zip'
set master2='..\IMAGE\S1A_IW_GRDH_1SDV_20191006T204325_20191006T204350_029343_035 5EA_1B68.zip'
set slave1='..\IMAGE\S1B_IW_GRDH_1SDV_20191012T204223_20191012T204248_018447_022C 0C_CBFE.zip'
set slave2='..\IMAGE\S1B_IW_GRDH_1SDV_20191012T204248_20191012T204313_018447_022C 0C_53C5.zip'
set subregion='POLYGON ((138.576 35.62,140.994 35.62,140.994 37.368,138.576 37.368,138.576 35.62))'
set graph=graph set project=output set sar_type=1 if %sar_type%==1 ( set sar_type=sentinel-1 set pol=VV
)
if %sar_type%==2 (
set sar_type=general Strip map SAR set pol=HH
)
gpt SliceAssembly -Ssource=%master1% %master2% -t
%project%\master.dim
gpt SliceAssembly -Ssource=%slave1% %slave2% -t
%project%\slave.dim
gpt subset Ssource=%project%\master.dim PgeoRegion=%subregion% -PcopyMetadata=true -t %project%\master.sub.dim
gpt subset Ssource=%project%\slave.dim PgeoRegion=%subregion% -PcopyMetadata=true -t %project%\slave.sub.dim
gpt %graph%\%orb_cal% -Ssource=%project%\master.sub.dim -t
%project%\master.sub.orb.cal.dim
gpt %graph%\%orb_cal% -Ssource=%project%\slave.sub.dim -t
%project%\slave.sub.orb.cal.dim
gpt Speckle-Filter -Ssource=%project%\master.sub.orb.cal.dim -t
%project%\master.sub.orb.cal.spk.dim -Pfilter="Refined Lee"
gpt Speckle-Filter -Ssource=%project%\slave.sub.orb.cal.dim -t
%project%\slave.sub.orb.cal.spk.dim -Pfilter="Refined Lee"
5.
計算雷達陰影區域 務。任完成後，請輸入 exit，並按下 Enter，回到主程式 runGPT*.bat。
10.
加入 HAND 作為判斷遮罩
gpt SAR-Simulation -Ssource=%project%\master.sub.orb.cal.spk.dim t %project%\master.sub.orb.cal.spk.OSM.dim
-PsaveLayoverShadowMask=true
gpt %graph%\%Coreg% e
PinFile1=%project%\master.sub.orb.cal.spk.OSM.dim -PinFile2=%project%\slave.sub.orb.cal.spk.dim -t
%project%\stack.dim
gpt %graph%\%dB% -Pfile=%project%\stack.dim -t
%project%\stack.db.dim
REM Using SRTM 3 sec
gpt Terrain-Correction -Ssource=%project%\stack.db.dim -t
%project%\stack.db.TC.dim -PnodataValueAtSea=true -PsaveDEM=true REM Using externalDEM
REM gpt Terrain-Correction -Ssource=%project%\stack.db.dim -t
%project%\stack.db.TC.dim PexternalDEMFile=5Mdem_wgs84.tif PexternalDEMNoDataValue=3.40282e+38 PnodataValueAtSea=true -PsaveDEM=true
REM gpt Terrain-Correction -Ssource=%project%\stack.db.dim -t
%project%\stack.db.TC.dim PexternalDEMFile=TW20mDEM_wgs84.tif PexternalDEMNoDataValue=32767 PnodataValueAtSea=true
-PsaveDEM=true
REM gpt Terrain-Correction -Ssource=%project%\stack.db.dim -t
%project%\stack.db.TC.dim PexternalDEMFile=ALOS30mDSM.tif PexternalDEMNoDataValue=32767 PnodataValueAtSea=true -PsaveDEM=true
echo ******************************************************
echo * Type "exit" after another cmd finish GRASS-GIS job * echo ******************************************************
copy %graph%\GRASS_JOB.bat GRASS_JOB.bat
copy %project%\stack.db.TC.data\elevation.* %tmp%
start /wait %graph%\runGRASSGIS.bat %tmp%\elevation.img copy HAND.* %project%\ /Y
del GRASS_JOB.bat HAND.* %tmp%\elevation.*
gpt Merge SmasterProduct=%project%\stack.db.TC.dim
-SsourceProducts=%project%\HAND.hdr -PgeographicError=NaN -t
%project%\stack.db.TC.HAND.dim
11.
以事件前後之雷達差異，將變異之淹水區以藍色標示。
以 2019/10/12 侵襲日本的哈吉貝颱風造成的淹水作為範例，將事件前設定 為藍色，事件後設定為綠色與紅色，由於事件後淹水區域沒有水體，雷達反射數 值較高，事件後淹水區沒有雷達訊號，將之設定為紅色與綠色，則淹水區會以藍 色系突顯出來，最後會儲存成有坐標系統的 GeoTiff 影像，我們可以將它匯入 QGIS 與其他地理資訊套疊，以進一步了解變異的區域（圖）。
圖 4-2.1 哈吉貝颱風在日本關東地區的雷達影像淹水判釋地圖
dir %project%\stack.db.TC.HAND.data\Sigma0_%pol%*mst*.hdr /b | gawk -F. "{print \"blue=\"$1}" >temp.rgb
dir %project%\stack.db.TC.HAND.data\Sigma0_%pol%*slv*.hdr /b | gawk -F. "{print \"green=\"$1}" >>temp.rgb
dir %project%\stack.db.TC.HAND.data\Sigma0_%pol%*slv*.hdr /b | gawk -F. "{print \"red=\"$1}" >>temp.rgb
pconvert.exe -f tif -p temp.rgb -o %project%
%project%\stack.db.TC.dim del temp.rgb
(二) 崩塌地災害偵測
⚫
01pre_sub_cal_mul_speckle_bat.xml – 前處理 xml 參數檔（一般 stripmap sar 使用）
⚫
01S1_pre_sub_cal_mul_speckle_bat.xml –前處理 xml 參數檔（Sentinel-1 使 用） 達影像為輸入、興趣範圍（W/E/S/N）為 141.86/142.12/42.62/42.89 與輸出產品名 稱為 ALOS2output 的範例：
設定好之後，打開 cmd，執行
>runGPT_dulPol.BAT
runGPT*.BAT 將執行以下步驟來完成山崩變異區域的標示
set master='ALSO2\0000236705_001001_ALOS2229840850-180825\VOL-ALOS2229840850-180825-HBQR1.1__A'
set slave='ALSO2\0000236706_001001_ALOS2231910850-180908\VOL-ALOS2231910850-180908-HBQR1.1__A'
set subregion='POLYGON ((141.86 42.65, 142.12 42.65, 142.12 42.89, 141.86 42.89, 141.86 42.65))'
set graph=graph set project=p1
REM sar_type = 1: sentinel-1 SAR(VV/VH); sar_type = 2: general stripmap SAR Type(HH/HV);
set sar_type=2
1.
包括裁切、輻射校正、多視處理、斑駁濾波等步驟的影像前處理，對於
gpt %graph%\%pre_process% -Pfile=%master% -Psubregion=%subregion%
-t %project%\master
gpt %graph%\%pre_process% Pfile=%slave% Psubregion=%subregion% -t %projec-t%\slave
gpt %graph%\02stack_bat.xml Pmaster=%project%\master.dim -Pslave=%project%\slave.dim -t %project%\stack
dir %project%\stack.data\Sigma0_%pol1%_mst*.hdr /b | gawk -F.
"{print \"set HH_MASTER=\"$1}" > tmp.cmd call tmp.cmd
dir %project%\stack.data\Sigma0_%pol2%_mst*.hdr /b | gawk -F.
"{print \"set HV_MASTER=\"$1}" > tmp.cmd call tmp.cmd
dir %project%\stack.data\Sigma0_%pol1%_slv*.hdr /b | gawk -F.
"{print \"set HH_SLAVE=\"$1}" > tmp.cmd call tmp.cmd
dir %project%\stack.data\Sigma0_%pol2%_slv*.hdr /b | gawk -F.
"{print \"set HV_SLAVE=\"$1}" > tmp.cmd call tmp.cmd
set exp1=20*log10(%HV_MASTER%/%HH_MASTER%) set exp2=20*log10(%HV_SLAVE%/%HH_SLAVE%)
gpt %graph%\03RVI_bat.xml -Pfile=%project%\stack.dim -PEXP1=%exp1%
-PEXP2=%exp2% -t %project%\stack_RVI
REM Using SRTM 3 sec
gpt Terrain-Correction -Ssource=%project%\stack_db.dim -t
%project%\stack_db_TC REM Using externalDEM
REM gpt Terrain-Correction -Ssource=%project%\stack_db.dim -t
%project%\stack_db.TC.dim PexternalDEMFile=5Mdem_wgs84.tif -PexternalDEMNoDataValue=-3.40282e+38 -PnodataValueAtSea=true
REM gpt Terrain-Correction -Ssource=%project%\stack_db.dim -t
%project%\stack_db_TC.dim PexternalDEMFile=TW20mDEM_wgs84.tif -PexternalDEMNoDataValue=-32767 -PnodataValueAtSea=true
REM gpt Terrain-Correction -Ssource=%project%\stack_db.dim -t
%project%\stack_db_TC.dim PexternalDEMFile=ALOS30mDSM.tif -PexternalDEMNoDataValue=32767 -PnodataValueAtSea=true
5. 以事件前後之雷達植被指標差異，將變異之山崩區以紅色標示。
本範例將事件前設定為紅色，事件後設定為藍色與綠色，由於事件後山崩區 域 RVI 數值的降低，藍色與綠色數值較低，因此會得到一張以紅色突顯山崩變 異的 Geotiff，我們可以將它匯入 QGIS 與其他地理資訊套疊，以進一步了解變異 的區域。
圖 4-2.2 2018 年 9 月 6 日日本北海道地震後的雷達影像崩塌地判釋地圖
dir %output%.data\RVI_SLAVE*.hdr /b | gawk -F. "{print
\"blue=\"$1}" >temp.rgb
dir %output%.data\RVI_SLAVE*.hdr /b | gawk -F. "{print
\"green=\"$1}" >>temp.rgb
dir %output%.data\RVI_MASTER*.hdr /b | gawk -F. "{print
\"red=\"$1}" >>temp.rgb
pconvert.exe -f tif -p temp.rgb %output%.dim del temp.rgb
伍、 緊急災害處置與教育訓練 (一) 緊急災害處置
由於西南風引進大量水氣， 2019/08/15 發生豪雨，包括嘉義縣茶山、高雄市 民權、高雄市那瑪夏國中、高雄市達卡努瓦、嘉義縣表湖、高雄市小林等地皆已 突破 200mm，中央氣象局 15 日早上已針對中南部 10 縣市發布大雨特報，現在 已擴大至 12 縣市，包含台中市、彰化縣、南投縣、雲林縣、嘉義縣市、台南市、
高雄市皆為「豪雨特報」等級(圖 5-1.1)。
圖 5-1.1 中央氣象局日降雨累積圖
使用 Sentinel-1 2019/8/9 與 8/15（18:00）的 SAR 影像進行變遷分析
離災情最近的時間 8/15（Local Tiime: 18:00）剛好有 C 波段 Sentinel-1 的升 軌雷達影像，16 日緊急下載該影像並與事件前 8/9 日之雷達影像進行分析(圖 5-1.2)。
圖 5-1.2 8/15 豪雨事件分析使用之 Sentinel-1 影像範圍。紅框：8/15，白框：
8/9
分析結果顯示，嘉義航空站、嘉義新營至白河的急水溪區域有較明顯的淹水 可能性，而台南善化至玉井區的曾文溪水位有明顯上漲，淹沒原來的高灘地(圖 5-1.3 與圖 5-1.4)。
圖 5-1.3 雷達影像變遷分析顯示嘉義地區可能淹水的區域，藍色區域代表可能 淹水之區域
圖 5-1.4 雷達影像變遷分析顯示曾文溪水體變化的區域，藍色區域代表可能被 水覆蓋之區域
8/9 8/15
8/9
8/15
另外在烏山頭水庫上游發現異常變異訊號，原來屬於水庫水域的區域，在 8/15 影像有異常大的雷達反射能量(圖 5-1.5)。這種情況推測可能會由以下情況 造成：
1. 原來的水體被其他地物覆蓋，造成較大的雷達反射訊號。
2. 局部性水面波浪或降雨造成水面粗糙度增加，造成較多的雷達反射訊號。
3. 局部密集的雷雨胞，使雷達波受到太厚的水氣干擾，在雲層就有反射訊 號。
圖 5-1.5 烏山頭水庫上游異常訊號
8/9 8/15
另外我們也使用雷達植被指標（RVI）進行分析，來看是否有較明顯因為山 崩造成植被變化的區域存在，目前的分析結果並沒有顯示明顯的山崩被偵測出 來(圖 5-1.6)。
圖 5-1.6 雷達影像變遷分析顯示可能山崩的區域，圖片由左至右分別為 8/9 的 影像、8/15 的影像與 RGB 多時變遷分析影像
8/9 8/15
(二) 教育訓練
教育訓練已於 108 年 10 月 28 日辦理完成，訓練時數共計 6 小時。內容包 括(1)雷達衛星影像之基礎理論；(2)崩塌地、堰塞湖等各類林地災害偵測原理及 處理流程；(3)雷達影像處理軟體安裝說明及教學及(4)雷達影像處理軟體時機操 作。
表 5-2.1 教育訓練簽到單
雷達衛星影像輔助林地災害偵測之研究(2/2)
【教育訓練】辦理實況 日期：108 年 10 月 28 日
地點：行政院農業委員會林務局農林航空測量所 201 會議室（上午）、 影像作業室（下午）
圖 5-2.1 教育訓練辦理實況-教育訓練開始
圖 5-2.2 教育訓練辦理實況-教育內容大綱
圖 5-2.3 教育訓練辦理實況-計畫緣起與目的
圖 5-2.4 教育訓練辦理實況-水體偵測原理及處理流程
圖 5-2.5 教育訓練辦理實況-崩塌地偵測原理及處理流程
圖 5-2.6 教育訓練辦理實況-回饋討論
圖 5-2.7 教育訓練辦理實況-雷達影像處理軟體實機操作
陸、 研究成果討論
去年度的研究成果說明只需事件前後兩張影像即可進行分析，這兩張影像 的拍攝參數（衛星、軌道、產品種類）差異越少，越能減少系統上的誤差。通常 會在影像搜尋網頁上選擇時間最接近與範圍最接近的兩張影像，降低隨時間變 化的影響，而產品的種類不論是 SLC 影像或 GRD 影像都可以，產品的等級需 要選擇未經過地形校正的等級，以方便使用更高解析度的 DEM 進行地形校正。
今年度的研究使用較長時間段的雷達影像進行分析，得到地物長期穩定的特徵 雷達訊號，除了可以在維持影像解析力進行斑駁雜訊的濾除，也可加強事件造成 的訊號差異性。但通常人工地物區域或森林茂密區域（如熱帶雨林）的雷達特徵 訊號在長時間是穩定的，可使用較長時間段的資料進行分析，但對於四季變化明 顯、植生變化明顯、人為變化明顯（農地灌溉與耕作）的區域，使用較長時間的 資料，也等同於加入了這段時間內更多的變化因子，這些變化因子在雷達影像造 成的影響可能會大過於本來研究想要觀測的事件，造成分析上的困難。如果可以 選擇的話，原則上盡量選擇離事件時間接近（約一個月內且時間間隔越短越好）
的影像進行分析，俾以降低時間降相關(temporal decorrelation)的干擾。
目前多時期分析的研究結果顯示，山崩偵測的效果較好，雖然需要使用 3 張 影像，但只需使用一般雷達衛星可以獲得的平行極化影像（HH 或 VV）即可進 行分析，可以使用的衛星影像上有更多的選擇。研究的結果也顯示高解析度 X 波段的雷達影像對於山崩後的雷達特性變化相當明顯，只需要使用單極化影像 的差異進行事件前後 RGB 套色處理，就能有效地把山崩區域偵測出來。本研究 也從偏移偵測法在雷達影像的應用發現位移量大與低相關性的區域可以比對到 山崩的區域，而且只需事件前後兩張影像即可分析，但由於相關性低的特性，顯 示其數值的可信度並不高，需要更進一步的分析，才能了解雷達影像利用偏移偵 測法在山崩偵測的能力。
而在淹水偵測方面，只需要使用 2 張雷達影像以 RGB 將水體變化造成的淹
水區以藍色呈現，即可快速標示淹水區域，不需要額外再多一張影像進行多時期
分析。
我們利用 Deutscher et al.（2017）方法分析火災造成的林地變遷，該方法主 要是考慮時間序列下高變異係數與背向散射係數趨勢負值的區域，再輔以相關 地理資訊作為遮罩，移除不相關區域的影響。但測試結果並不理想，臺灣臺中大 甲溪的火災區在變異係數與背向散射係數趨勢都沒有特別的差異，無法辨識出 火災的區域；美國加州火災的案例，雖然可以辨識出火災的區域，但也有許多非 火災區域被包括在其中。原因可能包括：
(1) 區域的差異：分析方法之文獻中的區域屬溼熱帶，植被終年茂密，火災所 造成的區域差異相當明顯，但是像臺灣或加州四季變化明顯的地方，因時 間造成的改變相當大，火災造成區域的變化不容易在雷達影像中分辨出 來。
(2) 地形的差異：文獻中的區域是剛果河盆地，地勢較平緩，而像台灣山區地 形起伏相當劇烈，就算經過地形平坦化的輻射校正，也無法完全移除地形 起伏對側視雷達系統造成的數值差異影響。
(3) 人為活動的影響：人為活動如灌溉、耕作也會造成變異，文獻中的方法有 利用額外的地理資訊作為遮罩，去除森林以外因人為活動造成的變異區 域。加州火災的案例就有許多類似耕地的區域被視為變異區，而本研究目 前沒有相關資訊以進行排除。
以上原因都會造成該方法在林地變異判釋能力的差異。加州火災分析結果 可以知道該方法的確可以實際應用在火災災跡偵測，雖然 2016 年 4 月的大甲火 災案例中沒有良好的結果，但由於目前在臺灣分析的案例有限，針對不同林相、
不同地形與不同得林地變異，可能都存在差異性，雷達影像在這一塊的研究，未
能仍需要更多實際案例的分析以了解其特性與可行性。
柒、 結論與未來方向 (一) 結論
雷達影像多時期的分析
1. 多時期濾波處理可以在保留影像空間解析度的情況下，有效降低斑駁雜 訊，提高影像特徵。
2. 多時期變遷分析在山崩有良好的偵測效果，雖然需要使用 3 張影像，但 只需要較易獲得的平行極化影像（HH 或 VV）即可進行分析。
3. 多時期分析可以應用在火災跡地的辨識，以美國加州火災為例，雷達影 像林地變遷偵測在火災跡地的辨識有一定的成果，但在臺灣仍需要更 多實際案例的分析以了解其特性與可行性。
4. 淹水區域判釋案例直接利用背向散射係數差異的效果比多時期分析的 成果要好。
偏移偵測法的應用
雷達影像在偏移偵測法的應用上，除了可以偵測斷層兩側的地表水平 變形，也能從低相關性與異常大變形量的區域反應山崩發生的位置，顯示 有能力作為山崩判釋的工具。
林地災害之雷達衛星影像分析判讀準則
本計畫兩年度的研究中使用的雷達衛星影像種類包括 L 波段 ALOS 的 SLC 影像、C 波段 Sentinel-1 的 GRD 影像與 X 波段 COSMO-SkyMed 的 SLC 影像，在使用這些影像做為案例分析的條件下顯示：
1. 雷達影像的背向散射係數的差異就足以作為洪水溢淹的判斷指標。
2. 具有雙極化影像時，可使用 RVI 差異可作為快速山崩判釋的指標。
3. 高解析度 X 波段的單極化影像有能力直接透過背向散射係數的強度差 異偵測山崩。
4. 多時變遷偵測與偏移偵測的方法都有能力做為山崩判釋的分析工具。
建立以開源軟體為主之雷達衛星影像林地災害偵測流程
根據過去林地災害偵測的研究結果進行自動化處理流程的建置，已完
成淹水災害偵測與崩塌地災害偵測，可提供緊急災害應變時的災害發生區
域偵測能量。
(二) 未來方向
在緊急應變方法已建立的情況下，災害發生時是否能快速取像作為分 析是未來需要考慮的方向。目前有新型合成孔徑雷達衛星包括加拿大的 RADARSAT Constellation Mission(C 波段)與芬蘭的 ICEYE(X 波段)有能力 在 24 小時內取得相同區域的影像，並且 ICEYE 未來預計要達到 3 小時的 再訪週期，是未來在防災應變使用影像上的良好選擇。
另外有鑑於遙感地球觀測數據可從許多自由訪問的開放性儲存庫中獲
得，但用戶需要執行一系列複雜的預處理步驟，否則無法實現對數據的有
效利用。如果能透過多維數據集(Data Cube, DC)的概念，建構可方便存儲
和分析大量柵格數據的地理數據基礎結構，則能降低由這些大數據挑戰引
起的障礙，並以分析就緒形式提供對大型時空數據的訪問。未來如果能建
立台灣地區的 Sentinel-1 SAR 與 Sentinel-2 光學影像的多維數據集，提供
經過輻射校正與地形校正的分析就緒數據，是有利於林地生態系統多樣性
的監測訊息(時間與空間)與防災緊急應變時的災害資訊提供。
英文專有名詞對照
Radar Vegetation Index, RVI：雷達植被指標。利用交叉極化與平行極化的差異 性，突顯植物在雷達回波訊號上的特徵。
Backscatter coefficient：背向散射係數。是雷達影像成像的紀錄值經過輻射校正 後之數值，代表目標物的雷達反應，使合成孔徑雷達可提供定量遙測的資 訊。校正後的資料可與不同時間、不同雷達、不同感測參數、不同地區的 目標訊息相互比較與匹配。
Coregistration：套疊。每張影像拍攝時的軌道與影像大小有差異，多時期影像 需經過影像強度的相關性計算，找出屬於相同地物的位置作為控制點，將 多張影像套疊至其中一幅主影像，形成影像範圍與網格大小一致的影像對。
Gamma naugh：背向散射係數的一種，更進一步校正地形起伏對雷達影像強度 的影響。
Interferometric Coherence：干涉同調性。兩張複數影像進行干涉處理時，其相 位相關的程度。
Intensity correlation：強度相關性。兩張影像強度相關的程度。
Multi-looking processing：多視處理。以多觀點平均來降低雜訊的方法。
Multi-temporal filtering：多時期濾波法。藉由多時期影像作為濾波處理的依據，
保留空間解析力來進行斑駁雜訊抑制的方法。
Normalized differences, ND：正規化差異指標。用來計算事件前後影像對差異性 變化的演算法。
Speckle noise：斑駁雜訊。雷達影像因其成像特性，地面像元內部分布大量的隨 機散射體，每個散射體的回波特性皆由地物結構、粗糙度、介電性質、雷 達波長、極化方式及雷達波入射角等綜合因素共同決定，在同一個解析力 單元內，由於雷達波波長與地表粗糙程度相近（公分級） ，造成的干涉現象 會使所有散射體回波相互的疊加或抵消，彼此之間相互干擾的結果則導致 雷達影像上出現隨機分布的亮點或暗點，此種雷達影像特性即稱為斑駁雜 訊。
Speckle filtering：以數學模式或統計模式演算法降低斑駁雜訊的方法。
參考文獻
Arciniegas, G. A., Bijker, W., Kerle, N., & Tolpekin, V. A. (2007). Coherence-and amplitude-based analysis of seismogenic damage in Bam, Iran, using ENVISAT ASAR data. IEEE Trans. Geosci. Remote, 45(6), 1571-1581, doi:
10.1109/TGRS.2006.883149.
Bruniquel, J., & Lopes, A. (1997). Multi-variate optimal speckle reduction in SAR imagery. Int. J. Remote Sens., 18(3), 603-627, doi: 10.1080/014311697218962.
Deutscher, J., Gutjahr, K., Perko, R., Raggam, H., Hirschmugl, M., & Schardt, M.
(2017). Humid tropical forest monitoring with multi-temporal L-, C-and X-band SAR data. In: Analysis of Multitemporal Remote Sensing Images (MultiTemp), 2017 9th International Workshop, 1-4, doi: 10.1109/Multi-Temp.2017.8035264.
Foumelis, M. (2015). ESA Sentinel-1 Toolbox Generation of SAR Backscattering Mosaics, Course Materials. In: Proceedings of the 6th ESA Adv., Bucharest, Romania.
Konishi, T., Suga, Y. (2018). Landslide detection using COSMO-SkyMed images: A case study of a landslide event on Kii Peninsula, Japan. Euro. J. Remote Sens., 51(1), 205-221, doi: 10.1080/22797254.2017.1418185.
Lee, J. S. (1981). Speckle analysis and smoothing of synthetic aperture radar images.
Comput. Graphics Image Processing, 17(1), 24-32, doi: 10.1016/S0146-664X(81)80005-6.
Lopes, A., Nezry, E., Touzi, R., & Laur, H. (1993). Structure detection and statistical adaptive speckle filtering in SAR images. Int. J. Remote Sens., 14(9), 1735-1758, doi: 10.1080/01431169308953999.
Löffler, E. (2013). Geographie und Fernerkundung: eine Einführung in die geographische Interpretation von Luftbildern und modernen Fernerkundung- sdaten. Springer-Verlag.
Matsuoka, M., & Yamazaki, F. (2000). Characteristics of satellite SAR images in the areas damaged by earthquakes. In: Geoscience and Remote Sensing Symposium, 2000. Proceedings. IGARSS 2000. IEEE 2000 International, 6, 2693-2696, doi:
10.1109/IGARSS.2000.859684.
Mirelva, P. R., & Nagasawa, R. (2017). Single and Multi-temporal Filtering Comparison on Synthetic Aperture Radar Data for Agriculture Area Classification.
In: Proceedings of the International Conference on Imaging, Signal Processing and Communication, 72-75, doi: 10.1145 / 3132300.3132316.
Plank, S. (2014). Rapid damage assessment by means of multi-temporal SAR—A comprehensive review and outlook to Sentinel-1. Remote Sens., 6(6), 4870-4906, doi: 10.3390/rs6064870.
Quegan, S., Le Toan, T., Yu, J. J., Ribbes, F., & Floury, N. (2000). Multitemporal ERS
SAR analysis applied to forest mapping. IEEE Trans. Geosci. Remote Sens., 38(2),
741-753, doi: 10.1109/36.842003.
Rosu, A. M., Pierrot-Deseilligny, M., Delorme, A., Binet, R., & Klinger, Y. (2015).
Measurement of ground displacement from optical satellite image correlation using the free open-source software MicMac. ISPRE J. PHOTOGRAMM, 100, 48-59.
Stewart, C. (2016) Exercise Sentinel-1 Processing, Course Materials. In: Proceedings of the 8th ESA Trng. Crse. Radar Optical Remote Sens., Cesis, Latvia.
Stramondo, S., Bignami, C., Chini, M., Pierdicca, N., & Tertulliani, A. (2006). Satellite radar and optical remote sensing for earthquake damage detection: results from different case studies. Int. J. Remote Sens., 27(20), 4433-4447, doi:
10.1080/01431160600675895.
Suga, Y., Takeuchi, S., Oguro, Y., Chen, A. J., Ogawa, M., Konishi, T., & Yonezawa, C. (2001). Application of ERS-2/SAR data for the 1999 Taiwan earthquake. Adv.
Space Res., 28(1), 155-163, doi: 10.1016/S0273-1177(01)00334-9.
Takeuchi, S., Suga, Y., Yonezawa, C., & Chen, A. J. (2000). Detection of urban disaster using InSAR. A case study for the 1999 Great Taiwan Earthquake. In: Geoscience and Remote Sensing Symposium, 2000. Proceedings. IGARSS 2000. IEEE 2000 International, 1, 339-341, doi: 10.1109/IGARSS.2000.860512.
Tzeng, Y. C., Chiu, S. H., Chen, D., & Chen, K. S. (2007). Change detections from SAR images for damage estimation based on a spatial chaotic model. In:
Geoscience and Remote Sensing Symposium, 2007. IGARSS 2007. IEEE International , 1926-1930, doi: 10.1109/IGARSS.2007.4423203.
Valkaniotis, S., Ganas, A., Tsironi, V. & Barberopoulou, A. (2018). A preliminary report on the M7.5 Palu earthquake co-seismic ruptures and landslides using image correlation techniques on optical satellite data, doi:
10.5281/zenodo.1467128.
Wang, C., Mao, X., & Wang, Q. (2016). Landslide Displacement Monitoring by a Fully Polarimetric SAR Offset Tracking Method. Remote Sens., 8(8), 624, doi:
10.3390/rs8080624.
Yaseen, M., & Anwar, S. (2013). Sensitivity analysis of sub-pixel correlation technique for measuring coseismic displacements using a pair of ASTER images.
J. Asian Earth Sci., 62, 349-362, doi: 10.1016/j.jseaes.2012.10.015.
Yonezawa, C., & Takeuchi, S. (1999). Detection of urban damage using interferometric SAR decorrelation. In: Geoscience and Remote Sensing Symposium, 1999. IGARSS'99 Proceedings. IEEE 1999 International, 2, 925-927, doi: 10.1109/IGARSS.1999.774487.
Yonezawa, C., & Takeuchi, S. (2001). Decorrelation of SAR data by urban damages caused by the 1995 Hyogoken-nanbu earthquake. Int. J. Remote Sens., 22(8), 1585-1600, doi: 10.1080/01431160118187.
王志添，陳錕山（2003） ，雷達影像多時處理研究， 航測及遙測學刊 ， 8，45-56。
邱俊穎，胡植慶與謝嘉聲（2019），合成孔徑雷達衛星影像於颱風豪雨後洪水溢
淹及堰塞湖之偵測研究，2019 臺灣地球科學聯合學術研討會 。
胡植慶，謝嘉聲，邱俊穎與李秀芳（2018），雷達衛星影像輔助林地災害偵測之 研究(1/2)， 行政院農業委員會林務局農林航空測量所委託研究計畫成果報 告書 ，共 163 頁。
極隼科技股份有限公司（2012） ，雷達遙測技術於林地災害判釋之研究， 農林航
空測驗所委託研究計畫成果報告書 。
委員審查意見
2. 利用SNAP、GRASS-GIS、R語言等開源軟體，
在工作計畫書並沒詳述這些開源軟體使用部
（ Grey Level Co-occurrence Matrix, GLCM）」，亦有中文翻譯為「灰階共生矩
期初審查意見 處理情形
期初審查意見 處理情形
（Positive Predictive Value，PPV）」兩者的 計算公式相同？建議一併寫出引用這些計 度（Incidence angle）。
12. 第4頁式 4-4「i：Moving Windows中第幾個 像元的指標」請補充說明擬採用的Window
16. 第6頁圖4-1.3的多時變遷RGB影像僅使用 兩個日期（2009/07/08及2009/08/23）的雷達
期初審查意見 處理情形 標（Coseismic Coherence Difference, CCD）」
的計算公式。
感謝委員建議，已補充。
22. 第 17 頁 「 次 像 元 相 關 法 （ Sbupixel correlation ）」 應訂 正 為「次 像元相 關法
（Subpixel correlation）」。
感謝委員建議，已更正。 window size計算得到的平均雷達植被指 標、平均RVI值呢？
期初審查意見 處理情形 究之案例。
3. 本計畫使用之Offset tracking有無分析其精 度及能力為何？之前所使用SNAP進行分
期初審查意見 處理情形
1. 本案除使用SNAP、GRASS-GIS、R語言等 開源軟體外，是否也有考量使用農航所目
2. 本所前於102年購買COSMO-skyMed雷達 衛星影像，惟無法從影像上直接分析出有
期中委員審查意見
2011-08-19 至 2012-05-09 期間的水平變形量，平均速率為1.6-10.9 mm/day(p. 35)。
5. 部分測試使用5m DEM進行地形校正（例如 第10頁），部分測試卻使用SRTM DEM進行 地形校正（例如第24頁），其原因為何？請 說明。
主要是因參考Deutscher et al.（2017）的研 究進行多時林地變遷偵測時，需要進行地
期中審查意見 處理情形 時間為2018年8/24與9/5（UTC+0）。」，所使 用之影像皆為災前，是否誤植，請確認。
10. 第43頁圖4-0.1及圖4-0.2兩者的圖標題太冗 長，建議宜修訂並將其冗長標題文字移到 TWD97及TWVD 2001坐標系統中？
如何計算式（2-1）的入射角？ 粗略 計算？還是嚴密的精確計算？例如：
期中審查意見 處理情形
「本案例分析結果與 Deutscher et al.（2017）的研究有不同的地
A. Plank(2014)多時期的災害偵測最少使 用3 張影像，所以採用離災害發生最近
期中審查意見 處理情形 Petroleum Survey Group，中文名稱為歐
洲 石 油 調 查 組 織
期中審查意見 處理情形 TerraSAR-X=96% 與 Sentinel-1=98%。」，如何求得準確率？
期中審查意見 處理情形 關法（Subpixel correlation）來進行分 析。」(p. 33)。
期中審查意見 處理情形
（2012）」、第31頁「Wang et al., 2016）」。
感謝委員建議，已更正。
2. 第8頁斑駁處理（Speckle noise）部分，因每 1張SAR影像range不同，不用正規化方式而
期中審查意見 處理情形 能與 MicMac 的差異，MicMac 在處理速 度與成果都比SNAP 內建功能要好，因此 目前考慮使用MicMac 來進行處理。
6. 第34頁有關圖2-4.4是否在SAR影像的角落 處，以致於東西方向的位移量偵測顯示較
期中審查意見 處理情形 林火災部分係以Deutscher et al.（2017）之 研究進行探討，建議於會議中詢問補充其
期中審查意見 處理情形 像對，應更明確定義「大小一致」之空間意
涵。
3. 有關林地災害之雷達衛星影像分析判讀準 則建立，於報告中有相當詳細的研究說明 與流程，惟因混列較不容易快速查找，建議 於該章節最後加上列表，依各類災害及分 析需求，列出建議之雷達衛星影像、波段及 極化組合、相關門檻值及指標，以及配合的 軟體程式功能、分析結果判讀方法等，並以 SOP手冊方式呈現，俾利於快速查找及操 作。
感謝委員建議。
李課長鴻德
1. 目前計畫內容有關衛SAR影像材料取用多 是以ALOS(L-band)及Sentinel-1(C-band)為 主，建議可增加不同頻段的SAR影像使用 分析說明，例如解析度較佳之SAR影像(X-band)。
感謝委員建議，已於山崩災害判釋列入 X band 影像分析案例(p.20~22)。
2. 建議手冊增加不同災害影像來源的範例說 明，使用衛星的參數也一併納入。
感謝委員建議。
期末委員審查意見
期末審查意見 處理情形 Range 與 Azimuth 為座標系統的雷達影像，
轉換為經緯度或分帶座標。以目前的案例
5. 第iv頁：遺漏第87至90頁的圖1至圖6的圖目 錄。且圖編號之格式應與各章節一致，如「圖
期末審查意見 處理情形 (classification error matrix, sometimes called a
confusion matrix or a contingency table)。
7. 第51頁「研究中嘗試找出可以用腳本操作來
「在一個解析力單元(resolution cell) 內」。
感謝委員建議，已修正。
(12) 「Speckle filtering：以數學模式演算 法降低斑駁雜訊的方法」：建議修訂
期末審查意見 處理情形 以使用的影像格式 UInt16，數值在 0-65535 之間。
期末審查意見 處理情形
本團隊使用 SNAP 內建功能 Coregistration 進行套疊。
19. 第12、17、36頁「地形校正：使用5m DEM進 行地形校正，將影像轉換為地理坐標（WGS 84）」，請補充說明研究團隊是如何進行地形 校正？以及為何不轉換為TWD97。
本 團 隊 使 用 SNAP 內 建 功 能 Terrain Correction 進行地形校正。
SNAP 使用 WGS84 座標格式進行處理。 響，轉換為gamma naught」：貴團隊是否考慮
感謝委員建議，這些差異有可能存在的，
但在本研究當中並沒有更適合的資料可以 使用。
期末審查意見 處理情形 (temporal decorrelation)的干擾」。
感謝委員建議，已修正。
期末審查意見 處理情形 抑制（Speckle Reduction）」。
感謝委員建議，已修正。
3. 有關airbone SAR或UAV SAR之可靠度、精度 為何，請補述之。
期末審查意見 處理情形
災害發生後要快速應變，短時間內也 只能取得一張影像，因此整體上是使 用災前2 張與災後 1 張影像。
2) 使用多時期影像分析需要考慮季節明 顯變化且植被有變化之區域是否會影 像判釋目標的變遷偵測。
2. 期末報告中大甲溪事業區森林火災案例不論 是在光學影像及雷達衛星影像都無法顯示火 災跡地，可否更換其他案例？
雷達影像林地變遷偵測在火災跡地的辨識 有一定的成果，但在臺灣仍需要更多實際 案例的分析以了解其特性與可行性。
(二) 許技正玉君
1. 期末報告較多以強度之差異性來分析討論，
另外有關干涉同調性是應用在哪方面？請補 充說明。
在都會區會有較高的干涉同調性，因此在 災害發生時的同調性變化，可作為變遷偵 測之依據。
2. 在淹水判釋部分以VV極化會得到較好的成 果，而在第47頁表3.1中的使用指標所列，可 否加入使用影像之優先順序？
感謝委員建議，已補充。
(三) 葉課長堃生
1. 第49頁所標示的圖2-2.4及圖2-4.7是否為誤 植？抑或引用前章節的圖來說明？
引用前章節的圖來說明
行政院農業委員會林務局 農林航空測量所
108年度委託專業服務採購案
雷達衛星影像處理操作手冊
「雷達衛星影像輔助林地災害偵測之研究 (2/2)」
全 程 計 畫：自107年3月至108年12月止
本年度計畫：自108年1月至108年12月止
執行單位：國立臺灣大學
中 華 民 國 108 年 12 月 10 日
目錄
一、 軟體安裝與設定 ... 1 (一) SNAP ... 1 (二) QGIS 與 GRASS-GIS ... 3 (三) R 與 R-Studio ... 5 (四) GAWK ... 6 二、 軟體基本操作說明 ... 7 (一) SNAP ... 7 (二) QGIS 與 GRASS-GIS ... 18 (三) R 與 R-Studio ... 20 (四) GAWK ... 22 三、 雷達影像災害偵測之處理腳本設置與操作 ... 23 (一) 淹水災害偵測 ... 23 (二) 崩塌地災害偵測 ... 28
一、 軟體安裝與設定 (一) SNAP
SNAP 是由歐洲太空總署開發的一款用於雷達影像資料處理的 開 源 軟 體 ， 目 前 版 本 為 7.0 ， 安 裝 檔 案 可 於 連 結 [ http://step.esa.int/
main/download/snap-download/ ]進行下載，下載後只要點擊下載檔案即可
進行安裝。
圖 1 SNAP 下載網頁畫面
安裝完成後，需要設定使用硬碟空間作為暫存，以防止計算時記
憶體不足的計算失敗情況：開啟 SNAP→Tools→Option→S1TBX→勾
選 Use FileCache in readers to conserve memory。
圖 2 SNAP FileCache 設定
SNAP 軟體有一些較常用的工具視窗在安裝完成時，並沒有打開，
需要手動開啟，包括同時比對不同影像的 Navigation、顯示 SAR 影像 於世界地圖位置的 World Map 與顯示指標上影像數值的 Pixel Info。
這些功能都在 View->Tool Windows 當中。
圖 3 SNAP 常用工具視窗設定
(二) QGIS 與 GRASS-GIS
QGIS 原稱 Quantum GIS 是一個開源的 GIS 軟體，目前版本為 3.8，安裝檔案於連結[https://qgis.org/en/site/forusers/ download.html]，
選擇單機版的版本進行下載。
圖 4 QGIS 下載網頁畫面
下載後只要點擊下載檔案即可進行安裝， QGIS 3.8 安裝時也會順 便安裝 GRASS-GIS 7.6.1，不需要另外安裝 GRASS-GIS。安裝完成後，
需先設定 GRASS-GIS：
⚫
開啟 GRASS GIS 7.6.1
⚫
選擇 GRASS 資料庫路徑[C:\GRASSDB]
⚫
選擇 GRASS 地點 [TW]
⚫
選擇 GRASS 圖集[PERMANENT]
⚫
以上設定完成時，會打開設定好 GRASS GIS 環境變數的命
令視窗（cmd） ，執行指令[g.extension r.stream.distance]新增
GRASS 模組
圖 5 GRASS-GIS 基本設定
圖 6 新增 GRASS-GIS 模組 r.stream.distance
(三) R 與 R-Studio
R 語言與 Python 同為目前開源數據分析社群當中相當受歡迎的 開源編程語言。 R-Studio 是為 R 語言設計的整合開發環境。安裝檔案 於 連 結 [ https ： //cran.r-project.org/bin/wind ows/base/, https ： //www.rstudio.com/products/rstudio/download/ ]，先安裝 R 語言，再安裝 R-Studio。
圖 7 R 語言與 R-Studio 下載網頁
(四) GAWK
GNU AWK（GAWK）是歷史悠久的 AWK 編程語言的開放源代 碼實現，可用於所有的 UNIX 系統。AWK 語言是一種 UNIX 備用 工具，它是一種功能強大的文本操作和模式匹配語言，特別適用於進 行信息檢索，這使得它非常適合用於當今的數據庫驅動的應用程序。
因為它集成於 UNIX 環境，所以可以設計、構建和快速地執行完整的 工作程序，並且立即就能得到結果。安裝檔案位於連結網頁當中 [http://gnuwin32.sourceforge.net/packages/gawk.htm] 。 選 擇 Complete package, except sources 的版本，執行安裝即可；或者選擇 Binaries 版 本，解壓縮下載後之檔案，將解壓縮後 bin 資料夾內容複製至 C:\Windows\System32 之內
圖 8 GAWK 下載網頁
二、 軟體基本操作說明
主要介紹 SNAP 如何進行雷達衛星影像的基本處理，另外說明輔 助的軟體 QGIS、GRASS-GIS、R、R-Studio 與 GAWK 如何運作。
(一) SNAP
本章節目標在提供有關 SNAP 處理 SAR 數據的基本操作。將學 習到幾本的產品的輸入與檢視、輻射校正（Calibrate）、多視處理
（ Multilooking）、斑駁雜訊抑制（Speckle Reduction）與地形校正
（Terrain Correction）。
Open a Product
Step 1 – Open a product: 使用 Open Product 按鈕，找到產品的 zip 檔，或 SNAP 已經處理過的.dim 檔案。或使用工具列上
File->Import->Sar Sensors->選擇衛星種類，再進入產品資料夾，即可選擇該產品 對應地開啟檔案。
圖 9 打開 Sentinel-1 zip 檔案
圖 10 輸入其他 SAR 衛星資料
在 Products View 你可以看見已經輸入 SNAP 的產品，每一個產 品都包括相關資訊在 metadata 與影像資料在 bands。
圖 11 Products View
雙擊 Amplitude_VH 可以顯示影像。
圖 12 Amplitude_VH Band
此產品為 Sentinel-1 Ground Range Detected（GRD）產品，已經 處理過聚焦、多視與地距投影，不再具有相位資料。影像 Y 方向微衛 星飛行的方向， X 方向為 Range 方向，也可以理解為雷達波射出的方 向。
Calibrating the Data
為了正確處理 SAR 數據，首先需要對資料進行輻射校正，尤其 在處理需要鑲嵌的資料。透過輻射校正，使影像的像素可以代表反射 物體的背向散射訊號，對於定量使用 SAR 資料非常重要。
Step 2 – Calibrate the product: 從 工 具 列 的
Radar->Radiometric->Calibrate
圖 13 Radiometric Menu
接下來要選擇輸入產品（Source）與輸出產品（Target Product） ， 再到參數（Processing Parameters）選擇要輸出的輻射校正種類，預設 為 Sigma0。
圖 14 Calibration Dialog
Multilooking
透過對 Range 方向與 Azimuth 方向的平均，可以降低雜訊，提高 輻射訊號的解析力，但空間解析度會降低。
圖 15 Multilooking an SLC Product
Step 3 – Multilook the Product: 從工具列的 Radar->SAR
Utilities->Multilooking
圖 16 Select Multilooking
設定 Range 與 Azimuth 方向多視的數量。
圖 17 Set Number of Looks
多視處理將使用 range 與 azimuth 方向 3 x 3 的平均來產生一個像 素，像素的平均 ground range 大小為 30.02 m。
當按下 Run 就會進行處理，完成後，可以在 Products View 看到 成果。
圖 18 多視處理前後的差異
Speckle Reduction
透過斑駁濾波，可降低雜訊，但代價是特徵模糊與空間解析度降 低。
Step 4 – Speckle Filtering: 從 工 具 列的 Radar->Speckle
Filtering->Single Product Speckle Filter
圖 19 Select Single Product Speckle Filtering
從 Speckle Filtering 視窗中，選擇 Refined Lee speckle filter。
圖 20 Select Refined Lee
按下 Run 進行處理。完成後，可以在 Products View 看到成果。
圖 21 Speckle Filter 前後的差異
Terrain Correction
地形校正將通過使用數值高程模型（DEM）校正 SAR 幾何失真 來對影像進行地理編碼，產生地圖投影的產品。
Step 5 – Terrain Correction: 從工具列的 Radar->Geometric->Terrain Correction->Range-Doppler Terrain Correction
圖 22 Select Range-Doppler Terrain Correction
預設使用 SRTM 3 sec DEM，可選擇自訂的數值高程，坐標系統
須為 WGS84，檔案格式為 GeoTiff。
圖 23 Terrain Correction Dialog
按下 Run 進行處理。完成後，可以在 Products View 看到成果。
圖 24 Terrain Correction 前後的差異
SAR 影像格式為 Int32，數值分布相當廣，且有相當多細節集中 於一小區段數值之間，不容易呈現其差異，可轉換為 dB 來呈現。在 Amplitude_VH 滑鼠右鍵->Linear to/from dB
圖 25 Select Linear to/from dB
新的影像將會是 10*log10（Amplitude_VH），完成後，雙擊產生
的新影像 Amplitude_VH_db 可以在 Products View 看到成果。
圖 26 dB 處理前後的差異
(二) QGIS 與 GRASS-GIS
本研究使用的是 QGIS 安裝後所附帶的 GRASS-GIS。GRASS GIS 為指令操作模式，使用 GRASS 必須了解各個指令的功能特性而選
系統資料管理模組(General Data Management) 網格資料處理模組(Raster Data Manipulation) 向量資料處理模組(Vector Data Manipulation) 區位點資料處理模組(Site Data Manipulation) 遙測影像處理模組(Imagery Data Manipulation) 外部資料處理模組(External Data Manipulation)
地圖展現模組(Display Graphics)
地圖印製模組(Paint Paper Graphics)
相關的技術文件可以從 GRASS 網站上查閱到
GRASS-GIS 指令運作範例
可將要處理的整個流程存成批次檔（BAT） ，透過 grass*.bat 來執 行整個流程。
GRASS-GIS 批次處理運作範例
:: 輸入 DEM
r.external input=%demfile% output=dem0 -o --overwrite :: 設定投影與範圍
g.proj epsg=4326 -p g.region raster=dem0 -p :: 移除小於 0
r.mapcalc "dem = if(dem0<0, null(), dem0)" --overwrite :: 計算 HAND raster
r.watershed elevation=dem drainage=drainage stream=streams threshold=%threshold% --overwrite
r.thin input=streams output=streams_thin
r.to.vect input=streams_thin output=streams_line type=line --overwrite
r.stream.distance stream_rast=streams direction=drainage elevation=dem method=downstream difference=above_stream --overwrite
r.mapcalc "HAND = if(above_stream<0,0,above_stream)" --overwrite :: 輸出
r.out.gdal input=HAND out=%output% --overwrite format=ENVI -f ::刪除不需要的資料 g.remove type=raster
name=dem0,dem,HAND,above_stream,drainage,streams,streams_thin,str eams_line -f g.remove type=vector name=streams -f
set GRASS76="C:\Program Files\QGIS 3.6\bin\grass76.bat"
set GISDBASE=C:\Users\choin0207\Documents\grassDB set tempLC=tw set demfile=%1
:: create tmp location for run
%grass76% %GISDBASE%/%tempLC%/PERMANENT --exec GRASS_JOB.bat
(三) R 與 R-Studio
R 語言與 R-Studio 安裝完成後，打開 R-Studio 軟體，主要分為四 大區。在 Console 區可輸入指令執行工作，也可以在 Workspace 開啟.r 檔案，將工作內容程式碼寫在裡面。 Environment 與 History 可以檢視 已經存在的變數與歷史命令。 Files、Plots、Packages、Help、Viewer 則 可檢視工作目錄、繪圖輸出、新增套件等功能。
圖 28 R-Studio 工作環境畫面
R 語言程式操作範例
#匯入函數庫
library(raster) library(rgdal)
#匯入自訂 function
source('X:\\R-code\\IntCorr.r')
#設定工作目錄
setwd(choose.dir('choose dir'))
#輸入影像
i1<-raster(choose.files('*.img','select img1')) i2<-raster(choose.files('*.img','select img2')) i3<-raster(choose.files('*.img','select img3'))
#計算 w<-19
c1<-intenCorr(i1,i2,w) c2<-intenCorr(i2,i3,w) ND<-(c1-c2)/(c1+c2)
#輸出 GTiff
writeRaster(c1,paste0('c1_W',w,'.tif'),'GTiff',overwrite=T,NAflag=
NaN)
writeRaster(c2,paste0('c2_W',w,'.tif'),'GTiff',overwrite=T,NAflag=
NaN)
writeRaster(ND,paste0('ND_W',w,'.tif'),'GTiff',overwrite=T,NAflag=
NaN)
(四) GAWK
gawk 為指令操作模式，可以用來處理檔案的資料，以下示範三 種處理的操作
1. 列出檔案裡的第一、第二與第五行 gawk "{ print $1, $2, $5 }" file.txt
2. 改變分欄符號為 , 並列出第一與第二欄 gawk -F, "{ print $1,$2 }" file.txt
3. 跳過 5 行檔頭，列出所有資料 gawk "NR>5 {print $0}" file.txt
本研究利用 gawk 來處理檔案內的文字。由於每次處理的 SAR 影 像產品檔案名稱略有差異，要進一步進行自動繪圖處理時，必須結合 dir 與 gawk 指令輸出專屬此 SAR 產品的繪圖設定檔。
gawk 程式操作範例
# [將 dir 指令找到的檔案與字串結合，輸出到新檔案]
dir %project%\stack_RVI_TC.data\RVI_SLAVE*.hdr /b | gawk -F.
"{print \"blue=\"$1}" >temp.rgb
dir %project%\stack_RVI_TC.data\RVI_SLAVE*.hdr /b | gawk -F.
"{print \"green=\"$1}" >>temp.rgb
dir %project%\stack_RVI_TC.data\RVI_MASTER*.hdr /b | gawk -F.
"{print \"red=\"$1}" >>temp.rgb
# [將產品 stack_RVI_TC.dim 根據 temp.rgb 之設定，輸出 RGB 影像成為 GeoTiff 格式]
pconvert.exe -f tif -p temp.rgb -o %project%
%project%\stack_RVI_TC.dim
三、 雷達影像災害偵測之處理腳本設置與操作 (一) 淹水災害偵測
為解決 SAR 影像應用於淹水實際操作之複雜性，目前建構一套 由以腳本操作來想對自動化的緊急災害偵側標準流程，需要輸入事件 前與事件後影像的檔名，處理的結果可產生一張由藍色系突顯可能淹 水變異區域的變遷偵測影像，做為緊急應變的依據。
整個處理程式包括以下檔案：
⚫
runGPT_General_1.BAT - 使用 Sentinel-1 以外 SAR 影像的串聯 程序之主程式
⚫
runGPT_S1_1.BAT -使用 Sentinel-1 單張 SAR 影像的串聯程序 的主程式
⚫
runGPT_S1_2.BAT -使用 Sentinel-1 兩張 SAR 影象的串聯程序 的主程式
⚫
orb.cal.xml - Sentinel-1 專用，給予精密軌道與輻射校正的 xml 參數檔
⚫
Coregistration.xml -事件前後影像套疊 xml 參數檔
⚫
runGRASSGIS.bat - GRASS-GIS 工作的啟動程式
⚫
GRASS_JOB - GRASS-GIS 計算 height above the nearest drainage（HAND）的內容
runGPT*.BAT 內需要設定一些參數，包括影像檔案位置（master*
and slave*） 、裁切範圍（subregion） 、前處理 xml 參數檔與 GRASS-GIS 程式所在目錄（graph） 、輸出資料夾（project）與影像種類（sar_type）
的選擇。以下是每個時間段兩張 Sentinel-1 影像為輸入、興趣範圍
（W/E/S/N）為 138.576/140.994/35.62/37.368 與輸出產品名稱為 output
的範例：
設定好之後，打開 cmd，執行
>runGPT_S1_2.BAT
runGPT*.BAT 將執行以下步驟來完成淹水變異區域的標示
1.
當興趣區間涵蓋一張以上，兩張以下的區域，需先將共軌道影象 合併。
2.
裁切至興趣區間來降低所需處理的資料量，並減少處理時間。
3.
給予精密軌道資訊與輻射校正
set master1='..\IMAGE\S1A_IW_GRDH_1SDV_20191006T204300_20191006T204325 _029343_0355EA_E1F1.zip'
set master2='..\IMAGE\S1A_IW_GRDH_1SDV_20191006T204325_20191006T204350 _029343_0355EA_1B68.zip'
set slave1='..\IMAGE\S1B_IW_GRDH_1SDV_20191012T204223_20191012T204248_
018447_022C0C_CBFE.zip'
set slave2='..\IMAGE\S1B_IW_GRDH_1SDV_20191012T204248_20191012T204313_
018447_022C0C_53C5.zip'
set subregion='POLYGON ((138.576 35.62,140.994 35.62,140.994 37.368,138.576 37.368,138.576 35.62))'
set graph=graph set project=output set sar_type=1 if %sar_type%==1 ( set sar_type=sentinel-1 set pol=VV
set graph=graph set project=output set sar_type=1 if %sar_type%==1 ( set sar_type=sentinel-1 set pol=VV