基于BFGS法融合InSAR和GPS技術監(jiān)測地表三維形變

胡 俊，李志偉＊，朱建軍，丁曉利，2，汪長城，馮光財，孫 倩

1 中南大學 地 球科學與信息物理學院，長沙 410083

2 香港理工大學土地測量與地理資訊學系，香港九龍

1 引 言

InSAR由于其覆蓋范圍大、空間分辨率高、測量精度高等優(yōu)勢，近年來已經被廣泛用于監(jiān)測地表形變[1］.但是其缺點是容易受到大氣[2］和失相關[3］的影響，時間分辨率低（幾天到幾十天，基于衛(wèi)星的重返周期），并且只能監(jiān)測地表在雷達視線方向（LOS）上的一維形變[4-7］.GPS技術是目前最常用的監(jiān)測地表形變的方法之一，通過區(qū)域布設的GPS網來連續(xù)測量，能夠提供高精度和高時間分辨率的地表三維運動信息[8］.但是GPS較高的安置和運行費用限制了GPS網的密度，因此空間分辨率很低.目前世界上布設密度最大的GPS網是美國南加州的SCIGN網和日本的GEONET網，分辨率也分別只有10～25km和25km左右.InSAR和GPS技術的特點正好形成互補，因此可以最大限度的挖掘地表形變等信息[9-12］.

針對InSAR和GPS技術的各自優(yōu)缺點，Gudmundsson等[13］在2002年首先提出將這兩種形變監(jiān)測技術進行最優(yōu)融合，將它們優(yōu)勢互補，獲取高空間和時間分辨率的地表三維形變.文中，Gudmundsson建立了InSAR和GPS獲取的地表形變量的能量方程，通過Markov隨機場理論和模擬退火法得到了冰島Reykjanes Peninsula地區(qū)的地表三維形變場.但模擬退火法計算復雜，有時甚至無法收斂.2006年，Samsonov等則對該方法進行了改進，簡化了Gibbs能量方程，并對其求偏導得到地表三維形變速率的解析解[14］.隨后，他利用該解析法獲取了美國南加州地區(qū)的地表三維形變速率場[15-16］.但該方法的缺點是計算時可能會出現(xiàn)數值不穩(wěn)定的現(xiàn)象.針對此缺點，羅海濱等[17］提出了一種直接分解法，將GPS插值所得的地表水平形變速率作為約束，直接將InSAR的LOS上的形變速率分解成垂直形變速率.該方法雖然簡單易操作，但過分依賴于GPS插值所得形變的精度，不能改善地表水平形變精度，對垂直形變精度的提高也有限.

鑒于上述問題，本文提出一種更穩(wěn)健更實用的綜合InSAR和GPS技術監(jiān)測地表三維形變的方法.該方法是基于文獻[14］中改進的綜合InSAR和GPS監(jiān)測地表形變速率的能量函數模型，并將其作為優(yōu)化反演的目標函數.首先，本文分析了解析法反演該能量函數模型存在的潛在數值不穩(wěn)定性（病態(tài)）問題，指出迭代法是解決數值不穩(wěn)定性的有效方法之一.其次，對該目標函數進行分析，從理論上證明了它是一個凸二次函數，保證了其局部最優(yōu)解即為全局最優(yōu)解.隨后，一種有效的局部最優(yōu)化迭代算法BFGS方法被用來求解地表的三維形變速率場的最優(yōu)估值.BFGS方法是一種擬Newton法，其特點是具有超線性收斂性和良好的自校正能力[18］.因此在本文研究中，該方法即能避免某些全局最優(yōu)化算法（如模擬退火法）計算復雜且難以收斂的問題，又能巧妙的避開解析法中數值不穩(wěn)定的現(xiàn)象.最后，本文通過模擬數據和南加州真實數據來驗證新方法的可行性和精度.

2 目標函數模型分析

升降軌InSAR干涉圖可以提供兩個不同LOS方向上的地表形變速率它們與地表在東西、南北和垂直三個正交方向上的形變速率[vxvyvz］T關系如下所示：

其中，α代表雷達方位角，即北方向與雷達的飛行方向的夾角（順時針），θ代表雷達入射角.

對于GPS而言，它僅提供若干個稀疏點上的地表三維形變速率，因此必須用某種插值方法（如克里金插值法[19］）提高其分辨率，得到與InSAR干涉圖相同分辨率下的地表三維形變速率：

可以看出，InSAR和GPS都提供了可用來重建地表真實三維形變的信息.根據Gibbs能量方程，可建立如下的能量函數[14］：

可以看出，（4）式中的I1和I2對應于升降軌InSAR，而Gx、Gy和Gz則對應于GPS.建立該能量函數的目的是將其作為參數優(yōu)化反演的目標函數，即當該函數的能量達到極小時，對應未知參數（[vxvyvz］T）的最優(yōu)結果.對于這個能量函數，Samsonov等提出解析法求解三維形變場，即對U（vx，vy，vz）中的未知量求偏導，并令其均為0，從而歸結為一個標準的線性方程組求逆的問題[14］：

值得注意的是，當只有升軌或降軌數據時（即I1=0或I2=0），該能量函數模型也同樣適用[15］.這種方法雖然簡單，但是在實際應用中，可能會帶來較大的誤差.文獻[20］中對這種直接矩陣逆法進行了誤差分析，得到其相對誤差限為

其中，kcond為矩陣A的條件數，‖Δx‖為解x的誤差，‖ΔA‖為系數矩陣A的誤差，‖Δb‖ 為常數項b的誤差.（6）式表明，x的相對誤差受到系數矩陣相對誤差、常數項相對誤差和條件數的影響.在系數矩陣誤差很小的情況下，x的相對誤差為b和A的相對誤差的kcond倍.因此，如果系數矩陣A是病態(tài)矩陣，其條件數kcond很大，那么即使‖Δb‖很小，也有可能會給解x帶來很大的影響.而在解析法中，系數矩陣完全有可能是病態(tài)的，這一點在后面數據驗證時會得證.因此，解析法存在潛在數值不穩(wěn)定性問題，在病態(tài)時將難以獲得準確的解.在這種情況下，需要尋找一種不受系數矩陣病態(tài)性影響的最優(yōu)化算法.

迭代法是一種常用的最優(yōu)化算法，其特點是通過循環(huán)計算來無限逼近問題的近似解，不直接對系數矩陣求逆，所以不受其病態(tài)性的影響[20］.因此，在本研究中，考慮引入迭代算法來解決病態(tài)系數矩陣帶來的數值不穩(wěn)定問題.首先對目標能量函數的特性進行進一步分析.分別給出（4）式相應的Hesse陣，以I1為例：

可求得2I1的特征值為

同理可得

可以看出2I1和2I2的特征值均≥0，因此它們在（-∞，+∞）上處處半正定.此外

也為 （- ∞，+ ∞）上的正定矩陣.因此，I1、I2、Gx、Gy和Gz均為（-∞，+∞）上的凸函數.由凸函數的性質可知，目標能量函數U（vx，vy，vz）也必為在（-∞，+∞）上的凸函數，其局部最優(yōu)解即為全局最優(yōu)解[21］.因此，無需設計復雜的全局最優(yōu)化迭代算法，較為簡單的局部最優(yōu)化迭代算法就可以達到形變速率[vxvyvz］T的全局最優(yōu)估值.

3 BFGS方法

本研究中，如何設計一種簡單有效且適用于前面建立的目標函數的局部最優(yōu)化迭代算法是一個重點.經過理論研究和實驗比較，最終選擇了一種目前最流行的、計算性能較好的擬Newton法——BFGS方法.擬Newton法的特點是既繼承了Newton法收斂速度快的優(yōu)點，又克服了Newton法計算量大、產生非下降方向的缺點.而由Broyden、Fletcher、Goldfard和Shanno在1970年提出的BFGS方法則是公認最為有效的擬Newton法，其優(yōu)勢在于在處理實際問題時具有超線性收斂性，并且具有非常好的自校正能力[18］.

對于上述改進的綜合InSAR和GPS監(jiān)測地表形變速率的能量函數模型反演優(yōu)化問題，可以寫成如下形式：

式中x∈R3是三維地表形變速率；f：R3→R為目標函數.

BFGS方法開始于以下的迭代過程[18］：

其中，k是迭代次數；αk為步長因子，取為某個滿足Wolfe條件的正常數；dk為一個單位矢量，決定下一次搜索的方向，其可由目標函數f（x）的梯度求得：

其中：

其中，I為單位矩陣.當給定初始點x0、誤差控制限ε和H0后，利用上式進行迭代運算，當‖Δ f（xk）‖ ≤ε時，迭代停止，輸出當前的xk即為最優(yōu)的地表三維形變速率.

一般當f（x）的Hesse陣的逆矩陣存在的情況下，取H0=Δ2f（x0）-1.但是值得注意的是，在本研究中，Δ2f（x0）-1的行列式值非常之小，因此可能會給結果帶來不確定性，事實上，這就是造成解析法不穩(wěn)定的原因[14］.幸運的是，由于BFGS方法良好的自校正能力，H0可取為單位矩陣，事實證明這并不會影響最終的結果，并且可以避免數值解算的不穩(wěn)定性[22］.

4 模擬數據驗證

本研究在100×100的網格點上模擬了地表10年的三維形變場.為了盡量貼合真實的地表形變情況，設計地表在東西向和南北向發(fā)生了勻速形變，而垂直向則作隨時間推移的周期性下沉和抬升.其地表三維形變速率場分別如圖1（a，b，c）所示.其中，Vx分量和Vy分量分別在東西向和南北向上從-6cm·a-1勻速增加至6cm·a-1；Vz分量則分為下沉、抬升和穩(wěn)定三個區(qū)域，總體跨度從-8cm·a-1到8cm·a-1.

圖1 在100×100的網格點上的模擬的地表三維形變速率場（a）東西向；（b）南北向；（c）垂直向.Fig.1 The simulated three-dimensional displacements velocity fields on 100×100grid（a）East-West direction；（b）North-South direction；（c）Up-Down direction.

圖2 加噪后的InSAR監(jiān)測的地表形變速率場（a）升軌；（b）降軌.Fig.2 The surface displacements velocity inferred by InSAR with Gauss noise（a）Ascending；（b）Descending.

在100×100的模擬地表三維形變速率場的網格點上隨機選取100個GPS點，提取出每個GPS點上三維形變速率.同樣，標準差為1cm和2.5cm的高斯白噪聲分別加到GPS點的水平和垂直形變上.利用普通克里金插值法[19］，該加噪后的稀疏GPS監(jiān)測數據被內插至整個100×100的網格點上，以達到和InSAR相同的空間分辨率.圖3（a，b，c）顯示的是100個GPS點內插所得的地表三維形變速率場[VxVyVz］.可以看出，GPS點在圖像上不均勻分布，在GPS點較為稀疏的地方，普通克里金方法的插值效果較差，如圖像的右下角.此外，由于水平方向上的噪聲較小，它們的插值結果也明顯好于垂直方向.

在研究中發(fā)現(xiàn)，如果升降軌InSAR的相干性非常之好，即其相應地表形變監(jiān)測標準差σins很小，而GPS觀測或插值標準差σgps比較大時，解析法中的法方程系數矩陣就很可能是出現(xiàn)嚴重病態(tài).例如，當InSAR的相干性達到0.99，而GPS觀測或插值標準偏差大于10cm時，σins／σgps＜10-4時，這時解析法的單精度解存在很大的誤差.圖3（g—i）顯示了此時解析法的解算結果.可以看出，三個方向上的形變結果都出現(xiàn)了不同程度的偏差，這是因為此時解析法中的系數矩陣為病態(tài)矩陣，條件數接近104，因此給結果帶來了很大的誤差.而本文所提出的新方法則可以很好的抵抗病態(tài)矩陣的不良影響，仍然可以得到穩(wěn)定的高精度結果（如圖3（d—f））.

圖3 （a）—（c）普通克里金插值法得到的地表三維形變速率場；（d）—（f）BFGS法得到的地表三維形變速率場；（g）—（i）在σins／σgps＜10-4情況下解析法得到的地表三維形變速率場.從左至右分別為東西向、南北向和垂直向結果.黑色點代表隨機選取的100個GPS觀測站.Fig.3 The three-dimensional displacements velocity fields estimated by（a）-（c）ordinary kriging，（d）-（f）BFGS method and （g）- （i）analytical optimization method withσins／σgps ＜ 10-4，respectively.The left row represent the East-West vectors，the middle row the North-South vectors and the right row the Up-Down vectors.The black dots are randomly picked locations of 100GPS stations.

表1 模擬實驗中形變速率的RMSE比較（單位：cm·a-1）Table 1 Comparison of the root mean square errors of the displacement velocities derived in the simulated experiments（Units：cm·a-1）

為了定量分析不同方法的結果，本文在表1中給出了它們的均方根誤差（RMSE）.可以看出，相對于GPS插值法，直接分解法無法改善水平方向上的地表形變監(jiān)測精度.而本文中所提出的新方法則可以提高東西向上平均39%左右的形變監(jiān)測精度；垂直向上的形變監(jiān)測精度提高的最為明顯，相比于GPS插值法和直接分解法則分別平均提高了大約96%和75%，這主要是因為升降軌InSAR對于垂直向的形變最為敏感（投影單位矢量為0.935）；而在南北向上，新方法難以改善其形變監(jiān)測精度，這則是由于升降軌InSAR對于南北向上的形變最不敏感造成的（投影單位矢量為±0.095）[23］.而在本次模擬實驗中，南北向的精度還略微有點降低，這可能是由于引入了InSAR監(jiān)測噪聲的原因，在實際應用中，該誤差是可以忽略的.而解析法的精度則明顯不令人滿意，在三個方向上都出現(xiàn)了較大的偏差，其結果是不可靠的.

5 真實數據驗證

在真實數據處理中，由于地面情況非常復雜，解析法存在的潛在數值不穩(wěn)定性（病態(tài)）問題可能出現(xiàn).例如遇到GPS站點少或分布不均的情況，GPS插值標準差就會急劇增大，使得解析法不一定可以得到合適的解，而BFGS法則可以避免這種情況的發(fā)生，保證了最終的解就是高精度、穩(wěn)定的全局最優(yōu)解.為了進一步說明，選取了ENVISAT衛(wèi)星在美國南加州地區(qū)獲取的二景降軌ASAR影像（見表2）和當地52個GPS連續(xù)觀測站提供的地表形變監(jiān)測資料作為研究對象.美國南加州是一個地質活動頻繁、地表形變復雜的地區(qū)，其地表不僅受到震間應力累積的影響而做水平運動，還由于地下水抽取和灌溉引起季節(jié)性的垂直運動[24］.

表2 本文所用到的降軌干涉影像對Table 2 The descending interferometric pair used in this study

首先對二景SAR影像進行二軌差分InSAR處理.其中地形相位的影響利用美國本土30m×30m分辨率的SRTM數據進行剔除；為了消除噪音對干涉圖的影響，首先采用了距離向2視、方位向10視的多視處理，隨后利用改進的Goldstein方法進行濾波[25］；為了避免相位纏繞，采用最小費用流法進行解纏處理；大氣的影響則可利用MERIS水汽數據減弱[26］.最后對干涉圖進行相位到形變的轉換，并地理編碼至WGS 84坐標系下.如圖4a所示，整個南加州地表都在向遠離衛(wèi)星的方向運動，其中最嚴重的區(qū)域超過了4cm·a-1.圖4b顯示的則是該地區(qū)的相干圖，很明顯，除了山區(qū)出現(xiàn)了失相關，大部分地區(qū)的相干性都非常好，這得益于兩景SAR影像之間較短的空間基線和南加州干燥的氣候.通過該相干圖和視數，可以得到該形變圖的標準偏差[15］.

本文研究所用到的52個GPS臺站的監(jiān)測資料由美國南加州綜合GPS網（SCIGN）提供.SCIGN網是世界上GPS臺站密度最高的網之一，但是在其它的大部分地區(qū)，GPS站點的數量都較少且位置分布不均.因此為了使得InSAR和GPS融合監(jiān)測地表三維形變的方法具有適用性，在本次實驗中，只利用其中47個GPS臺站（位置見圖5中黑色小三角形所示）進行普通克里金插值，而將另外5個GPS站點（位置見圖4和圖5中黑色小正方形所示）留作最后的檢核.從圖中可以看出，這5個檢驗點均位于右下區(qū)域，且該區(qū)域沒有插值點，從而能夠檢驗優(yōu)化方法在GPS站點分布不均勻的情況下的性能.插值后的東西向、南北向和垂直向形變場分別如圖5（a，b，c）所示.其標準偏差可通過 GPS測量標準差和插值標準差計算得到[15］.

圖4 （a）南加州地區(qū)ASAR降軌視線方向上的地表形變速率圖；（b）南加州地區(qū)的相干圖其中底圖顯示的是SRTM數據提供的當地地形，圖中三角形和正方形代表此次研究中所用的52個GPS站的位置.Fig.4 （a）The interferometric displacements velocity map and（b）coherence map of Southern California The shaded relief map of Southern California from SRTM is used as base map.The triangles and squares represent the locations of 52GPS sites used in this study.

圖5 （a，b，c）南加州地區(qū)GPS插值所得的地表三維形變速率場；（d，e，f）BFGS法所得的地表三維形變速率場；（g，h，i）解析法所得的地表三維形變場從左至右分別為東西向、南北向和垂直向結果.黑色小三角形代表插值用的47個GPS觀測站，黑色小正方形代表驗證用的5個GPS觀測站.P和SFS分別代表Pomona和Santa Fe Springs地區(qū).Fig.5 The three-dimensional displacements velocity fields of Southern California respectively by（a，b，c）ordinary kriging，（d，e，f）BFGS and（g，h，i）Analytical method The left row represent the East-West vectors，the middle row the North-South vectors and the right row the Up-Down vectors.The 47 GPS stations with small black triangles are used for interpolation，while the 5station with small black squares are used for validation.P and SFS represent the Pomona and the Santa Fe Springs areas，respectively.

隨后，分別采用BFGS法和解析法對InSAR和GPS形變結果進行最優(yōu)融合，得到優(yōu)化后的地表三維形變速率場.在計算中，為了能夠得到更為精確的地表三維形變結果，利用了InSAR形變結果中每個像素的雷達方位角和入射角信息來獲取該像素的雷達視線方向的單位投影矢量.圖5（d，e，f）給出的是BFGS法的結果.從圖中可以看出，與GPS插值法相比，BFGS法雖然對水平向上的插值結果改進不大，但在垂直向上則有明顯的改善，許多在插值結果中捕捉不到的小區(qū)域形變被展示出來.如圖5f中的紅色圓圈所示，在Pomona地區(qū)，大約3cm的地面沉降被BFGS法探測出來；而在Santa Fe Springs地區(qū)，則發(fā)現(xiàn)了1cm左右地表抬升.這些結果與已有的研究很吻合[27］，前者是由于地下水抽取造成的，而后者則是受到油氣回灌的作用.然而，GPS插值法卻無法反映出來（圖5c）.這主要是由于ENVISAT衛(wèi)星入射角較?。s23°），使得InSAR監(jiān)測對垂直形變更加敏感.此外，缺乏升軌數據也可能是水平結果改進不大的原因之一.圖5（g，h，i）顯示的則是解析法的結果，可以看出三個方向的形變圖在右下角位置都出現(xiàn)了明顯的不連續(xù)的異常形變（如藍色矩形框所示），與已有的研究完全不符[26-27］，因此是錯誤的.這是由于在這個區(qū)域GPS插值標準差過大引起的系數矩陣病態(tài)問題，從而導致該區(qū)域形變結果不穩(wěn)定.然而，BFGS法的結果則非常穩(wěn)定（圖5（d，e，f））.

由于無法獲取其它手段得到的地表形變場來驗證最終的結果，因此本文采用BLSA等5個站點上的GPS測量結果與InSAR／GPS融合結果進行比較.如表3所示，解析法在東西、南北和垂直向上的形變都與GPS測量結果出現(xiàn)了很大程度的偏差，其結果是不可靠的，這與模擬數據得出的結論是一致的；而新方法在三個方向上的結果都與GPS測量值非常吻合，且比GPS插值法和直接分解法的精度都有了一定的提高.但是，由于被用來驗證的5個GPS點本身的形變量比較?。ú辉谛巫兠黠@的區(qū)域），而且GPS觀測值本身也存在一定誤差，特別是垂直向的結果精度較差，一般為水平向結果精度的2～3倍，因此相比其它方法，新方法的改善程度沒有模擬實驗中那么顯著.

表3 南加州地區(qū)實驗中不同方法得到的形變速率的RMSE比較（單位：cm·a-1）Table 3 Comparison of the root mean square errors of the displacements velocities derived from different methods in the Southern California experiments（Units：cm·a-1）

6 結 論

為了提高地表三維形變監(jiān)測的時空分辨率和精度，國內外許多學者將InSAR技術和GPS技術進行融合，以求優(yōu)勢互補.本文在已有研究的基礎上，對綜合InSAR和GPS監(jiān)測地表形變速率的能量函數模型進行深入分析和研究，針對目前常用的解析法存在的潛在數值不穩(wěn)定性（病態(tài)）問題，指出迭代法更適合用于模型的反演；并從數學上證明了一種局部最優(yōu)化算法就能夠得到該能量函數模型全局最優(yōu)解，因此無需引入復雜的全局最優(yōu)化算法.隨后，本文引入BFGS局部最優(yōu)算法對該目標函數模型進行優(yōu)化反演.

相比于復雜的全局最優(yōu)化算法，BFGS方法計算簡單且具有超線性收斂性，僅需要幾次迭代就可以達到全局最優(yōu)解；另一方面，BFGS方法良好的自校正能力使其在計算中可利用單位矩陣替代Hesse陣的逆矩陣，從而避免了解析法在計算過程中由于矩陣近似奇異而帶來的數值不穩(wěn)定現(xiàn)象.模擬和真實試驗結果表明，在GPS站點較少或分布不均造成GPS插值標準差過大的情況，解析法往往會帶來不穩(wěn)定性，從而大大降低了解的精度；而BFGS方法則不受GPS數量和分布的影響，能夠得到最優(yōu)的地表三維形變估值，且精度優(yōu)于插值法和直接分解法.致 謝 本次研究所用的ENVISAT ASAR數據為歐洲空間局（ESA）提供（AO-4458，4914）；GPS數據由美國南加州綜合GPS網（SCIGN）所提供（http：／／www.scign.org／），在此一并表示感謝.

（References）

[1]Massonnet D，F(xiàn)eigl K J.Radar interferometry and its application to Changes in the Earth′s surface.Reviews of Geophysics，1998，36（4）：441-500.

[2]Li Z W，Ding X L，Huang C，et al.Atmospheric effects on repeat-pass InSAR measurements over Shanghai region.Journal of Atmospheric Solar-Terrestrial Physics，2007，69（12）：1344-1356.

[3]Zebker H A，Villasenor J.Decorrelation in interferometric radar echoes.IEEE Transactions on Geoscience and Remote Sensing，1992，30（5）：950-959.

[4]Ding X L，Liu G X，Li Z W，et al.Ground Subsidence Monitoring in Hong Kong with Satellite SAR Interferometry.Photogrammetric Engineering and Remote Sensing，2004，70（10）：1151-1156.

[5]Lu Z，Masterlark T，Dzusisin D，et al.Magma supply dynamics at westdahl volcano，alaska，modeled from satellite radar interferometry.Journal of Geophysical Research，2003，108（B7）：2354，doi：10.1029／2002JB002311.

[6]王超，張紅，劉智等.基于D-InSAR的1993—1995年蘇州市地面沉降監(jiān)測.地球物理學報，2002，45（增刊I）：244-253.Wang C，Zhang H，Liu Z，et al.Subsidence of Suzhou city（China）from 1993to 1995detected by ERS differential SAR interferometry.Chinese J.Geophys.（in Chinese），2002，45（Suppl.I）：244-253.

[7]Hu J，Li Z W，Ding X L，et al.Two-dimensional Co-Seismic surface displacements field of the Chi-Chi earthquake inferred from SAR image matching.Sensors，2008，8（10）：6484-6495.

[8]Zhu J J，Santerre R，Chang X W.A bayesian method for linear，inequality-constrained adjustment and its application to GPS positioning.Journal of Geodesy，2005，78（9）：528-534.

[9]Sudhaus H，Jónsson S.Improved source modelling through combined use of InSAR and GPS under consideration of correlated data errors：application to the June 2000Kleifarvatn earthquake，Iceland.Geophysical Journal International，2009，176（2）：389-404.

[10]萬永革，沈正康，王敏等.根據GPS和InSAR數據反演2001年昆侖山口西地震同震破裂分布.地球物理學報，2008，51（4）：1074-1084.Wan Y G，Shen Z K，Wang M，et al.Coseismic slip distribution of the 2001Kunlun mountain pass west earthquake constrained using GPS and InSAR data.Chinese J.Geophys.（in Chinese），2008，51（4）：1074-1084.

[11]張勤，趙超英，丁曉利等.利用GPS與InSAR研究西安現(xiàn)今地面沉降與地裂縫時空演化特征.地球物理學報，2009，52（5）：1214-1222.Zhang Q，Zhao C Y，Ding X L，et al.Research on recent characteristics of spatio-temporal evolution and mechanism of Xi＇an land subsidence and ground fissure by using GPS and InSAR techniques.Chinese J.Geophys.（in Chinese），2009，52（5）：1214-1222.

[12]陳國滸，單新建，Wooil M等.基于InSAR、GPS形變場的長白山地區(qū)火山巖漿囊參數模擬研究.地球物理學報，2008，51（4）：1085-1092.Chen G H，Shan X J，Wooil M，et al.A modeling of the magma chamber beneath the Changbai Mountains volcanic area constrained by InSAR and GPS derived deformation.Chinese J.Geophys.（in Chinese）：2008，51（4）：1085-1092.

[13]Gudmundsson S，Sigmundsson F，Carstensen J M.Three-Dimensional surface motion maps estimated from combined interferometric synthetic aperture radar and GPS data.Journal of Geophysical Research，2002，107（B10）：2250，doi：10.1029／2001JB000283.

[14]Samsonov S V，Tiampo K F.Analytical optimization of a DInSAR and GPS dataset for derivation of three-dimensional surface motion.IEEE Geoscience and Remote Sensing Letters，2006，3（1）：107-111.

[15]Samsonov S V，Tiampo K F，Rundle J B，et al.Application of DInSAR-GPS optimization for derivation of fine-scale surface motion maps of southern California.IEEE Transcations on Geosciences and Remote Sensing，2007，45（2）：512-521.

[16]Samsonov S V，Tiampo K F，Rundle J B.Application of DInSAR-GPS optimization for derivation of three-dimensional surface motion of the southern California region along the san andreas fault.Computers &Geosciences，2008，34（5）：503-514.

[17]羅海濱，何秀鳳，劉焱雄.利用DInSAR和GPS綜合方法估計地表3維形變速率.測繪學報，2008，37（2）：960-963.Luo H B，He X F，Liu Y X.Estimation of three-dimensional surface motion velocities using integration of DInSAR and GPS.Acta Geodaetica et Cartographica Sinica （in Chinese），2008，37（2）：960-963.

[18]董云達.數值優(yōu)化引論.鄭州：黃河水利出版社，2007.Dong Y D.Introduction to the Numerical Optimization（in Chinese）.Zhengzhou：Yellow River Conservancy Press，2007.

[19]張仁鐸.空間變異理論及應用.北京：科學出版社，2005.Zhang R Z.Spatial Variability Theory and its Application（in Chinese）.Beijing：Science Press，2005.

[20]韓旭里，萬中.數值分析與實驗.北京：科學出版社，2006.Han X L，Wan Z.Numerical Analysis and Experiments（in Chinese）.Beijing：Science Press，2006.

[21]何堅勇.最優(yōu)化方法.北京：清華大學出版社，2007.He J Y.Optimization Method （in Chinese）.Beijing：Tsinghua University Press，2007.

[22]Xu P L.A hybrid global optimization method：the multidimensional case.Journal of Computational and Applied Mathematics，2003，155（2）：423-446.

[23]查顯杰，傅容珊，戴志陽.DInSAR技術對不同方位形變的敏感性研究.測繪學報，2006，35（2）：133-137.Zha X J，F(xiàn)u R S，Dai Z Y.The sensitivity of DInSAR to surface deformation in different direction.Acta Geodaetica et Cartographica Sinica （in Chinese），2006，35（2）：133-137.

[24]Bawden G W，Thatcher W，Stein R S，et al.Tectonic contraction across los angeles after removal of groundwater pumping effects.Nature，2001，412（6849）：812-815.

[25]孫倩，朱建軍，李志偉等.基于信噪比的InSAR干涉圖自適應濾波.測繪學報，2009，38（5）：437-442.Sun Q，Zhu J J，Li Z W，et al.A new adaptive INSAR interferogram filter based on SNR.Acta Geodaetica et Cartographica Sinica （in Chinese），2009，38（5）：437-442.

[26]許文斌，李志偉，丁曉利等.利用 MERIS水汽數據改正ASAR干涉圖中的大氣影響.地球物理學報，2010，53（5）：1073-1084.Xu W B，Li Z W，Ding X L，et al.Correcting atmospheric effects in ASAR Interferogram with MERIS integrated water vapor data.Chinese J.Geophys.（in Chinese），2010，53（5）：1073-1084.

[27]Lanari R，Lundgren P，Manzo M，et al.Satellite radar interferometry time series analysis of surface deformation for the metropolitan los angeles and San Francisco，California，areas.Geophysical Research Letters，2004，31（L23613），doi：10.1029／2004gl021294.