杜自豪，孫哲，劉智敏，武鵬達(dá)，戴昭鑫

(1.山東科技大學(xué) 測繪與空間信息學(xué)院，山東 青島 266590； 2.中國測繪科學(xué)研究院，北京 100039)

1 引 言

地表沉降如何解決是在全球范圍內(nèi)的主要工程地質(zhì)問題之一，它是由非人類或人類相關(guān)導(dǎo)致的地表以下松散層的鞏固和壓實(shí)[1～3]。近數(shù)十年來，由于城市工業(yè)化進(jìn)程和大型城鎮(zhèn)化工程的迅速實(shí)施，城市沉降現(xiàn)象越來越頻繁，其會(huì)造成地表塌陷、城市市政基礎(chǔ)設(shè)施被嚴(yán)重破壞，沿海城市還面臨著海水倒灌的威脅，嚴(yán)重威脅到社會(huì)的可持續(xù)發(fā)展[4～6]。地表沉降嚴(yán)重程度和因素隨著地域的差異而產(chǎn)生巨大變化，所以，通過探究地表沉降的危害機(jī)理就可以為地表沉降有效防止提出科學(xué)合理的決策根據(jù)。

目前，關(guān)于地表沉降時(shí)空演化規(guī)律和其形成機(jī)制的研究成果，主要集中在兩方面：①對地表下沉及其影響因子的定性研究，如張金芝[7]等人研究發(fā)現(xiàn)，黃河三角洲地表沉降與油層開采和地表荷載有明顯的相關(guān)性；Le[8]等人深入研究了巖溶坍塌、建筑荷載、土地開墾和利用，以及人類活動(dòng)對地表下沉的影響；Tei je[9]等通過箱線圖解析緬甸仰光地區(qū)的地表沉降發(fā)現(xiàn)，地下水抽取是導(dǎo)致該區(qū)域地表下沉的決定性原因；在剖析成都市地面變形和發(fā)生機(jī)制的過程中，薛廉[10]等人提出正在加快的城市城鎮(zhèn)化建設(shè)和軌道交通快速發(fā)展，是造成城市地面沉降的主因。肖巍峰[11]在分析北京市地表沉降時(shí)，發(fā)現(xiàn)土地利用變化對地表沉降是有影響的；程霞[12]等研究黃河三角洲地表形變時(shí)發(fā)現(xiàn)鹽化工廠開采地下鹵水對地表沉降帶來的影響是很大。②各因素之間的“橫向”關(guān)系定量研究，如Merad[13]等人利用多因子評判模型，對礦井塌陷的影響開展了深入研究；Zhou[14]等人使用梯度提升決策樹模型，定量分析了地下水位、可壓縮沉積物厚度和地表負(fù)荷變化對沉降的影響；Omid[15]等人通過最大熵值法和規(guī)則集的遺傳算法，研究了地表下沉特征與地質(zhì)環(huán)境影響因子的相互關(guān)系?？偟膩碚f，地表沉降的原因是復(fù)雜的。盡管已有許多研究討論過地面沉降的誘因，但目前仍面臨著若干問題：對傳統(tǒng)的研究中沉降影響各種因素(如地殼運(yùn)動(dòng)、地下水開采、地下水水位、斷裂帶等)都無法掌握，且資料收集時(shí)間較長，在適用程度上不能廣泛的應(yīng)用分析，作出合理有效的防治措施；關(guān)于影響地表沉降的因素，多是各種因素的“橫向”對比分析，并沒有對影響因素進(jìn)行“縱向”分析，而這些研究偏重分析各影響因素對地表下沉的獨(dú)立影響與意義，對影響因素相互影響則有相對較少的研究。

