石鵬遠，余 潔，朱 琳，王彥兵

(1.首都師范大學城市環(huán)境過程與數(shù)字模擬國家重點實驗室培育基地，北京 100048；2.首都師范大學資源環(huán)境與旅游學院, 北京 100048；3.首都師范大學三維信息獲取與應用教育部重點實驗室，北京 100048)

0 引言

地面沉降作為一種在自然因素和人為因素下形成的地表垂直下降現(xiàn)象[1]，能夠對城市基礎設施和建筑物(如農(nóng)田、住宅、道路和輸電線路等)造成破壞，使政府的財政和社會成本增加[2-3]。進行地面沉降危險性評估可以為城市規(guī)劃、防災減災以及地質環(huán)境的恢復治理提供基礎資料和科學依據(jù)，以避免或減輕地面沉降造成的損失[4]。

對地面沉降危險性進行評估的模型有多種。一些學者[5-8]采用土水模型法，從水位模型和土體力學模型著手，對地面沉降危險性進行評估，土水模型法的優(yōu)點是能夠基于精確的水文、土壤數(shù)據(jù)和地面沉降成因機理模擬預測地面沉降，但該方法需要大量的實地勘測數(shù)據(jù)以使土水模型的模擬條件與實際的水文、土壤狀態(tài)吻合。另一些學者采用加權綜合評價法[9-12]，綜合考慮各因子對地面沉降的影響程度，把各指標的優(yōu)劣綜合起來，用一個無量綱的評估分值加以集中，以表示地面沉降危險性，其中部分學者[13-16]采用因子疊加模型將與地面沉降相關的因子進行指標打分并進行多個柵格數(shù)據(jù)層的數(shù)學運算，按照人工賦予的疊加函數(shù)計算得到評估結果；另一部分學者[17-21]采用模糊層次分析模型，建立模糊評價矩陣為各項評估指標賦予不同的權重，計算各指標的評分加權總和并進行評估。加權綜合評價法的優(yōu)點在于靈活性和綜合性強，能夠根據(jù)評價目標有目的地調整參與評估的因子組，適合于對技術、決策或方案進行綜合分析評價和優(yōu)選，是一種較為常用的方法[22-24]，然而在因子選取和指標權重確定的過程中，人為因素的判斷會影響評估結果。在一些研究中，指標的權重由特爾斐法或學者直接確定[25-27]；而采用模糊層次分析模型計算權重時，因子之間的重要性比較也常依靠特爾斐法或學者的專業(yè)判斷[28-29]。這使因子的指標項權重的確定受到較強的人工干預，具有不確定性，進而影響評估模型的準確性。

為了減少指標權重確定過程中人工干預導致的不確定性，本文將地理探測器方法應用于加權綜合評價法構建的地面沉降危險性評估模型中，以因子疊加模型和模糊層次分析模型為例進行改進。改進后的模型將根據(jù)各評估因子在地理位置上的分布與地面沉降分布的關聯(lián)程度，從統(tǒng)計學的角度出發(fā)，定量地確定各項評估指標的權重，以期提高評估模型的準確性。

1 地面沉降危險性評估模型及因子指標

1.1 評估模型

本文共選取三種較常被使用的、便于操作的評估模型，考慮到因子疊加模型有多種運算函數(shù)，選擇了較為常見的最大值因子疊加模型和加權因子疊加模型。第三種模型是模糊層次分析模型。對三種模型分別應用地理探測器進行評估指標權重確定環(huán)節(jié)的改進。

(1)最大值因子疊加模型

