張青鎖，于松暉，徐郅杰

1.河南省地質(zhì)環(huán)境監(jiān)測(cè)院 地質(zhì)災(zāi)害防治技術(shù)中心，鄭州 450016；2.河南省地質(zhì)災(zāi)害防治重點(diǎn)實(shí)驗(yàn)室，鄭州 450016

0 引言

中國(guó)上海最早在20世紀(jì)20年代發(fā)現(xiàn)地面沉降，50年代末天津市、西安市等地發(fā)現(xiàn)地面沉降；到了70年代，長(zhǎng)江三角洲主要城市、河北滄州市、山西太原市、安徽阜陽(yáng)市等相繼出現(xiàn)地面沉降；截至目前為止,在全國(guó)發(fā)現(xiàn)地面沉降現(xiàn)象的城市和地區(qū)至少有96個(gè)[1--5]。河南省地面沉降發(fā)現(xiàn)于20世紀(jì)70年代，主要分布在許昌、開封、濮陽(yáng)、商丘等城市。進(jìn)入80年代以后，鄭州市、洛陽(yáng)市以及河南北部平原區(qū)的安陽(yáng)市、鶴壁市、濮陽(yáng)市、新鄉(xiāng)市等均有不同程度地面沉降[6--11]。

建設(shè)水準(zhǔn)點(diǎn)和基巖標(biāo)及分層標(biāo)等觀測(cè)點(diǎn)并進(jìn)行長(zhǎng)期監(jiān)測(cè)是傳統(tǒng)測(cè)量地面沉降的主要工作方法。這些監(jiān)測(cè)方法需要對(duì)觀測(cè)點(diǎn)進(jìn)行長(zhǎng)期維護(hù)。測(cè)量過(guò)程需要較長(zhǎng)的野外作業(yè)時(shí)間，并且需要耗費(fèi)數(shù)量較大的人財(cái)物資源。且只能監(jiān)測(cè)有限離散點(diǎn)，很難達(dá)到在宏觀上掌握區(qū)域地面沉降規(guī)律的效果[12--13]。

合成孔徑雷達(dá)干涉測(cè)量(Interferometric Synthetic Aperture Radar，InSAR)技術(shù)是空間大地測(cè)量技術(shù)中發(fā)展最快的技術(shù)之一，在各行各業(yè)得到了廣泛的應(yīng)用。該技術(shù)實(shí)現(xiàn)了大范圍、全天候、高精度地表形變監(jiān)測(cè)，是傳統(tǒng)測(cè)量方法的有益補(bǔ)充，并且具有非接觸的特點(diǎn)，在監(jiān)測(cè)地面形變方面具有很大優(yōu)勢(shì)，目前已成功應(yīng)用于地面沉降觀測(cè)中[14--16]。短基線集(Small Baseline Subset，SBAS)干涉測(cè)量技術(shù)采用多主影像的干涉對(duì)，以高相干點(diǎn)為基礎(chǔ)恢復(fù)地表時(shí)間序列形變信息。SBAS--InSAR技術(shù)降低了SAR數(shù)據(jù)需求量，提高了運(yùn)算效率[17--19]。當(dāng)前SBAS--InSAR方法應(yīng)用中人工交互較多，操作人員的經(jīng)驗(yàn)直接影響形變監(jiān)測(cè)的效果[20]。

本次研究利用2014—2016年間RADASAT--2雷達(dá)衛(wèi)星數(shù)據(jù)，采用短基線集技術(shù)對(duì)河南北部平原區(qū)進(jìn)行了大時(shí)空維度上的地表形變監(jiān)測(cè)。對(duì)研究區(qū)形變監(jiān)測(cè)成果進(jìn)行分析，劃分地面沉降嚴(yán)重程度，圈定地面沉降重點(diǎn)區(qū)，分析沉降成因，為制定地面沉降的防控措施、指導(dǎo)后續(xù)監(jiān)測(cè)提供了基礎(chǔ)數(shù)據(jù)依據(jù)。

1 數(shù)據(jù)處理