關(guān)于地表沉降監(jiān)測，常規(guī)的監(jiān)測手段有水準(zhǔn)測量、全球?qū)Ш叫l(wèi)星系統(tǒng)(GNSS)等科學(xué)技術(shù)雖可提升較大的準(zhǔn)確度，但具有監(jiān)測范圍較小、工作量大和監(jiān)測成本高以及站點(diǎn)較少等缺點(diǎn)。SBAS-InSAR技術(shù)采用的是多主影像的沉降觀測，降低了主從影像發(fā)生時(shí)間和空間上的失相干現(xiàn)象的可能，且節(jié)約人力物力，大范圍的沉降監(jiān)測更加容易，同時(shí)結(jié)果展示也更加直觀。為此，本文通過運(yùn)用SBAS-InSAR技術(shù)獲得濱州市地表沉降點(diǎn)數(shù)據(jù)，并利用地理探測器模型分析自然與人類活動(dòng)影響因素對地表沉降空間分異性的影響，從而確定出地表沉降影響原因以及交互作用機(jī)理，為濱州市地面沉降科學(xué)治理提出了重要的理論依據(jù)。

2 研究區(qū)概況

濱州市(36°41′～38°16′N，117°15′～118°37′E)位于黃河三角洲西北部邊緣，渤海灣西南岸，沿海濕地十分廣闊，能源礦產(chǎn)和非金屬礦產(chǎn)資源豐富，濱州地勢由南向東北略有坡度。地表坡度下降1/6000～1/8000，整體海拔約為 -10 m～4 m[16]。2010年以來，濱州市市政府批復(fù)了《濱州北海經(jīng)濟(jì)開發(fā)區(qū)城市總體規(guī)劃(2011-2030)》，使得土地?cái)U(kuò)張，城鎮(zhèn)化建設(shè)加快，加快了地面沉降。選取了濱州市部分區(qū)縣(惠民縣、陽信縣、濱城區(qū)、無棣縣、沾化區(qū))作為研究區(qū)，位置如圖1所示。

圖1 研究區(qū)地理位置示意圖

3 數(shù)據(jù)與方法

3.1 研究數(shù)據(jù)

本文數(shù)據(jù)中主要包含地面沉降點(diǎn)數(shù)據(jù)和影響因子數(shù)據(jù)。地表沉降點(diǎn)數(shù)據(jù)是基于SBAS-InSAR方法利用2016-1-9—2020-9-20獲取的覆蓋濱州市部分地區(qū)的18幅上升Sentinel-TOPS合成孔徑雷達(dá)圖像(C波段)和輔助數(shù)據(jù)SRTM-1(shuttle radar topography mission)DEM數(shù)據(jù)及ESA提供的相應(yīng)精密軌道數(shù)據(jù)進(jìn)行地表沉降速率估算得到的。為了驗(yàn)證SBAS干涉測量結(jié)果，使用濱州市2016-2019年水準(zhǔn)點(diǎn)數(shù)據(jù)進(jìn)行校驗(yàn)。

自然因素和人類活動(dòng)是影響地表沉降的兩個(gè)因素，二者共同決定了地表沉降發(fā)生的情況。所以在選擇沉降影響因子時(shí)必須兼顧以上這兩個(gè)方面的因素外，還需要兼顧可獲取性與適用性，這里選定3類共16種影響因子，包括：土地利用類型。這里采用的是第三次國土調(diào)查數(shù)據(jù)，按照 1 km×1 km分別統(tǒng)計(jì)格網(wǎng)內(nèi)一級類土地利用類型。土壤質(zhì)地，其根據(jù) 1∶100萬土壤類型圖和第二次土壤普查數(shù)據(jù)編制而成，包括砂土、粉砂土、黏土。POI數(shù)據(jù)從高德開放平臺抓取獲得，這里篩選出鹽化工產(chǎn)業(yè)和化學(xué)化工類企業(yè)，對其進(jìn)行核密度分析和重分類處理。

3.2 研究方法

3.2.1 SBAS-InSAR技術(shù)原理

SBAS-InSAR技術(shù)旨在提高變形觀測的質(zhì)量，通過使用多景SAR數(shù)據(jù)來準(zhǔn)確觀測長時(shí)間的小變形現(xiàn)象，并使用新形成的具有相對較小時(shí)間和空間基線特征的干涉圖序列來最好地補(bǔ)償過去相關(guān)的影響。

假設(shè)測區(qū)N+1景雷達(dá)影像所在(t0，…tN)的時(shí)間段里，根據(jù)干涉關(guān)系的自由組合，可以得到M幅干涉圖，并滿足：