最大值因子疊加模型(Maximum Value Factor Overlapping, MVFO)在進行地面沉降危險性評估時，僅使用沉降特征因子指標(如沉降速率，累計沉降量)。對沉降特征因子進行危險性指標打分，并按柵格進行疊加運算，柵格中危險性分值最高的因子指標分值視為最終的沉降危險性估值。這一模型的優(yōu)勢在于參與評估因子的數(shù)量少，因子數(shù)據(jù)易于收集，評估過程簡單快捷，基于沉降特征的評估較為符合地區(qū)近期的沉降狀態(tài)，對高危險性估值的地區(qū)漏檢率低；不足之處在于該模型不考慮推動沉降發(fā)展的誘發(fā)因子，缺少對沉降未來的發(fā)展趨勢的估計，且對于沉降危險性較低的地區(qū)容易產(chǎn)生誤判。危險性評估值計算公式為：

H危險性=max{N1,N2,…,Nn}(1)

式中：H危險性——沉降危險性最終估值；

Nn——第n項因子的指標評分；

n——因子總數(shù)。

(2)加權因子疊加模型

加權因子疊加模型(Weighted Factor Overlapping, WFO)在評估中既考慮沉降特征因子對沉降狀態(tài)的表達刻畫，又考慮沉降誘發(fā)因子對沉降發(fā)展趨勢的預示作用。由于不同的地區(qū)沉降誘發(fā)因子不同，且誘發(fā)因子打分的指標也不同，因此該模型不取指標評分最大值為沉降危險性的最終估值，而是將每一個指標都視為具有單獨的權重，在簡化計算的前提下，各指標被認為權重相等。各因子進行指標打分并單獨成圖，按照柵格進行疊加運算，柵格的危險性最終估值為所有指標的平均分值。該評估模型的優(yōu)勢在于即包含能反映沉降狀態(tài)的沉降特征因子又包含能揭示沉降未來趨勢的沉降誘發(fā)因子，且操作簡單便捷；不足之處在于模型忽視了不同因子與沉降關聯(lián)的差異性，在評估中將各因子視為同等重要，對評估精度產(chǎn)生影響。危險性評估值計算公式為：

(2)

式中：H危險性——沉降危險性最終估值；

Ni——第i項因子的指標評分；

n——因子總數(shù)。

(3)模糊層次分析模型

模糊層次分析模型(Fuzzy Analytic Hierarchy Process, FAHP)是根據(jù)模糊數(shù)學方法與層次分析法建立的評估模型?？紤]到不同因子與地面沉降的關聯(lián)性具有差異，該模型通過不同因子之間的兩兩比較構建各因子對評估對象的隸屬度函數(shù)，量化確定評估中各因子的指標分值所占權重。FAHP模型既解決了判斷矩陣的一致性問題，也解決了解的收斂速度和精度問題，最終求得與實際相符的指標權重排序向量[18]。該模型的優(yōu)勢在于半定量地體現(xiàn)了不同因子作用的差異性，使評估結果更加合理有效；不足之處在于，在確定指標權重的過程中，對各因子進行的人工比較存在著不確定性，對評估精度產(chǎn)生影響。

該方法的模型建立過程如下[30]：

a.通過對因子的兩兩比較，建立優(yōu)先關系矩陣F。矩陣標度取值范圍為0.1～0.9。

(3)

式中：fij——矩陣F中因子i與因子j重要性判斷標度；

s(i)——因子i的重要程度；

s(j)——因子j的重要程度。

fij值越高，說明因子i在與因子j的對比中顯得更為重要。fij為0.5時，兩個比較的因子同等重要。矩陣F是一個模糊互補矩陣，即fij=1-fji。

b.計算模糊一致性判斷矩陣R。

c.通過和行歸一法計算排序向量W0。

W0=(w1,w2,…,wn)T=

式中：W0——計算得到的排序向量；

wn——向量W0在第n行的值；

n——因子總數(shù)；

rij——矩陣R在i行j列的值。

e.以排序向量W0作為迭代初值V0，進一步求解精度較高的排序向量Vk。

(5)

式中：Vk+1,n——特征向量Vk+1的第n行值。所得向量Wk=Vk+1，即為排序向量，迭代結束。

f.結合加權綜合評價法模型計算最終的評估結果：

(6)

式中：H——最終評估指數(shù)；

