基于SBAS-InSAR的濱州市地表沉降及其對影響因素分析

杜自豪，孫哲，劉智敏，武鵬達，戴昭鑫

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

1 引 言

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

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

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

2 研究區概況

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

圖1 研究區地理位置示意圖

3 數據與方法

3.1 研究數據

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

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

3.2 研究方法

3.2.1 SBAS-InSAR技術原理

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

假設測區N+1景雷達影像所在(t0，…tN)的時間段里，根據干涉關系的自由組合，可以得到M幅干涉圖，并滿足：

(1)

首先，從tA，tB兩景的圖像中去除不同時間的地形相位，根據假設tA

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

(2)

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

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

(3)

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

(4)

將其改寫為矩陣形式：

Bv=δφ

(5)

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

3.2.2 精度驗證

為檢驗使用SBAS-InSAR技術分析成果的真實性，根據《CHT6006-2018時間序列InSAR地表形變監測數據處理規范》[17]，在這里使用二等水準實測數值對InSAR地表變形監測結果進行精度分析。為確保精度分析的客觀與精確，采用鄰近點原則，即對水準點數據設置 80 m緩沖區，并比較水準點數據與InSAR監測結果的差值，精度評價指標為：

①誤差絕對值，宜優于15 mm/a；

②誤差標準差，宜≤8 mm/a。

3.2.3 地理探測器

地理探測器是王勁峰[18]發明的一種統計學方法，用于檢測空間結構異質性現象并揭示其背后的驅動力，它被廣泛用于分析自然環境和人類經濟機制的空間變化的驅動力。地理探測器模型由四個模塊組成：因子探測、交互探測、風險探測和生態探測。地理探測器不僅可以分析每個影響因素之間的相互作用，還可以分析兩個影響因素中每個因素之間的相互作用。本文將分別使用單因素探測和交互探測模塊，解析單因素與不同影響因素之間的交互作用及其對地面沉降的影響。

(1)因子探測器

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

(6)

(2)交互探測

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

因子之間交互作用分別表現為增強、雙因子增強、非線性增強、減弱、非線性減弱、獨立，判斷方式如表1所示：

交互作用類型 表1

4 結果與討論

4.1 基于SABS-InSAR的地表沉降監測結果及精度驗證

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

圖2 基于SBAS-InSAR技術監測的2016-1—2020-9研究區沉降速率分布結果

為了驗證SBAS-InSAR技術對地表監測的可靠性，現用21個二等水準點數據建立80米緩沖區與緩沖區內的沉降數據進行比較，得到結果如表2，可以得到緩沖區內的沉降點有259個，與水準點進行比較，可以看出，誤差與標準差都符合《規范》[17]，其中差值絕對值的 ≤6 mm/a占134個，差值絕對值分布如圖3所示，證明InSAR監測結果具有較高的可靠性，可用于研究和分析。

水準數據與InSAR監測結果比較 表2

圖3 差值絕對值分布

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

圖4 研究區沉降面積統計

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

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

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

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

從表3中可以看出，工礦用地對地表沉降速率的作用強度最大，工礦用地包含鹽田、工業用地、采礦用地。鹽田區的曬鹽及地下鹵水的開采，導致地下水位下降，從而引起地表沉降。

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

圖5 實驗數據

POI類型對地表沉降速率的驅動力 表4

土壤質地對地表沉降速率的驅動力 表5

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

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

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

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

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

根據表6，表7的結果，其可能解釋如下：

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

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

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

量30億立方米。因此，地下鹵水的長時間開采，再加上砂土質地的液化，導致地下水位下降，加劇地表沉降。

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

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

5 結 語

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

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

(2)單因子探測結果。對土地利用數據、設施類型POI數據、土壤質地數據中對地面沉降影響力較大的因素依次是工礦用地、鹽化工產業和砂土土質。

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

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