河南北部平原區(qū)位于河南省的東北部，西邊緊靠太行山，南部以黃河為界，是華北平原的組成部分。范圍涉及31個(gè)縣(市、區(qū))，面積約21 641.00 km2。區(qū)域氣候干旱，降雨量由西向東呈帶狀分布，西部較大，東部較小。地表水系分屬黃河流域和海河流域，地貌類型分為殘丘崗地、沖積--洪積平原、沖積平原3類。在區(qū)域上屬華北地層區(qū)，區(qū)內(nèi)廣泛分布第四系，前第四系僅在局部出露。該區(qū)位于中朝準(zhǔn)地臺(tái)南部，區(qū)域新構(gòu)造活動(dòng)主要表現(xiàn)為升降運(yùn)動(dòng)，局部表現(xiàn)為斷裂活動(dòng)。第四紀(jì)以來(lái)構(gòu)造活動(dòng)在喜山運(yùn)動(dòng)影響下不斷增強(qiáng)，沉降幅度由山前至平原呈階梯狀增大，太行山前小于東部平原，凸起區(qū)小于凹陷區(qū)。新構(gòu)造時(shí)期斷裂通常發(fā)育在盆地邊緣的新近系中，多呈高角度斷層，斷距一般<100 m，個(gè)別達(dá)150 m±。構(gòu)造線以北北東、北東向?yàn)橹鳎渲斜睎|向斷裂呈階梯狀分布。為研究河南北部平原區(qū)地面沉降特征，利用短基線集技術(shù)獲取地表形變數(shù)據(jù)是有效手段。

1.1 數(shù)據(jù)源

本次數(shù)據(jù)源選擇RADASAT--2數(shù)據(jù)，覆蓋研究區(qū)大部分范圍需要3個(gè)條帶(條帶1、條帶2、條帶3)(圖1)，實(shí)際覆蓋面積21 192.80 km2。其中每個(gè)條帶獲取2014年2月至2016年11月之間的圖像20景，共計(jì)60景圖像。

圖1 SAR數(shù)據(jù)覆蓋范圍

1.2 方法原理

短基線集(SBAS)技術(shù)利用獲取的合成孔徑雷達(dá)影像數(shù)據(jù)集，設(shè)置時(shí)間和空間基線閾值，經(jīng)過(guò)自由組合以后，生成若干個(gè)子集合。在各個(gè)子集合內(nèi)部通過(guò)最小二乘法獲取子集合的形變數(shù)據(jù)，在各個(gè)子集合之間，使用奇異值分解(SVD)方法聯(lián)合求解，獲得區(qū)域的地表形變信息數(shù)據(jù)。

假定獲取同一區(qū)域在(t0，t1，… ，tN)時(shí)刻的N+1幅不同的SAR圖像，依據(jù)不同的干涉條件，經(jīng)過(guò)自由組合形成M組干涉圖，可知

(N+1)/2≤M≤[N(N+1)]/2

(1)

如果想要得到第j幅干涉圖，需要獲得tA，tB兩個(gè)時(shí)刻的SAR圖像，消除平地效應(yīng)和地形相位以后，在x處的差分干涉相位為:

(2)

式中：λ為雷達(dá)的工作波長(zhǎng)；d(tB，x)為tB時(shí)刻累計(jì)形變；d(tA，x)為tA時(shí)刻累計(jì)形變；d(t0，x)=0；φ(tB，x)為d(tB，x)引起的形變相位；φ(tA，x)是d(tA，x)引起的形變相位。運(yùn)用線性模型估算所取得的N幅圖像的形變量為:

AΦ=ΔΦ

(3)

式中：A為系數(shù)矩陣，其中行是干涉圖，共M行，列是時(shí)刻SAR圖像，共N列；用+1表示主圖像所在列，用-1表示輔圖像所在列，用0表示其余列；Φ是相位矩陣，表示目標(biāo)點(diǎn)在N個(gè)時(shí)刻M幅圖像上的形變相位；ΔΦ為相位值矩陣，表示目標(biāo)點(diǎn)在N個(gè)時(shí)刻在M幅差分干涉圖上相位值。假如M≥N，而且A秩為N，通過(guò)最小二乘法可以得到形變相位

Φ=(ATA)-1AT·ΔΦ

(4)

若要想得到上述結(jié)果，需要滿足獲得的所有SAR圖像都必須在1個(gè)短基線集范圍內(nèi)。但在實(shí)際中能真正滿足這樣條件的可能性極小，往往由于單個(gè)集合內(nèi)時(shí)間樣品不足的原因，造成ATA是個(gè)奇異矩陣，使得方程存在無(wú)數(shù)解。但是可以通過(guò)連接多個(gè)短基線集選擇SVD的求解方法，得到最小范數(shù)解。

1.3 處理過(guò)程

(1)對(duì)每個(gè)條帶的20景SAR數(shù)據(jù)按照式(1)計(jì)算，獲得每個(gè)條帶干涉像對(duì)配對(duì)數(shù)量范圍為11～190對(duì)。根據(jù)3個(gè)條帶獲得SAR數(shù)據(jù)時(shí)間基線分布不同，結(jié)合每個(gè)條帶的地形和植被發(fā)育狀況，分別設(shè)置3個(gè)條帶的時(shí)間基線閾值(表1)。空間基線閾值按照SAR數(shù)據(jù)空間臨界基線的10%取250 m。經(jīng)過(guò)最優(yōu)組合方式配對(duì)，獲得3個(gè)條帶的干涉相干組合(圖2)。

