基于地形的大氣水汽插值方法比較

俞曉瑩，許文斌，楊亞夫

（1.中南大學 地球科學與信息物理學院，湖南 長沙，410083；2.長沙理工大學 土木與建筑學院，湖南 長沙，410076）

衛星合成孔徑雷達干涉（Synthetic aperture radarnterferometry，INSAR）是目前正在發展中的極具潛力的微波遙感新技術[1]。它對于大范圍測圖和地表形變監測有著極其重要的應用價值。然而，由于雷達成像過程中信號傳播透過大氣層時會產生相位延遲，因此，其監測精度受大氣尤其是大氣中水汽的影響非常大[2]。Massonnet等[3]最早于1994年識別出InSAR中的大氣影響。1997年Zebker等發現當大氣相對濕度時空變化達20%時，將引起10～14 cm的形變測量誤差或80～290 m的高程誤差[4]。因此，減少干涉圖中的大氣影響對提高 INSAR觀測精度有著極其重要的意義。目前，大氣改正方法主要分為2個方面：（1）堆棧法；（2）外部數據校正法。堆棧法主要指利用多個干涉對或者相干圖分析INSAR中的大氣影響，主要包括干涉圖平均法[4]、相關分析法[5?6]、干涉圖比較法[7]和永久散射體法[8?10]。干涉圖平均法和永久散射體法對于非線性形變較明顯的區域的結果將會產生較大影響。干涉對比較法只能比較大氣影響，無法準確計算出大氣影響的大小。外部數據校正法主要有利用地面氣象數據[11]、GPS數據[12]、高分辨率的氣象模型[13]、衛星數據、MODIS水汽數據[14?15]和 MERIS水汽數據[16]等。在利用外部數據改正大氣影響時，由于可獲得的外部數據的分辨率遠低于ASAR數據的分辨率，需要采用一定的插值手段實現有效的空間插值，才有可能去除影響形變的大氣因素。Emardson等[17]提出了利用連接GPS監測網獲取的可降水汽量（PWV）來計算天頂濕延遲，并利用計算出來的天頂濕延遲校正INSAR圖上的大氣影響。該方法計算出來的監測點的天頂濕延遲精度很高，但由于GPS空間分辨率很低，因此，利用空間插值仍無法精確計算整個區域大氣影響的情況。MERIS和ASAR是搭載在同一顆衛星上的覆蓋同一地區的近紅外產品，這為獲取與ASAR數據實時的水汽數據提供了很好的條件，同時，MERIS空間分辨率（300 m×300 m）比GPS高很多，且利用無云天氣情況下的MERIS水汽產品可以大大提高INSAR圖像上的大氣影像監測精度。在前人研究的基礎上[11?16,18?20]，本文作者提出基于區域殘差的 ONN[18]地形模型及von Karman Kriging法[20]的空間組合插值方法，并將該方法與目前主流的幾種大氣空間插值方法進行比較，以便為ASAR數據的大氣改正提供可靠的水汽數據。

1 大氣空間插值方法

空間插值方法有很多且具有代表性的有趨勢面法、反向距離插值法（IDW）和克里金插值法。由于大氣中水汽分布與地形密不可分，本研究在利用水汽與高程的地形模型——ONN模型的基礎上，提出了對普通克里金插值進行改進的插值方法——基于區域殘差的普通克里金法和基于區域殘差von Karman Kringing法，并對這3種方法的結果精度進行比較分析。

1.1 基于ONN地形模型插值方法

利用 ONN提出的地形與水汽的回歸模型，得到擬合的SAR圖像分辨率下的不同高程下的水汽值。

式中：Δl（h）為高程為h處的天頂濕延遲；Δlmin為高程最高處的天頂濕延遲；C為水平面的天頂濕延遲；α為垂直水汽剖面延遲率；其中C和α由回歸確定。

1.2 ONN地形模型+基于區域殘差的普通克里金法的空間組合插值方法

克里格法（Kriging）是Kriging在1951年首次將統計學技術運用到地礦評估而命名的。克里格法的核心是利用半變異函數模型代表空間中隨距離變化的函數，再以無偏估計與最小估計變異函數為準則，決定各采樣點的權重系數，最后，以各采樣點與已求得的權重線性組合，求空間任意點或塊的內插估計值。

