劉計(jì)洪，胡 俊，李志偉，朱建軍

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

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

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

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

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

1 SM-VCE方法

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

1.1 基于SM建立觀測方程

dk=H·Δk+d0

(1)

(2)

文獻(xiàn)[30]將矩陣H表示為一個對稱矩陣和一個非對稱矩陣之和

H=S+R

(3)

(4)

(5)

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

(6)

(7)

(8)

式中

(9)

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

(10)

式中

若假設(shè)點(diǎn)P0周圍有Kj個第j類InSAR觀測值，則最終可得

L=B·l

(11)

式中

此時，建立了綜合考慮InSAR觀測值與三維形變之間的幾何關(guān)系，以及地表鄰近點(diǎn)三維形變之間的力學(xué)關(guān)系的InSAR三維形變估計(jì)函數(shù)模型。

1.2 利用VCE實(shí)現(xiàn)不同類觀測值之間的精確定權(quán)

(12)

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

(13)

根據(jù)方差分量估計(jì)算法可知，當(dāng)各類觀測值單位權(quán)中誤差近似相等時，即

(14)

(15)

當(dāng)上述迭代過程結(jié)束時，可根據(jù)各個單位權(quán)中誤差及其相應(yīng)權(quán)陣計(jì)算各類觀測值對應(yīng)的方差矩陣為

(16)

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

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

(17)

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

(18)

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

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

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

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

結(jié)合已有的地表應(yīng)力應(yīng)變模型在火山[35]、地震[31,40]方面的研究發(fā)現(xiàn)，參與SM建模的像素點(diǎn)個數(shù)大于200之后，即可得到較為穩(wěn)健的未知參數(shù)平差值?；诖?，本文設(shè)置窗口大小的初值S=15像素，像素點(diǎn)個數(shù)閾值Kthr=200。為了確定Sthr的取值大小，本文基于InSAR觀測值進(jìn)行變異函數(shù)擬合[40,42]

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

(19)

式中，C(D)表示距離為D的像素之間的協(xié)方差；C0代表區(qū)域內(nèi)的方差；D0即為相關(guān)距離。事實(shí)上，當(dāng)兩點(diǎn)間距離大于D0時，根據(jù)相關(guān)距離及像素空間分辨率即可確定不同InSAR觀測值對應(yīng)的閾值Sthr,j。最終的閾值Sthr為

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

(20)

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

1.4 基于迭代加權(quán)最小二乘的窗口內(nèi)同一類觀測值自適應(yīng)定權(quán)方法

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

(21)

為了簡單起見，這里省略了下標(biāo)j。顧及周圍的K個點(diǎn)，則可得

L=Bsm,L·lL

(22)

式中

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

(23)

(24)

(25)

式中，°代表哈達(dá)瑪運(yùn)算；G為對角矩陣，其對角線元素可認(rèn)為是降權(quán)因子。G的第k個對角線元素gk可基于以下權(quán)函數(shù)確定

(26)

(27)

(28)

2 模擬試驗(yàn)

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

表1 真實(shí)試驗(yàn)中所用到的哨兵1號衛(wèi)星SAR數(shù)據(jù)基本信息

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

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

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

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

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

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

注：觀測值的權(quán)重是無量綱。圖4 行列號(454, 357)處15×15窗口內(nèi)(a)—(f)6類觀測值數(shù)據(jù)及(g)—(l)IWLS得到的每一類觀測值內(nèi)部的相對權(quán)重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)

另外，基于地表應(yīng)力應(yīng)變模型構(gòu)建形變在空間上的關(guān)聯(lián)關(guān)系在一定程度上可以認(rèn)為是一個濾波過程，進(jìn)而會導(dǎo)致SM-VCE方法對一些空間上的局部高頻形變信號不敏感。但是，地表應(yīng)力應(yīng)變模型作為一個力學(xué)模型，用于三維形變求解時，肯定比純數(shù)學(xué)的濾波效果更好。此外，傳統(tǒng)逐像素求解的WLS方法，由于缺少多余觀測數(shù)據(jù)，難以利用VCE算法進(jìn)行驗(yàn)后精確定權(quán)，對于SM-VCE方法而言，正是因?yàn)轭櫦暗乇響?yīng)力應(yīng)變模型建立的函數(shù)模型具有較多的多余觀測數(shù)據(jù)，才能夠利用方差分量估計(jì)算法進(jìn)行觀測值精確定權(quán)。

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

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

注：均方根誤差的單位是m。圖6 不同粗差量級(即不同的M取值)和不同粗差比例(即不同的p取值)情況下，利用本文SM-VCE方法獲取的模擬試驗(yàn)三維形變結(jié)果的精度對比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

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

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

3.1 研究區(qū)域及所用數(shù)據(jù)

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

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

圖7 研究區(qū)域與數(shù)據(jù)覆蓋情況Fig.7 Shaded relief map of the study area

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

3.2 InSAR三維地表形變場

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

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

注：正值代表地面朝向衛(wèi)星運(yùn)動，或沿衛(wèi)星飛行方向運(yùn)動,紫色曲線代表斷層線。圖8 不同數(shù)據(jù)不同方法獲取的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)

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

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

4 總 結(jié)

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

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

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

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