表1 干涉組合統(tǒng)計(jì)表

(2)在干涉的處理過(guò)程中采用8×2的濾波窗口，對(duì)SAR圖像進(jìn)行地理坐標(biāo)匹配，進(jìn)行兩幅復(fù)共軛相乘，獲得干涉條紋圖。為獲得較好的解纏相位效果，消除處理方位向上平地效應(yīng)造成的條紋(圖3)。

(3)選擇使用最小費(fèi)用流的方法開展相位解纏；為了降低噪聲誤差對(duì)解纏結(jié)果的影響，基于本次圖像的分辨率和當(dāng)?shù)氐匚锏念愋吞卣?，只?duì)相干系數(shù)>0.4的柵格點(diǎn)相位解纏(圖4)。

a.條帶1;b.條帶2;c.條帶3。

圖3 去除平地效應(yīng)后的干涉條紋圖(a)和相干圖(b)

(4)相干圖經(jīng)過(guò)平均處理以后，獲取平均相干系數(shù)圖，再設(shè)定合適的相干系數(shù)閾值就可以獲得高相干的目標(biāo)點(diǎn)。本次結(jié)合區(qū)域范圍大小、時(shí)間基線長(zhǎng)短和地物類型特征，相干系數(shù)的閾值設(shè)定為0.6；獲得可以用于形變反演的高相干目標(biāo)點(diǎn)約207萬(wàn)個(gè)(圖5)。

(5)根據(jù)選取的高相干點(diǎn)，解纏后獲得觀測(cè)范圍內(nèi)的干涉條紋圖數(shù)據(jù)，建立對(duì)應(yīng)的相位模型，選擇SVD方法求解就能夠獲得形變速率(圖6)。

(6)通過(guò)控制點(diǎn)把不同圖像的形變數(shù)據(jù)統(tǒng)一到相同的坐標(biāo)系統(tǒng)?？刂泣c(diǎn)的地理位置信息存儲(chǔ)于圖像頭文件中，根據(jù)控制點(diǎn)坐標(biāo)信息對(duì)形變速率圖進(jìn)行校正處理，拼接3個(gè)條帶解譯結(jié)果并進(jìn)行裁剪，獲得研究區(qū)地理編碼后的形變速率圖(圖7)。

a.條帶1;b.條帶2;c.條帶3。

a.條帶1;b.條帶2;c.條帶3。

a.條帶1;b.條帶2;c.條帶3。

圖7 地理編碼后的形變速率圖

2 地面沉降分析

2.1 地面沉降空間分布特征

研究區(qū)采用InSAR解譯地面沉降，監(jiān)測(cè)面積21 192.80 km2。利用地理編碼后的形變速率圖制作地面沉降速率等值線圖(圖8)，2014—2016年研究區(qū)均發(fā)生了地面沉降，大部分區(qū)域沉降速率10～30 mm/a，最大的沉降速率達(dá)114.85 mm/a。參照地面沉降嚴(yán)重程度的劃分標(biāo)準(zhǔn)[21]，依據(jù)地面沉降速率等值線劃分研究區(qū)地面沉降嚴(yán)重程度(表2)。分區(qū)結(jié)果顯示：區(qū)內(nèi)僅有個(gè)別點(diǎn)沉降速率>80 mm/a，沒(méi)有圈出地面沉降嚴(yán)重程度高區(qū)；分布最為廣泛的是地面沉降嚴(yán)重程度較低區(qū)，面積19 609.34 km2，占監(jiān)測(cè)區(qū)92.53%。結(jié)合不同嚴(yán)重區(qū)域大小、人類活動(dòng)的分布特征、集中連片等原則，綜合確定8個(gè)重點(diǎn)沉降區(qū)，面積約3 006 km2(表3)。

圖8 地面沉降速率等值線圖及重點(diǎn)沉降區(qū)分布

表2 地面沉降嚴(yán)重程度劃分表

表3 重點(diǎn)沉降區(qū)特征表

2.2 地面沉降與地質(zhì)環(huán)境因素的關(guān)系分析

(1)研究區(qū)廣泛分布全新統(tǒng)粉質(zhì)黏土及砂層，全新統(tǒng)風(fēng)積層在黃河故道附近呈零星砂丘分布，由淡黃色、黃白色細(xì)砂、粉細(xì)砂組成，中更新統(tǒng)較單一的粉質(zhì)黏土主要出露于西部太行山山前及北部山前崗地。統(tǒng)計(jì)表明，研究區(qū)內(nèi)重點(diǎn)地面沉降區(qū)大部分處于砂及粉細(xì)砂和粉質(zhì)砂土分布區(qū)。

