InSAR三維同震地表形變監測
——窗口優化的SM-VCE算法

劉計洪，胡 俊，李志偉，朱建軍

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

干涉合成孔徑雷達(interferometric synthetic aperture radar,InSAR)技術具有高空間分辨率，高精度等優點，已被廣泛地應用于地形測量[1-3]和地表形變監測領域[4-9]。然而，傳統差分InSAR(differential InSAR,DInSAR)技術僅可以獲得真實地表形變沿視線向(line of sight,LOS)的一維形變投影，在地質災害解譯過程中具有極大的局限性[10-11]。為了克服這一缺點，文獻[12—14]分別提出了像素偏移量跟蹤(pixel offset-tracking,POT)、多孔徑雷達(multi-aperture InSAR，MAI)技術和哨兵數據方位向子帶重疊區域干涉(burst overlap InSAR，BOI)技術來獲取地表真實形變沿雷達衛星方位向(azimuth,AZI)的投影形變。結合升降軌DInSAR和MAI/POT技術，利用加權最小二乘(weighted least square,WLS)即可解算地震、火山、冰川、滑坡等導致的真實三維地表形變場[5,15-24]。此外，全球導航衛星系統(global navigation satellite system，GNSS)觀測數據與InSAR數據在空間分辨率和測量精度方面優勢互補，二者融合也可實現高精度、高空間分辨率三維地表形變場的獲取[25-29]。

然而，上述研究均基于單個點進行三維形變解算，并未考慮臨近點三維形變之間的力學關系。文獻[30]基于彈性理論提出了一種SISTEM(simultaneous and integrated strain tensor estima-tion from geodetic and satellite deformation measure- ments)方法，首次將應力應變模型(strain model,SM)引入InSAR與GNSS數據相結合的地表三維形變求解當中。但是，現實情況往往很難獲取大范圍的GNSS觀測，文獻[31]對SISTEM方法進行了擴展，僅需InSAR觀測值即可進行高精度三維形變解算。為了方便起見，本文將DInSAR、MAI、POT和BOI等技術獲取的形變觀測數據統稱為InSAR觀測值。然而，在解算真實三維地表形變時，必須要融合不同方向的多源異質觀測數據，因此精確確定不同數據之間的權重比例具有重要意義[6,32]。為了實現此目的，科研人員研究利用一個固定的滑動窗口求解窗口內數據的標準差作為先驗信息進行定權[17]，這種定權方式是基于數據各態歷經性的假設，即假設觀測數據在某一數值上下一定范圍內是隨機分布的，但是對于地震等導致的大梯度形變而言，往往難以滿足。也有相關學者在假設遠場區域不包含形變信號的前提下，將遠場一定區域數據的均方根值作為該類觀測值整個區域的先驗方差[33-34]。這種方法忽略了InSAR觀測噪聲在空間上的差異性，因此也會影響三維地表形變的求解精度。文獻[25]在求解時序數據的平均形變速率時引入了方差分量估計(variance component estimation,VCE)進行定權，驗證了VCE在多源InSAR觀測數據定權方面的優越性。但是，對于瞬時形變場的求解，由于缺少多余觀測，應用VCE進行定權受到了極大限制。

基于此，文獻[35]提出了一種基于地表應力應變模型和方差分量估計的InSAR三維地表形變監測方法(SM-VCE)。該方法基于一定窗口范圍內不同像素點三維形變之間的力學關系(SM)建立函數模型，顯著增加了冗余觀測個數，進而可利用VCE實現不同類InSAR觀測值之間的精確定權。研究表明，基于窗口的SM-VCE方法比傳統單點解算的加權最小二乘方法得到的三維地表形變精度更高[36]。在原始SM-VCE方法中，求解不同目標點的三維形變時，窗口大小是固定不變的，并且窗口內的InSAR觀測值被認為是等精度的。對于同震形變場而言，在地震破裂帶附近往往難以獲得有效的InSAR觀測值，使得固定窗口大小的SM-VCE方法難以獲取完整的三維同震地表形變。即使在斷層附近有InSAR觀測值的情況下，窗口內往往會包含斷層兩側的異質InSAR觀測值，使得原始SM-VCE方法難以得到可靠的近場三維形變。而近場地表形變結果對于約束斷層淺部滑動又是十分重要的[23,37-38]。因此，有必要重建地震破裂帶附近三維地表形變場。同時，InSAR觀測值極易受到失相干等噪聲影響，窗口內同一類觀測值中不同像素點的形變測量值精度往往各不相同，因此，忽略窗口內InSAR觀測值精度差異的做法將會降低三維地表形變精度。

