高 晨 馬 棟 曹 筠 屈 曼 徐 強(qiáng)
1 河北紅山巨厚沉積與地震災(zāi)害國家野外科學(xué)觀測研究站,河北省邢臺(tái)市隆堯縣山口鎮(zhèn),055350 2 河北省地震局,石家莊市槐中路262號(hào),050021
張家口-宣化盆地位于張渤帶和山西斷陷帶交匯部位,盆地內(nèi)地震活動(dòng)受張家口斷裂和宣化盆地南緣斷裂控制,區(qū)內(nèi)曾發(fā)生一系列中強(qiáng)以上地震,但目前該盆地及周邊區(qū)域地震活動(dòng)特征表現(xiàn)為中等地震活動(dòng)超常平靜。國內(nèi)學(xué)者利用地質(zhì)學(xué)方法、水準(zhǔn)和GPS觀測數(shù)據(jù)等對張家口斷裂和宣化盆地南緣斷裂的滑動(dòng)速率和活動(dòng)性開展相關(guān)研究[1-5],但未給出2條斷裂沿?cái)鄬用娴膬A滑速率和閉鎖深度,無法對張家口-宣化盆地中-長期地震危險(xiǎn)性分析提供參考。
InSAR技術(shù)作為一種新興的大地測量手段,具有全天候、全天時(shí)、覆蓋廣、高空間分辨率和高測量精度等優(yōu)勢,可獲取連續(xù)的空間地表形變場[6]。本文使用2018~2020年升降軌模式下3條軌道的Sentinel-1數(shù)據(jù),利用PS-InSAR技術(shù)結(jié)合先驗(yàn)條件的最小二乘迭代逼近法,獲取張家口-宣化盆地地表三維形變速率場,并以此為基礎(chǔ)反演張家口斷裂和宣化盆地南緣斷裂的傾滑速率和閉鎖深度,為區(qū)域構(gòu)造模式研究和中-長期地震危險(xiǎn)性評價(jià)提供依據(jù)。
本文選用歐空局3條軌道的Sentinel-1A和Sentinel-1B雷達(dá)影像數(shù)據(jù),成像波段為C波段,TOP模式為干涉寬幅(IW)模式,時(shí)間跨度為2018~2020年,影像參數(shù)見表1。47軌道、120軌道和142軌道分別選取2019-03-10、2018-11-27和2019-06-03的SLC影像作為主影像,空間基線見圖1。采用ISCE軟件中topsStack數(shù)據(jù)處理模塊進(jìn)行影像配準(zhǔn)和差分干涉處理,其中影像配準(zhǔn)采用綜合考慮主輔影像脈沖響應(yīng)函數(shù)和重疊區(qū)域偏差的增強(qiáng)譜分集方法(ESD),距離向和方位向多視比為4∶1,利用StaMPS軟件實(shí)現(xiàn)干涉圖的時(shí)序分析[7]。InSAR外界誤差源主要包括大氣誤差、軌道誤差和地形誤差等,本文采用PS-InSAR方法提取雷達(dá)波后向散射較強(qiáng)且時(shí)序上相對穩(wěn)定的地物目標(biāo),數(shù)據(jù)信噪比較高,大氣誤差較小;C波段雷達(dá)衛(wèi)星發(fā)射和接收電磁波主要受對流層影響而產(chǎn)生相位延遲,在PS點(diǎn)相位解纏之前使用GACOS大氣延遲相位模型削弱對流層引起的干涉誤差;分別使用精密軌道星歷數(shù)據(jù)AUX_POEORB和30 m分辨率的STRM高程模型數(shù)據(jù)改進(jìn)軌道精度和模擬地形相位,以提高永久散射體的形變速率精度[7]。最終得到3條軌道的LOS向最優(yōu)形變速率。
圖1 數(shù)據(jù)覆蓋范圍及干涉數(shù)據(jù)時(shí)空基線分布Fig.1 Data coverage area and spatio-temporal baselines of interferograms
表1 研究區(qū)Sentinel-1衛(wèi)星影像參數(shù)Tab.1 Sentinel-1 image parameters of the study area
由于Sentinel-1衛(wèi)星升降軌使用同一衛(wèi)星平臺(tái),成像時(shí)各物理參數(shù)幾乎未發(fā)生變化,因此影像質(zhì)量一致性較好。升軌和降軌SAR衛(wèi)星成像的最大差異為雷達(dá)電磁波束入射角和衛(wèi)星飛行方向方位角不同,根據(jù)該特征,構(gòu)建地表形變監(jiān)測結(jié)果的三維分解模型[8]。InSAR觀測結(jié)果為地表真實(shí)三維形變在LOS向的投影,同時(shí)LOS向形變又可被分解到NS向、EW向和UD向,其幾何關(guān)系式為:
dLOS=DUcosθ-DNsinθcos(α-3π/2)-
DEsinθsin(α-3π/2)
(1)
式中,dLOS為LOS向觀測值;DN、DE和DU分別為NS向、EW向和UD向形變值(約定U向、N向和E向?yàn)檎?;θ為雷達(dá)側(cè)視角;α為衛(wèi)星軌道方位角。
理論上,融合3個(gè)及更多平臺(tái)或軌道的InSAR形變觀測值就可對地表三維形變場進(jìn)行構(gòu)建。本文使用結(jié)合先驗(yàn)條件的最小二乘迭代逼近法,該方法利用已知條件或假設(shè)可信已知條件作為多余觀測,迭代使用最小二乘方法求解矛盾方程,以獲取變量近似解[8]。模型方程及約束條件為:
Fn×3×D3×1=dn×1
(2)
Bn×3×D3×1+Wn×1=0
(3)
式中,dn×1為n組視線向觀測值組成的矩陣;Fn×3為投影系數(shù)矩陣,由SAR成像幾何條件確定;Wn×1為約束條件;Bn×3為條件系數(shù)矩陣。
現(xiàn)有SAR衛(wèi)星多采用近極地軌道飛行和側(cè)視成像模式,LOS向形變觀測值在三維分解過程中對NS向形變信息非常不敏感,對UD向形變敏感度最高??蓪S向位移值向極小逼近作為初始約束條件,結(jié)合式(2)進(jìn)行最小二乘解算,根據(jù)LOS向形變觀測值中三維分量貢獻(xiàn)程度和解算精度的高低,依次選擇可信度最高的分量作為約束條件,使用最小二乘解算得到研究區(qū)真實(shí)可靠的三維形變場。
本文利用PS-InSAR技術(shù)和結(jié)合先驗(yàn)條件的最小二乘迭代逼近法獲取張家口-宣化盆地及周邊區(qū)域三維形變速率場。
使用文獻(xiàn)[9-10]中2009~2015年水平向和垂直向GPS速率與基于InSAR觀測數(shù)據(jù)并結(jié)合先驗(yàn)條件的最小二乘迭代逼近法獲取的三維形變速率進(jìn)行對比,結(jié)果如表2所示。由表可知,二者差值大多在2 mm/a以內(nèi),這可能是由局部地表變形在不同監(jiān)測時(shí)段的運(yùn)動(dòng)速率差所致。統(tǒng)計(jì)結(jié)果表明,NS向、EW向和UD向的均方根誤差分別為0.76 mm/a、1.64 mm/a和1.51 mm/a,說明本文三維解算結(jié)果總體上具有可靠性。
圖2(a)為張家口-宣化盆地NS向形變速率場,由圖可知,形變速率整體為負(fù),但絕對值較小,盆地所在區(qū)域表現(xiàn)為整體向南的微弱運(yùn)動(dòng)趨勢,且盆地內(nèi)部由西向東N向運(yùn)動(dòng)速率逐漸變大。圖2(b)為張家口-宣化盆地EW向形變速率場,其中張家口斷裂北側(cè)和宣化盆地南緣斷裂南側(cè)區(qū)域整體為正,僅在張家口市區(qū)北部、張家口斷裂與宣化盆地南緣斷裂之間區(qū)域?yàn)樨?fù),葛峪堡村、南灘村附近區(qū)域?yàn)檎医^對值較大,經(jīng)野外實(shí)地調(diào)查發(fā)現(xiàn)均為人為干擾。圖2(c)為張家口-宣化盆地UD向形變速率場,其中張家口斷裂中段北側(cè)區(qū)域整體呈微弱抬升趨勢,斷裂東段北側(cè)區(qū)域整體為沉降運(yùn)動(dòng),沉降量最大區(qū)域?yàn)槌缍Y區(qū)東坪黃金礦區(qū)附近;張家口斷裂和宣化盆地南緣斷裂之間區(qū)域整體為正,且盆地東部形變速率值大于西部,但在陳家堡村和葛峪堡村分別受地下水開采和土方開采等人為干擾導(dǎo)致局部發(fā)生顯著沉降,南灘村附近山體上大面積安裝太陽能板使地表形變表現(xiàn)為顯著抬升;宣化盆地南緣斷裂南部區(qū)域整體表現(xiàn)為沉降趨勢。
張家口斷裂中段為該斷裂主體部分且構(gòu)造行跡清晰。為進(jìn)一步分析張家口-宣化盆地、張家口斷裂中段及宣化盆地南緣斷裂的形變速率變化特征,在跨斷裂兩側(cè)12 km近場區(qū)域選取3條垂直于斷裂走向的速率剖面,剖面寬1 km,均跨越活動(dòng)斷裂且盡可能遠(yuǎn)離人為干擾較大區(qū)域。如圖2所示,AA1剖面位于張家口斷裂中段(西太平山),長21 km;BB1剖面線位于張家口斷裂中段(人頭山村),長24 km;CC1剖線位于宣化盆地南緣斷裂(龍門坡村東),長24 km。
根據(jù)InSAR形變速率剖面,提取跨斷層兩側(cè)的形變速率均值,估算得到張家口斷裂中段和宣化盆地南緣斷裂LOS向形變速率(圖3~5)。由圖可知,張家口斷裂中段兩側(cè)LOS向跨斷裂形變速率差異為0.19~0.66 mm/a,宣化盆地南緣斷裂兩側(cè)速率差異為0.01~0.66 mm/a,在斷層近場區(qū)域均表現(xiàn)出較小的運(yùn)動(dòng)差異。根據(jù)反演得到的張家口-宣化盆地三維形變速率剖面,分別提取跨斷層兩側(cè)NS向、EW向和UD向平均形變速率,由于斷層附近形變點(diǎn)位數(shù)量較大,部分點(diǎn)位因誤差過大并不能真實(shí)反映斷層實(shí)際運(yùn)動(dòng)狀態(tài),因此僅保留3倍標(biāo)準(zhǔn)差以內(nèi)的點(diǎn)位。可以看出,3條剖面在斷層近場處均有微小差異,其中張家口斷裂中段和宣化盆地南緣斷裂跨斷層水平向速率剖面顯示,斷層兩側(cè)地表存在微小形變差異。由NS向跨斷層形變速率最佳擬合曲線(圖6(a)~8(a))可知,盆地南部區(qū)域和張家口斷裂中段北部區(qū)域表現(xiàn)為微弱拉張運(yùn)動(dòng),這與張家口-宣化盆地NNW-SSE向拉張性質(zhì)水平構(gòu)造應(yīng)力場特征基本一致[11];由EW向跨斷層形變速率最佳擬合曲線(圖6(b)~8(b))可知,張家口斷裂中段和宣化盆地南緣斷裂兩側(cè)地表相對運(yùn)動(dòng)存在微小差異;將斷層兩側(cè)NS向和EW向形變速率投影到平行斷層走向方向,最佳擬合曲線(圖6(c)~8(c))結(jié)果顯示,張家口斷裂中段呈右旋走滑運(yùn)動(dòng)特征,速率為0.29~0.57 mm/a,宣化盆地南緣斷裂呈左旋走滑運(yùn)動(dòng)特征,速率為1.01 mm/a。圖9為3條剖面跨斷層UD向形變速率最佳擬合曲線,可以看出,張家口-宣化盆地內(nèi)部和張家口斷裂中段北側(cè)區(qū)域構(gòu)造活動(dòng)整體表現(xiàn)為抬升,宣化盆地南緣斷裂南側(cè)區(qū)域則表現(xiàn)為沉降;圖9(a)顯示張家口斷裂中段在該處兩盤均呈上升趨勢,且南盤相對北盤下降,表現(xiàn)為正斷特征;圖9(b)顯示張家口斷裂中段在該處南盤相對北盤呈微弱上升趨勢,表現(xiàn)為逆斷特征;圖9(c)顯示宣化盆地南緣斷裂北盤上升,南盤下降,表現(xiàn)為正斷特征。
圖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 地殼垂直運(yùn)動(dòng)速率剖面Fig.9 Crustal vertical motion rate profiles
張家口斷裂和宣化盆地南緣斷裂構(gòu)造運(yùn)動(dòng)以垂向?yàn)橹鱗11],本文基于跨斷層剖面提取的跨斷層兩側(cè)垂向平均形變速率,采用由刃型位錯(cuò)提出的傾滑斷層震間形變數(shù)學(xué)表達(dá)式進(jìn)行剖面反演[12]。反演模型數(shù)學(xué)表達(dá)式為:
(4)
(5)
式中,ux為水平面垂直斷層方向位移,uz為垂直地面方向位移,s為沿?cái)鄬用鎯A滑位移,d為斷層閉鎖深度,δ為斷層傾角,x為觀測點(diǎn)到斷層的垂直距離。
由于式(5)高度非線性,無法確定參數(shù)的近似值,直接使用最小二乘法進(jìn)行反演容易使方程產(chǎn)生奇異并無解。為解決該問題,本文使用擬牛頓法和通用全局優(yōu)化法來確定傾滑速率、閉鎖深度和斷層傾角3個(gè)參數(shù)的最優(yōu)估計(jì)結(jié)果。在反演過程中,根據(jù)區(qū)域歷史地震的震源深度,將斷層閉鎖深度范圍設(shè)置為1~30 km;根據(jù)已有的斷層結(jié)構(gòu)研究結(jié)果[1],將斷層傾角范圍設(shè)置為55°~65°;研究區(qū)內(nèi)斷裂平均垂直滑動(dòng)速率小于1 mm/a,因此將斷層傾滑速率搜索范圍設(shè)置為0.1~1 mm/a[2]。利用上述方法反演得到的張家口斷裂中段和宣化盆地南緣斷裂的斷層參數(shù)如表3所示,斷層閉鎖深度分布如圖10所示。
圖10 斷層閉鎖深度分布Fig.10 Distribution of fault locking depth
基于由刃型位錯(cuò)提出的傾滑斷層震間形變數(shù)學(xué)表達(dá)式,利用2018~2020年獲取的InSAR垂直形變速率場數(shù)據(jù)進(jìn)行擬合,反演得到張家口斷裂中段和宣化盆地南緣斷裂的傾滑速率分別為0.42~0.79 mm/a和0.51 mm/a,與前人采用地質(zhì)學(xué)方法和水準(zhǔn)觀測得到的斷層垂直滑動(dòng)速率結(jié)果[1-3]基本一致。
使用剖面投影法分別獲得斷層兩側(cè)水平向運(yùn)動(dòng)差異特征,其中張家口斷裂中段和宣化盆地南緣斷裂的水平向滑動(dòng)速率分別為0.29~0.57 mm/a和1.01 mm/a,與陳長云等[4]利用1999~2007年GPS觀測數(shù)據(jù)計(jì)算得到的張家口斷裂左旋走滑速率結(jié)果基本一致,而本文計(jì)算的張家口斷裂中段水平向走滑運(yùn)動(dòng)特征與已有研究結(jié)果相反[2,4]。張家口-宣化盆地所處的華北地區(qū)在印度洋板塊和太平洋板塊的共同作用下,現(xiàn)代構(gòu)造應(yīng)力場表現(xiàn)為以NEE-EW向擠壓為主,兼NNW向拉張,使得NW向和近EW向斷裂表現(xiàn)為左旋錯(cuò)動(dòng),但區(qū)域構(gòu)造應(yīng)力場在演化過程中會(huì)出現(xiàn)應(yīng)力調(diào)整作用[13]。
洞體應(yīng)變觀測是監(jiān)測地殼運(yùn)動(dòng)和變形的重要手段之一,主要用來測量地殼表面兩點(diǎn)之間的應(yīng)變量(觀測曲線上升為拉張,下降為擠壓),進(jìn)而揭示彈性應(yīng)變能在局部區(qū)域逐步積累和應(yīng)力、應(yīng)變逐步演變特征。本文利用位于張家口斷裂中段北側(cè)附近的張家口地震臺(tái)洞體應(yīng)變儀2017~2020年觀測數(shù)據(jù),該儀器NS分量觀測數(shù)據(jù)于2019-07-13~2020-04-01出現(xiàn)與往年同期不同的區(qū)域應(yīng)力異?,F(xiàn)象,且持續(xù)時(shí)間較長,局部區(qū)域構(gòu)造應(yīng)力場在NS向呈擠壓變化趨勢,并持續(xù)至2020-05-30(圖11)。張家口臺(tái)洞體應(yīng)變儀NS分量觀測數(shù)據(jù)反映出張家口-宣化盆地北部局部區(qū)域構(gòu)造應(yīng)力場在2019~2020年發(fā)生變化,導(dǎo)致NWW向張家口斷裂中段在2018~2020年沿?cái)鄬幼呦蚍较虮憩F(xiàn)出微弱的右旋走滑特征。
圖11 洞體應(yīng)變觀測NS分量年變對比Fig.11 Comparison of the annual variation of NS component of cave strain
張家口斷裂中段和宣化盆地南緣斷裂的閉鎖深度分別為9.05~12.98 km和8.38 km,與張希等[14]使用負(fù)位錯(cuò)模型反演得到的鄂爾多斯塊體及其周緣斷層的閉鎖深度結(jié)果基本一致,說明2條斷裂均具有發(fā)生中強(qiáng)以上地震的應(yīng)變能積累背景條件。
本文利用PS-InSAR技術(shù)分別選用3條軌道(升軌142、降軌47和降軌120)的Sentinel-1數(shù)據(jù)獲取張家口-宣化盆地及周邊區(qū)域地表形變速率場,并利用結(jié)合先驗(yàn)條件的最小二乘迭代逼近法計(jì)算盆地及周邊區(qū)域三維地表形變場。分別使用剖面投影法和由刃型位錯(cuò)提出的傾滑斷層震間形變數(shù)學(xué)表達(dá)式計(jì)算和分析張家口斷裂中段及宣化盆地南緣斷裂現(xiàn)今水平向滑動(dòng)速率、傾滑速率及斷層閉鎖深度,得出以下結(jié)論:
1)張家口斷裂中段和宣化盆地南緣斷裂的水平向滑動(dòng)速率分別為0.29~0.57 mm/a和1.01 mm/a,傾滑速率分別為0.42~0.79 mm/a和0.51 mm/a;張家口斷裂中段表現(xiàn)為右旋走滑兼正(逆)斷傾滑活動(dòng)特征,宣化盆地南緣斷裂表現(xiàn)為正斷傾滑兼左旋走滑活動(dòng)特征。
2)張家口斷裂中段和宣化盆地南緣斷裂的斷層閉鎖深度分別為9.05~12.98 km和8.38 km,二者相差較小,均具有發(fā)生中強(qiáng)以上地震的應(yīng)變能積累背景條件。