多視角InSAR數據解算2017兩伊地震三維同震形變場

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

(長安大學地質工程與測繪學院，陜西 西安 710054)

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

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

1 地質構造背景

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

圖1 研究區震中與地質構造背景

2 三維同震形變場解算方法

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

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

(1)

圖2 衛星軌道方向

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

(2)

令

(3)

3 視線向同震形變場提取與分析

3.1 數據選取與處理流程

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

表1 Sentinel1-A(來源ESA)衛星數據基本信息

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

3.2 視線向同震形變場提取分析

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

圖3 升軌-A72視線向形變

圖4 降軌-D79視線向形變

圖5 降軌-D6視線向形變

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

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

4 三維形變結果分析

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

圖7 三維同震形變場

5 結 語

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