應用短基線集技術監(jiān)測河南北部平原區(qū)地面沉降

張青鎖，于松暉，徐郅杰

1.河南省地質(zhì)環(huán)境監(jiān)測院 地質(zhì)災害防治技術中心，鄭州 450016；2.河南省地質(zhì)災害防治重點實驗室，鄭州 450016

0 引言

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

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

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

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

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

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

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

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

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

1.2 方法原理

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

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

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

(1)

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

(2)

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

AΦ=ΔΦ

(3)

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

Φ=(ATA)-1AT·ΔΦ

(4)

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

1.3 處理過程

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

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

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

(3)選擇使用最小費用流的方法開展相位解纏；為了降低噪聲誤差對解纏結(jié)果的影響，基于本次圖像的分辨率和當?shù)氐匚锏念愋吞卣鳎粚ο喔上禂?shù)>0.4的柵格點相位解纏(圖4)。

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

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

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

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

(6)通過控制點把不同圖像的形變數(shù)據(jù)統(tǒng)一到相同的坐標系統(tǒng)。控制點的地理位置信息存儲于圖像頭文件中，根據(jù)控制點坐標信息對形變速率圖進行校正處理，拼接3個條帶解譯結(jié)果并進行裁剪，獲得研究區(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)測面積21 192.80 km2。利用地理編碼后的形變速率圖制作地面沉降速率等值線圖(圖8)，2014—2016年研究區(qū)均發(fā)生了地面沉降，大部分區(qū)域沉降速率10～30 mm/a，最大的沉降速率達114.85 mm/a。參照地面沉降嚴重程度的劃分標準[21]，依據(jù)地面沉降速率等值線劃分研究區(qū)地面沉降嚴重程度(表2)。分區(qū)結(jié)果顯示：區(qū)內(nèi)僅有個別點沉降速率>80 mm/a，沒有圈出地面沉降嚴重程度高區(qū)；分布最為廣泛的是地面沉降嚴重程度較低區(qū)，面積19 609.34 km2，占監(jiān)測區(qū)92.53%。結(jié)合不同嚴重區(qū)域大小、人類活動的分布特征、集中連片等原則，綜合確定8個重點沉降區(qū)，面積約3 006 km2(表3)。

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

表2 地面沉降嚴重程度劃分表

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

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

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

(2)研究區(qū)東北部斷裂方向以北東向為主，西南部以東西向、北西向為主。從空間上分布顯示，研究區(qū)的大部分重點沉降區(qū)與區(qū)域構(gòu)造、主要活動斷裂構(gòu)造線具有較強的一致性。在研究區(qū)東北部的重點沉降區(qū)也呈北東方向展布，在南部呈東西向展布。表明研究區(qū)構(gòu)造運動與地面沉降密切相關(圖9)。

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

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

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

3 結(jié)論

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

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