張家口-宣化盆地現今地殼變形特征與斷層滑動速率反演

高 晨 馬 棟 曹 筠 屈 曼 徐 強

1 河北紅山巨厚沉積與地震災害國家野外科學觀測研究站,河北省邢臺市隆堯縣山口鎮(zhèn),055350 2 河北省地震局,石家莊市槐中路262號,050021

張家口-宣化盆地位于張渤帶和山西斷陷帶交匯部位,盆地內地震活動受張家口斷裂和宣化盆地南緣斷裂控制,區(qū)內曾發(fā)生一系列中強以上地震,但目前該盆地及周邊區(qū)域地震活動特征表現為中等地震活動超常平靜。國內學者利用地質學方法、水準和GPS觀測數據等對張家口斷裂和宣化盆地南緣斷裂的滑動速率和活動性開展相關研究[1-5],但未給出2條斷裂沿斷層面的傾滑速率和閉鎖深度,無法對張家口-宣化盆地中-長期地震危險性分析提供參考。

InSAR技術作為一種新興的大地測量手段,具有全天候、全天時、覆蓋廣、高空間分辨率和高測量精度等優(yōu)勢,可獲取連續(xù)的空間地表形變場[6]。本文使用2018～2020年升降軌模式下3條軌道的Sentinel-1數據,利用PS-InSAR技術結合先驗條件的最小二乘迭代逼近法,獲取張家口-宣化盆地地表三維形變速率場,并以此為基礎反演張家口斷裂和宣化盆地南緣斷裂的傾滑速率和閉鎖深度,為區(qū)域構造模式研究和中-長期地震危險性評價提供依據。

1 數據與處理方法

1.1 InSAR數據處理

本文選用歐空局3條軌道的Sentinel-1A和Sentinel-1B雷達影像數據,成像波段為C波段,TOP模式為干涉寬幅(IW)模式,時間跨度為2018～2020年,影像參數見表1。47軌道、120軌道和142軌道分別選取2019-03-10、2018-11-27和2019-06-03的SLC影像作為主影像,空間基線見圖1。采用ISCE軟件中topsStack數據處理模塊進行影像配準和差分干涉處理,其中影像配準采用綜合考慮主輔影像脈沖響應函數和重疊區(qū)域偏差的增強譜分集方法(ESD),距離向和方位向多視比為4∶1,利用StaMPS軟件實現干涉圖的時序分析[7]。InSAR外界誤差源主要包括大氣誤差、軌道誤差和地形誤差等,本文采用PS-InSAR方法提取雷達波后向散射較強且時序上相對穩(wěn)定的地物目標,數據信噪比較高,大氣誤差較小;C波段雷達衛(wèi)星發(fā)射和接收電磁波主要受對流層影響而產生相位延遲,在PS點相位解纏之前使用GACOS大氣延遲相位模型削弱對流層引起的干涉誤差;分別使用精密軌道星歷數據AUX_POEORB和30 m分辨率的STRM高程模型數據改進軌道精度和模擬地形相位,以提高永久散射體的形變速率精度[7]。最終得到3條軌道的LOS向最優(yōu)形變速率。

圖1 數據覆蓋范圍及干涉數據時空基線分布Fig.1 Data coverage area and spatio-temporal baselines of interferograms

表1 研究區(qū)Sentinel-1衛(wèi)星影像參數Tab.1 Sentinel-1 image parameters of the study area

1.2 三維形變場構建

由于Sentinel-1衛(wèi)星升降軌使用同一衛(wèi)星平臺,成像時各物理參數幾乎未發(fā)生變化,因此影像質量一致性較好。升軌和降軌SAR衛(wèi)星成像的最大差異為雷達電磁波束入射角和衛(wèi)星飛行方向方位角不同,根據該特征,構建地表形變監(jiān)測結果的三維分解模型[8]。InSAR觀測結果為地表真實三維形變在LOS向的投影,同時LOS向形變又可被分解到NS向、EW向和UD向,其幾何關系式為:

dLOS=DUcosθ-DNsinθcos(α-3π/2)-

DEsinθsin(α-3π/2)

(1)

式中,dLOS為LOS向觀測值;DN、DE和DU分別為NS向、EW向和UD向形變值(約定U向、N向和E向為正);θ為雷達側視角;α為衛(wèi)星軌道方位角。

