韓 鳴，張永志，程 冬，尹 鵬

(長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院，陜西 西安 710054)

傳統(tǒng)的大地測(cè)量技術(shù)包括GPS、水準(zhǔn)測(cè)量等采集的數(shù)據(jù)都是離散分布的，并且受時(shí)空條件的影響，難以準(zhǔn)確和及時(shí)獲取大面積地形測(cè)量數(shù)據(jù)[1]。差分合成孔徑雷達(dá)干涉測(cè)量(DInSAR)技術(shù)是一種快速發(fā)展的大地測(cè)量技術(shù)，能夠全天候獲取地面高精度、大面積、低成本的地表變形信息，在地質(zhì)災(zāi)害監(jiān)測(cè)和評(píng)估方面已經(jīng)取得一系列成果。

在各種地質(zhì)災(zāi)害中，地震和火山運(yùn)動(dòng)及地表沉降和山體滑坡等均為需要利用測(cè)繪技術(shù)研究的問(wèn)題[2]。DInSAR技術(shù)通過(guò)處理兩幅或兩幅以上的雷達(dá)影像，提取地表形變相位差值，計(jì)算得到厘米級(jí)甚至毫米級(jí)的地表沿視線方向(line-of-sight,LOS)的形變信息。利用覆蓋發(fā)震區(qū)域的震前和震后兩幅雷達(dá)影像，通過(guò)DInSAR結(jié)合數(shù)字高程模型數(shù)據(jù)DEM處理就能獲取視線向的同震形變場(chǎng)。在伊拉克和伊朗邊境發(fā)生的Mw7.3級(jí)(USGS發(fā)布)地震，美國(guó)地質(zhì)調(diào)查局(USGS)研究認(rèn)為，該地震的震源機(jī)制主要是逆沖兼具輕微右旋走滑地震，為了更深入地研究該地震，本文利用InSAR技術(shù)及升降軌3個(gè)視角的雷達(dá)影像數(shù)據(jù)，重建了此次地震的三維同震形變場(chǎng)，為更深入研究和認(rèn)識(shí)該地震發(fā)震機(jī)制提供必要的基礎(chǔ)。

1 地質(zhì)構(gòu)造背景

2017年11月12日在伊拉克的蘇萊曼尼亞省和伊朗的克爾曼沙阿省之間發(fā)生了Mw7.3級(jí)地震，造成了至少530人死亡，超過(guò)7200人受傷。震中位置為如圖1所示的五角星處(34.911°N，45.959°E，USGS)[3]。研究區(qū)屬于伊朗高原，伊朗高原由許多地形起伏較大的陸地塊體邊界組成，阿拉伯板塊以相對(duì)于歐亞板塊20 mm/a的速度向北移動(dòng)，使阿拉伯板塊和歐亞板塊持續(xù)碰撞，形成了扎格羅斯造山帶，如圖1所示。圖1中的曲線是阿拉伯板塊與歐亞板塊的邊界線[4]。碰撞帶邊界的斜向應(yīng)力聚集引起了此次地震，發(fā)震區(qū)域處于阿拉伯-歐亞板塊碰撞帶的中心，緊鄰扎格羅斯山，USGS公布的震中位于扎格羅斯山前斷層(ZMFF)附近，這是一條逆沖斷層,圖1中的小圓圈為余震分布，大部分余震也都發(fā)生在扎格羅斯山前斷層(ZMFF)附近[5]。

圖1 研究區(qū)震中與地質(zhì)構(gòu)造背景

2 三維同震形變場(chǎng)解算方法

視線向形變場(chǎng)由兩幅SAR影像干涉處理得到，同一地區(qū)由于衛(wèi)星的飛行方向和雷達(dá)照射方向的不同，處理得到的視線向形變場(chǎng)也是不同的，甚至是相反的[6]。衛(wèi)星雷達(dá)傳感器成像模式分為升軌和降軌，如圖2所示。三維形變分別為北向形變，東向形變和垂直向形變，任何變形都可以通過(guò)三維形變來(lái)合成。但是3個(gè)方向的形變對(duì)視線向形變的貢獻(xiàn)各不相同[7]，Sentinel1-A衛(wèi)星的視線方向?yàn)榇怪庇谛l(wèi)星飛行方向的右下方，因此雷達(dá)傳感器對(duì)垂直向運(yùn)動(dòng)最敏感；而衛(wèi)星的飛行方向是近南北向的，如果地震斷層引起的地表運(yùn)動(dòng)是近NS向的，則將對(duì)InSAR觀測(cè)結(jié)果非常不利，甚至觀測(cè)不到任何形變信息，比如巴姆地震和美國(guó)的Hector Mine地震[8-9]。根據(jù)雷達(dá)成像的幾何關(guān)系，影像干涉得到的視線向形變用dLOS表示，約定地表接近衛(wèi)星的dLOS值為正，遠(yuǎn)離衛(wèi)星的dLOS值為負(fù)，三維形變(垂向dU,北向dN,東向dE)與視線向形變dLOS之間的幾何關(guān)系用公式表達(dá)為

