高宏偉，史先琳，陳 晨，尹 勇，戴可人,2

(1.成都理工大學(xué) 地球科學(xué)學(xué)院，四川 成都 610059； 2.成都理工大學(xué) 地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國家重點(diǎn)實(shí)驗室,四川 成都 610059)

0 引 言

據(jù)中國地震網(wǎng)臺測定，2021年5月21日21時48分，云南大理州漾濞縣(北緯25.67°，東經(jīng)99.87°)發(fā)生6.4級地震，震源深度為8 km.此次地震造成房屋、公路等基礎(chǔ)設(shè)施受損，人民生命安全受到威脅，生態(tài)環(huán)境遭到破壞.及時掌握地震同震形變信息，對于災(zāi)后重建有著重要意義.

水準(zhǔn)測量、GNSS測量等傳統(tǒng)的地表沉降測量受到地理環(huán)境影響，難以及時、全面地監(jiān)測整個震域[1-3].合成孔徑雷達(dá)差分干涉測量技術(shù)(Differential Interferometric Synthetic Aperture Radar, D-InSAR)具有全天時、全天候、覆蓋范圍廣、空間分辨率高等優(yōu)點(diǎn)[4-7]，彌補(bǔ)了傳統(tǒng)監(jiān)測技術(shù)的不足，已經(jīng)成為大范圍沉降監(jiān)測和地表形變的一種重要監(jiān)測手段[8].由于合成孔徑雷達(dá)差分干涉測量技術(shù)獲取的是沿雷達(dá)視線方向(line of sight，LOS)的形變，并不能完整地表達(dá)地表形變信息[9-11]，需要采用多個軌道SAR數(shù)據(jù)得到的雷達(dá)視線方向(LOS)的形變聯(lián)合解算垂直方向上的形變、水平東西方向上的形變和水平南北方向上的形變[12-13]，聯(lián)合升降軌解算使用雷達(dá)衛(wèi)星的平均入射角，分析如何獲取精準(zhǔn)雷達(dá)衛(wèi)星入射角度和提高解算形變量的精度.

圖1 研究區(qū)域地形圖 Fig.1 The topographic map of study area

本文利用合成孔徑雷達(dá)差分干涉測量技術(shù)對漾濞地震區(qū)域進(jìn)行地表形變監(jiān)測，根據(jù)SAR衛(wèi)星飛行幾何關(guān)系，精確雷達(dá)衛(wèi)星的入射角，提高垂直方向和水平東西方向的形變解算精度.

圖2 雙軌差分干涉原理圖Fig.2 Two-track differential interference schematic

1 研究區(qū)域及數(shù)據(jù)處理

1.1 研究區(qū)概況

漾濞縣位于大理州中部，東與大理市相鄰，西與永平縣接壤，南與長寧縣相連，北與洱海縣相交，屬橫斷山滇西高山峽谷區(qū)，地勢由北向南逐漸降低，該區(qū)域附近地形起伏較大，斷裂帶較多，容易發(fā)生地震，影響地表基礎(chǔ)設(shè)施的安全.研究區(qū)域如圖1所示.

1.2 InSAR原理

如圖2所示，SAR天線S1、S2距離地物的斜距為R1、R2，S1、S2表示同一衛(wèi)星在不同時間不同位置對同一地面目標(biāo)進(jìn)行觀測[14-16]，S1、S2因為空間位置不同，把S1、S2的相距距離稱為空間基線，設(shè)定S1為主影像，S2為從影像，入射角為θ，H表示S1相對于地面參考面的高度.干涉相位表達(dá)公式為：

φ=φref+φtop+φdef+φnoise

(1)

式中：φ代表干涉相位，φref代表參考橢球相位，φtop代表地形相位，φdef代表形變相位[17]，φnoise代表干涉相位中噪聲引起的相位誤差.為了得到地表的真實(shí)形變，需要將參考橢球相位、地形相位和干涉相位中噪聲引起的相位誤差去掉[18].

所以形變相位可以表達(dá)為：

φdef=φ-φref-φtop-φnoise

(1)

在入射角、地面高程空間基線已知的情況下，參考橢球相位和地形相位為：

(3)

采用幾何成像參考橢球面相位模型去除參考橢球相位，通過外部DEM去除地形相位[19]，使用濾波減弱噪聲引起的相位誤差[20-21].

1.3 InSAR數(shù)據(jù)處理

本文采用雷達(dá)數(shù)據(jù)為Sentinel-1的干涉寬幅模式的VV極化方式影像數(shù)據(jù)，Sentinel-1是歐洲空間局(European Space Agency，ESA)發(fā)射的對地觀測衛(wèi)星[22]，采用C波段成像，波長為5.6 cm，重訪周期為 12 d[23-24]，具體參數(shù)見表1.分別使用降軌軌道號為135和升軌軌道號為99的Sentinel-1A數(shù)據(jù)對地震區(qū)域進(jìn)行沉降監(jiān)測，使用30 m空間分辨率的AW3D30的DEM，利用Goldstein濾波方法對干涉圖進(jìn)行濾波處理，采用最小費(fèi)用流方法(Minimum Cost Flow，MCF)對相位進(jìn)行解纏[25-26]，再經(jīng)過軌道精煉、重去平和地理編碼，得到地理坐標(biāo)系的形變圖[27].

表1 Sentinel-1參數(shù)

2 地表形變分析