理論上,融合3個及更多平臺或軌道的InSAR形變觀測值就可對地表三維形變場進行構建。本文使用結合先驗條件的最小二乘迭代逼近法,該方法利用已知條件或假設可信已知條件作為多余觀測,迭代使用最小二乘方法求解矛盾方程,以獲取變量近似解[8]。模型方程及約束條件為:

Fn×3×D3×1=dn×1

(2)

Bn×3×D3×1+Wn×1=0

(3)

式中,dn×1為n組視線向觀測值組成的矩陣;Fn×3為投影系數矩陣,由SAR成像幾何條件確定;Wn×1為約束條件;Bn×3為條件系數矩陣。

現有SAR衛(wèi)星多采用近極地軌道飛行和側視成像模式,LOS向形變觀測值在三維分解過程中對NS向形變信息非常不敏感,對UD向形變敏感度最高。可將NS向位移值向極小逼近作為初始約束條件,結合式(2)進行最小二乘解算,根據LOS向形變觀測值中三維分量貢獻程度和解算精度的高低,依次選擇可信度最高的分量作為約束條件,使用最小二乘解算得到研究區(qū)真實可靠的三維形變場。

2 InSAR三維形變場分析

本文利用PS-InSAR技術和結合先驗條件的最小二乘迭代逼近法獲取張家口-宣化盆地及周邊區(qū)域三維形變速率場。

2.1 三維形變結果與檢驗

使用文獻[9-10]中2009～2015年水平向和垂直向GPS速率與基于InSAR觀測數據并結合先驗條件的最小二乘迭代逼近法獲取的三維形變速率進行對比,結果如表2所示。由表可知,二者差值大多在2 mm/a以內,這可能是由局部地表變形在不同監(jiān)測時段的運動速率差所致。統(tǒng)計結果表明,NS向、EW向和UD向的均方根誤差分別為0.76 mm/a、1.64 mm/a和1.51 mm/a,說明本文三維解算結果總體上具有可靠性。

2.2 三維形變場空間變化特征

圖2(a)為張家口-宣化盆地NS向形變速率場,由圖可知,形變速率整體為負,但絕對值較小,盆地所在區(qū)域表現為整體向南的微弱運動趨勢,且盆地內部由西向東N向運動速率逐漸變大。圖2(b)為張家口-宣化盆地EW向形變速率場,其中張家口斷裂北側和宣化盆地南緣斷裂南側區(qū)域整體為正,僅在張家口市區(qū)北部、張家口斷裂與宣化盆地南緣斷裂之間區(qū)域為負,葛峪堡村、南灘村附近區(qū)域為正且絕對值較大,經野外實地調查發(fā)現均為人為干擾。圖2(c)為張家口-宣化盆地UD向形變速率場,其中張家口斷裂中段北側區(qū)域整體呈微弱抬升趨勢,斷裂東段北側區(qū)域整體為沉降運動,沉降量最大區(qū)域為崇禮區(qū)東坪黃金礦區(qū)附近;張家口斷裂和宣化盆地南緣斷裂之間區(qū)域整體為正,且盆地東部形變速率值大于西部,但在陳家堡村和葛峪堡村分別受地下水開采和土方開采等人為干擾導致局部發(fā)生顯著沉降,南灘村附近山體上大面積安裝太陽能板使地表形變表現為顯著抬升;宣化盆地南緣斷裂南部區(qū)域整體表現為沉降趨勢。

張家口斷裂中段為該斷裂主體部分且構造行跡清晰。為進一步分析張家口-宣化盆地、張家口斷裂中段及宣化盆地南緣斷裂的形變速率變化特征,在跨斷裂兩側12 km近場區(qū)域選取3條垂直于斷裂走向的速率剖面,剖面寬1 km,均跨越活動斷裂且盡可能遠離人為干擾較大區(qū)域。如圖2所示,AA1剖面位于張家口斷裂中段(西太平山),長21 km;BB1剖面線位于張家口斷裂中段(人頭山村),長24 km;CC1剖線位于宣化盆地南緣斷裂(龍門坡村東),長24 km。