dLOS=dUcosθ-dNsinθcosα-dEsinθsinα

(1)

圖2 衛(wèi)星軌道方向

式中,dLOS為InSAR視線向觀測(cè)值；dU為地表垂向形變；dN為地表南北向形變；dE為地表東西向形變；θ為雷達(dá)入射角；α為方位視線方向，即距離向與北向的夾角，該式同時(shí)適用于升降軌的投影計(jì)算。由式(1)可知，要獲取真實(shí)的地表三維形變,至少需要3個(gè)不同視線向的觀測(cè)量。發(fā)生地震的區(qū)域有3對(duì)不同成像視角的雷達(dá)影像時(shí)，即可利用這3對(duì)影像數(shù)據(jù)得到的視線向形變場(chǎng)進(jìn)行直接解算地表三維形變[10]，根據(jù)式(1)有

(2)

令

(3)

3 視線向同震形變場(chǎng)提取與分析

3.1 數(shù)據(jù)選取與處理流程

Sentinel1-A是由歐空局2014年4月發(fā)射的對(duì)地觀測(cè)衛(wèi)星，在近地太陽(yáng)同步軌道上運(yùn)行，重訪周期為12 d，選取Sentinel1-A衛(wèi)星的干涉寬幅(interferometric wide swath,IW)模式，對(duì)地表進(jìn)行遞進(jìn)掃描得到3個(gè)字條帶，該成像模式采用中等分辨率(5 m×20 m)，對(duì)地面數(shù)據(jù)提取的范圍達(dá)到240 km[11]。本文選取了升降軌不同的3對(duì)影像，共6景數(shù)據(jù)，覆蓋范圍如圖1中黑色虛線框所示，基本信息見表1。選取的DEM數(shù)據(jù)為SRTM航天飛機(jī)雷達(dá)地形測(cè)量任務(wù)數(shù)據(jù)，地面分辨率為90 m。

表1 Sentinel1-A(來(lái)源ESA)衛(wèi)星數(shù)據(jù)基本信息

常規(guī)二通差分DInSAR處理流程為影像配準(zhǔn)、影像干涉、去平地相位、干涉圖濾波、相位解纏、地理編碼[12]。Sentinel1-A數(shù)據(jù)覆蓋范圍較寬，本文裁剪了比震區(qū)范圍稍大的部分影像進(jìn)行處理，在處理過(guò)程中使用地面分辨率為90 m的SRTM數(shù)據(jù)去除干涉圖中的地形相位分量，研究區(qū)處于植被覆蓋率低的中東地區(qū)，忽略由地表植被引起的干涉失相干影響[13]。為了盡可能消除干涉圖的噪音，先對(duì)干涉影像進(jìn)行距離向4視數(shù)，方位向1視數(shù)的多視處理，得到距離向分辨率為16.621 8 m，方位向分辨率為13.968 m，利用Goldstein濾波方法進(jìn)行濾波,其濾波器是可變的，提高了干涉條紋的清晰度，減少了由空間基線或時(shí)間基線引起的失相干的噪聲[14-15]。采用最小費(fèi)用流法(MCF)進(jìn)行相位解纏，設(shè)定相干性閾值，對(duì)相干性低于0.2的區(qū)域進(jìn)行掩膜處理。

3.2 視線向同震形變場(chǎng)提取分析

利用最終的解纏結(jié)果，轉(zhuǎn)換為視線向的形變信息，獲取升降軌3種視角的干涉圖和視線向形變信息。如圖3—圖5所示，分別為差分干涉條紋圖和視線向形變圖，圖中F1—F2虛線為USGS推測(cè)出的發(fā)震斷層頂部，3幅影像的重疊區(qū)沒(méi)能完全覆蓋形變區(qū)，但仍涵蓋大部分形變區(qū)域，黑色虛線框?yàn)檫x取的三維形變解算區(qū)。升軌-A72的同震形變場(chǎng)視線向隆升中心集中在西南部，最大隆升量為84.7 cm，視線向沉降值最大為-26 cm；降軌-D79的同震形變場(chǎng)呈兩個(gè)扇形沿NW向?qū)ΨQ分布，西南部為視線向隆升區(qū)，最大隆升量為58.55 cm，東北部為視線向沉降區(qū),最大沉降量為-41.88 cm；降軌-D6的同震形變場(chǎng)的形變分布范圍和D79一致，隆升量與沉降量都比降軌-D79的結(jié)果稍小。