基于合成孔徑雷達(dá)差分干涉測量技術(shù)得到了漾濞地震形變場見圖3、圖4.在降軌形變圖中可知，此次地震形成一條北西-南東走向的分界線，最后的結(jié)算結(jié)果顯示分界線的東北邊是抬升，西南邊是沉降，而當(dāng)前降軌形變場沿雷達(dá)視線方向最大沉降量為 9.00 cm，最大的抬升量為9.17 cm;升軌形變場最大沉降量為9.39 cm,最大抬升量為6.26 cm.

圖3 升軌形變圖Fig.3 Rail lifting deformation diagram

圖4 降軌形變圖Fig.4 Descending deformation diagram

如圖3升軌剖線分布圖顯示，在沉降區(qū)域繪制AA′和BB′剖面線，在圖5(a)中剖線走向如同一個“V”字形，在7 km處，沉降量達(dá)到了最大，最大的沉降量為8.60 cm.在圖5(b)中，在8～14 km處剖線走向也如同“V”字形，在9 km處沉降量達(dá)到最大，沉降量最大可達(dá)5.80 cm.如圖4降軌剖線分布圖顯示，繪制CC′和DD′剖面線橫跨抬升和沉降區(qū)域.由圖5(c)可知在5～8 km處出現(xiàn)了由沉降到抬升的急劇變化，在 7 km 前沉降量最大為6.40 cm，在7 km后抬升量最高為8.20 cm.由圖5(d)可知在5～8 km處同樣出現(xiàn)了由沉降到抬升的急劇變化，在8 km前沉降量最大為5.80 cm，在8 km后抬升量最高為7.30 cm.由剖線圖5(a)、5(c)和5(d)得出，在5～8 km處均發(fā)生急劇的變化，可以初步判斷此處為形變較嚴(yán)重的點(diǎn).在圖5(c)和5(d)在6～14 km處的走勢相同，而圖5(b)沉降量比圖5(a)的沉降量小一些，走勢有一定差異，可以說明剖線BB′在6～14 km處受到的影響比剖線AA′在6～14 km處受到的影響小.

圖5 剖線圖Fig.5 Sectional drawing

3 地表二維形變分析

D-InSAR技術(shù)是獲取地表一維形變，即真實(shí)地表形變在LOS的投影[28]，并不能反映地表真實(shí)形變，單一軌道LOS形變值，不容易通過差分干涉形變場分求得地表真實(shí)的形變場，利用三個以上不同軌道的LOS形變值，可聯(lián)合求解沿水平東西方向、水平南北方向和垂直方向的形變值[12]，由于現(xiàn)在的在軌衛(wèi)星屬于近極地飛行特點(diǎn)，監(jiān)測南北方向的形變不敏感，因此忽略水平南北方向的形變.根據(jù)整幅影像的行列關(guān)系(如圖6所示)，確定發(fā)生地表形變的位置所屬的行列，得到形變位置是屬于近斜距、中斜距和遠(yuǎn)斜距，由此得到雷達(dá)衛(wèi)星精確的入射角度，這樣就可以提高垂直方向和水平東西向的形變以求解精度.聯(lián)合升、降軌LOS方向的觀測值，可得到形變場的水平東西方向和垂直方向的地表形變場[29-31].

在圖6中，Dlos是水平東西方向、水平南北方向和垂直方向三個分矢量在視線方向上的投影之和，水平東西方向的形變?yōu)閐E，水平南北方向為dN，垂直方向的形變?yōu)閐U[32]，根據(jù)雷達(dá)衛(wèi)星的空間幾何關(guān)系可以表示為：

(4)

由于SAR影像監(jiān)測南北方向能力弱，所以忽略南北方向的形變，得到下式

(5)

式中：θ為衛(wèi)星入射角，α代表衛(wèi)星飛行方向的方位角[33].

將升、降軌的形變數(shù)據(jù)帶入式子中，下角標(biāo)A代表升軌，下角標(biāo)D代表降軌.GR代表地距方向.

圖6 SAR影像幾何關(guān)系Fig.6 Geometric relationship of SAR image

(6)

由式(6)可得：

(7)

根據(jù)二維形變解算可以得到水平東西方向和垂直方向的形變，在垂直方向上如圖7(a)所示，抬升區(qū)域最大的形變值為 5.52 cm，沉降區(qū)域最大的形變值為 -6.96 cm，在水平方向上如圖7(b)所示，抬升區(qū)域最大的形變值為 7.60 cm，沉降區(qū)域最大的形變值為 -9.76 cm.由表2可知，在垂直方向上的平均形變量為 0.74 cm，在水平方向上的平均形變量為 -1.27 cm，從平均值形變量中得出漾濞縣此次地震在垂直方向上的形變量級比水平方向的形變量級小.

(a) 垂直形變圖 (b) 水平形變圖圖7 形變分解圖Fig.7 Deformation decomposition diagram

表2 形變量統(tǒng)計表

4 結(jié) 論

本文采用了歐空局提供的Sentinel-1A升、降軌數(shù)據(jù)，通過合成孔徑雷達(dá)差分干涉測量技術(shù)，對漾濞縣地震區(qū)域進(jìn)行了形變監(jiān)測，得到了升、降軌沿雷達(dá)視線方向(LOS)的地表形變量，由升軌形變圖可得，形變場抬升最大量為6.26 cm，沉降最大量為9.39 cm.由降軌形變圖可知，形變場抬升最大量為9.17 cm，沉降最大量為9.00 cm.根據(jù)發(fā)生地表形變在影像上的位置，來確定雷達(dá)衛(wèi)星入射角.然后聯(lián)合升、降軌對形變場進(jìn)行地表形變分解，得到垂直方向上的形變平均值為-0.74 cm，水平方向上的形變平均值為 -1.27 cm，結(jié)果說明在垂直方向上的形變影響要比水平東西方向上的形變小.