地面沉降災害是一個全球性問題,僅在美國,已經有遍及45個州的土地受到地面沉降的影響.傳統基於水準測量等手段的地面沉降監測技術需要大量的外業工作,強度大、作業周期長、費用高,覆蓋範圍和採樣密度均受到很大的限制,難以適用於大範圍的地表形變調查和監測.進入21世紀以來,隨著雷達衛星觀測能力的增強和合成孔徑雷達干涉測量技術(InSAR)的發展,星載InSAR技術已被廣泛應用於城市地表沉降監測(王艷等,2007;張勤等,2009;Huetal.,2018).
本文研究區域美國南加州地區地殼運動活躍,地表形變模式複雜多樣,局部地區地表沉降、抬升和季節性形變等現象並存.眾多學者利用InSAR技術對該地區地表形變開展了大量研究工作,Ferretti等(2000)將永久散射體InSAR技術(PS-InSAR)用於監測Pomona-Ontario盆地形變;Bawden等(2001)監測到SantaAna盆地季節性形變;Zhang等(2012)提出臨時相干點InSAR技術(TCP-InSAR),利用1995年10月至2000年12月的32景ERS-1/2數據得到的監測結果與GPS結果較為一致;Hu等(2013)提出加權最小二乘方法,利用18景ENVISATASAR數據得到2003年至2007年間洛杉磯地區的地表形變.上述應用均採用單一平台衛星數據,但受其重訪周期和一維形變測量能力的限制,很難反映真實的地表形變特徵及其演化規律.隨著越來越多的SAR衛星系統發射運行,使得相同時間段內覆蓋同一區域的多源多軌道InSAR數據的融合成為可能,不僅可以提高形變時間序列的連續性,還可以重建三維形變場(Wangetal.,2018).然而,要融合多源InSAR數據,除了要建立不同平台觀測數據準確的函數模型之外,還需要確定觀測數據的隨機模型,進而根據觀測數據質量確定各數據源對應的權重.但是,由於時空失相干、地形誤差、軌道誤差、解纏誤差、大氣延遲誤差等因素的影響,導致不同平台通常具有不同誤差統計特性,所以很難甚至不可能準確獲取先驗方差.鄧琳等(2016)建立約束模型並結合多源數據提取南加州地區長時間序列形變結果,但該方法沒有考慮到不同平台之間誤差分布的差異性,而且忽略水平方向形變,結果僅提取出了垂直方向的形變.
本文針對時間和空間上均重疊的多源多軌道序列InSAR數據,根據先驗形變趨勢信息改進小基線集技術(SBAS)形變信號數學模型,並採用方差分量估計(VarianceComponentEstimation,VCE)方法(張勤等,2011)對多源InSAR觀測值進行分組,通過解算出來的單位權中誤差對觀測值的方差(或權重)進行疊代估計更新,最終得到形變參數的最優估值.實驗結果表明,本文方法得到的自適應融合形變測量結果在垂直方向的形變時間序列和水平方向的形變速率與位於研究區域內的9個GPS觀測站點的監測資料均較為一致.
1、融合多源InSAR數據獲取時間序列和三維形變場
1.1 改進的時間序列形變信號數學模型
首先回顧單一平台獲取時間序列形變信號的方法,並提出改進的形變信號數學模型.
假設研究區域在時間序列t0到tN上以重複軌道觀測方式獲取了N+1幅SAR影像,按照SBAS方法可組合生成M幅差分干涉圖,所有干涉圖均已經過解纏處理,對於某一像素,其解纏干涉相位δφT=[δφ1,…,δφM]與對應的SAR影像相位φT=[φ1,…,φN](定義t0時刻的相位φ0=0)之間的關係可以表示為(Berardinoetal.,2002):
其中,A為M×N維繫數矩陣,每行對應一幅干涉圖,其主影像和輔影像對應係數分別為1和-1,其他均為0.如果M≥N,A的秩等於N,由最小二乘原則:φ=(ATA)-1ATδφ.
為了得到一個物理意義上平滑連續的解,將公式(1)中的未知參數用平均形變相位速率來替代(Berardinoetal.,2002),即:
公式(1)經過變換可以表示為:
其中,B是M×N維繫數矩陣,其位置(j,k)的數值可以表示為B(j,k)=tk+1-tk,ISj+1≤k≤IEj,j=1,…,M,IE表示主影像,IS表示輔影像;否則B(j,k)=0.通常,採用奇異值分解(SingularValueDecomposition,SVD)方法計算B的偽逆矩陣(Berardinoetal.,2002).
公式(1)還可以表示為
其中,d(n)表示第n次觀測對應的視線方向累計形變量,n=1,…,N,R為傳感器到目標的斜距,θ為入射角,Bk⊥為垂直基線,Δh為高程誤差,Δφk為噪聲相位.
該式與公式(3)聯合可以寫成:
上述形變速率為視線向(Line-of-Sight,LOS)的一維形變反演結果,通過分析南加州的GPS監測資料可發現這些站點的位移在垂直向上波動較為明顯,且具有一定周期性,而在水平方向(東西向和南北向)上僅呈現線性趨勢.為此建立相應的形變信號數學模型,同時引入東西向和南北向形變速率參數,以獲取三維形變測量結果.假設[SeSnSu]為LOS向在東西向、南北向和垂直向上的投影矢量,則公式(5)可以表示成:
其中,v″是維數為(N+3)×1的未知參數矢量,viU,vE,vN分別為以相位記錄的垂直向、東西向和南北向形變速率,i=1,…,N,D是維數為M×(N+3)的係數矩陣,δφ是維數為M×1的觀測值矢量,tmb為第m幅干涉圖的時間基線,m=1,…,M.
根據上述模型,假設觀測值的方差(權重P)已知,採用最小二乘估計可得:
直接利用公式(13)進行單一衛星平台數據求解三維形變參數時會導致嚴重的秩虧問題,因此需要融合多源InSAR數據進行多維時序形變估計.
1.2 驗後方差分量估計
要融合多源InSAR數據,不僅需要建立不同平台觀測數據的函數模型,還需要確定觀測數據的隨機模型.然而如前所述,InSAR觀測值的方差是難以精確獲取的.因此本文利用方差分量估計對觀測值的方差或權陣進行驗後估計.方差分量估計的基本理論是:首先對各觀測量初始定權,並採用最小二乘方法進行預平差;然後依據一定原則,利用平差得到的觀測值改正數來疊代估計觀測量的方差,直至各觀測量的單位權中誤差相等(張勤等,2011).多類觀測值要根據它們的統計特性分成不同的組,由於不同的InSAR平台的觀測值相互獨立,我們將一類InSAR數據觀測量分為一組.本文採用的實驗數據為三種InSAR平台獲取的觀測數據,所以將觀測值分為三組.設每一組為Li,其權重為Pi,其中i=1,2,3.首先利用最小二乘方法得到第一次估值:
就是最後方差分量估計的各組觀測值的權陣,將其代入公式(14)即可得到未知參數最優估值.
1.3 數據處理流程
由於不同的SAR衛星系統具有不同的成像幾何,導致其獲取形變測量結果的坐標系和參考基準不一致,因此需要對雷達坐標系下的解纏相位進行地理編碼,使得多源平台的觀測值位於統一的地理坐標系下.為了實現參考基準的統一,不同平台的參考點均選在同一位置.在重疊區域內,選取坐標相近的測量點作為同名點,提取其解纏相位,建立1.1節所述形變信號數學模型,按照1.2節提出的驗後方差分量估計方法進行解算融合處理.數據處理流程如圖1所示.
圖1數據處理流程圖
2、研究區域與數據
2.1 研究區域概況
美國南加州地區由於受到太平洋板塊和北美板塊的南北向擠壓,成為世界上構造運動最為活躍的區域之一.該區域地震多發,有文獻研究記錄的強震包括:1933年MW6.4的LongBeach地震、1971年MW6.7的SanFernando地震、1987年MW6.0的WhittierNarrows地震、1994年MW6.7的Northridg地震和1999年MW7.1的HectorMine地震.同時,該地區地質構造複雜,分布了許多斷層,如Newport-Inglewood斷層、SanAndreas斷層、Chino斷層、Raymond斷層等(Wasonetal.,2002).
該區域人口密集,地下水開採強度大,如Pomona-Ontari盆地(Ferrettietal.,2000)、SanBernardino盆地(LuandDanskin,2001)為典型的地下水開採區,同時礦產資源豐富,油田眾多,如Inglewood油田、SantaFeSprings油田和Wilmington油田(許文斌等,2012).由於地下水抽取與回灌和石油的開採,同時受到斷層影響,導致局部地區同時出現地表沉降、抬升以及季節周期性形變現象.位於該區域內的南加州綜合GPS觀測網(SCIGN)是世界上站點分布最密集的GPS觀測網之一,但其較低的空間採樣率(≥10km)使得無法獲取廣域範圍高空間解析度的連續地表形變場特徵信息.而InSAR技術因具有監測範圍大、測量精度高、不受天氣影響等特點,使其在地面沉降監測領域具有獨特的優勢.
2.2 數據集
本實驗研究共收集到22景L波段ALOSPALSAR升軌數據(時間跨度為2006年12月至2010年10月),47景C波段ENVISATASAR升軌數據(時間跨度2005年1月至2010年9月),以及51景ASAR降軌數據(時間跨度為2005年5月至2010年9月),三組數據集的參數如表1所示,覆蓋範圍如圖2所示.從表1和圖2中我們可以看出,三個數據集在時間和空間上均有一定的重疊度,這為提高時間序列連續性和反演三維形變場提供了數據保證.
我們首先採用StaMPS/SBAS方法對三組數據分別進行處理(時間基線和空間基線分布如圖3所示),30m解析度的SRTMDEM用於去除地形相位和平地效應,參考點均選擇位於GPS網的參考基準站ELSC站處,獲取其各自的形變速率估計和干涉圖解纏結果並進行地理編碼,然後在重疊區域內,在地理坐標系下識別同名點,並提取同名點在不同數據集得到所有干涉圖中的解纏相位值,建立形變觀測模型.針對該模型,利用方差分量估計方法得到驗後方差,進行疊代處理,直到獲取最優形變參數.位於三個數據集重疊區域內的9個GPS監測站數據從南加州綜合GPS網(SCIGN,http://www.scign.org)下載,用以和多源InSAR數據解算結果進行比對驗證.
圖2數據覆蓋範圍
黑色實線矩形框為PALSAR升軌數據；黑色虛線矩形框為ASAR升軌數據；黑色點線矩形框為ASAR降軌數據．背景圖為30m解析度SRTMDEM.
圖3時空基線分布圖
3、結果與分析
3.1 年平均形變速率
由三組時間序列InSAR數據分別解算得到各自觀測時間段內的雷達視線方向(LOS)年平均形變速率如圖4所示.
從三個數據集的解算結果可以看出,三種數據可識別到的形變區域基本一致,形變量級大致相同:在位於油氣田和城市區域的主要沉降區,由石油開採和地下水抽取和回灌引起的年平均形變量最大可達20mm;Newport-Inglewood斷層兩側表現為明顯不均勻形變(如圖4中標註所示).其中ASAR升降軌數據表現出較好的一致性,因為其具有相同的波長和相近的入射角;而PALSAR數據與ASAR數據存在較小的差異,由於ASAR數據相對於PALSAR數據波長較短,對於地表植被穿透能力較弱,因此ASAR數據對時間失相干更為敏感,導致在山區等植被茂密的地區均無法識別到足夠數量的測量點目標,另外ASAR升降軌數據集獲取時間較為接近,個別數據時間間隔僅為一天,而PALSAR數據觀測時間間隔較大,且起始自2006年,也是導致其形變速率存在微小差異的原因.
3.2 時間序列和三維形變反演結果
上述年形變速率測量結果均為雷達視線方向,我們可以根據此結果獲知形變區的位置和形變量級,但受到單一平台重訪周期等因素的限制,很難得到時間解析度較高的時間序列結果,而且僅僅從一維形變也很難確定該研究區域的形變特徵,針對以上難點,我們採用驗後方差分量估計的方法將三個序列的InSAR形變測量結果進行融合.為了驗證本文驗後方差分量估計方法的優勢,採用等權估計方法作為對照組.等權估計方法不考慮多源觀測數據之間的不同誤差統計特性,對各類觀測值賦以相同的權重,解算過程不需要疊代.方差分量估計方法和等權估計方法得到的垂直向上InSAR時間序列形變解算結果與位於重疊區域內的9個GPS監測站觀測的時間序列形變進行對比,結果如圖5所示.
從圖5中可以看出,等權估計方法明顯受到不同平台數據的誤差影響,解算的形變時間序列結果波動較大,且與方差分量估計方法存在偏差,說明本文採用的方差分量估計方法可以有效自適應調節不同平台觀測數據的權重,經過疊代處理,最終獲取形變參數的最優估值.多源InSAR數據融合後,時間節點增加至120個,提高了形變時間序列的連續性,與GPS監測數據對比來看,較為準確地反映了地表形變波動特徵,且周期性也比較一致,說明本文採用的形變信號數學模型對於重建垂直向形變時間序列具有較好效果.
圖6a—b為東西向和南北向形變速率結果,為了便於比較,GPS水平方向上監測數據經過擬合得到的形變速率在圖中用帶有顏色的三角形表示,從多源InSAR數據解算結果以及GPS監測結果可以看出,該區域存在水平形變,而且不可忽略,東西向主要表現為向西形變(向東為正值,向西為負值),年平均速率為20～40mm·a-1,南北向主要表現為向北形變(向北為正值,向南為負值),年平均速率為10～30mm·a-1,但在SantaFeSprings油田區域,存在方向相反的水平形變,這與該油田的開採導致形成漏斗狀的形變有關(圖6中剖線A-A′提取形變如圖7a—b所示).另外在Newport-Inglewood斷層兩側,也存在明顯不均勻水平形變(圖6中B-B′剖線提取形變如圖7c—d所示).水平方向上形變速率方向和量級分布如圖8所示,可以看到該區域水平形變以西北方向為主,說明該區域受到太平洋板塊和北美洲板塊的相互作用,產生了向西北方向形變的趨勢.
圖4三個數據集分別解算得到的年平均形變速率
多源InSAR數據融合解算結果與GPS監測資料在水平方向上的差異統計結果列於表2,東西向和南北向的均方根誤差分別為7.567mm·a-1和10.294mm·a-1,可以看到兩個方向的形變測量精度大致相當.雖然InSAR對南北向形變不敏感,但由於本文中使用的數據包含升降軌數據,而且水平形變具有明顯線性特徵,因此建立的形變信號數學模型對南北向形變也起到約束作用,有助於改善南北向的形變測量精度.
4、討論
由於受到SAR衛星成像幾何的影響,InSAR測量對南北方向形變不敏感,而要得到三維形變結果,至少需要包括升降軌在內的三個獨立的InSAR數據集,然而目前普遍採用的多源InSAR數據融合方法均為針對大尺度形變監測設計,不適用於緩慢形變情況,或者忽略南北向形變甚至水平形變,只解算垂直向上的形變.因此在建立準確的函數模型和隨機模型的基礎上發展更為有效的多源InSAR數據自適應融合方法是本文的研究重點.
圖5方差分量估計方法(藍色虛線)與等權估計方法(綠色虛線)解算垂直向形變時間序列與GPS監測數據(灰色實線)對比結果.各子圖左下角給出了對應GPS觀測站點名稱
圖6水平方向年形變速率結果
首先我們改進小基線集技術時序InSAR分析採用的形變反演模型,加入東西向和南北向的水平形變速率參數,從而得到垂直向上形變時間序列結果和水平向上形變速率結果,雖然該形變信號數學模型不具有普遍性,但可以利用研究區域的先驗信息,發掘研究目標的形變特徵,建立模型時做適當取捨,這樣就會避免忽略有用形變信息,造成誤判.另外,多源數據融合要考慮到不同平台數據之間誤差分布的差異性,採用方差分量估計方法將多源InSAR數據進行分組,對驗後方差進行疊代計算,用以精確確定權函數,最終得到形變參數的最優估值,實現多源InSAR數據融合.
5、總結
本文獲取覆蓋南加州地區的22景升軌ALOSPALSAR數據,47景升軌和51景降軌ENVISATASAR數據,結合南加州綜合GPS網觀測數據,發現該區域在垂直向上呈現明顯形變波動,且帶有一定周期性,而在東西向和南北向則表現為線性形變.因此,本文借鑑小基線集技術獲取時間序列的方法,根據GPS觀測數據提供的先驗信息,建立相應的形變信號數學模型,以得到垂直向上的形變時間序列和水平向上的形變速率結果.另外,要融合多源InSAR數據,除了要建立不同平台之間準確的函數模型,還需要確定觀測數據的隨機模型,進而用以獲取精確的權函數,但是考慮到多源InSAR數據通常具有不同的誤差統計特性,很難確定其先驗方差,因此本文採用驗後方差分量估計方法,對觀測值的方差(或權重)進行疊代估計,最終獲得形變參數最優估值,從而實現提高時間序列的連續性,重建三維形變場.與位於重疊區域的9個GPS站的監測數據對比來看:垂直向上,經過多源InSAR數據融合之後,時間節點增加至120個,較為準確地反映了地表形變波動,且周期性也比較一致,說明本文採用的形變信號數學模型對於恢復垂直向形變時間序列具有較好效果;水平方向上的形變速率結果顯示,該區域存在水平形變,而且不可忽略,東西向表現為向西形變,年平均速率為20～40mm·a-1,南北向表現為向北形變,年平均速率為10～30mm·a-1,東西向和南北向的形變速率與GPS數據擬合得到的形變速率均方根誤差分別為7.567mm·a-1和10.294mm·a-1,也證明本文使用的包含升降軌在內的多源InSAR數據以及建立的形變信號數學模型對南北向形變起到一定約束作用,有助於改善南北向的形變測量精度.
圖7沿圖6中剖線位置提取的形變速率結果
(a)—(b）沿剖線A-A′提取SantaFeSprings油田東西向和南北向的形變速率結果；（c)—(d）沿剖線B-B′提取Newport-Inglewood斷層兩側東西向和南北向的形變速率結果．虛線位置為斷層所在位置．
圖8水平形變的速率和方向統計分布