0 引言
地面沉降是在自然與人類活動交互影響下，由於地下鬆散地層固結壓縮導致區域性地面高程降低的一種不可逆的地質災害現象，具有發展緩慢、不易引起人們察覺、波及範圍廣、持續時間長等特點。近幾十年來，隨著工業化進程和大規模城市化建設的快速推進，沉降現象越來越普遍，嚴重的沉降會導致地面塌陷、建築物開裂、市政設施被破壞等一系列災害問題，影響人類社會的可持續發展。截至2011年底，中國50多個城市出現地面沉降，其中，華北平原、長三角地區和汾渭盆地已發展成了沉降重災區。地面沉降通常呈現區域性特徵，即沉降量及成因隨區域的不同而發生改變，因此，探測區域地面沉降的影響機制能夠為地面沉降有效防治防控提供科學的決策依據。
目前，關於地面沉降影響機制的研究主要集中於以下兩方面：①地面沉降與影響因素的定性研究，如Teije等人利用箱線圖等分析緬甸仰光地區的地面沉降，發現地下水抽取是引起該地區地面沉降的決定性因素；薛廉等人在分析成都市地面變形及發生機制的過程中，認為快速的城市化進程和軌道交通建設的迅猛發展是導致市中心區域地面沉降的主因；荊創利等人揭示了地殼運動、地下水開採、地表建設等自然因素和人為因素所導致的地面沉降規律。②採用基於機器學習或空間統計方法分析各影響因素對地面沉降的作用，如Zhou等人基於梯度提升決策樹模型分別定量分析了地下水位變化、可壓縮沉積層厚度、地面荷載對沉降的影響程度；Omid等人採用最大熵值法和規則集的遺傳算法研究了地面沉降特徵與地質環境影響因素的關係。雖然學者們對地面沉降的影響機制進行了一系列研究並取得了較為豐碩的成果，但大多數研究側重於分析各影響因素對地面沉降的獨立影響作用，對因素間交互作用鮮有探討。而實際上，多種因子之間的複雜相互作用共同決定了地面沉降的空間格局。探測地面沉降分布的主要影響因子以及影響因子間交互作用有助於更全面地掌握區域地面沉降分布機制。
傳統的地面沉降監測手段（如地面監測站、高精度水準測量網、全球導航衛星系統（ G l o b a l Navigation Satellite System，GNSS）測量網）能提供高精度的變形信息，但具有監測成本高、監測範圍小等缺點。永久散射體幹涉測量（Persistent Scatterer Interferometric Synthetic Aperture Radar，PSInSAR）技術利用穩定且小於像元尺寸的永久散射體，能有效地去除大氣效應貢獻值，獲得高精度的地表形變值，具有監測範圍廣、成本低、精度高等特點，被廣泛地用於大範圍的地表變形監測。因此，本文利用PS-InSAR技術獲得北京市地面沉降點（Persistent Scatterer，PS）數據，藉助於地理探測器模型分析自然與人為要素對地面沉降分異的影響，識別出地面沉降影響因素及其交互作用機制，為理解北京市地面沉降的影響機制以及科學防治提供重要的依據。
1 研究區概況及方法
1.1 研究區與數據來源
北京，是我國的首都，也是世界上著名的國際化大都市，位於115°25′～117°30′E，39°26′～41°30′N，是國內地面沉降嚴重的城市之一。選取北京市主城區（海澱區、朝陽區、東城區、西城區、石景山區、豐臺區、順義區、通州區）為研究區，如圖1所示，該區域位於北京市東南部，是人口較為集中、社會經濟發展水平較高的區域，也是沉降嚴重區域。據《2018年北京統計年鑑》記載，截至2018年底，該區面積是3 303.25 km²，常住人口達1 440.6萬，2018年實現地區生產總值20 267億元。研究區地形以平原為主，西部與東北部有少量山地丘陵，平原面積約佔總面積的98%，屬於暖溫帶半溼潤大陸性季風氣候，夏季高溫多雨，冬季寒冷乾燥，降水季節分布極不均勻，全年約80%的降水集中在夏季6—8月。
本文數據主要包括地面沉降PS點集數據和影響因子數據。PS點集數據是利用GAMMA軟體處理24景空間解析度為30 m的SAR影像，獲取了2003年6月至2010年9月的23個時間段總計230 028個PS點數據，包括經緯度坐標、時間、年均沉降速率、累積沉降量等屬性信息，獲取時間間隔不定，由2～10個月不等，PS點集密度如圖1所示。
圖1 研究區及PS點密度
Fig.1 Study area and PS point density
地面沉降發生受自然因素和人類活動影響，兩者的作用方向和強度共同決定了地面沉降發生的狀況，因此在選取沉降影響因子時需要綜合考慮這兩方面的因素，考慮到研究區的實際情況和數據的可獲取性，選定下列7類共14種影響因子，包括：①粘性土可壓縮層厚度。可壓縮層粘性土受地面荷載的等外力作用容易壓縮，進而引起地面沉降，研究區內可壓縮層厚度取值範圍是51.2～236.8 m，其中大於200 m部分佔很小比例，部分區域數據缺失，為此，將其分為5個等級（50～100 m、100～150 m、150～200 m、大於200 m以及缺失類）。②土壤類型。通州區、順義區主要以潮土為主，西城區、東城區、朝陽區、海澱區、豐臺區等建成區以城市混凝土為主。③地形因子。包括高程、坡度和坡向。④地下水位高度。以統一高程作為基準面，地下水水位高度數值越小，表示地下水水位越低，地下水含量越少，地下水被超採。考慮到研究區地下水水位高度範圍，以5 m和10 m為等高距，將其分為12個等級（小於-10 m、-10～-5 m、-5～0 m、0～5 m、5～10 m、10～15 m、15～20 m、20～25 m、25～30 m、30～40 m、40～50 m、50～60 m）。⑤地質類影響因子。包括距斷裂帶距離、距5級以下地震距離(簡稱小地震源)和距5級以上地震距離(簡稱大地震源)。⑥土地利用類型。分別為草地、林地、地表水、農用地、建設用地和未利用地。⑦交通可達類影響因子。包括距國道距離、距省道距離、距鐵路距離、距河流距離。影響因子數據來源於地理空間數據云平臺、國家地球系統科學數據中心等，詳細數據來源見表1。
表1 影響因子數據來源
Tab.1 Data source of associated factors
利用ArcGIS軟體將點、線類型的因子數據轉換成到格網的距離值，在此基礎上，通過Python語言對各因子屬性數據表格進行處理，基於面積佔優原理對各影響因子數據實施格網化，並與PS點網格一一對應，實現屬性關聯。
1.2 研究方法
1.2.1 基於空間自相關的PS點格網化
由於地面沉降數據是離散PS點，空間上呈離散分布，為了便於後續的空間關聯分析，需對離散PS點進行不同尺寸的網格化。空間自相關性是空間數據所特有的性質之一，網格化後的數據應以最大限度保留原始數據的空間自相關性為前提，常用的測度方法有莫蘭指數（Moran's I），故本文採用空間自相關的Moran's I係數差法來確定最優空間格網，其計算公式如下：
式中，Ig和Ib分別為格網化後數據和原始沉降均值數據的Moran's I係數；S為格網化數據與原始沉降均值數據的Moran's I係數差值；對莫蘭指數差值取絕對值，S值越小，表明該尺寸下格網化後數據與原始數據的自相關係數越相近，網格採樣後的數據較好地保留了原始數據的自相關結構，則可將該尺寸作為最優的空間格網化參數。
1.2.2 地理探測器
地理探測器是探測空間分異性並揭示其背後驅動力的一組統計學方法，廣泛應用於地理要素空間變異的驅動力分析、精準扶貧、共享單車騎行時空分布特徵等諸多領域。地理探測器模型包括因子探測、交互探測、風險探測和生態探測4個模塊，這些模塊既可以分析單影響因子的作用，也可以分析兩影響因子交互的作用。本文分別採用因子探測與交互探測兩個模塊解釋單因子和不同影響因子之間的交互作用對地面沉降影響力的作用。
1）因子探測。其核心思想是比較某一地理環境因子（自變量X）與地理事物（因變量Y）的變化在空間上是否具有顯著的一致性，如地理環境因子與地理事物的變化具有一致性，則說明該地理環境因子對地理事物的發生和發展具有決定意義，用q值來度量，其計算公式為：
式中，h=1,2,…,L為地面沉降發生風險各指標的分層；Nh和N分別為h層和全區的單元數；σh²和σ²分別為h層和全區因變量Y值的方差；SSW和SST分別為層內方差之和和全區方差；q值的範圍為[0,1]，q值越大說明因子X對地面沉降分布的空間分布影響力越強。進一步，通過非中心的F分布的顯著性檢驗結果（用p值表示）對q統計量進行顯著性檢驗，若p值小於預設的顯著性水平（0.05），則拒絕原假設（分層條件下Y變量呈現隨機分布），而接受備選假設（分層條件下Y變量呈現空間聚集分布），否則，q值不顯著，沒有充分的證據拒絕原假設。
2）交互探測。主要用於研究兩個不同影響因子之間的交互對地面沉降格局的作用關係，如對於影響地面沉降格局影響因子X1和X2，通過空間疊加X1和X2形成新圖層X3，X3的屬性由X1和X2共同決定。通過比較X1和X2層的影響力和X3的影響力，可以判斷兩個因子交互作用對地面沉降空間格局與單個影響因子對地面沉降空間格局的影響力是強化還是弱化。
交互探測由下列表達式組成：若q (X1∩X2 )＜min(q (X1 ) ,q (X2 ))，則表示因子X1和X2交互後非線性減弱；若min(q (X1 ) ,q (X2 ))＜q (X1∩X2 )＜max(q(X1),q(X2))，則表示因子X1和X2交互後單因子非線性減弱；若q(X1∩X2)＞max(q(X1),q(X2))且q(X1∩X2) ＜q(X1)+q(X2)，則表示因子X1和X2交互後雙因子增強；若q(X1∩X2)＞q(X1)+q(X2)，則表示因子X1和X2交互後非線性增強；若q(X1∩X2)=q(X1)+q(X2)，則說明X1和X2相互獨立。本文計算結果均使用軟體GeoDetector計算獲得。
2 結果分析
2.1 最優網格尺寸
採用ArcGIS軟體分別計算PS原始點集與不同格網尺寸下的Moran's I係數，並求得其係數差。由圖2可知，當格網尺寸為1 700 m時，Moran's I係數差值最小為0.062，因此，本文基於1 700 m×1 700 m網格對原始沉降PS點實施重採樣，得到地面沉降柵格數據，作為後續關聯分析和地理探測的數據源。
圖2 格網數據Moran's I係數差
Fig.2 Moran's I difffference of grid data
2.2 影響因子離散數據分級
以時序沉降數據獲得的多年年平均沉降速率定義為地理探測器中因變量Y，其分布格局如圖3所示，上文提到的14種影響因子作為地理探測器中的自變量X1， X2，…，X14，其中，因變量多年年平均速率和自變量距斷裂帶距離、距河流距離等10個影響因子均為連續型變量。由於地理探測器中的影響因子變量應為離散型，因此，首先需對連續型的自變量與因變量進行基於自然斷點法的離散化。自然斷點法是一種通過減少類內距離而增大類間距離的自然聚類方法，即通過迭代比較每個分組和分組中元素的均值和觀測值之間的平方差之和來確定值分組中的最佳排列，分類後數據具有強度大、表達效果良好的特點，被廣泛應用於離散化數值型數據。對因變量多年年平均速率和自變量距斷裂帶距離、距河流距離等10個連續型影響因子進行分級，當類間距離與類內距離之差最大，則表明該分類級別是最優的。其中，多年平均沉降速率的分級如圖4所示，總的分類結果見表2（表中：Y為多年年平均速率，X1為距斷裂帶距離，X2為距河流距離，X3為距國道距離，X4為距省道距離，X5為距鐵路距離，X6為高程，X7為坡度，X8為坡向，X9為距5級以上地震源距離，X10為距5級以下地震源距離）。如離散指標Y，當分級數量為11時，類間距離與類內距離之差為24.81達到最大，故對指標Y分類的最優分類數目為11。類似地，計算其他指標因子的類間距離和類內距離，以兩者差值最大為原則來確定其相應的最優分類數目。
圖3 多年平均沉降速率分布格局
Fig.3 Distribution pattern of annual average settlement rate
圖4 多年平均沉降速率分級
Fig.4 Classifification of annual average settlement rate
表2 離散數據自然斷點法分級
Tab.2 Classifification of discrete data based on natural break method
2.3 地面沉降影響單因子探測
從圖3可以發現，研究區內的沉降速率具有明顯的空間分異性，因此有必要採用地理探測器模型進一步定量分析地面沉降空間分異的影響機制。通過地理探測器計算各影響因子對地面沉降速率分布的影響力（q值），結果表明只有8種地理環境要素影響了北京市地面沉降速率分布的空間格局，其作用強度由大到小分別為：地下水位(q=0.41,p=0.00)、可壓縮層厚度(q=0.33,p=0.00)、坡度(q=0.08,p=0.00)、土地利用類型(q=0.07,p=0.00)、土壤類型(q=0.04,p=0.00)、高程(q =0.03,p =0.00)、距5級以上地震源距離(q=0.02,p=0.01)和距鐵路距離(q=0.01,p=0.01)，其他6種因子沒明顯的影響力。
探測結果表明人類活動引起的地下水位下降對研究區的地面沉降作用強度遠遠大於其他因子。研究區北京是一個嚴重缺水的特大型城市，城區2/3以上的供水來源於地下水，地下水開採引起的地下水位下降是北京地面沉降形成與發展的關鍵外在因素。在滿足顯著性水平0.05的前提下，可壓縮層厚度、坡度、土地利用類型、土壤類型、高程、距5級以上地震源距離這6種因子的作用強度存在較大差異，其中，可壓縮層厚度因子的影響力遠遠大於其他因子，與地下含水層所對應的可壓縮層組的分布、結構、厚度是形成地面沉降的主要內在因素。此外，距鐵路距離對整個研究區沉降的作用強度不大，但對鐵路沿線地區的沉降還是有一定作用強度。
2.4 地面沉降影響因子交互探測
交互作用探測結果顯示見表3，大部分不同因素交互作用對地面沉降的影響力均大於單獨作用的影響力，各影響因素間的交互作用類型有非線性增強型和雙因子增強型兩種（表中A↑B表示A和B非線性增強；A↑↑B表 示A和B雙因子增強）。其中，地下水分別與距大地震源距離、距小地震源距離、土壤類型、距省道距離、距河流距離、距國道距離、距斷裂帶距離、距鐵路距離、坡向在影響地面沉降方面是非線性增強的；而地下水分別與可壓縮層厚度、土地利用、坡度具有交互協同的作用，即兩因子在影響地面沉降方面是雙因子增強的，效果沒有非線性增強型顯著。可壓縮層厚度分別與土壤類型、距大地震源距離、距省道距離、距小地震源距離在影響地面沉降方面是非線性增強的。
表3 交互因子對地面沉降的影響(q>0.4)
Tab.3 Inflfluence of interaction factors on land subsidence (q>0.4)
各類因素交互探測的結果顯示，地下水位和可壓縮層厚度作為對地面沉降最大影響的單因子，在與其它較小影響力的因子疊加時，對地面沉降的綜合作用是增強的。根據表3的結果，其可能解釋如下：
1）單因子影響力最大的地下水位和可壓縮層厚度指標，交互作用後的影響力並不是最大，說明兩種因子相互疊加後在解釋地面沉降過程中可能存在作用重複，地下水水位與可壓縮層厚度相互影響，地下水大量抽取導致地下水水位下降，可壓縮層內缺少填充物而發生壓縮，外來補充水源無法充分吸收補充地下水，兩者相互促進，地面沉降現狀不斷惡化。
2）土壤類型與地下水位存在密切關係，降水是地下水補給的主要來源，降水形成的土壤水分入滲將穿過非飽和地帶形成對地下水的補給。由於不同土壤類型的土壤顆粒間縫隙不同，其滲水性能不一致，同時土壤蒸發鋒面也不同，導致降水形成地下水補給的時間也不同。降水對地下水補給與地下水位埋深具有密切聯繫，如地下水位埋深較淺時，降雨形成的土壤水分入滲將在較短時間內穿過非飽和地帶，形成對地下水的補給，雨停後發生的蒸發，其土壤蒸發鋒面將緩慢下移最終抵達地下水面，形成潛水蒸發。此外，土壤類型作為地面沉降的下墊面，在影響地下水分布的同時，導致土壤空隙減小，有可能在一定程度上影響可壓縮層狀況。
3）土地利用類型變化主要是由於人類活動導致的，如成片的高大建築物和鐵路工程建設，大大增加了地面荷載，導致地面高程損失而直接形成地面沉降。同時，由於工程建設導致的土地利用類型變化而形成的不透水地面，嚴重阻礙了降水的地面滲水率，不能快速有效地形成地下水補給，間接影響了地下水位的提升，加劇了地面沉降的發生。不同類型的地面植被，需要不同的灌溉條件，也會對地下水位的變化產生影響。
4）在北京區域地下水超採不嚴重的地區引起的水位趨勢性變化中，包含了地震前兆信息，也就意味著地下水位降低有可能會誘發地震的發生。同時，離震源越近的地區地下水位變化的可能性越大，可能的情況是地震的發生與地下水位變化具有密切聯繫，而地下水位變化又與地面沉降的發生密切相關。
3 結束語
以北京市城中心的8個區為例，基於PS-InSAR技術獲取了時間跨度為2003年6月至2010年9月的23個時間段PS點地面沉降數據，以及包括地下水、可壓縮層厚度等多源數據，建立了以格網為單元的地理探測器模型，精細地探測了地面沉降分布的影響因素，並深入分析了各影響因素間交互作用產生的影響力。研究發現：
1）單因子探測結果。對地面沉降影響力較大的因子依次是地下水位、可壓縮層厚度、坡度、土地利用類型、土壤類型、高程、距大地震源距離、距鐵路距離，其中，除了地下水位和可壓縮層厚度兩個因子對地面沉降產生顯著的作用外，其餘要素產生的作用強度顯著減弱。
2）因子交互探測結果。作用強度較弱的單因子與地下水位或可壓縮層厚度分別相互疊加後，其綜合作用強度均表現為增強作用，大部分表現為非線性增強，極少部分表現為雙因子增強。
隨著地下水的過度開採與城市化進程的迅猛發展，中國的許多城市和地區出現了嚴重的地面沉降現象，帶來了一系列的地質環境問題。本文主要聚焦地面沉降分布空間分布特徵與其影響因素差異，其研究結果可為城市地面沉降防控提高科學的參考依據。本研究是將空間統計學方法應用於PS-InSAR大數據的一次有益的嘗試，但目前還處於探索階段，在後續研究中將從地面沉降時空變化特徵來探測其影響因子的交互作用變化特徵。
致謝：感謝中國測繪科學研究院張永紅研究員提供的研究區PS-InSAR數據以及深圳大學黃金彩博士對研究區PS-InSAR數據的整理與預處理。
作者簡介：
肖巍峰（1977-），男，湖南湘潭人，講師，博士，主要從事三維點雲數據處理及表面重構、地理空間數據分析等研究工作。
E-mail:472364382@qq.com
通訊作者：楊文濤（1988-），男，湖北孝感人，講師，博士，碩士生導師，主要從事地理空間數據分析與挖掘理論、方法與技術等研究工作。
E-mail:yangwentao8868@126.com