(1)

首先，從tA，tB兩景的圖像中去除不同時(shí)間的地形相位，根據(jù)假設(shè)tA

δφi(x，y)=φ(tB，x，y)-φ(tA，x，y)

(2)

其中φ(tA，x，y)和φ(tB，x，y)分別表示與tA，tB時(shí)間相對應(yīng)的相位，λ表示雷達(dá)波長，d(tA，x，y)和d(tB，x，y)分別代表在tA和tB各自時(shí)間每個(gè)像元(x，y)在LOS(line of sight)方向上的地表形變，相對于參考時(shí)刻t0。

方程(2)中的相位φ表示為平均速度v和任何兩個(gè)測量時(shí)間t之間的時(shí)間的乘積，被定為具有物理含義的沉降序列，表達(dá)式為：

(3)

故，第i幅干涉圖的相位可表示為：

(4)

將其改寫為矩陣形式：

Bv=δφ

(5)

式中：B—M×N矩陣。最后，用SVD方法可以得到速率矢量v的廣義逆解，然后通過對每個(gè)時(shí)間段的速率在一個(gè)時(shí)間區(qū)間內(nèi)進(jìn)行微分，可以找到每個(gè)時(shí)間段的形變量。

3.2.2 精度驗(yàn)證

為檢驗(yàn)使用SBAS-InSAR技術(shù)分析成果的真實(shí)性，根據(jù)《CHT6006-2018時(shí)間序列InSAR地表形變監(jiān)測數(shù)據(jù)處理規(guī)范》[17]，在這里使用二等水準(zhǔn)實(shí)測數(shù)值對InSAR地表變形監(jiān)測結(jié)果進(jìn)行精度分析。為確保精度分析的客觀與精確，采用鄰近點(diǎn)原則，即對水準(zhǔn)點(diǎn)數(shù)據(jù)設(shè)置 80 m緩沖區(qū)，并比較水準(zhǔn)點(diǎn)數(shù)據(jù)與InSAR監(jiān)測結(jié)果的差值，精度評價(jià)指標(biāo)為：

①誤差絕對值，宜優(yōu)于15 mm/a；

②誤差標(biāo)準(zhǔn)差，宜≤8 mm/a。

3.2.3 地理探測器

地理探測器是王勁峰[18]發(fā)明的一種統(tǒng)計(jì)學(xué)方法，用于檢測空間結(jié)構(gòu)異質(zhì)性現(xiàn)象并揭示其背后的驅(qū)動(dòng)力，它被廣泛用于分析自然環(huán)境和人類經(jīng)濟(jì)機(jī)制的空間變化的驅(qū)動(dòng)力。地理探測器模型由四個(gè)模塊組成：因子探測、交互探測、風(fēng)險(xiǎn)探測和生態(tài)探測。地理探測器不僅可以分析每個(gè)影響因素之間的相互作用，還可以分析兩個(gè)影響因素中每個(gè)因素之間的相互作用。本文將分別使用單因素探測和交互探測模塊，解析單因素與不同影響因素之間的交互作用及其對地面沉降的影響。

(1)因子探測器

因子探測器用來分析自變量(不同驅(qū)動(dòng)因素X)對因變量(Y)空間分布形成的影響程度，識別不同驅(qū)動(dòng)因素對因變量空間的驅(qū)動(dòng)力大小(q值度量)，表達(dá)式為：

(6)

(2)交互探測

交互探測主要是為了確定各種影響因子間的交互多對地面沉降格局的影響關(guān)系，即各種影響因子在聯(lián)合相互作用時(shí)，是否會(huì)增加或降低對因變量Y空間分布的影響程度。通過分別計(jì)算值q(X1)、q(X2)值及兩因子圖層疊加后的q(X1∩X2)數(shù)值，并對比前三者數(shù)值大小，從而確定了兩因子之間是否具有相互作用，以及對地面沉降空間格局的影響力是加強(qiáng)還是削弱，并由此推廣到三因子的相互影響。