其中：Z（si）為已測得的第i個位置的值；λi為在第i個位置上測得值的未知的權重；s0為預測的位置；N為已知點（已測得值的點）的數目。在 OK插值中，權重不僅建立在已知點和預測點位置間距離的基礎上，而且要依據已知點的位置和已知點的值整體的空間分布和排列。應用權重的空間排列，空間自相關必須量化。因此，運用普通克里格插值（Ordinary Kriging），權重λi取決于已知點的擬合模型、距預測位置的距離和預測點周圍的已知點間的空間關系。

確定ONN方程參數后，利用ONN模型插值，將MERIS的觀測水汽值與ONN模型估計出來的水汽值進行差分，得到殘差 r。利用 OK對區域內水汽的殘差進行估值，即為區域殘差的OK法。對估計出來的殘差值與ONN方程的估計值相加即得到ASAR的水汽估計值。但OK插值的結果存在明顯的平滑效應，即較小的值常常被高估，而較大的值往往被低估[20]。因此，本研究利用Yamamoto提出的對OK估值的處理辦法，即將區域內殘差變化較大的殘差估計值剔除。將殘差r的ERMS控制在一定范圍內進行OK插值：

其中：h0為所求點的高程值；?z（hi）為滿足一定的條件下的殘差；在這里作為輸入變量。在進行OK插值過程中，變程a反映變量的影響范圍，其設置對處理結果影響較大：a設置過大，可能會因采用了空間相關較小的點而使估值產生偽差；a過小，則有部分已知數據沒有利用。本研究利用最小二乘法按一定的半徑搜索最近的相關點，得到估計的殘差水汽值。在ONN回歸中正確的點中找到的離 h0處最近的幾個點作為輸入值，將ONN模型估計值 z*（h0）與對應位置的殘差修正值?z*（h0）相加，即得到對應位置的估計值z**（h0）：

1.3 ONN地形模型+基于區域殘差von Karman克里金法的空間組合插值方法

式中：h0為所求點的高程；a為變程，反映變量受影響的范圍。變程太大，則有過多的點被用，當數據出現聚態時，導致插值失真；變程太小，則導致學習樣本不足。? z*（hi）為在回歸正確的點中找到的離 h0處最近的幾個點作為輸入值，通過von Karman Kriging法得到調整的水汽值z*（h0）。

參照Yamamoto提出的對OK估值的處理辦法進行基于區域殘差von Karman Kriging插值。利用ONN模型插值確定ONN方程后，將MERIS的觀測水汽值與ONN模型估計的水汽值進行差分，得到殘差r。將殘差r滿足固有假定條件下的估計值作為殘差估值，將ONN估值與殘差估值相加，即得到原始ASAR水汽的估計值。根據滿足固有假定條件的殘差 r，利用最小二乘擬合來尋找最優的變程a。對殘差r利用von Karman平滑的自協方差函數模型來改進 Kriging插值，對于von Karman模型非常重要。而v的求取又沒有經驗公式可以利用。本研究采用遍歷搜索法計算v，在 v的取值范圍內，以一定步長求取每個 v對應的Kriging插值函數的ERMS，從而求得v。

2 實例分析

2.1 研究區域概況

本實驗的研究區域位于太平洋沿岸（西經118.428 333°—117.095 833°，北緯 33.413 333°—34.582 5°）。研究范圍為100 m×100 m。該區域是美國重大的油氣生產基地和地下水回灌區域。由于常年人工作業，導致該區域的地表沉降備受關注。然而，由于該區域常年受海洋氣候影響，空氣中的水汽成分較多，因此，產生的大氣延遲嚴重影響對該區域地表沉降的監測。

2.2 數據說明

本研究的數據來自歐洲空間局Category 1/AO合作項目。通過本項目獲取ENVISAR衛星上的2004年8月7日18:01:36的MERIS水汽數據（MER_RR_2P）和對應時間的ASAR數據，為開展利用MERIS數據改正大氣提供了數據。為了便于插值計算，ASAR數據處理成4∶20的多視圖像，并將MERIS得到的水汽數據投影到ASAR視線空間上。

2.3 MERIS水汽數據處理