Wi——計算出的第i項指標權重；

Ni——第i項因子的指標評分；

n——因子總數(shù)。

1.2 地理探測器對模型的改進

1.2.1地理探測器

地理探測器(Geographical Detector, GD)是由王勁峰等[31]所開發(fā)的一種新的統(tǒng)計學方法，該方法能夠探測某一地理屬性的空間分異性，并揭示其背后的驅動因子，例如探測某一地理屬性Y和其解釋因子X之間的關系。其核心思想基于的假設是：若某個因變量受到某一自變量的重要影響，那么因變量和自變量的空間分布應該具有相似性[32]。該方法的優(yōu)勢在于沒有過多的假設條件與約束，能有效地克服傳統(tǒng)統(tǒng)計分析方法處理類別變量時的局限性[33]。

地理探測器由四個探測器組成：因子探測器、交互作用探測器、風險區(qū)探測器、生態(tài)探測器，分別用于探測因子X對地理屬性Y分布的影響程度、兩個不同的因子X1和X2共同作用下對Y分布的影響、因子X的兩個子區(qū)域屬性值均值差異是否顯著、X1和X2對Y的空間分布影響的差異是否顯著。

本文使用因子探測器參與評估模型中因子指標的權重確定。因子探測器計算結果在本文中稱為“因子驅動力”，其得到的計算結果用q值度量，q值越高，說明自變量X對因變量Y的影響越強；反之，則說明自變量X對因變量Y的影響較弱。q值的計算公式為：

(7)

式中：h=1,2,…,L——因變量Y或自變量X內部的分類或分區(qū)；

Nh和N——分別為層h和全區(qū)的單元數(shù)；

SSW和SST——分別為層內方差和與全區(qū)總方差；

q——值域為[0,1]。

1.2.2模型的改進

(1)地理探測器——最大值因子疊加模型

地理探測器——最大值因子疊加模型(Geographical Detector-Maximum Value Factor Overlapping, GD-MVFO)對MVFO模型沉降特征因子的指標分值比較過程進行了修改。根據(jù)地理探測器探測得到的各因子q值計算得到因子的指標權重。在選取危險性最高的指標時，各因子的比較分值等于因子的指標分值與指標權重的乘積。比較分值最高的因子，其指標分值將作為沉降危險性最終估值。

與MVFO模型相比，改進后的GD-MVFO模型在評估中不僅考慮沉降特征因子的指標分值大小，還考慮不同因子的影響差異，而不是將各因子視為同等重要。GD-MVFO模型危險性評估值計算公式為：

(8)

Cm=max{N1w1,N2w2,…,Nnwn}(9)

式中：wi——第i項因子的指標權重；

qi——第i項因子的因子驅動力；

n——因子總數(shù)；

Nn——第n項因子的指標分值；

Cm——因子m的比較分值；

wm——因子m的指標權重；

H危險性——沉降危險性最終估值。

(2)地理探測器——加權因子疊加模型

地理探測器——加權因子疊加模型(Geographical Detector-Weighted Factor Overlapping, GD-WFO)使用探測到的q值計算得到各因子的指標權重。與WFO模型相比，改進后的GD-WFO模型不僅綜合了各因子的指標打分值，還考慮了不同因子的影響能力差異，不再將各因子視為同等重要。沉降危險性最終估值等于各因子對應的指標分值與指標權重的乘積之和。GD-WFO模型的公式如下：

(11)

式中：wi——第i項因子的指標權重；

qi——第i項因子的因子驅動力；

n——因子總數(shù)；

Nj——因子j的指標分值；

wj——因子j的指標權重；

H危險性——沉降危險性最終估值。

(3)地理探測器——模糊層次分析模型

地理探測器——模糊層次分析模型(Geographical Detector-Fuzzy Analytic Hierarchy Process,GD-FAHP)在原有的FAHP模型的基礎上，對生成因子優(yōu)先關系矩陣的環(huán)節(jié)進行改進。根據(jù)對因子驅動力的探測結果而不是對因子的人工對比估值計算產(chǎn)生優(yōu)先關系矩陣F。