因子之間交互作用分別表現(xiàn)為增強(qiáng)、雙因子增強(qiáng)、非線性增強(qiáng)、減弱、非線性減弱、獨(dú)立，判斷方式如表1所示：

交互作用類型 表1

4 結(jié)果與討論

4.1 基于SABS-InSAR的地表沉降監(jiān)測結(jié)果及精度驗(yàn)證

通過SARscape軟件對2016-1-9—2020-9-2018景Sentinel-1A數(shù)據(jù)進(jìn)行SBAS-InSAR處理，獲得地表沉降速率矢量點(diǎn)，經(jīng)過克里金插值得到研究區(qū)的沉降速率圖，如圖2所示，正值代表地表抬升，負(fù)值代表地表下降。其中，研究時(shí)段內(nèi)地表沉降平均沉降速率范圍為 -62.8 mm/a～34.6 mm/a，沉降較嚴(yán)重區(qū)域集中分布在濱州市無棣縣和沾化區(qū)的北部沿海地區(qū)。

圖2 基于SBAS-InSAR技術(shù)監(jiān)測的2016-1—2020-9研究區(qū)沉降速率分布結(jié)果

為了驗(yàn)證SBAS-InSAR技術(shù)對地表監(jiān)測的可靠性，現(xiàn)用21個(gè)二等水準(zhǔn)點(diǎn)數(shù)據(jù)建立80米緩沖區(qū)與緩沖區(qū)內(nèi)的沉降數(shù)據(jù)進(jìn)行比較，得到結(jié)果如表2，可以得到緩沖區(qū)內(nèi)的沉降點(diǎn)有259個(gè)，與水準(zhǔn)點(diǎn)進(jìn)行比較，可以看出，誤差與標(biāo)準(zhǔn)差都符合《規(guī)范》[17]，其中差值絕對值的 ≤6 mm/a占134個(gè)，差值絕對值分布如圖3所示，證明InSAR監(jiān)測結(jié)果具有較高的可靠性，可用于研究和分析。

水準(zhǔn)數(shù)據(jù)與InSAR監(jiān)測結(jié)果比較 表2

圖3 差值絕對值分布

按照《地質(zhì)災(zāi)害危險(xiǎn)性評估規(guī)范》對地表沉降速率劃分等級，探究地表沉降空間分布情況(圖2)和沉降面積統(tǒng)計(jì)(圖4)。由圖4可知：研究區(qū)以較低沉降為主，占有 3 351.17 km2，其主要分布在鹽田區(qū)域的埕鎮(zhèn)和馬山子鎮(zhèn)北部，以及濱城鎮(zhèn)；沉降區(qū)域最嚴(yán)重僅有 4.66 km2，分布在濱州市北部與渤海相接區(qū)域。研究區(qū)總體以較低沉降為主，呈北部沉降，南部沉降較小分布狀態(tài)。

圖4 研究區(qū)沉降面積統(tǒng)計(jì)

4.2 基于地理探測器的地表沉降影響因素分析

4.2.1 地表沉降影響因子單因子探測

圖2顯示，研究區(qū)域內(nèi)的沉降速率存在明顯空間異質(zhì)性現(xiàn)象，因此，可以利用地理探測器對其進(jìn)行定量分析，研究地面沉降-空間分布的影響機(jī)理。在本研究中，以2016-1-9—2020-9-20地表沉降速率為因變量Y，以土地利用類型分類面積為自變量X，進(jìn)行地理探測器探測，計(jì)算其類內(nèi)影響因子對地表沉降速率分布的驅(qū)動(dòng)力值(q值)，其強(qiáng)度大于0.01的驅(qū)動(dòng)力如表3所示。

土地利用類型對地表沉降速率的影響力 表3

從表3中可以看出，工礦用地對地表沉降速率的作用強(qiáng)度最大，工礦用地包含鹽田、工業(yè)用地、采礦用地。鹽田區(qū)的曬鹽及地下鹵水的開采，導(dǎo)致地下水位下降，從而引起地表沉降。