利用DORIS精軌數據和SRTM-5°的DEM數據來糾正ASAR的軌道誤差和地形相位的去除。對應同一時間的 MERIS數據用來提取大氣可降解水汽量（Precipitable water vapour,PWV）。由PWV轉化成天頂濕延遲（ZWD）：

Π為無量綱的轉換因子；ρ為水的密度；Rv為水汽氣體常數；Tm為加權平均溫度；k1，k2和k3為大氣折射常量；w=0.062 2。根據研究區域的氣象數據和經驗值，這里的Π取為6.4。

3 大氣插值算法比較研究

3.1 各種模型的插值結果分析

3.1.1 基于地形模擬的ONN插值模型

由于MERIS數據的分辨率為1 200 m×1 200 m，而 ASAR數據經過 4∶20的多視處理之后，分辨率為80 m×80 m左右，這就需要將MERIS水汽數據插值到ASAR對應的像素空間中，實現利用MERIS水汽數據去除ASAR圖像上的大氣延遲信息，從而計算出可靠的形變量。首先利用外部SRTM數據對ASAR數據進行編碼獲得該區域地圖坐標系下各個點的高程h，然后，計算對應的高程下MERIS格網數據的大氣水汽含量IWV（Integrated water vapor），利用ONN模型回歸水汽和高程的關系，尋得 ONN模型下的水平面的天頂濕延遲C和垂直水汽剖面延遲率參數α。根據得到的參數計算出Zmin，再用ONN模型計算出所有的高程h處的水汽值。圖1所示為投影到SAR視線方向上的ONN插值的水汽圖。

圖1 投影到SAR視線方向（LOS）上的ONN插值水汽分布圖Fig.1 Water vapor maps based on ONN interpolation in LOS project

從圖1可以發現：近太平洋海域的水汽值達130 mm以上，這說明這個區域的形變特征受大氣影響很大，需要先考慮去掉大氣影響導致的相位延遲效應。ONN插值運用表征估計值可能的誤差范圍的平均絕對誤差（EMAE）、反映插值相對精確性的平均相對誤差（EMRE）、反映利用數據估值的靈敏度和極值效應的均方根誤差（ERMSE）3個參數評估插值方法效果，其結果會在很多局部地方是不太正確的。這時，需要利用局部自相關函數模型進行調整。

3.1.2 ONN地形模型+基于區域殘差的普通克里金法（OK）的空間組合插值

OK法根據變差函數的不同，相應的半變異函數模型有球狀模型、指數模型、高斯模型和冪函數模型。根據前人的研究成果發現球形模擬在對氣象要素場插值時結果最好。本研究在比較分析各種變差函數的適用性發現球狀模型的半變異函數曲線與真實的測量數據最擬合，因此，選擇球狀模型作為本實驗的半變異函數模型。

根據學習樣本數據建立實驗半變異函數，求出的距離為h的任意點對的實驗半變異函數值γ（h），圖2（a）所示為擬合出的半變異函數曲線，其方程為：

圖2 投影到SAR視線方向（LOS）上的ONN+OK空間插值的大氣分布圖Fig.2 Water vapor maps based on ONN + OK interpolation in LOS project

根據式（10）中的半變異函數，求解出任意點對的距離h的半變異函數矩陣K，再根據距離矩陣，求出權重λi，從而得到各插值點上的水汽值。

圖 2（b）所示是在 ONN插值的基礎上采用了基于區域殘差的 OK法進行局部自相關的插值處理結果圖。從圖2（b）可以看出：離太平洋較遠的區域的水汽值明顯降低，這充分體現了距海洋遠近直接影響空中水汽含量。由于大氣環流的任何改變，必然通過流場和風速場的改變而影響到水汽輸送的方向、路徑和強度。從平滑度和局部變異性2方面對照圖1和圖2，單獨用ONN法插值的水汽圖不平滑，ONN+區域殘差的 OK空間組合插值的效果具有連續性和光滑性，這與自然界中水汽分布呈流動狀比較吻合。從圖2（b）左上角發現：用 ONN+OK空間組合插值表現為較大的水汽值被低估。這是由于區域殘差的OK空間插值仍然過于考慮空間局部自相關的特點，這與 MERIS本身的水汽分布是不相符的。