根據InSAR形變速率剖面,提取跨斷層兩側的形變速率均值,估算得到張家口斷裂中段和宣化盆地南緣斷裂LOS向形變速率(圖3～5)。由圖可知,張家口斷裂中段兩側LOS向跨斷裂形變速率差異為0.19～0.66 mm/a,宣化盆地南緣斷裂兩側速率差異為0.01～0.66 mm/a,在斷層近場區(qū)域均表現出較小的運動差異。根據反演得到的張家口-宣化盆地三維形變速率剖面,分別提取跨斷層兩側NS向、EW向和UD向平均形變速率,由于斷層附近形變點位數量較大,部分點位因誤差過大并不能真實反映斷層實際運動狀態(tài),因此僅保留3倍標準差以內的點位。可以看出,3條剖面在斷層近場處均有微小差異,其中張家口斷裂中段和宣化盆地南緣斷裂跨斷層水平向速率剖面顯示,斷層兩側地表存在微小形變差異。由NS向跨斷層形變速率最佳擬合曲線(圖6(a)～8(a))可知,盆地南部區(qū)域和張家口斷裂中段北部區(qū)域表現為微弱拉張運動,這與張家口-宣化盆地NNW-SSE向拉張性質水平構造應力場特征基本一致[11];由EW向跨斷層形變速率最佳擬合曲線(圖6(b)～8(b))可知,張家口斷裂中段和宣化盆地南緣斷裂兩側地表相對運動存在微小差異;將斷層兩側NS向和EW向形變速率投影到平行斷層走向方向,最佳擬合曲線(圖6(c)～8(c))結果顯示,張家口斷裂中段呈右旋走滑運動特征,速率為0.29～0.57 mm/a,宣化盆地南緣斷裂呈左旋走滑運動特征,速率為1.01 mm/a。圖9為3條剖面跨斷層UD向形變速率最佳擬合曲線,可以看出,張家口-宣化盆地內部和張家口斷裂中段北側區(qū)域構造活動整體表現為抬升,宣化盆地南緣斷裂南側區(qū)域則表現為沉降;圖9(a)顯示張家口斷裂中段在該處兩盤均呈上升趨勢,且南盤相對北盤下降,表現為正斷特征;圖9(b)顯示張家口斷裂中段在該處南盤相對北盤呈微弱上升趨勢,表現為逆斷特征;圖9(c)顯示宣化盆地南緣斷裂北盤上升,南盤下降,表現為正斷特征。

圖3 跨張家口斷裂中段InSAR速率剖面AA1Fig.3 InSAR velocity profiles AA1 across the middle section of Zhangjiakou fault

圖4 跨張家口斷裂中段InSAR速率剖面BB1Fig.4 InSAR velocity profiles BB1 across the middle section of Zhangjiakou fault

圖5 跨宣化盆地南緣斷裂InSAR速率剖面CC1Fig.5 InSAR velocity profiles CC1 across the southern margin fault of Xuanhua basin

圖6 跨張家口斷裂中段形變速率剖面AA1Fig.6 Deformation rate profiles AA1 across the middle section of Zhangjiakou fault

圖7 跨張家口斷裂中段形變速率剖面BB1Fig.7 Deformation rate profiles BB1 across the middle section of Zhangjiakou fault

圖8 跨宣化盆地南緣斷裂西段形變速率剖面CC1Fig.8 Deformation rate profiles CC1 across the west section of southern margin fault of Xuanhua basin

圖9 地殼垂直運動速率剖面Fig.9 Crustal vertical motion rate profiles

3 斷裂滑動速率和閉鎖深度反演

張家口斷裂和宣化盆地南緣斷裂構造運動以垂向為主[11],本文基于跨斷層剖面提取的跨斷層兩側垂向平均形變速率,采用由刃型位錯提出的傾滑斷層震間形變數學表達式進行剖面反演[12]。反演模型數學表達式為:

(4)

(5)

式中,ux為水平面垂直斷層方向位移,uz為垂直地面方向位移,s為沿斷層面傾滑位移,d為斷層閉鎖深度,δ為斷層傾角,x為觀測點到斷層的垂直距離。

由于式(5)高度非線性,無法確定參數的近似值,直接使用最小二乘法進行反演容易使方程產生奇異并無解。為解決該問題,本文使用擬牛頓法和通用全局優(yōu)化法來確定傾滑速率、閉鎖深度和斷層傾角3個參數的最優(yōu)估計結果。在反演過程中,根據區(qū)域歷史地震的震源深度,將斷層閉鎖深度范圍設置為1～30 km;根據已有的斷層結構研究結果[1],將斷層傾角范圍設置為55°～65°;研究區(qū)內斷裂平均垂直滑動速率小于1 mm/a,因此將斷層傾滑速率搜索范圍設置為0.1～1 mm/a[2]。利用上述方法反演得到的張家口斷裂中段和宣化盆地南緣斷裂的斷層參數如表3所示,斷層閉鎖深度分布如圖10所示。

