InSAR三維同震地表形變監(jiān)測(cè)
——窗口優(yōu)化的SM-VCE算法

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

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

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

然而，上述研究均基于單個(gè)點(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觀測(cè)，文獻(xiàn)[31]對(duì)SISTEM方法進(jìn)行了擴(kuò)展，僅需InSAR觀測(cè)值即可進(jìn)行高精度三維形變解算。為了方便起見(jiàn)，本文將DInSAR、MAI、POT和BOI等技術(shù)獲取的形變觀測(cè)數(shù)據(jù)統(tǒng)稱為InSAR觀測(cè)值。然而，在解算真實(shí)三維地表形變時(shí)，必須要融合不同方向的多源異質(zhì)觀測(cè)數(shù)據(jù)，因此精確確定不同數(shù)據(jù)之間的權(quán)重比例具有重要意義[6,32]。為了實(shí)現(xiàn)此目的，科研人員研究利用一個(gè)固定的滑動(dòng)窗口求解窗口內(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è)觀測(cè)數(shù)據(jù)在某一數(shù)值上下一定范圍內(nèi)是隨機(jī)分布的，但是對(duì)于地震等導(dǎo)致的大梯度形變而言，往往難以滿足。也有相關(guān)學(xué)者在假設(shè)遠(yuǎn)場(chǎng)區(qū)域不包含形變信號(hào)的前提下，將遠(yuǎn)場(chǎng)一定區(qū)域數(shù)據(jù)的均方根值作為該類觀測(cè)值整個(gè)區(qū)域的先驗(yàn)方差[33-34]。這種方法忽略了InSAR觀測(cè)噪聲在空間上的差異性，因此也會(huì)影響三維地表形變的求解精度。文獻(xiàn)[25]在求解時(shí)序數(shù)據(jù)的平均形變速率時(shí)引入了方差分量估計(jì)(variance component estimation,VCE)進(jìn)行定權(quán)，驗(yàn)證了VCE在多源InSAR觀測(cè)數(shù)據(jù)定權(quán)方面的優(yōu)越性。但是，對(duì)于瞬時(shí)形變場(chǎng)的求解，由于缺少多余觀測(cè)，應(yīng)用VCE進(jìn)行定權(quán)受到了極大限制。

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

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

1 SM-VCE方法

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

1.1 基于SM建立觀測(cè)方程

dk=H·Δk+d0

(1)

(2)

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

H=S+R

(3)

(4)

(5)

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

(6)

(7)

(8)

式中

(9)

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

(10)

式中

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

L=B·l

(11)

式中

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

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

(12)

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

(13)

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

(14)

(15)

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

(16)

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

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

(17)

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

(18)

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

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

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

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

結(jié)合已有的地表應(yīng)力應(yīng)變模型在火山[35]、地震[31,40]方面的研究發(fā)現(xiàn)，參與SM建模的像素點(diǎn)個(gè)數(shù)大于200之后，即可得到較為穩(wěn)健的未知參數(shù)平差值?；诖?，本文設(shè)置窗口大小的初值S=15像素，像素點(diǎn)個(gè)數(shù)閾值Kthr=200。為了確定Sthr的取值大小，本文基于InSAR觀測(cè)值進(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時(shí)，根據(jù)相關(guān)距離及像素空間分辨率即可確定不同InSAR觀測(cè)值對(duì)應(yīng)的閾值Sthr,j。最終的閾值Sthr為

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

(20)

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

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

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

(21)

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

L=Bsm,L·lL

(22)

式中

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

(23)

(24)

(25)

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

(26)

(27)

(28)

2 模擬試驗(yàn)

本文利用2019年Ridgecrest地震的已有的斷層幾何和滑動(dòng)分布數(shù)據(jù)[45]正演了三維地表形變場(chǎng)，整個(gè)形變場(chǎng)的大小為890×1120像素，像素分辨率約為100 m×100 m。根據(jù)真實(shí)數(shù)據(jù)的衛(wèi)星成像幾何參數(shù)(表1)計(jì)算得到InSAR視線向與方位向的形變信息。在此基礎(chǔ)上，加入不同量級(jí)的高斯噪聲(升軌DInSAR視線向、升軌POT方位向、升軌POT視線向、降軌DInSAR視線向、降軌POT方位向、降軌POT視線向的高斯噪聲均方差分別為5、300、100、5、300和100 mm)，同時(shí)在不同的觀測(cè)值上隨機(jī)挑選了5%的像素點(diǎn)加入了粗差。需要說(shuō)明的是，模擬試驗(yàn)中粗差的模擬方法為：將整個(gè)影像上隨機(jī)選取一定比例p=5%的像素作為粗差點(diǎn)，粗差點(diǎn)上的數(shù)值為整個(gè)形變場(chǎng)中最大形變值的±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é)開(kāi)發(fā)的Matlab工具包中的分形函數(shù)fracsurf進(jìn)行模擬的，分形維度為2.2，大氣噪聲的范圍為[-2.8 cm,2.8 cm]，在C波段哨兵數(shù)據(jù)的干涉圖上約為2個(gè)條紋。通過(guò)以上步驟即可得到用于模擬實(shí)驗(yàn)的InSAR觀測(cè)數(shù)據(jù)(圖2)。

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

注：紅色加號(hào)的行列號(hào)為(454, 357)，是圖4中15×15窗口內(nèi)中心像素的位置。黑色折線為斷裂線。圖2 模擬試驗(yàn)中所用的InSAR觀測(cè)值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é)果過(guò)于平滑、殘差值較大。本文的SM-VCE方法則通過(guò)窗口大小自適應(yīng)變化使得高精度的InSAR形變觀測(cè)值也可以參與解算，同時(shí)通過(guò)顧及斷裂線位置剔除了窗口內(nèi)的異質(zhì)點(diǎn)，最終得到的斷層破裂區(qū)域的三維形變結(jié)果殘差更小、精度更高。通過(guò)兩種方法在斷層破裂區(qū)域的對(duì)比，說(shuō)明了本文提出的顧及斷裂線和自適應(yīng)變化窗口的觀測(cè)值選取策略的優(yōu)越性。

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

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

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

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