鑒于此，本文在SM-VCE方法的框架下，提出一種顧及斷層破裂帶及自適應變化窗口的觀測值選取策略，進而基于鄰近的InSAR有效觀測值即可重建地震破裂帶附近的三維地表形變場，同時，引入迭代加權最小二乘(iterative weighted least squares，IWLS)方法[39]實現窗口內同一類InSAR觀測值的內部自適應定權。本文首先利用模擬試驗驗證了本文中窗口優化SM-VCE方法的可靠性和優勢，進而基于升降軌哨兵1號衛星數據，獲取了2019年7.1級Ridgecrest地震的高精度三維地表形變場。

1 SM-VCE方法

文獻[35]較為系統地介紹了SM-VCE方法，本文在此基礎上，提出了顧及斷裂線和自適應變化窗口的觀測值選取策略和基于IWLS方法實現窗口內同一類觀測值內部的自適應定權方法，旨在基于SM-VCE方法的框架，獲取精度更高、更為完整的三維同震地表形變場。

1.1 基于SM建立觀測方程

dk=H·Δk+d0

(1)

(2)

文獻[30]將矩陣H表示為一個對稱矩陣和一個非對稱矩陣之和

H=S+R

(3)

(4)

(5)

式中，ξ和ω代表地表應力應變模型中的應變參數和旋轉參數。在原始SM-VCE方法中，H用式(3)表示，進而可將式(1)寫成

(6)

(7)

(8)

式中

(9)

綜合式(6)、式(8)，可得

(10)

式中

若假設點P0周圍有Kj個第j類InSAR觀測值，則最終可得

L=B·l

(11)

式中

此時，建立了綜合考慮InSAR觀測值與三維形變之間的幾何關系，以及地表鄰近點三維形變之間的力學關系的InSAR三維形變估計函數模型。

1.2 利用VCE實現不同類觀測值之間的精確定權

(12)

式中，Wj代表第j類InSAR觀測的初始權矩陣，一般取值為單位矩陣。進而根據方差分量估計算法[41]可得各類InSAR觀測值的單位權中誤差向量為

(13)

根據方差分量估計算法可知，當各類觀測值單位權中誤差近似相等時，即

(14)

(15)

當上述迭代過程結束時，可根據各個單位權中誤差及其相應權陣計算各類觀測值對應的方差矩陣為

(16)

同時，通過式(12)即可得到三維地表形變的最優估值。

在觀測值方差已知的前提下，根據誤差傳播定律可以計算得到三維形變的方差。在使用誤差傳播定律時，應盡可能地保證觀測值之間相互獨立。因此，在計算三維地表形變的方差時，本文將觀測方程式(11)簡化為

(17)

即僅考慮了觀測值與三維地表形變之間的幾何關系。此時，三維地表形變的方差矩陣為

(18)

針對不同像元重復利用本章介紹的算法，直至所有像素點的三維形變解算完成。

1.3 顧及斷裂線和自適應變化窗口的觀測值選取策略

在利用SM-VCE方法進行三維地表形變求解時，其中的關鍵一步是選取一定大小窗口范圍內的像素點建立觀測方程。原始SM-VCE方法選取固定大小窗口(如15×15像素)內的所有點來建立觀測方程[35]。這種選取方法在地表形變未發生跳變或者原始形變觀測值相干性較高時可以得到較高精度的三維形變[36]。然而，地震導致的地表形變在震中區域往往會發生跳變，且精度較高相位觀測值(如DInSAR獲取的形變觀測值)會發生失相干現象。鑒于此，本文提出顧及斷裂線和自適應變化窗口的觀測值選取策略，旨在獲取可靠的近場三維地表形變結果。

根據光學影像或者POT獲取的形變數據人工勾勒出地震發生后導致的主要斷裂線，基于該斷裂線即可排除窗口內與中心點不在斷層同一側的像素點。在設置初始窗口大小為S×S像素的情況下，當解算震中某像素點的三維形變時，S×S像素大小的窗口內的失相干像素將會被掩膜，以保證解算結果的精度和可靠性。但是，在這種情況下，精度較高的相位觀測值的像素點個數可能少于某一閾值Kthr，進而無法得到可靠的三維形變結果。本文提出的窗口大小自適應變化的觀測值選取策略將會在Kj