圖10 斷層閉鎖深度分布Fig.10 Distribution of fault locking depth

4 討論與分析

4.1 與地質、水準和GPS結果進行對比

基于由刃型位錯提出的傾滑斷層震間形變數學表達式,利用2018～2020年獲取的InSAR垂直形變速率場數據進行擬合,反演得到張家口斷裂中段和宣化盆地南緣斷裂的傾滑速率分別為0.42～0.79 mm/a和0.51 mm/a,與前人采用地質學方法和水準觀測得到的斷層垂直滑動速率結果[1-3]基本一致。

使用剖面投影法分別獲得斷層兩側水平向運動差異特征,其中張家口斷裂中段和宣化盆地南緣斷裂的水平向滑動速率分別為0.29～0.57 mm/a和1.01 mm/a,與陳長云等[4]利用1999～2007年GPS觀測數據計算得到的張家口斷裂左旋走滑速率結果基本一致,而本文計算的張家口斷裂中段水平向走滑運動特征與已有研究結果相反[2,4]。張家口-宣化盆地所處的華北地區(qū)在印度洋板塊和太平洋板塊的共同作用下,現代構造應力場表現為以NEE-EW向擠壓為主,兼NNW向拉張,使得NW向和近EW向斷裂表現為左旋錯動,但區(qū)域構造應力場在演化過程中會出現應力調整作用[13]。

洞體應變觀測是監(jiān)測地殼運動和變形的重要手段之一,主要用來測量地殼表面兩點之間的應變量(觀測曲線上升為拉張,下降為擠壓),進而揭示彈性應變能在局部區(qū)域逐步積累和應力、應變逐步演變特征。本文利用位于張家口斷裂中段北側附近的張家口地震臺洞體應變儀2017～2020年觀測數據,該儀器NS分量觀測數據于2019-07-13～2020-04-01出現與往年同期不同的區(qū)域應力異常現象,且持續(xù)時間較長,局部區(qū)域構造應力場在NS向呈擠壓變化趨勢,并持續(xù)至2020-05-30(圖11)。張家口臺洞體應變儀NS分量觀測數據反映出張家口-宣化盆地北部局部區(qū)域構造應力場在2019～2020年發(fā)生變化,導致NWW向張家口斷裂中段在2018～2020年沿斷層走向方向表現出微弱的右旋走滑特征。

圖11 洞體應變觀測NS分量年變對比Fig.11 Comparison of the annual variation of NS component of cave strain

4.2 斷裂滑動速率、閉鎖深度與地震危險性

張家口斷裂中段和宣化盆地南緣斷裂的閉鎖深度分別為9.05～12.98 km和8.38 km,與張希等[14]使用負位錯模型反演得到的鄂爾多斯塊體及其周緣斷層的閉鎖深度結果基本一致,說明2條斷裂均具有發(fā)生中強以上地震的應變能積累背景條件。

5 結 語

本文利用PS-InSAR技術分別選用3條軌道(升軌142、降軌47和降軌120)的Sentinel-1數據獲取張家口-宣化盆地及周邊區(qū)域地表形變速率場,并利用結合先驗條件的最小二乘迭代逼近法計算盆地及周邊區(qū)域三維地表形變場。分別使用剖面投影法和由刃型位錯提出的傾滑斷層震間形變數學表達式計算和分析張家口斷裂中段及宣化盆地南緣斷裂現今水平向滑動速率、傾滑速率及斷層閉鎖深度,得出以下結論:

1)張家口斷裂中段和宣化盆地南緣斷裂的水平向滑動速率分別為0.29～0.57 mm/a和1.01 mm/a,傾滑速率分別為0.42～0.79 mm/a和0.51 mm/a;張家口斷裂中段表現為右旋走滑兼正(逆)斷傾滑活動特征,宣化盆地南緣斷裂表現為正斷傾滑兼左旋走滑活動特征。

2)張家口斷裂中段和宣化盆地南緣斷裂的斷層閉鎖深度分別為9.05～12.98 km和8.38 km,二者相差較小,均具有發(fā)生中強以上地震的應變能積累背景條件。