GD-FAHP模型的建立過程如下：

a.首先根據(jù)探測出的q值計算生成因子驅動力比較矩陣Q，公式如下：

(13)

式中：Qij——因子i與因子j的因子驅動力比較值；

q——因子驅動力。

b.根據(jù)因子力比較矩陣Q計算生成優(yōu)先關系矩陣F：

(14)

式中：Fij——因子i與因子j重要性判斷標度；

Qij——因子i與因子j的因子驅動力比較值；

Qji——因子j與因子i的因子驅動力比較值。

c.將矩陣F帶入模糊層次模型中，求得權重向量Wk。

1.3 因子指標

1.3.1評估因子選取

根據(jù)北京市地面沉降發(fā)生機理[34-35]，綜合考慮地面沉降災害對社會經(jīng)濟可持續(xù)發(fā)展帶來的嚴重影響，依照少而精、實用易操作、數(shù)據(jù)易獲取等原則，參照地質災害評估的行業(yè)規(guī)范、沉降因素研究文獻和以往的地面沉降危險性評估文獻，本文在進行地面沉降危險性評估時將地面沉降特征因子、地面沉降誘發(fā)因子作為評估因子。

(1)沉降特征因子

參照中華人民共和國國土資源部發(fā)布的《地質災害危險性評估規(guī)范》[36]對地面沉降調查的內容規(guī)定及地面沉降發(fā)育程度分級表，本文選擇地面沉降速率和地面累計沉降量兩個因子作為沉降特征因子參與地面沉降危險性評估。

本文的地面沉降速率，采用基于永久散射體合成孔徑雷達干涉測量技術和TerraSAR-X的微波遙感影像計算得到研究區(qū)2010年年均地面沉降速率。地面累計沉降量的選取，則根據(jù)其它文獻中[37]已發(fā)布的北京市歷史累計沉降量數(shù)據(jù)，采用1955～2010年的沉降值。沉降速率及累計沉降量分布見圖1(a)和(b)。

(2)沉降誘發(fā)因子

有學者[38]的研究顯示，北京市地下水開采對地面沉降影響最大，沉降中心與地下水位降落漏斗區(qū)高度吻合，是地面沉降的主要貢獻因素。

北京地區(qū)的地面沉降分層沉降量被認為與對應層位上的黏性土占比呈正比例關系，其空間分布特征及變化趨勢與平原區(qū)的地層結構及可壓縮黏性土層厚度具有很好的一致性，整體由北西向的單一結構區(qū)向南東方向的多層結構區(qū)擴張。沉降速率大于50 mm/a的沉降區(qū)大多分布在黏性土層厚度大于100 m的地區(qū)，幾大沉降中心與黏性土層厚度較大地區(qū)吻合較好。

還有學者[39]認為降水會通過對地下水補給和土壤含水量的變化間接影響到地面沉降發(fā)育狀態(tài)。

根據(jù)以上原因，本文選取地下水水位變化速率、地下水水位標高、可壓縮層厚度和年均降水量四個因子作為沉降誘發(fā)因子參與評估，見圖1(c)～(f)。

1.3.2指標評分準則

目前在地面沉降危險性評估工作中，不同地區(qū)對各因子的指標評分標準存在差異。本文參照國家發(fā)布的地質災害危險性評估規(guī)范文件[36]、危險性評估指導專著[1]及北京地區(qū)地面沉降危險性評估文獻[4,15,18,21]，結合研究區(qū)各因子的分布情況，制定指標評分準則(表1)。根據(jù)準則對各評估因子賦予危險性指標分值，并參與后續(xù)計算。

2 實驗及驗證

2.1 實驗區(qū)概況