結合已有的地表應力應變模型在火山[35]、地震[31,40]方面的研究發現，參與SM建模的像素點個數大于200之后，即可得到較為穩健的未知參數平差值。基于此，本文設置窗口大小的初值S=15像素，像素點個數閾值Kthr=200。為了確定Sthr的取值大小，本文基于InSAR觀測值進行變異函數擬合[40,42]

C(D)=C0·e-D/D0

(19)

式中，C(D)表示距離為D的像素之間的協方差；C0代表區域內的方差；D0即為相關距離。事實上，當兩點間距離大于D0時，根據相關距離及像素空間分辨率即可確定不同InSAR觀測值對應的閾值Sthr,j。最終的閾值Sthr為

Sthr=min([Sthr,1Sthr,2…Sthr,j…Sthr,J])

(20)

式中，min(■)代表取最小值。本文取一系列Sthr,j的最小值作為最終閾值Sthr的取值，從而保證每一類觀測中距離為Sthr的兩點對應的觀測值都是相關的。在本文中，利用式(19)、式(20)對真實數據進行擬合，計算得到Sthr=63像素。

1.4 基于迭代加權最小二乘的窗口內同一類觀測值自適應定權方法

Lk=[?L/?xe,?L/?xn,?L/?xu]·Δk+L0

(21)

為了簡單起見，這里省略了下標j。顧及周圍的K個點，則可得

L=Bsm,L·lL

(22)

式中

lL=[?L/?xe?L/?xn?L/?xuL0]T

(23)

(24)

(25)

式中，°代表哈達瑪運算；G為對角矩陣，其對角線元素可認為是降權因子。G的第k個對角線元素gk可基于以下權函數確定

(26)

(27)

(28)

2 模擬試驗

本文利用2019年Ridgecrest地震的已有的斷層幾何和滑動分布數據[45]正演了三維地表形變場，整個形變場的大小為890×1120像素，像素分辨率約為100 m×100 m。根據真實數據的衛星成像幾何參數(表1)計算得到InSAR視線向與方位向的形變信息。在此基礎上，加入不同量級的高斯噪聲(升軌DInSAR視線向、升軌POT方位向、升軌POT視線向、降軌DInSAR視線向、降軌POT方位向、降軌POT視線向的高斯噪聲均方差分別為5、300、100、5、300和100 mm)，同時在不同的觀測值上隨機挑選了5%的像素點加入了粗差。需要說明的是，模擬試驗中粗差的模擬方法為：將整個影像上隨機選取一定比例p=5%的像素作為粗差點，粗差點上的數值為整個形變場中最大形變值的±m倍，其中m的取值與預設的M值有關，不同粗差點處m的取值不一定相同，但所有粗差點處m的取值服從[M-5,M+5]區間內的均勻分布，并且不同粗差點處m取值的正負也是隨機的。在此次模擬試驗中M=5。另外，為了更加符合真實數據中的InSAR噪聲源，本文在升降軌InSAR視線向的模擬數據中加入了模擬的大氣噪聲，其中大氣噪聲是利用代爾夫特理工大學開發的Matlab工具包中的分形函數fracsurf進行模擬的，分形維度為2.2，大氣噪聲的范圍為[-2.8 cm,2.8 cm]，在C波段哨兵數據的干涉圖上約為2個條紋。通過以上步驟即可得到用于模擬實驗的InSAR觀測數據(圖2)。

表1 真實試驗中所用到的哨兵1號衛星SAR數據基本信息

注：紅色加號的行列號為(454, 357)，是圖4中15×15窗口內中心像素的位置。黑色折線為斷裂線。圖2 模擬試驗中所用的InSAR觀測值Fig.2 The InSAR measurements in the simulation