3.1.3 ONN地形模型+基于區域殘差的 von Karman Kriging法的空間組合插值

在調整球狀模型的半變異函數的基礎上，利用von Karman平滑的自協方差函數模型來改進 OK插值。對比圖1、圖2和圖3，圖1中的插值效果不如圖2和圖3中的插值效果，表現為大氣分布圖不具備連續性和光滑性。ONN+von Karman Kriging克服了OK方法的局部自相關的缺點，整個圖表現為平滑而且均質（見圖3），它與大氣的實際分布特征和MEIRS插值前的水汽分布也基本吻合。

3.2 大氣空間插值方法精度比較分析

本研究采用交叉驗證法對不同空間插值方法插值效果進行檢驗，即將由ONN法、基于區域殘差的ONN地形模型及OK法的空間組合插值和基于區域殘差的ONN地形模型及von Karman Kriging法的空間組合插值進行比較分析。在最優參數模型條件下通過隨機選擇的高程大于零的10000個點所得的擬合值和MERIS水汽實際值進行比較。一般來說，運用表征估計值可能的誤差范圍的平均絕對誤差（EMAE）、反映插值相對精確性的平均相對誤差（EMRE）、反映利用數據估值的靈敏度和極值效應的均方根誤差（ERMSE）和相位標準偏差4個參數評估插值方法效果。3種插值方法交叉驗證的結果如表1所示。

式中：Zx,i是第i個點的實際觀測值；Zm,i是第i個點的觀測插值估計值；n為樣本點個數。

圖3 ONN+von Karman Kriging空間組合插值的大氣分布圖Fig.3 Water vapor maps based on ONN +von Karman Kriging interpolation in LOS

由表1可知：3種插值方法中，基于ONN+基于區域殘差的von Karman Kriging方法的4項交叉驗證結果均優于其他2種插值結果。均方根誤差能較好地反映插值結果中極值的情況，這說明基于地形的局部空間插值方法是有效的且是必要的。平均絕對誤差是所有單個觀測值與算術平均值的偏差的絕對值的平均。方法3的平均絕對誤差值最小，這說明方法3的局部平滑效果最佳，方法3的平均相對誤差和相位標準偏差均遠比方法1和方法2的低，也均比ONN插值法的低。

通過水汽的空間插值方法進行交叉驗證比較分析發現：ONN+基于區域殘差的von Karman Kriging插值法的精度遠高于ONN插值法和ONN+基于區域殘差的OK法的精度，而且von Karman Kriging插值法能有效地克服ONN+基于區域殘差的OK方法的局部插值過度平滑的問題，因此，相對于水汽的插值而言，ONN+基于區域殘差的von Karman Kriging是最優的。

表1 交叉驗證結果對比Table 1 Comparison of cross-validation

4 結論

（1）利用MERIS水汽產品分析同一時間的ASAR數據的大氣延遲效應。將基于地形模擬的 ONN水汽插值模型分別結合基于殘差分析的OK和von Karman Kriging模型進行插值水汽估計。即在利用地形插值模型進行空間插值計算的基礎上，對局部區域內殘差也進行空間插值。為了提高殘差的局部插值精度，利用對普通克里金估計值光滑效應進行后處理的方法，修正OK及von Karman Kriging的平滑效應的影響。

（2）這種結合地形和局部自相關的插值方法可以達到全局最優，同時也極大地提高了局部插值的精度，它更符合空氣中水汽分布的特征，這為后續利用ASAR干涉去除大氣延遲提供了可靠思路。

（3）目前該方法存在的主要問題是：利用Kriging法尋找合適的變異函數和參數選擇過程需要花費大量的時間，而且需要反復驗證修改。

[1]Zebker H A,Goldstein R M.Topographic mapping from interferometric synthetic aperture radar observations[J].Journal of Geophysical Research,1986,91（B5）:4993?4999.

[2]Hanssen Romon F.Radar interferometry:Data Interpretation and Error Analysis[M].Dordretch:Kluwer Academic Publishers,2002:22?24.

[3]Massonnet D,Feigl K L,Rossi M,et al.Radar interferometry mapping of deformation in the year after the Landers earthquake[J].Nature,1994,369（6477）:227?230.