注：觀測(cè)值的權(quán)重是無(wú)量綱。圖4 行列號(hào)(454, 357)處15×15窗口內(nèi)(a)—(f)6類觀測(cè)值數(shù)據(jù)及(g)—(l)IWLS得到的每一類觀測(cè)值內(nèi)部的相對(duì)權(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)建形變?cè)诳臻g上的關(guān)聯(lián)關(guān)系在一定程度上可以認(rèn)為是一個(gè)濾波過(guò)程，進(jìn)而會(huì)導(dǎo)致SM-VCE方法對(duì)一些空間上的局部高頻形變信號(hào)不敏感。但是，地表應(yīng)力應(yīng)變模型作為一個(gè)力學(xué)模型，用于三維形變求解時(shí)，肯定比純數(shù)學(xué)的濾波效果更好。此外，傳統(tǒng)逐像素求解的WLS方法，由于缺少多余觀測(cè)數(shù)據(jù)，難以利用VCE算法進(jìn)行驗(yàn)后精確定權(quán)，對(duì)于SM-VCE方法而言，正是因?yàn)轭櫦暗乇響?yīng)力應(yīng)變模型建立的函數(shù)模型具有較多的多余觀測(cè)數(shù)據(jù)，才能夠利用方差分量估計(jì)算法進(jìn)行觀測(cè)值精確定權(quán)。

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

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

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

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

3 2019年Ridgecrest地震的InSAR三維地表形變場(chǎng)

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

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

針對(duì)2019年Ridgecrest地震，本次研究使用了升降軌哨兵1號(hào)衛(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處理過(guò)程中，利用基于最小二乘的濾波方法對(duì)多視后的干涉圖進(jìn)行濾波處理[47]，基于最小費(fèi)用流方法對(duì)濾波后的干涉圖進(jìn)行解纏[48]。在POT的數(shù)據(jù)處理過(guò)程中，哨兵數(shù)據(jù)的匹配窗口為128×128像素，過(guò)采樣因子為2。

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

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

3.2 InSAR三維地表形變場(chǎng)

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

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

注：正值代表地面朝向衛(wèi)星運(yùn)動(dòng)，或沿衛(wèi)星飛行方向運(yùn)動(dòng),紫色曲線代表斷層線。圖8 不同數(shù)據(jù)不同方法獲取的Ridgecrest地震InSAR形變觀測(cè)值Fig.8 The InSAR measurements of the Ridgecrest earthquake

圖9 在使用升軌DInSAR視線向、升軌POT方位向、升軌POT視線向、降軌DInSAR視線向、降軌POT方位向、降軌POT視線向等6個(gè)觀測(cè)值情況下，不同方法得到的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個(gè)觀測(cè)值情況下，不同方法得到的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)

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

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

4 總 結(jié)

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

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

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

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