圖3 升軌-A72視線向形變

圖4 降軌-D79視線向形變

圖5 降軌-D6視線向形變

圖3—圖5中設(shè)置線條AB為視線向形變值剖面線，剖面提取信息如圖6所示,升軌-A72的視線向隆升量和沉降量與兩個(gè)降軌相比差別較大，是由于成像視角不同造成的，但3種視角下的視線向上升量與沉降量總體趨勢(shì)基本一致。3對(duì)影像的震后成像時(shí)間差值最大不超過(guò)5 d，可以忽略震后余滑的影響，認(rèn)為是對(duì)同一同震形變場(chǎng)的3種不同的觀測(cè)成果。

圖6 沿AB剖線的升降軌視線向形變

4 三維形變結(jié)果分析

利用3種視角的視線向觀測(cè)值，通過(guò)式(3)重建了研究區(qū)三維同震形變場(chǎng)，計(jì)算過(guò)程中選用表1中所列的平均入射角和平均視線方位角，計(jì)算結(jié)果如圖7所示。圖7中主要的形變區(qū)分為兩個(gè)部分，西南部的主要形變區(qū)域稱為C區(qū)，東北部的主要形變區(qū)域稱為D區(qū)。圖7(a)顯示，C區(qū)向上隆升，最大上升量達(dá)到了84 cm，D區(qū)沉降，最大沉降量為-34.5 cm；圖7(b)顯示，C區(qū)和D區(qū)都朝西運(yùn)動(dòng)，最大形變量為-50 cm;圖7(c)顯示，形變區(qū)整體向南運(yùn)動(dòng)，最大形變量為-91 cm，運(yùn)動(dòng)主要集中在形變區(qū)南部。結(jié)合三維形變場(chǎng)特征，得出整個(gè)形變區(qū)呈向上抬升且在水平方向上向西南方向運(yùn)動(dòng)的特點(diǎn)。通過(guò)上文所述，本文可以通過(guò)幾條線索來(lái)確認(rèn)發(fā)震主斷層：第一，視線向同震形變場(chǎng)和三維同震形變場(chǎng)的形變邊界約束位置緊鄰扎格羅斯山前斷層(ZMFF)；第二，主震和余震都發(fā)生在扎格羅斯山前斷層(ZMFF)周圍；第三，主震具有低傾角及純逆沖運(yùn)動(dòng)的特征，與扎格羅斯山前斷層的運(yùn)動(dòng)特征一致。結(jié)合地質(zhì)構(gòu)造和形變分析，推測(cè)此次地震的發(fā)震斷層為扎格羅斯山前斷層(ZMFF)。

圖7 三維同震形變場(chǎng)

5 結(jié) 語(yǔ)

本文利用Sentinel1-A衛(wèi)星的3對(duì)升降軌SAR數(shù)據(jù)，通過(guò)二軌差分DInSAR技術(shù)處理得到2017年兩伊地震區(qū)域的3種視角視線向同震形變場(chǎng)，將3幅干涉圖重疊區(qū)域裁剪為三維解算區(qū)，利用直接解算法得到研究區(qū)的三維形變場(chǎng)。通過(guò)試驗(yàn)得到以下結(jié)論：①形變區(qū)域的視線向同震形變場(chǎng)中，升軌的最大隆升和最大沉降分別為84.7和-26 cm，降軌的最大隆升和最大沉降分別為58.55和-41.88 cm，且3種視線向的形變場(chǎng)總體形變趨勢(shì)基本一致；②解算得到的三維形變場(chǎng)影響范圍大致為75 km×80 km，分為C區(qū)和D區(qū)兩個(gè)主要形變區(qū)，最大的垂直向、東西向、南北向的形變值分別為84、50、-91 cm，重建的三維形變場(chǎng)的顯示研究區(qū)呈現(xiàn)整體抬升并朝西南向運(yùn)動(dòng)的特點(diǎn)；③綜合地質(zhì)構(gòu)造背景和視線向同震形變場(chǎng)以及三維形變場(chǎng)的形變邊界和趨勢(shì)等因素，本文推測(cè)發(fā)震斷層為扎格羅斯山前斷層。