隨后，分別利用原始SM-VCE方法和本文的SM-VCE方法計算了模擬數據的三維形變，如圖3(d)—(f)和(j)—(l)所示。整體而言，兩種方法得到的三維形變與模擬數據圖3(a)—(c)較為一致。圖3(g)—(i)和圖3(m)—(o)分別給出了兩種方法得到的三維形變殘差圖。可以看出，在斷層破裂區域，原始SM-VCE方法僅能利用精度較低的POT方法進行三維形變解算，且未區分窗口內斷層兩側的異質點，使得該區域的三維形變結果過于平滑、殘差值較大。本文的SM-VCE方法則通過窗口大小自適應變化使得高精度的InSAR形變觀測值也可以參與解算，同時通過顧及斷裂線位置剔除了窗口內的異質點，最終得到的斷層破裂區域的三維形變結果殘差更小、精度更高。通過兩種方法在斷層破裂區域的對比，說明了本文提出的顧及斷裂線和自適應變化窗口的觀測值選取策略的優越性。

注：黑色折線為斷裂線。圖3 模擬試驗中不同方法得到的三維形變場Fig.3 The 3D deformation obtained by different methods in the simulation

表2 模擬試驗中不同方法得到的三維地表形變均方根誤差對比

此外，鑒于IWLS主要是為了抑制粗差觀測值對三維形變結果的影響，本文進一步利用模擬試驗研究了不同粗差量級(即不同的M取值)和不同粗差個數(即不同的p取值)情況下，利用本文的IWLS方法獲取三維形變結果的精度(圖6)。為了避免大氣噪聲的影響，此處的模擬試驗中未添加大氣噪聲。不難發現，隨著觀測值中粗差比例的增加，得到的三維形變精度在不斷下降。當粗差比例小于等于40%時，不同粗差量級的模擬試驗均可得到較為可靠的形變結果。而且，當粗差比例小于等于40%時，粗差量級越大得到的三維形變結果精度越高。這是因為粗差量級越大，IWLS方法識別出粗差的正確率就會越高。但是，當粗差觀測值比例過大時(例如大于40%)，本文提及的IWLS方法也無法十分有效地抑制粗差觀測值對三維形變結果的影響。此外，南北向和垂直向形變結果均方根誤差呈現先減小后增大的趨勢(拐點出現在粗差比例約為20%時)。這是因為南北向形變主要是由POT方位向觀測值計算而來，并且在計算三維形變時，垂直向形變結果受南北向形變影響較大[46]，而POT方位向觀測值本身的方差就比較大，當粗差比例較小時，IWLS方法對粗差的敏感度并不高，隨著粗差比例的增大，使得觀測值噪聲的分布更加非高斯性，進而IWLS可以更為準確地識別出粗差、提高結果精度。而當粗差比例過大時，形變結果的精度會不斷降低。為了驗證上述觀點，本文將模擬試驗中DInSAR和POT觀測值中的高斯噪聲方差均設置為5 mm，得到的南北向和垂直向均方根誤差變化曲線和東西向變化曲線較為一致，未發現局部極小值的現象。

需要說明的是，SM-VCE方法在獲取三維同震地表形變時的優勢主要有兩點：①基于地表應力應變模型構建了鄰近點三維形變之間的模型關系，有效抑制了空間高頻噪聲，使得引入力學模型的InSAR三維地表形變估計函數模型更加符合實際情況；②顧及地表應力應變模型建立的InSAR三維地表形變估計函數模型顯著增加了多余觀測個數，使得利用方差分量估計進行驗后迭代定權成為可能，相比于傳統的先驗定權方法，方差分量估計算法可以顯著提高多源觀測數據的定權精度，進而可以獲得更加精確、可靠的三維地表形變場。相對于傳統僅利用觀測值與三維形變之間幾何關系進行求解的方法[17]而言，SM-VCE方法通過顧及形變在空間上的關聯關系以及觀測值之間的準確定權，可以得到更為精確的三維地表形變結果。

注：觀測值的權重是無量綱。圖4 行列號(454, 357)處15×15窗口內(a)—(f)6類觀測值數據及(g)—(l)IWLS得到的每一類觀測值內部的相對權重Fig.4 (a)—(f) The six kinds of measurements and (g)—(l) the corresponding weights from the IWLS method in the 15×15 window of the pixel (454, 357)

另外，基于地表應力應變模型構建形變在空間上的關聯關系在一定程度上可以認為是一個濾波過程，進而會導致SM-VCE方法對一些空間上的局部高頻形變信號不敏感。但是，地表應力應變模型作為一個力學模型，用于三維形變求解時，肯定比純數學的濾波效果更好。此外，傳統逐像素求解的WLS方法，由于缺少多余觀測數據，難以利用VCE算法進行驗后精確定權，對于SM-VCE方法而言，正是因為顧及地表應力應變模型建立的函數模型具有較多的多余觀測數據，才能夠利用方差分量估計算法進行觀測值精確定權。