[4]Zebker H A,Rosen P A,Hensley S.Atmospheric effects in interferometric synthetic aperture radar surface deformation and topographic maps[J].Journal of Geophysical Research,1997,102（B4）:7547?7563.

[5]Fruneau B,Sarti F.Detection of ground subsidence in the city of Paris using radar interferometry:isolation of deformation from atmospheric artifact using correlation[J].Geophysical Research Letters,2000,27（24）:3981?3984.

[6]Remy D,Bonvalot S,Briole P,et al.Accurate measurements of tropospheric effects in volcanic areas from SAR interferometry data:Application to Sakurajima volcano（Japan）[J].Earth and Planetary Science Letters,2003,213（3/4）:299?310.

[7]Massonnet D,Feigl K.Discrimination of geophysical phenomena in satellite radar interferograms[J].Geophysical Research letters,1995,22（12）:1537?1540.

[8]Ferretti A,Prati C,Rocca F.Nonlinear subsidence rate estimation using permanent scatters in differential SAR interferometry[J].IEEE Transactions on Geosciences and Remote Sensing,2000,38（5）:2202?2212.

[9]Ferretti A,Prati C,Rocca F.Permanent scatters in SAR interferometry[J].IEEE Transactions on Geoscience and Remote Sensing,2001,39（1）:8?20.

[10]Ferretti A,Savio G,Barzaghi R,et al.Submillimeter accuracy of InSAR time series:Experimental validation[J].IEEE Transactions on Geoscience and Remote Sensing,2007,45（5）:1142?1153.

[11]Li Z W,Ding X L,Liu G X.Modeling atmospheric effects on InSAR with meteorological and continuous GPS observations:Algorithms and some test results[J].Journal of Atmospheric and Solar-Terrestrial Physics,2004,66（11）:907?917.

[12]宋小剛,李德仁,廖明生.基于GPS觀測量的INSAR干涉圖中對流層改正方法及其論證[J].武漢大學學報:信息科學版,2008,33（3）:233?236.SONG Xiao-gang,LI De-ren,LIAO Ming-sheng,et al.A method to correct tropospheric delay in SAR interferometry from GPS observations[J].Geomatics and Information Science of Wuhan University,2008,33（3）:233?236

[13]Webley P W,Wadge G,James I N.Determining radio wave delay by non-hydrostatic atmospheric modeling of water vapour over mountains[J].Physics and Chemistry of the Earth,2004,29（2/3）:139?148.

[14]Li Z H,Muller J P,Cross P,et al.Interferometric synthetic aperture radar（InSAR）atmospheric correction:GPS,moderate resolution imaging spectroradiometer（MODIS）,and InSAR integration[J].Journal of Geophysical Research,2005,110（3）:B03410.

[15]Li Z H,Fielding E J,Cross P,et al.Interferometric synthetic aperture radar atmospheric correction:Medium resolution imaging spectrometer and advanced synthetic aperture radar integration[J].Geophysical Research Letters,2006,33:L06816.

[16]許文斌,李志偉,丁曉利,等.利用 MERIS水汽數據改正ASAR干涉圖中的大氣影響[J].地球物理學報,2010,53（5）:1073?1084.XU Wen-bin,LI Zhi-wei,DING Xiao-li.Correcting atmospheric effects in ASAR interferogram with MERIS integrated water vapor data[J]Chinese Journal of Geophysics,2010,53（5）:1073?1084.

[17]Emardson T R,Simons M,Webb F H.Neutral atmospheric delay in interferometric synthetic aperture radar applications:Statistical description and mitigation[J].Journal of Geophysical Research,2003,108（B5）:2231.

[18]Onn F,Zebker H A.Correction for interferometric synthetic aperture radar atmospheric phase artifacts using time series of zenith wet delay observations from a GPS network[J].J Geophys,2006,11:B09102.

[19]Rolf Sidler.Kriging and conditional geostatistical simulation based on scale-invariant covariance models[D].Zurich:Swiss Federal Institute of Technology Zurich,2003:12?16.

[20]Yamamoto J K.Correcting the smooth effect of ordinary Kriging estimates[J].Mathematical Geology,2005,37（11）:69?94.