行政院國家科學委員會補助專題研究計畫
□成果報告
期中進度報告
雷達干涉處理之大氣效應修正研究(1/2)
計畫類別:
個別型計畫 □ 整合型計畫
計畫編號:NSC 98-2116-M -151 -001-MY2
執行期間:98 年 01 月 01 日至 99 年 07 月 31 日
計畫主持人:謝嘉聲
共同主持人:
計畫參與人員:朱良治、曾昱峰、王啟增
成果報告類型(依經費核定清單規定繳交):精簡報告 □完整報告
本成果報告包括以下應繳交之附件:
□赴國外出差或研習心得報告一份
□赴大陸地區出差或研習心得報告一份
□出席國際學術會議心得報告及發表之論文各一份
□國際合作研究計畫國外研究報告書一份
處理方式:除產學合作研究計畫、提升產業技術及人才培育研究計畫、
列管計畫及下列情形者外,得立即公開查詢
□涉及專利或其他智慧財產權,□一年□二年後可公開查詢
執行單位:
中 華 民 國 98 年 05 月 31 日
摘要
雷達波在大氣層中受水蒸汽影響而產生延遲誤差,因此雷達干涉量測技術產生之干涉 圖同樣受水汽的影響而有大氣效應雜訊,影響偵測地表變形之精度。本研究目前已完成 GPS 資料的蒐集及計算,配合屏東地區雷達干涉圖及變形圖的處理,擬利用 GPS 接收儀中天頂 訊號的延遲量來估計雷達干涉圖中同樣位置的相位延遲,並以二維的多項式模式進行內插 計算,期望能有效的干涉圖中移除大氣效應。 關鍵詞:雷達干涉量測技術、變形量、大氣效應Abstract
Radio wave propagation delay introduced by the water vapor content in the Earth's atmosphere is the major error source in Interferometric Synthetic Aperture Rarar (InSAR) technique. In this project we used the zenith wet delay data obtained from GPS receivers to have estimated the atmosphere phase delay at corresponding locations in the radar interferogram. Using low order 2-D polynomial model to represent the phase delay in the whole image, we managed to largely remove the low frequency atmosphere noise from the interferogram.
應用 GPS 資料移除雷達干涉圖中大氣效應影響量研究
Atmosphere Delay Measurement using GPS data and Removal of Atmosphere Effect from InSAR Images
一、
簡介
雷達干涉技術係利用不同時間或不同位置的雷達天線所獲得的兩幅或多幅複數影 像,可以產生地表 DEM [1,2],可以精確的量測地表的位移[3]及冰川、洋流的移動等[4]。在 相關的應用中,皆藉由相位值的差異量來獲取地表的三維資訊。 在雷達干涉處理模式中,假設雷達波的行進是不受影響的,引起相位差的來源只有距 離差值。然而,在雷達的成像過程中,雷達波從衛星天線發射出來至與地面交互作用後返 回接收器之間,訊號的路徑需經過大氣層,接收的雷達訊號一定會受大氣層的影響,所以 探討高精度的地表變形研究,就要考慮大氣效應的影響[5,6]。在重複軌道干涉量測地表變 形的應用中,兩幅雷達影像必須在不同的時間觀測,大氣條件的差異會引起不同的相位延 遲,對干涉結果產生影響,因此要降低大氣效應引起的相位差,才能獲得較精確的量測結 果[7]。 本研究的目的在了解及分析大氣層中水蒸汽的影響量,因為水蒸汽的效應,使得量測 的雷達波產生路徑延遲,進而影響相位干涉圖。因此,本研究擬藉由同為電磁波的量測工 具 GPS 及相關的氣象資料,以同一時間、同一點位有相同影響量的特性,來計算干涉圖中 的大氣延遲量,並修正該延遲量,以產生較精確的量測結果。二、
大氣效應影響
在大氣層中主要造成訊號延遲誤差,一般可分為對流層折射誤差與電離層折射誤差兩 部分。此折射誤差造成的主要原因是由於大氣並非由單一均勻介質所構成,因此電磁波在 傳遞時便會受到大氣折射影響而造成延遲,使得傳遞速度以小於真空光速的群速度在傳 遞,另外傳遞路徑也將產生彎曲。上述兩種情形即稱為電磁波之大氣折射效應,此影響將 使所量測距離增大,所造成的距離變化量便稱為遲延量。因此,為了得到觀測衛星至測站 間真實的直線距離,就必須考慮大氣折射修正。 電離層的延遲效應所產生的路徑增長可以用下式估計Idelay=40.3×TEC/f2 (meters) (1)
其中 TEC(Total Electronic Content) 為立方米單位中總電子含量,f 為無線電磁波的頻 率。 電離層為色散介質(Dispersive Medium),所以延遲量與波長相關。當 GPS 接收儀內配 置有雙頻機收儀時,就可以量測電離層誤差量[8]。因為電離層的延遲與空間的關係沒有相 關性[7],對干涉圖的影響僅是一個常數的相位延遲量,甚至在不同的時間段,其影響量變 化都不明顯,所以電離層延遲對干涉圖結果的影響就不明顯,也就不是本研究探討的重點。 事實上,雷達干涉圖本身對於絕對相位値就不敏感,常態相位値通常可以利用固定點或 GPS 的量測等附加資訊來移除。
對 於 頻 率 小 於 50GHz 訊 號傳 遞而 言, 對流 層 乃屬 於非 色散 介質 (non-Dispersive Medium),即電磁波訊號在對流層傳遞速度與訊號頻率是無關的。故此層所造成之遲延量, 是無法以雙頻 L1、L2 之線性組合來消除,所以在探討此層遲延量時,大多是以空氣折射 率與對流層折射度兩者來著手。由於雷達波訊號在對流層的折射度,與訊號傳遞路徑上之 空氣折射率是有關的,而空氣折射率則與大氣溫度、溼度及壓力有著一定函數關係存在。 且當大氣因時間、季節及氣候環境而改變時,空氣折射率也會隨之而產生變化。 延遲乃是沿著傳遞路徑上之折射率(Refractive Index)所引起的,不同介質折射率會隨著 行進路線上不同溫度與壓力而改變,因此可將折射率視為傳遞路徑的函數。進而才可利用 訊號通過中性大氣分子以及訊號在真空傳遞時,兩者到達接收器時間的不同來做比較,或 是以訊號偏折後所進行的傳遞路徑與原直線傳遞距離之差值,來直接表示成傳遞路徑延遲。 對流層延遲可以用雷達波傳遞的延長路徑長來描述,公式如下
∫
− = 2 1 ( 1) l l delay n dl T (meter) (2) 其中 n 為折射率,該値在傳播路徑中可能不是常數。 對流層的折射率與大氣壓力、溫度及溼度關係密切,但因該層的對流作用很強,且大 氣壓力、溫度及溼度等因素變化複雜,所以對於大氣折射率的變化及影響,目前尚未準確 的模型化。為便於分析及研究,通常將對流層折射對觀測量的影響可分為乾分量和濕分量 兩部份,表示為 N=Nd+Nw (3) 其中,乾分量與大氣壓力和溫度有關,而濕分量主要與信號傳播路徑上的大氣狀況有 關。 對於濕分量而言,主要為對流層中位置較低的水蒸汽所產生的影響,雖然濕分量的影 響僅佔對流層 10%的影響量,但是因其不可預估,且與空間分布有關,所以濕分量為干涉 圖中誤差量的主要來源。因此,我們研究大氣效應的影響就集中於大氣中水蒸汽的分布及 影響量。三、
利用 GPS 資料移除大氣效應原理
因為 GPS 的訊號也經過大氣層,同樣受大氣效應的影響,所以利用同時間、同地區的 GPS 資料,可以推求各項的大氣模式參數,利用這些參數資料,可以進行干涉處理的改正 資料。利用 GPS 資料進行雷達干涉技術大氣效應修正的流程如圖 1。 點位 斜距延遲量 影像修正量 對映函數 內插處理 點位 天頂延遲量 原始干涉圖 大氣修正 干涉圖 圖 1 以 GPS 資料修正大氣效應流程圖對流層的天頂方向延遲量可以利用固定站的 GPS 觀測量獲得,但是影響雷達干涉技術 的大氣效應則在斜距方向,因此須將天頂影響量轉換至斜距方向,其關係如圖 2。 圖 2 天頂方向延遲與斜距延遲關係圖[修改自[9]] 將天頂的延遲量轉換至雷達視距方向延遲量的方法可用對應函數的方法計算,簡單的 對應關係可以用三角函數的方式計算 θ δ δ cos z = (4) 式中,δ 為天頂方向的延遲量,δ 為斜距方向延遲量,θ 為點位之雷達視角。 z 斜向遲延量經映射後並不能完全表示為實際天頂向遲延量,因此映射函數只能將斜向 遲延量相較於天頂向而有所降低,至於是否能降低接近實際天頂向遲延量的大小,則是依 所使用映射函數之準確性而定。所以在利用 GPS 訊號求解對流層遲延量時,若考慮每一個 訊號方向的大氣遲延量時,終將使未知數比觀測量多而無法求解。因此在解算時才需映射 函數,以降低未知數的數量。此外映射函數大都為仰角的函數,一般以(4)式表示,其中乾、 濕映射函數在仰角 10~15 度以上相當類似,又 GPS 觀測低於 10 度仰角情況是非常少的, 故可將乾、濕遲延量使用同一個映射函數一起計算[10]。
四、
研究成果
4-1. GPS資料處理 本文所使用 GPS 計算軟體為瑞士伯恩大學(University of Berne)以研究為目的而開發 之 Bernese 4.2 軟體,其它轉換程式(如:天頂向遲延量轉換等)則是自行開發程式語言撰寫。 Bernese 4.2 提供了許多不同功能可供選擇,本文以 GPS 觀測資料處理過程為例,說明如 下: 1.原始觀測資料轉換:Bernese 4.2 軟體內含有多種轉換程式,可提供不同種類接收器的使 用者,將GPS 原始觀測資料(raw data)轉換成標準交換格式RINEX(Receiver Independent Exchange format)檔案。接著由RINEX檔案格式轉成Bernese 4.2 軟體所使用的二進位檔案 (包含電碼觀測及載波相位觀測檔案)。2.衛星軌道資料處理:使用衛星即時下傳的廣播星曆,或使用精密星曆建立衛星軌道資料 與衛星時鐘誤差資料,以提供任何時刻衛星之正確位置。 3.GPS 觀測資料解算: (1) 資料預先處理:首先利用電碼觀測量,個別對參考站與待測站做單點定位,目的是為 了計算出接收器之時鐘偏差量,以存入載波相位觀測檔中。而後再透過三次差分載波 相位觀測方程式偵測週波脫落,以進行週波脫落之補償。 (2) 週波模稜求定:以二次差分載波相位觀測量為平差計算之基本觀測量,利用QIF(Quasi Ionospere-Free)法[11]逐一求解每一條基線之週波模稜值(ambiguity)。 (3) 對流層參數估計:利用對流層附加參數法以最小二乘平差法求解天頂向對流層遲延之 改正量,此改正量加上先前以Saastamoinen模式計算的天頂向總遲延量,即為GPS 所 估計之大氣總遲延量。 4-2. GPS量測天頂延遲量 我們首先蒐集台灣地區 GPS 的點位資料,點位資料如表 1 所示,點位分別位於台灣地 區不同的位置,如經內插計算,可以涵蓋台灣大部分地區。 表 1 本研究搜集 GPS 點位資料 WGS84-X WGS84-Y WGS84-Z 地區 地點 KDNM -3028999.5453 5084820.7446 2369241.3629 台灣 墾丁 PKGM -2951343.6691 5049506.1173 2535725.1924 台灣 北港 TMAM -3034565.2496 5048871.0218 2437550.6415 台灣 太麻里 YMSM -3024807.9053 4921747.0277 2696033.5289 台灣 陽明山 KMNM -2761837.7491 5110347.1955 2625150.9904 台灣 金門 MZUM -2858574.9481 4964558.2545 2794722.2513 台灣 馬祖 CAOT -2975992.8218 5014387.9109 2576255.1022 台灣 南投草屯 YILN -3049510.9117 4928790.1320 2653583.3111 台灣 宜蘭大學 SINY -2997053.8388 5017005.9756 2547753.9858 台灣 南投信義 TACH -2955322.3628 5010114.0768 2607692.7107 台灣 台中梧棲 SHJU -2975958.6807 4968143.4450 2663522.4513 台灣 新竹南寮 GS10 -3015272.5731 4927897.7470 2693761.0937 台灣 桃園復興鄉 SOFN -3057733.9355 4970632.0715 2565173.5715 台灣 花蓮壽峰鄉 LANY -3095826.2720 5040448.8950 2378373.1672 台灣 蘭嶼 為方便與干涉圖比較,先計算較早成像的時間 96 年 5 月 16 日上午 10 時的大氣延遲量, 經二維內插後,所得結果如圖 3 所示。從圖中可以看出,在該時間段台灣地區有兩個大氣 水蒸汽含量較少的中心點,一個在北台灣,另一個在台灣中部山區部份,以此中心往外輻 射擴散,所以在海岸及平原的部份大氣的影響量反而較大。
圖 3 台灣地區 96 年 5 月 16 日 10 時大氣延遲量。 因此次計算的干涉圖位於屏東地區,所以就選擇位於屏東南部的點位 KDMN,並顯示 當天不同時間段的影響量如圖 4。從圖中可以看出在清晨及白天時間 0-16 時間延遲量較少, 爾後就逐步上升至 19 點為最大值,其後再逐步下降。此分布與其他站比較,並無一定規則, 完全為區域局部的特性,也就驗證水蒸汽影響量不規則的特性。 KDNM 2.6 2.65 2.7 2.75 0 2 4 6 8 10 12 14 16 18 20 22 24 (Time:H) Z W D (mm) 圖 4 屏東墾丁 96 年 5 月 16 日全天大氣延遲量。 4-3. 干涉圖的大氣延遲量 本研究先處理屏東地區干涉影像,結果如圖 5 所示,影像像對時間分別為 96 年 1 月 31 日至同年 5 月 6 日,第二對為同年 4 月 10 日至 5 月 6 日,理論上第一對的變形量涵蓋第二 對的變形量,因本圖未經變形量計算,故無法看出實際變形結果,僅能辨識相對變形量。 圖 5 中差分干涉圖之干涉環(fringe)表示視衛星方向的變形量,假若顏色變化由紅到藍表示 變化量為+2.8cm,反之,若顏色變化由藍到紅,則變化量為-2.8cm。圖左下方之圖例分別 代表像對日期、垂直及水平基線長、干涉環比例尺。
圖 5 屏東地區 96 年 11 月 12 日至 99 年 5 月 6 日之干涉圖 為了解屏東地區實際的變形結果,我們擷取變形量較大的林邊溪出口附近干涉圖進行 相位回復處理,並藉由 GPS 的水平位移量求得實際上的高程變形量。從變形圖中可以明顯 看出,前一像對的變形量遠大於後面像對,與預估結果相符。但僅從兩圖中無法了解大氣 的影響量。 圖 6 屏東地區 98 年 11 月 12 日至 99 年 5 月 6 日之變形量圖 為有效評估大氣效應的影響,進行另一時間段的干涉處理,此次分別進行 98 年 11 月 12 日至 99 年 5 月 6 日三個像對的處理,干涉圖結果如圖 7,變形量如圖 8 所示。 圖 7 屏東地區 98 年 11 月 12 日至 99 年 5 月 6 日之干涉圖
圖 8 屏東地區 98 年 11 月 12 日至 99 年 5 月 6 日之變形量圖 從處理結果中可以看到,依變形的時間判斷,981112-990121 及 990112-990506 的變形 量相加應等於 981112-990121 的變形量,但從圖中可以明顯的看出,990112-990506 像對的 變形量大於整體時間段的變形量,其原因可能為該時間段為屏東平原的乾季,超抽地下水 嚴重,導致地層下陷量較大外,另有可能原因即為大氣效應影響,如無大氣效應影響,前 兩時間段的變形量應等於第三像對的變形量,但結果明顯顯示並非如此,證明在處理的時 間內,大氣效應的確影響結果。
五、
結論
大氣效應延遲是雷達干涉圖進行地表偵測時主要的誤差來源之ㄧ,從研究過程中可以明 顯發現,干涉圖中有大氣效應的影響,因此移除大氣效應的影像量,才能獲得精確的地表 變形量。干涉圖偵測地表變形量測的結果,可藉由 GPS 資料的處理來進行大氣效應的延遲 量模擬,進而有效的移除大氣效應影響。從大氣影響的原理可知 GPS 點位的資料與雷達干 涉圖同一點位的天頂(ZWD)影響量相當類似,因此只要進行二維的內插處理,即可獲得干 涉圖全區的大氣影響量,藉以移除干涉圖的大氣效應。 本研究目前已完成屏東地區雷達干涉圖及變形量的計算,並蒐集部分時間段的 GPS 資 料及計算延遲量。為確實掌握大氣效應的實際影響層面,應有更密集的點位資料及其他的 氣象資料,配合更多的干涉圖及計算量,進一步的了解在台灣地區大氣效應對干涉圖的影 響。六、
參考文獻
[1] H. A. Zebker, C. L. Werner, P. A. Rosen, S. Hensley, Accuracy of Topographic Maps Derived from ERS-1 Interferometric Radar, IEEE Transactions on Geoscience and Remote Sensing, vol.32, no.4, p. 823-36, July 1994.
[2] H. A. Zebker, R. P. Salazar, T. H. Dixon, Mapping the World's Topography Using Radar Interferometry: The TOPSAT Mission, Proceedings of the IEEE, Vol.82, No.12, p. 1774-86, Dec. 1994.
[3] H. A. Zebker, P. A. Rosen, R. M. Goldstein, A. Gabriel, C. L. Werner, On the Derivation of Coseismic Displacement Fields Using Differential Radar Interferometry: The Landers Earthquake, J. of Geophysical Research, vol.99, no.B10, p. 19617-34, Oct. 1994.
[4] R. M. Goldstein, T. P. Barnett, H. A. Zebker, Remote Sensing of Ocean Currents, Science, vol.246, pp. 1282-1285, Dec. 1989.
[5] H. A. Zebker, P. A. Rosen, Scott Hensley, Atmospheric Effects in Interferometric Synthetic Aperture Radar Surface Deformation and Topographic Maps, J. Geophysical Research, Vol. 102, No. B4, pp. 7547-7563, April 1997.
[6] J. L. Davis, T. A. Herring, I. I. Shapiro, A. E. E. Rogers, G. Elgered, Geodesy by Radio Interferometry: Effects of Atmospheric Modeling Errors on Estimates of Baseline Length,
Radio Science, Vol. 20, No. 7, pp. 1593-1607, Nov.- Dec. 1985.
[7] R. F. Hanssen, T. M. Weckwerth, H. A. Zebker, R. Klees, High-Resolution Water Vapor Mapping from Interferometric Radar Measurements, Science, Vol. 283, pp. 1297-1299, Feb. 1999.
[8] B. Hofmann-Wellenhof, H. Lichtenegger, J. Collins, GPS Theory and Practice, Springer-Verlag Wien, 1994.
[9] R. Hanssen, Radar Interferometry : Data Interpretation and Error Analysis, Kluwer Academic Publishers, Boston, 2001.
[10] D. M. Tralli, S. M. Lichten, Stochastic Estimation of Troposphereic Path Delay in Global Postitioning System Geodetic Measurements, Bull. Geod., Vol.64, pp. 127-159, 1990.
[11] M. Rothacher, L. Mervart, G. Beutler, E. Brockmann, S. Fankhauser, W. Gurtner, J. Johnson, S. Schaer, T. Springer, R. Weber, The Bernese GPS Software Version 4.0, Astronomical Institute, University of Berne, 1996
七、
計畫成果自評
本計畫執行第一年將近結束(98 年 1 月 1 日至 98 年 7 月 31 日),研究內容與原計畫大致 相符,預期目標除部份項目待加強執行外,其餘已達成,執行進度請參閱附表一。目前已 整理包含台南地區及屏東地區可進行變形量計算的雷達影像像對資料,並完成計算屏東地 區地表變形干涉結果;同時並蒐集研究區可用之 GPS 資料,並完成 GPS 資料大氣延遲量的 計算。 本年度內執行期間至 98 年 7 月 31 日止,期間內仍將繼續進行地區氣象資料及 MODIS 資料的搜集及處理,以進行大氣效應影響的評估及比較。研究的成果可對台灣地區大氣效 應對雷達干涉的影響量有明確的評估及了解,進而探討合適的處理方式以移除大氣效應雜 訊,產生更精確的地表變形量。附表一、 預計完成進度 時間 完成 (1).整理可應用於計算變形研究之ERS衛星雷達影像。 第一季 已完成 (2).計算出初步的地表變形干涉成果。 第一季 已完成 (3).整理研究區之GPS資料、氣象資料等 第二季 部份完成 (4).計算模擬GPS、地區氣象資料及MODIS資料產生的大 氣影響量 第二季 部份完成 進度 (5).產生變形量的干涉圖,可以得知影像範圍內全區含大 氣效應的地表變化情形。 第二季 部份完成