同時，SM-VCE方法中的窗口大小是一個必不可少的參數。理想情況下，窗口范圍內的形變應滿足地表應力應變模型，即形變梯度應為一個常數。因此，窗口的范圍應遠小于感興趣的形變場尺度。對于高空間分辨率影像，或者大尺度形變場而言，本文中所述的15×15像素，甚至63×63像素均可以得到可靠的三維地表形變。對于本文提出的顧及斷裂線和自適應變化窗口的觀測值選取策略，當面對淺部發生的中小型地震時，由于斷裂帶小范圍產生的復雜形變特征造成窗口內的觀測值呈現異質性，此時應減小窗口閾值Sthr以保證近場三維地表形變結果的可靠性。

注：σj的單位為m。圖5 行列號(454, 357)處方差分量估計迭代過程中的單位權中誤差收斂Fig.5 Convergence of the variance component estimation at the pixel (454, 357)

注：均方根誤差的單位是m。圖6 不同粗差量級(即不同的M取值)和不同粗差比例(即不同的p取值)情況下，利用本文SM-VCE方法獲取的模擬試驗三維形變結果的精度對比Fig.6 The root mean squares errors of the three-dimensional deformations estimated by the enhanced SM-VCE method in this paper with the simulated experiments under different magnitude and different proportion of the gross error

在觀測值噪聲服從理想的高斯分布時，方差分量估計可以較為準確地得到各類觀測值的驗后權重。但是，真實的形變觀測值往往會包含各種復雜的誤差源信號(例如大氣噪聲、地形殘差、軌道誤差等)，在個別情況下，多源誤差信號整體上表現為嚴重的非高斯性，進而使得方差分量估計難以收斂，甚至出現負方差，因此，需要針對個別情況進行人為干預的先驗定權，或者引入更為穩健的算法來抑制相關誤差源的影響。

3 2019年Ridgecrest地震的InSAR三維地表形變場

3.1 研究區域及所用數據

2019年7月6日(UTC),美國加州Ridgecrest地區(震中35.770°N，117.599°W)發生了M7.1級地震，震源深度約為8.0 km。Ridgecrest地震的發震斷層位于東加利福尼亞剪切帶上(east California shear zone，ECSZ)。該區域是圣安德烈亞斯斷層(San Andreas fault，SAF)南段以東的一個地震活躍區域，與莫哈韋沙漠大致重合，具有多個與SAF平行的右旋走滑斷層。因此，對加州地震展開地表形變監測可為研究該地震，ECSZ及SAF的構造運動提供基礎資料。

針對2019年Ridgecrest地震，本次研究使用了升降軌哨兵1號衛星SAR數據(如表1和圖7)。基于瑞士的GAMMA軟件，分別利用DInSAR和POT技術進行了數據處理。其中，SRTM(shuttle radar topography mission) 30 m分辨率的DEM數據被用于去除地形相位及地形起伏引起的配準誤差。為了抑制失相干噪聲、提高計算效率，采用了32×8的多視因子(距離向×方位向)。在DInSAR處理過程中，利用基于最小二乘的濾波方法對多視后的干涉圖進行濾波處理[47]，基于最小費用流方法對濾波后的干涉圖進行解纏[48]。在POT的數據處理過程中，哨兵數據的匹配窗口為128×128像素，過采樣因子為2。

圖7 研究區域與數據覆蓋情況Fig.7 Shaded relief map of the study area

圖8為本文獲取的6個InSAR形變觀測數據。可以看出，Ridgecrest地震震中的形變可達米級。由于震中地表破壞較為嚴重，DInSAR獲取的形變觀測結果在斷裂帶附近失相干現象較為嚴重。相反，POT技術可以較好地抵制失相干現象，獲得了較為完整的地表形變場。但是，由于POT技術是根據強度匹配獲取地表形變，該技術往往比DInSAR技術獲取的形變精度低。因此，盡管融合升降軌POT數據即可獲得較為完整的三維形變，本文仍有必要通過窗口大小自適應變化的觀測值選取策略來使得更高精度的DInSAR形變觀測數據參與斷裂帶附近的三維形變解算。此外，圖8(b)、(c)、(e)和(f)顯示，基于POT獲取的地表形變觀測數據噪聲較大。圖8(b)、(c)和(e)中的形變場左側還出現了明顯的條帶信號，通過對比光學影像發現，條帶信號與光學影像中的公路較為相似，并且公路走向與衛星飛行方向越相近，POT方位向觀測值中的條帶信號越強，進而推測此條帶信號是因為公路對電磁波反射信號過強，導致影像匹配時公路上的強反射信號沿公路走向較為相似，進而引起了配準錯誤。本文則基于IWLS方法，通過測量平差數據處理過程中的權重因子來抑制相關噪聲(包括圖8(b)、(c)和(e)中的條帶信號)對SM-VCE方法解算三維形變的影響。