針對土地利用類型對地表沉降速率的影響，得到工礦用地是影響地表沉降的主要因素，在研究區(qū)內(nèi)，工礦用地占比最大的是鹽田，如圖5(c)所示，因此，這里對研究區(qū)內(nèi)的POI類型中的鹽化工產(chǎn)業(yè)及化學(xué)化工類企業(yè)進(jìn)行核密度分析，對其按照自然間斷法重分類，作為自變量X輸入地理探測器，得到其驅(qū)動(dòng)力如表4所示。從結(jié)果表中可以得到鹽化工產(chǎn)業(yè)的q值最大，鹽化工產(chǎn)業(yè)分布如圖5(a)所示，從圖中可以看到，沉降區(qū)內(nèi)的鹽化工產(chǎn)業(yè)分布相對較多，而鹽化廠是生產(chǎn)鹽和化工原料的基地，如鹵水開采，溴素提取。如果是長期或長期不合理利用地下鹵水將會(huì)導(dǎo)致地下鹵水水位的降低，從而產(chǎn)生地面沉降。

圖5 實(shí)驗(yàn)數(shù)據(jù)

POI類型對地表沉降速率的驅(qū)動(dòng)力 表4

土壤質(zhì)地對地表沉降速率的驅(qū)動(dòng)力 表5

從圖5(a)鹽化工產(chǎn)業(yè)POI的分布可以看出，除了沉降區(qū)分布有鹽化工廠外，在濱州的東部和西南部均有分布，但對比發(fā)現(xiàn)，這些地方未發(fā)生嚴(yán)重沉降，初步猜想，可能是工礦用地區(qū)域與其他因子結(jié)合，對沉降區(qū)域的影響有所增強(qiáng)，因此，引入土壤質(zhì)地類型數(shù)據(jù)，對地表沉降進(jìn)行地理探測器分析，進(jìn)一步探究沉降分布機(jī)制，得到土壤質(zhì)地對地表沉降的驅(qū)動(dòng)力如表5所示。由此可以得出，砂土的驅(qū)動(dòng)力遠(yuǎn)大于粉砂土與黏土對地表沉降的影響，砂土液化導(dǎo)致地表沉降。

4.2.2 地表沉降影響因子交互探測

土地利用類型、設(shè)施類型POI數(shù)據(jù)、土壤質(zhì)地中對地表沉降影響的最強(qiáng)的是工礦用地(q=0.294，p=0.000)、鹽化工產(chǎn)業(yè)(q=0.063，p=0.000)、砂土(q=0.373，p=0.000)，對研究區(qū)進(jìn)行雙因子交互和三因子交互，其作用探測顯示如表6、表7，在表中可發(fā)現(xiàn)，不同因素相互影響對地表沉降的影響力均超過單一作用的影響力，并體現(xiàn)為雙因子相互增強(qiáng)，其中，砂土分別和鹽化工產(chǎn)業(yè)、工礦用地進(jìn)行交互作用，對地表沉降的影響力均達(dá)到了 0.415 9、0.402 9，三因子交互作用影響力達(dá)到了0.624。

雙因子交互對地表沉降的影響 表6

三因子交互對地表沉降的影響 表7

根據(jù)表6，表7的結(jié)果，其可能解釋如下：

(1)濱州市位于斷裂帶之上[19]又位于渤海灣南岸，歷史上曾出現(xiàn)多次地震[20]，使得濱州無棣縣和沾化區(qū)北部沿海地區(qū)沙土地質(zhì)受到外力后，使得大面積區(qū)域的砂土發(fā)生液化，其位置區(qū)域如圖6所示。當(dāng)沙土變成液體時(shí)，孔隙水在超孔壓力下上下移動(dòng)，當(dāng)沙土上不再有不透水的覆蓋物時(shí)，地下水就會(huì)自動(dòng)排出土壤；如果沙土層上有透水性低的黏土，超孔水的壓力超過了覆蓋物的強(qiáng)度，地下水就會(huì)通過覆蓋物輸送沙粒，沙土的保水能力很差。因此導(dǎo)致地表不均勻沉降，再加上鹽化工產(chǎn)業(yè)主要是以開采地下鹵水為主要原料制取溴化合物和鹽，而過量利用地下鹵水都會(huì)導(dǎo)致地下鹵水水位下降，兩者相互促進(jìn)，從而導(dǎo)致地面沉降狀況越來越嚴(yán)峻。