(2)研究區(qū)東北部斷裂方向以北東向?yàn)橹?，西南部以東西向、北西向?yàn)橹鳌目臻g上分布顯示，研究區(qū)的大部分重點(diǎn)沉降區(qū)與區(qū)域構(gòu)造、主要活動(dòng)斷裂構(gòu)造線具有較強(qiáng)的一致性。在研究區(qū)東北部的重點(diǎn)沉降區(qū)也呈北東方向展布，在南部呈東西向展布。表明研究區(qū)構(gòu)造運(yùn)動(dòng)與地面沉降密切相關(guān)(圖9)。

圖9 研究區(qū)重點(diǎn)沉降區(qū)與地質(zhì)構(gòu)造及地下水空間關(guān)系圖

(3)研究區(qū)長(zhǎng)期地下水開采加之近年來(lái)氣候干旱，降雨量偏少，使得研究區(qū)地下水位呈大面積下降勢(shì)態(tài)，淺層地下水已形成滑縣—濮陽(yáng)—清豐—南樂(lè)漏斗和溫縣—孟州漏斗2個(gè)區(qū)域性地下水位降落漏斗。其中滑縣—清豐—南樂(lè)淺層地下水區(qū)域降落漏斗為一開放式的降落漏斗，北部和河北、山東淺層地下水降落漏斗相連，20世紀(jì)70、80年代漏斗中心水位埋深4～7 m，80年代到90年代，漏斗中心水位埋深達(dá)14～20 m，年平均降幅高達(dá)0.9～1.3 m/a；進(jìn)入21世紀(jì)以來(lái)漏斗中心水位埋深維持在30 m±。溫縣—孟州漏斗形成于20世紀(jì)70年代，漏斗中心水位埋深12 m，1982年漏斗中心水位埋深16 m；1999年漏斗中心水位埋深22.28 m；2005年漏斗中心水位埋深33.50 m；2010年漏斗中心水位埋深減小為26.18 m；2014年漏斗中心水位埋深12.63 m。根據(jù)河南省超采區(qū)評(píng)價(jià)資料，淺層地下水在濮陽(yáng)—封丘—原陽(yáng)—新鄉(xiāng)—鶴壁—安陽(yáng)、孟州—溫縣—武陟、濟(jì)源、淇縣等區(qū)域?qū)僖话愠蓞^(qū)；中深層承壓水在濮陽(yáng)—清豐—南樂(lè)及以東區(qū)域?qū)俪蓞^(qū)。除了輝縣沉降區(qū)與礦山開采有關(guān)之外，其他重點(diǎn)沉降區(qū)均處于地下水超采區(qū)范圍；東部的重點(diǎn)沉降區(qū)還位于滑縣—濮陽(yáng)—清豐—南樂(lè)漏斗范圍或邊緣(圖9)。重點(diǎn)沉降區(qū)多分布于超采區(qū)和地下水漏斗范圍，地下水超量開采是發(fā)生地面沉降的重要因素。

(4)研究區(qū)地?zé)衢_發(fā)形式主要為井采，開采的熱儲(chǔ)層主要為新近系、古近系松散巖熱儲(chǔ)層，個(gè)別為其他基巖裂隙熱儲(chǔ)層，以零星開采為主，對(duì)地面沉降具有一定影響。隨著近30年來(lái)大面積的石油開采，高強(qiáng)度注水采油和超采地下水造成地下水位下降、地面沉降和地裂縫等一系列生態(tài)環(huán)境問(wèn)題。城鎮(zhèn)化的發(fā)展，城市建成區(qū)規(guī)模不斷擴(kuò)大，城市人口不斷增加，截止到2016年研究區(qū)城市化率已經(jīng)達(dá)到50%以上，城市基礎(chǔ)設(shè)施的建設(shè)在一定程度上也引起了城市局部沉降。

3 結(jié)論

(1)2014—2016年研究區(qū)全部是地面沉降區(qū)，大部分區(qū)域沉降速率10～30 mm/a，最大的沉降速率達(dá)114.85 mm/a。在安陽(yáng)、濮陽(yáng)、新鄉(xiāng)和鶴壁的局部區(qū)域地面沉降較嚴(yán)重，大部分區(qū)域?qū)儆诘孛娉两祰?yán)重程度較低區(qū)。綜合確定了8個(gè)重點(diǎn)沉降區(qū)，面積約3 006 km2。

(2)地下水超量開采是研究區(qū)地面沉降的最重要原因，區(qū)域構(gòu)造活動(dòng)、軟弱土層、城鎮(zhèn)化發(fā)展、石油和地?zé)岬荣Y源開發(fā)也與地面沉降的發(fā)生具有相關(guān)性。