3.2 InSAR三維地表形變場

與模擬試驗類似，本文分別利用原始SM-VCE方法和本文SM-VCE方法獲取了該地震的三維地表形變場(圖9)。可以看出，不同方法得到的三維地表形變具有高度一致性。其中，水平向和垂直向的最大形變分別約為2.1 m和0.38 m。然而，由于原始SM-VCE方法在斷層附近無法兼顧高精度的相位觀測數據，且未區分窗口內斷層兩側的異質點，圖9(d)—(f)中，由該方法得到的三維地表形變在斷層附近的三維形變場過于平滑，嚴重偏離真實情況。從圖9(g)—(i)，本文SM-VCE方法和原始SM-VCE方法得到的三維形變場差異圖中可以看出，兩種方法得到的三維形變在大部分區域是基本相同的，在斷裂帶附近差異最大，說明了本文SM-VCE方法獲取可靠斷裂帶附近三維形變的優越性。同時，在遠離斷裂帶區域(如圖9(h)中主斷裂帶的西南側)，本文SM-VCE方法可以較好地抑制一些局部誤差信號對三維地表形變結果的影響，說明了IWLS在抑制粗差觀測值方面的優勢。

此外，為了更加突出窗口自適應選點方法的優越性，本文僅基于DInSAR獲取的視線向和POT獲取的方位向4個形變觀測值，利用本文的SM-VCE方法獲取了該地震的三維地表形變結果，如圖10(a)—(c)所示。與圖9(a)—(c)基于6個觀測值的結果對比發現，二者的差異幾乎可以忽略，如圖10(d)—(f)，尤其在斷層破裂帶區域，在DInSAR完全失相干的情況下，本文的SM-VCE方法也可得到較為可靠的三維地表形變結果。然而，在DInSAR失相干區域，原始SM-VCE方法在15×15像素固定窗口內僅有兩個觀測值(即升降軌獲取的POT方位向的形變觀測值)，導致最終的三維地表形變場圖10(g)—(i)在DInSAR失相干區域必定會缺失，并且在斷裂帶附近，由于15×15像素固定窗口內的像素點個數較少，三維形變結果會出現較多的跳變信號。盡管在窗口大小為49×49像素情況下，原始SM-VCE方法可以得到完整的三維形變場，如圖10(j)—(l)，但是在斷裂帶附近，由于窗口內包含了較多的異質點，原始SM-VCE方法得到的三維形變結果會出現較多的錯誤形變信號，進一步說明了本文方法在獲取高精度、完整三維同震形變場的優勢。

注：正值代表地面朝向衛星運動，或沿衛星飛行方向運動,紫色曲線代表斷層線。圖8 不同數據不同方法獲取的Ridgecrest地震InSAR形變觀測值Fig.8 The InSAR measurements of the Ridgecrest earthquake

圖9 在使用升軌DInSAR視線向、升軌POT方位向、升軌POT視線向、降軌DInSAR視線向、降軌POT方位向、降軌POT視線向等6個觀測值情況下，不同方法得到的Ridgecrest地震三維地表形變Fig.9 The 3D deformation of the Ridgecrest earthquake obtained by different methods based on the six measurements in Fig.8

圖10 在使用升軌DInSAR視線向、升軌POT方位向、降軌DInSAR視線向、降軌POT方位向等4個觀測值情況下，不同方法得到的Ridgecrest地震三維地表形變Fig.10 The 3D deformation of the Ridgecrest earthquake obtained by different methods base on the four measurements in Fig.8(a), (b), (d), and (e)