實驗區(qū)位于北京市平原地區(qū)，地處東經(jīng)116度20分至116度40分，北緯39度45分至40度5分之間(圖2)。實驗區(qū)面積約為385.51 km2，覆蓋北京市東城區(qū)大部、朝陽區(qū)南部，并與西城區(qū)、豐臺區(qū)、大興區(qū)和通州區(qū)有少許交疊。實驗區(qū)西側的西城區(qū)和東城區(qū)地面沉降發(fā)育程度較低，西單到東單地區(qū)被認為是北京市最早發(fā)現(xiàn)沉降的地區(qū)[37](20世紀30年代中期)。東側涵蓋了朝陽來廣營—金盞和東八里莊—通州城區(qū)—黑莊戶兩個沉降帶的一部分，其中東郊八里莊大郊亭和來廣營兩個地區(qū)的沉降產(chǎn)生于20世紀70年代初到80年代初期工業(yè)發(fā)展需求導致的對地下水的大量開采；通州沉降區(qū)形成于80年代初至20世紀末，并在21世紀初邁入快速發(fā)展階段[40]。到目前為止，實驗區(qū)東側部分地區(qū)平均年沉降速率超過50 mm/year，局部地區(qū)沉降速率超過100 mm/year。

表1 指標評分表

圖1 因子分布圖Fig.1 Distribution of factors

2.2 評估實驗

2.2.1因子驅動力探測

在研究區(qū)中以500 m為間距生成1 545個采樣點(圖3)，提取各因子位于采樣點的值。評估實驗的假設條件是基于2010年以前的因子分布數(shù)據(jù)對研究區(qū)2010年至2015年沉降發(fā)展做出估計。從評估者的角度出發(fā)，2010年以前的因子分布數(shù)據(jù)為已知量，可以作為自變量參與探測計算；2010至2015年的沉降量為未知量，不能作為因變量進行探測。根據(jù)沉降速率在時間上的連續(xù)性，將2010年沉降速率作為替代的因變量Y，所有因子作為自變量X進行因子驅動力(q值)計算(表2)。表中的p值為顯著性檢驗值。

圖3 采樣點分布Fig.3 Distribution of sampling points

探測因子q值p值沉降速率0.999 70.000 0累計沉降0.013 80.003 6地下水水位變化速率0.261 00.000 0地下水水位標高0.025 40.000 0可壓縮層厚度0.018 90.000 0年均降水0.010 00.032 6

2.2.2實驗結果

(1)未改進的評估模型

根據(jù)模型公式和指標評分表，MVFO、WFO、FAHP模型的評估結果見圖4地面沉降危險性評估結果(a)～(c)。其中，F(xiàn)AHP模型的計算矩陣和最終權重如下所示：

a.優(yōu)先關系矩陣F：

b.模糊一致性矩陣R：

c.排序向量W0：

d.互反型判斷矩陣E：

e. 權重向量Wk：

(2)改進后的模型

根據(jù)模型公式和指標評分表，GD-MVFO、GD-WFO、GD-FAHP模型的評估結果見圖4地面沉降危險性評估結果(d)～(e)。

GD-MVFO模型中評估因子的指標權重為：

[0.986 4 0.013 6]T

GD-WFO模型中評估因子的指標權重為：

GD-FAHP模型計算矩陣和最終權重如下所示：

a.因子驅動力比較矩陣Q：

b.優(yōu)先關系矩陣F：

c.將矩陣F帶入模糊層次模型后，求得最終權重向量Wk：

圖4 地面沉降危險性評估結果Fig.4 Assessment results of land subsidence hazard

2.3 結果驗證

2.3.1InSAR監(jiān)測的沉降危險性

圖5 InSAR監(jiān)測沉降速率及沉降危險性分布Fig.5 Actual subsidence rate and subsidence hazard distribution