圖6 濱州市無棣縣、沾化區(qū)砂土液化程度分區(qū)圖[16]

(2)從圖5(c)可以看出，無棣、沾化和北海新區(qū)工礦區(qū)域分布最多的是鹽田，宜鹽面積可達(dá)9.6萬公頃，現(xiàn)已有鹽田100萬畝，年產(chǎn)均200萬噸，根據(jù)2016-2019年濱州《統(tǒng)計(jì)年鑒》可得原鹽產(chǎn)量如表8所示，“井灘曬鹽”的主要采用的是地下鹵水，已探明靜態(tài)儲

量30億立方米。因此，地下鹵水的長時(shí)間開采，再加上砂土質(zhì)地的液化，導(dǎo)致地下水位下降，加劇地表沉降。

2016-2019濱州原鹽產(chǎn)量(萬噸) 表8

(3)土地利用類型因子與設(shè)施類型POI因子中影響最大的工礦用地和鹽化工產(chǎn)業(yè)，雖然這種交互作用后的影響力是雙因子增強(qiáng)，但并不像砂土分別與鹽化工產(chǎn)業(yè)和工礦用地交互后，對地表沉降作用強(qiáng)度大，表明這兩個(gè)因素在影響地表沉降過程中可能有重疊，鹽田制鹽和鹽化工產(chǎn)業(yè)的生產(chǎn)都需要對地下鹵水進(jìn)行開采，導(dǎo)致地表沉降。

5 結(jié) 語

以濱州市5個(gè)區(qū)縣為主要研究區(qū)，利用SBAS-InSAR技術(shù)獲取了2016-1-9—2020-9-20的18個(gè)時(shí)間段地表沉降數(shù)據(jù)，并與土地利用類型數(shù)據(jù)、土壤質(zhì)地?cái)?shù)據(jù)和設(shè)施類型POI等多源數(shù)據(jù)相結(jié)合，基于地理探測器模型，定量監(jiān)測了影響地表沉降的各種因素，并深入分析了影響因素之間交互作用對地表沉降的影響力。該研究表明：

(1)研究區(qū)地表沉降以中等沉降為主，沉降區(qū)域主要分布于濱州市無棣縣和沾化區(qū)北部，南部沉降區(qū)域較少。

(2)單因子探測結(jié)果。對土地利用數(shù)據(jù)、設(shè)施類型POI數(shù)據(jù)、土壤質(zhì)地?cái)?shù)據(jù)中對地面沉降影響力較大的因素依次是工礦用地、鹽化工產(chǎn)業(yè)和砂土土質(zhì)。

(3)因子交互探測結(jié)果。單因子影響力最強(qiáng)的三個(gè)因子雙交互作用后都表現(xiàn)為雙因子增強(qiáng)，其中砂土分別與工礦用地、鹽化工產(chǎn)業(yè)相互疊加后，綜合影響強(qiáng)度體現(xiàn)得更為突出；三因子作用后對地表沉降的影響則尤為突出。

由于濱州市北部鹽田產(chǎn)量的逐漸增加和《濱州北海經(jīng)濟(jì)開發(fā)區(qū)城市總體規(guī)劃(2011-2030)》，隨著地下鹽鹵過量開發(fā)以及城市化飛躍型發(fā)展趨勢，目前濱州市及其以北區(qū)域已發(fā)生了強(qiáng)烈的地面沉降現(xiàn)象，建議在以后的開發(fā)建設(shè)過程中，合理利用資源，減緩地表沉降。本研究在考慮研究區(qū)位置和人文情況下，僅選用了土地利用和土壤質(zhì)地因子，并引入大數(shù)據(jù)方面的POI因子，并將空間統(tǒng)計(jì)學(xué)方法應(yīng)用于地表沉降，定量分析影響因子對地表沉降的作用強(qiáng)度，今后可收集地下鹵水或地下水位數(shù)據(jù)作為補(bǔ)充數(shù)據(jù)研究地表沉降，并探測其與他影響因子交互對地表沉降的作用。