需要說明的是，本文所提出的自適應選點方法主要是針對位于發震斷層處的地震破裂帶區域。一般情況下，結合歷史資料、實地調研及大地測量等數據可以較為準確地獲取發震斷層斷裂線數據，進而利用本文的自適應選點方法可以獲取可靠的近場三維同震形變場。實際地震中，部分非發震斷層在受到應力擠壓、拉張等作用下，也可能表現出形變場不連續的現象。與發震斷層處不同的是，InSAR技術在這些區域往往可以獲取完整的地表形變場，進而基于斷裂線的人工識別即可得到可靠的三維地表形變場。對于Ridgecrest地震而言，真實地表破裂較多[49]，在本文中為了簡單起見，僅使用了兩個主要的發震斷層破裂線計算三維同震形變場來說明本文方法的優越性。

此外，本文利用Nevada Geodetic Laboratory(http:∥geodesy.unr.edu/)免費公開的5個近場GNSS站點數據(站點分布如圖7所示)對兩種方法得到的三維地表形變進行了精度評估。原始SM-VCE方法在東西向、南北向和垂直向形變精度分別為7.8、6.5和2.6 cm，而本文SM-VCE方法則分別為7.9、5.6和2.6 cm。可以發現，本文方法相比于原始方法僅在南北向分量有一定的改善效果。與模擬試驗對比發現，真實觀測數據中的粗差量級和個數相對較小，僅在POT方位向形變觀測值中具有較為明顯的粗差觀測值信號，因此，在與GNSS對比中，本文SM-VCE僅在南北向分量改進較為明顯。此外，InSAR數據和GNSS數據的時空分辨率均不相同，GNSS數據本身就包含誤差，因此，利用GNSS進行InSAR精度評估并不完全可靠[16]。

4 總 結

融合多源InSAR形變觀測數據是實現InSAR三維同震地表形變監測的基本途徑，其中三維同震地表形變求解過程中建立的函數模型與隨機模型對三維形變結果的精度影響至關重要。本文詳細介紹了一種基于地表應力應變模型(SM)和方差分量估計(VCE)的InSAR三維地表形變監測方法(SM-VCE)。該方法通過顧及一定窗口范圍內不同像素點三維形變之間的力學關系(即SM)建立函數模型，進而利用VCE算法實現多源InSAR形變觀測數據隨機模型的準確確定。相比于傳統的單點求解方法，SM-VCE方法可得到更為準確的三維地表形變結果。

為了使得該方法可以得到更為精確的、完整的三維同震地表形變場，本文進一步提出了顧及斷裂線位置和窗口大小自適應變化的觀測值選取策略和一種基于迭代加權最小二乘的窗口內同一類觀測值相對定權方法。其中，觀測值選取策略是為了解決斷裂帶附近InSAR形變觀測數據易缺失且不連續導致的SM-VCE方法精度較低和形變結果缺失問題，而同一類觀測值內部相對定權方法是通過迭代加權最小二乘方法抑制窗口內的異常形變觀測值，進而提高三維形變的解算精度。

模擬試驗和基于升降軌哨兵數據2019年Ridgecrest地震三維形變的研究表明，本文提及的方法可得到更為精確的三維同震形變結果，尤其在斷層破裂帶附近，窗口優化的SM-VCE方法可在數據缺失的情況下，利用鄰近的有效InSAR形變觀測值恢復得到可靠的三維地表形變場。同時，本文通過模擬試驗驗證了觀測值中不同的粗差比例對三維形變結果的影響。結果表明，當粗差比例小于等于40%時，窗口優化的SM-VCE方法通過迭代加權最小二乘可顯著抑制粗差觀測值對結果的影響。

本文介紹的SM-VCE方法也可聯合星載、地基InSAR面觀測數據和GNSS、水準等傳統大地測量點觀測數據實現多種地質災害(例如火山噴發、滑坡等)的三維地表形變準確測量。同時，在獨立觀測數據較少的情況下，例如只有升降軌的DInSAR觀測數據，也可利用SM-VCE方法實現二維地表形變測量，或者，通過引入可靠的先驗信息約束(例如坡度約束、地下開采沉陷模型約束)實現三維地表形變測量。此外，本文介紹的窗口大小自適應變化的準則是窗口內的有效觀測值個數，只是從數值解算的穩定性方面進行考慮，未來可發展基于數據本身或者局部實際地表形變梯度的窗口大小自適應確定方法，以提高三維地表形變的精度及可靠性。