根據(jù)InSAR監(jiān)測的2010年至2015年沉降值計算得到年均沉降(圖5(a))，并參照危險性的指標評分表中對沉降速率的危險性評語，得到作為驗證標準的沉降危險性分布(圖5(b))，該圖將作為衡量各評估模型利用2010年因子數(shù)據(jù)對2011～2015年(相對于因子數(shù)據(jù)的未來五年)的沉降危險性進行評估的結果準確性參照標準。

以250 m為間距獲取驗證結果的采樣，根據(jù)驗證采樣點提取模型評估結果和作為對照的沉降危險性。

2.3.2模型評估結果與目標評語的匹配系數(shù)

根據(jù)驗證采樣點提取到的評估數(shù)據(jù)計算得到模型評估結果與目標評語的匹配準確度，驗證結果如表3所示，在未改進的模型中，MVFO模型從“輕度”到“重度”五項評語準確度及總體準確度分別為30.15%、40.37%、44.31%、78.67%、100%和60.51%；WFO模型為40.25%、40.91%、44.31%、50.95%、74.26%和53.73%；FAHP模型為63.36%、43.85%、50.00%、49.05%、75.81%和65.05%。在應用地理探測器改進的模型中，GD-MVFO模型從“輕度”到“重度”五項評估準確度及總體準確度分別為49.48%、93.32%、89.82%、83.89%、100.00%和75.23%；GD-WFO為91.67%、63.37%、68.56%、70.85%、97.88%和89.56%；GD-FAHP模型為72.95%、54.28%、69.46、53.79%、89.45%和76.38%。

表3 模型評估結果準確度

2.3.3Kappa系數(shù)

在通常情況下，總體準確度能很好地表示評估結果與目標評語的匹配準確度，但是當某一類目標評估的驗證采樣點占據(jù)多數(shù)時，總體準確度受到該類別評估的影響較大，可能出現(xiàn)評估模型的偶然性高準確度。為進一步驗證各模型的評估效果，使用Kappa系數(shù)進行分析。Kappa系數(shù)計算結果通常在0到1之間，系數(shù)越大表示一致性越高。

計算結果表明，MVFO模型、WFO模型、FAHP模型的Kappa系數(shù)分別為0.485 1、0.405 2、0.525 5。GD-MVFO模型、GD-WFO模型、GD-FAHP模型的Kappa系數(shù)分別為0.663 7、0.842 9、0.665 1。

2.3.4F1-measure函數(shù)

匹配準確度從驗證采樣點目標真實評估值的角度出發(fā)驗證模型評估效果，此外還可以從評估結果的角度出發(fā)驗證評估模型的某一類危險性評語對目標真實評語值的估計能力。為兼顧兩種驗證角度，進一步驗證各模型的評估效果，使用F1-measure函數(shù)進行計算，F(xiàn)1-measure值越高說明評估模型越有效。

計算結果表明，MVFO模型、WFO模型、FAHP模型的F1-measure平均值分別為55.58%、49.43%、53.82%。GD-MVFO模型、GD-WFO模型、GD-FAHP模型的Kappa系數(shù)分別為75.39%、76.27%、63.74%。

3 結論

本文應用地理探測器，分別對MVFO模型、WFO模型和FAHP模型進行了改進，將改進模型的地面沉降危險性評估與三種未改進的評估模型進行了結果驗證和對比分析，得到如下結論：

(1)改進后的評估模型在評估因子的指標權重確定環(huán)節(jié)中使用地理探測器進行因子驅動力的探測，根據(jù)因子與未來發(fā)育的地面沉降在空間分布上的相似關系，定量地計算并確定各項因子指標的權重值，替代了人工打分在指標權重確定環(huán)節(jié)中的作用，減少了評估工作因人工干預作用的差異而產(chǎn)生不同結果的可能性。

(2)與三種未改進的評估模型相比，應用地理探測器方法改進的GD-MVFO模型、GD-WFO模型、GD-FAHP模型與用作驗證的地面沉降分布擬合度更好，表明將地理探測器應用于地面沉降的危險性評估是一種可行的方法，提升了評估模型的性能。