郭小萌，宮輝力，朱 鋒，郭 琳，沈媛媛

(1．首都師范大學(xué)資源環(huán)境與旅游學(xué)院，北京 100048;2．北京市水文地質(zhì)工程地質(zhì)大隊(duì)，北京 100195)

GM(1，1)模型是鄧聚龍教授于20世紀(jì)80年代前期提出的用于控制和預(yù)測(cè)的新理論、新技術(shù)，目前已廣泛地應(yīng)用于各個(gè)學(xué)科領(lǐng)域，并取得了顯著成就[1]。GM(1，1)模型通過鑒別系統(tǒng)因素之間發(fā)展趨勢(shì)的相異程度，即進(jìn)行關(guān)聯(lián)分析，對(duì)原始數(shù)據(jù)進(jìn)行處理以尋找系統(tǒng)的變化規(guī)律，生成規(guī)律性強(qiáng)的數(shù)據(jù)序列，進(jìn)而建立微分方程，預(yù)測(cè)事物未來發(fā)展趨勢(shì)。

地面沉降是指在自然和人為因素作用下，由于地殼表層土體壓縮而導(dǎo)致區(qū)域性地面標(biāo)高降低的一種環(huán)境地質(zhì)現(xiàn)象，是地質(zhì)環(huán)境系統(tǒng)破壞所導(dǎo)致的惡果，是一種不可補(bǔ)償?shù)挠谰眯原h(huán)境和資源損失，對(duì)資源利用、環(huán)境保護(hù)、經(jīng)濟(jì)發(fā)展、城市建設(shè)和人民生活構(gòu)成了嚴(yán)重威脅[2]。

地面沉降影響因素復(fù)雜多樣，主要包括地下水開采、動(dòng)靜載荷、新構(gòu)造運(yùn)動(dòng)等，是一個(gè)包含土層-地質(zhì)作用-人類活動(dòng)的復(fù)雜系統(tǒng)，其中既包含已知信息又包含未知信息，因此可以將其概化為一個(gè)不具有物理原型的灰色系統(tǒng)，進(jìn)而運(yùn)用GM(1，1)模型對(duì)其進(jìn)行建模和預(yù)測(cè)[3]。

1 基于矩陣的GM(1，1)模型簡(jiǎn)介

傳統(tǒng)的GM(1，1)沉降預(yù)測(cè)模型都是基于單點(diǎn)的運(yùn)算，無法獲取整個(gè)區(qū)域的變化趨勢(shì)。本文利用Matlab對(duì)GM(1，1)預(yù)測(cè)模型進(jìn)行改進(jìn)，將原有的單點(diǎn)模型擴(kuò)展為矩陣模型，從而掌握區(qū)域地面沉降的發(fā)展規(guī)律，預(yù)測(cè)地面沉降的發(fā)展趨勢(shì)。

針對(duì)矩陣各個(gè)元素，按以下步驟建立模型[4]:

(1)對(duì)矩陣內(nèi)各元素對(duì)應(yīng)的序列進(jìn)行級(jí)比平滑檢驗(yàn)，由P(k)=X(k-1)/X(k)，若P(k)?(0．1343，7．389)，則原始序列是平滑的，可用作預(yù)測(cè)，再對(duì)序列作級(jí)比界區(qū)檢驗(yàn)，若P(k)?(e-2/(n+1)，e2/(n+1))(n為序列長(zhǎng)度)，則原始序列可以獲得精度較高的GM(1，1)模型。

(2)采用累加的方式，把原始時(shí)間序列轉(zhuǎn)換為一個(gè)隨時(shí)間變化趨勢(shì)明顯的時(shí)間序列，設(shè)原始時(shí)間序列:

其累加生成的序列為:

則有:

按X(1)序列建立時(shí)間函數(shù)為:

(4)將參數(shù)a、u帶入時(shí)間函數(shù)，然后對(duì)X(1)求導(dǎo)還原得到:

(5)后驗(yàn)差檢驗(yàn)。即先計(jì)算實(shí)測(cè)數(shù)據(jù)離差:

及殘差的離差:

再計(jì)算后驗(yàn)差之比c=s2/s1以及小誤差概率p=

(6)基于矩陣原理，按上述5個(gè)步驟對(duì)矩陣內(nèi)所有元素建立模型。

2 實(shí)例分析

2.1 研究區(qū)概況

北京市位于華北平原西北邊緣，全市總面積16807．8km2，包括 18 個(gè)區(qū)縣，其中平原區(qū) 6236km2，占總面積的38%。地勢(shì)西北高、東南低，由山前向南、東南傾斜，地貌分為西部山區(qū)、北部山區(qū)和東南平原三塊。東南部是由大清河、永定河、北運(yùn)河、潮白河、薊運(yùn)河沖洪積作用形成的平原。

近年來，北京市水資源需求不斷增加，許多地區(qū)由于過量開采地下水，造成地下水位大幅下降，并直接導(dǎo)致地面沉降范圍不斷擴(kuò)大。截至2009年底，最大年沉降量達(dá)到137．51mm，最大累計(jì)沉降量達(dá)到1163mm[5]。北京市平原早期沉降區(qū)有東八里莊—大郊亭沉降區(qū)，朝陽區(qū)來廣營(yíng)沉降區(qū)，昌平沙河—八仙莊沉降區(qū)，順義平各莊沉降區(qū)，大興區(qū)南部榆墊—禮賢沉降區(qū)，近年來還形成了多個(gè)新沉降區(qū)[6]。地面沉降不僅引起樓房墻壁開裂、地基下沉、地下管道工程損壞，而且降低了建筑物抗震能力、并使得測(cè)量水準(zhǔn)點(diǎn)失準(zhǔn)[7]。

2.2 模型構(gòu)建及驗(yàn)證

本文采用的原始數(shù)據(jù)為北京市平原區(qū)23景ENVISAT ASAR影像(2003～2009年)及30個(gè)水準(zhǔn)測(cè)量點(diǎn)數(shù)據(jù)，利用PS-InSAR技術(shù)獲取得到該地區(qū)的各年累積沉降量，在 Matlab平臺(tái)建立基于矩陣的GM(1，1)模型，采用水準(zhǔn)測(cè)量數(shù)據(jù)進(jìn)行精度驗(yàn)證，并進(jìn)行區(qū)域地面沉降趨勢(shì)的預(yù)測(cè)分析。

級(jí)比平滑檢驗(yàn)表明，針對(duì)矩陣內(nèi)所有元素建立的模型，P(k)?(0．8744，1．132)?(0．1343，7．389)，級(jí)比界區(qū)檢驗(yàn)表明，P(k)的界區(qū) =(0．7788，1．284)，由P(k)?(0．8744，1．132)?(0．7788，1．284)，說明本次實(shí)驗(yàn)數(shù)據(jù)可以獲得精度較高的GM(1，1)模型。

在后驗(yàn)差檢驗(yàn)中，需去除模型精度不達(dá)標(biāo)的元素[8]。根據(jù)后驗(yàn)比c和小誤差概率p，得到模型精度的指標(biāo)見表1。

表1 模型精度指標(biāo)Table1 Precision index of the model

本文將精度為“合格”以上的模型都用作預(yù)測(cè)。統(tǒng)計(jì)表明，對(duì)于矩陣內(nèi)所有元素建立的模型，小誤差概率p最大值為1，最小值為0．667，后驗(yàn)差之比c最大值為0．408，最小值為0．111，不合格的模型有4個(gè)，占總數(shù)的0．11‰(表2)。

利用模型計(jì)算2003～2009年北京平原區(qū)的累積沉降值(以2002年為基準(zhǔn)開始累積)及多年平均沉降量，并與實(shí)際監(jiān)測(cè)值進(jìn)行比較。其分布特征見圖1～2。

圖1～2表明:實(shí)測(cè)結(jié)果與模擬結(jié)果在空間上具有高度的相關(guān)性，最大沉降量相對(duì)誤差不超過12%，沉降區(qū)空間分布范圍基本一致。利用ArcMap在圖中隨機(jī)選取100個(gè)點(diǎn)，對(duì)這些點(diǎn)的多年平均實(shí)測(cè)值與模擬值進(jìn)行獨(dú)立樣本t檢驗(yàn)，檢驗(yàn)結(jié)果表明，當(dāng)置信度為95%時(shí)，P=0．413，說明實(shí)測(cè)值與模擬結(jié)果無顯著性差異，兩者的擬合曲線見圖3。

為進(jìn)一步驗(yàn)證模型精度，根據(jù)區(qū)域內(nèi)水準(zhǔn)測(cè)量點(diǎn)數(shù)據(jù)，計(jì)算得到各點(diǎn)處2010年的累積沉降量，并與模型預(yù)測(cè)結(jié)果進(jìn)行對(duì)比分析(圖4)。統(tǒng)計(jì)表明，30個(gè)水準(zhǔn)測(cè)量點(diǎn)中，有28個(gè)點(diǎn)的實(shí)測(cè)結(jié)果與預(yù)測(cè)結(jié)果誤差不超過5mm，對(duì)二者建立回歸分析，相關(guān)系數(shù)r=0．99，說明實(shí)測(cè)結(jié)果與預(yù)測(cè)結(jié)果存在顯著性相關(guān)。

表2 模型運(yùn)算結(jié)果表(僅列部分元素)Table2 Calculate results of the model(part of the elements only)

圖1 研究區(qū)地面沉降實(shí)測(cè)結(jié)果Fig.1 Measured results of land subsidence of the study area

圖2 研究區(qū)地面沉降模擬結(jié)果Fig.2 Similation results of land subsidence of the study area

綜上所述，應(yīng)用基于矩陣的GM(1，1)模型可以分析區(qū)域地面沉降趨勢(shì)，且具有較高精度。

2.3 模型預(yù)測(cè)

利用基于矩陣的GM(1，1)模型，計(jì)算得到北京市平原區(qū)2011～2015年地面沉降預(yù)測(cè)結(jié)果(圖5)。

結(jié)合圖1、圖2及圖5，可以得出以下結(jié)論:

圖3 實(shí)測(cè)值與模擬值擬合曲線Fig.3 Match of measured value and simulation values

圖4 水準(zhǔn)測(cè)量結(jié)果與模型預(yù)測(cè)結(jié)果對(duì)比圖Fig.4 Comparison of measured values and simulation values

圖5 研究區(qū)地面沉降預(yù)測(cè)結(jié)果(2011年，2015年)Fig.5 Prediction results of land subsidence of the study area(in 2011 and 2015)

(1)地面沉降范圍逐漸擴(kuò)大。平原區(qū)從中部幾個(gè)小沉降區(qū)“擴(kuò)散”為一個(gè)范圍大、沉降嚴(yán)重的集中沉降區(qū)，2004年北京市平原區(qū)已有3218km2地區(qū)呈現(xiàn)地面沉降，到 2011年該數(shù)據(jù)達(dá) 4730km2，到 2015年有4900km2的地區(qū)呈現(xiàn)地面沉降，占平原區(qū)總面積的93．3%，其中，年均沉降5～10mm的地區(qū)最廣，占沉降區(qū)總面積的45．1%;年均沉降0～5mm和10～20mm的分別占27%和24．8%;年均沉降20～30mm的占2．9%;沉降大于30mm的占0．2%。據(jù)北京市地面沉降危害分區(qū)參考標(biāo)準(zhǔn)，城市建設(shè)發(fā)展相對(duì)較緩的郊區(qū)縣，年沉降速率30～50mm/a區(qū)域?yàn)榈孛娉两蛋l(fā)育較嚴(yán)重區(qū)，大于50mm/a為嚴(yán)重區(qū);市區(qū)內(nèi)及重點(diǎn)城鎮(zhèn)區(qū)，年沉降10～30mm/a區(qū)域?yàn)榈孛娉两蛋l(fā)育較嚴(yán)重區(qū)，大于30mm/a區(qū)域?yàn)閲?yán)重區(qū)。由此可見，北京市平原區(qū)地面沉降較嚴(yán)重區(qū)和嚴(yán)重區(qū)面積不斷在擴(kuò)大，已由2004年的93km2擴(kuò)大到2011年的571km2和2015年的1281km2，地面沉降問題日益嚴(yán)重，覆蓋面涉及海淀、朝陽、昌平、順義、通州、大興等區(qū)縣。

(2)局部地區(qū)地面沉降程度日益嚴(yán)重。2004年全區(qū)最大沉降量?jī)H 66．80mm，2011 年為 248．14mm，2015年達(dá)421．37mm，其中通州區(qū)徐辛莊附近沉降最嚴(yán)重，2011年、2015年單年沉降量分別36．2mm和37．9mm，多年平均沉降量達(dá)35．1mm/a。

(3)地面沉降趨勢(shì)與地下水位下降趨勢(shì)呈現(xiàn)正相關(guān)。已有研究表明，超量開采地下水是引起北京市地面沉降的主要成因[9]，近年來北京市每年開采25×108～27×108m3地下水，年超采約 2×108m3，加上1999年以來的連年干旱，地下水資源消耗過度、補(bǔ)給不足，致使水位持續(xù)下降，地面沉降形勢(shì)愈來愈嚴(yán)峻，在經(jīng)濟(jì)發(fā)展迅速的城區(qū)和近郊，表現(xiàn)尤為顯著。以朝陽東郊牛場(chǎng)625孔、昌平立水橋水文公司C205-A孔為例(圖6)，地面沉降變化與地下水位動(dòng)態(tài)變化趨勢(shì)基本一致:從2004到2011年，各孔地下水位基本處于持續(xù)下降狀態(tài)，地面沉降也以較快的速度發(fā)展，累積沉降量不斷增加。

圖6 625孔、C205-A孔處累積沉降量與水位變化關(guān)系Fig.6 Relation between accumulated subsidence and water level in observation wells 625 and C205-A

3 結(jié)論與討論

本文以北京市平原區(qū)為研究區(qū)，利用基于矩陣的GM(1，1)模型對(duì)區(qū)域地面沉降進(jìn)行預(yù)測(cè)，模擬結(jié)果顯示，地面沉降的空間分布與PS-InSAR處理結(jié)果一致;采用獨(dú)立樣本t檢驗(yàn)方法分析PS-InSAR處理結(jié)果與模型實(shí)測(cè)值，二者沒有顯著性差異;回歸分析表明，地面水準(zhǔn)點(diǎn)實(shí)測(cè)值與模型預(yù)測(cè)值具有密切的相關(guān)性。因此，利用基于矩陣的GM(1，1)模型進(jìn)行區(qū)域地面沉降預(yù)測(cè)的精度較高，具有較好的可行性。

地面沉降的誘發(fā)因素有天然的及人為的，其形式及性質(zhì)多種多樣，在時(shí)間和空間上存在一定的不確定性，如果針對(duì)大空間尺度的長(zhǎng)期預(yù)測(cè)，使用基于矩陣的GM(1，1)模型可能會(huì)存在一定的誤差[10]。欲解決此問題，需要結(jié)合其他地理學(xué)方法對(duì)預(yù)測(cè)結(jié)果進(jìn)行校正，以提高模型的預(yù)測(cè)精度。

[1] 劉賢趙，張安定，李嘉竹．地理學(xué)數(shù)學(xué)方法[M]．北京:科學(xué)出版社，2009．[LIU X Z，ZHANG A D，LI J Z．Math Method in Geography[M]．Beijing:China Science Press，2009．(in Chinese)]

[2] 閻世駿，劉長(zhǎng)禮．城市地面沉降研究現(xiàn)狀與展望[J]．地學(xué)前緣，1996，3(1):93-98．[YAN S J，LIU C L．Status And Prospect of Urban Land Subsidence[J]．Earth Science Frontiers，1996，3(1):93-98．(in Chinese)]

[3] 宮相霖，董榮鑫．灰色系統(tǒng)在地面沉降分析中的應(yīng)用[J]．上海地質(zhì)，2003(3):16-21．[GONG X L，DONG R X．The Application of Grey System in the Land Subsidence Analysis[J]．Shanghai Geology，2003(3):16-21．(in Chinese)]

[4] 周振民，趙明亮，李玲．GM(1，1)模型在灤河下游地區(qū)地下水位預(yù)測(cè)中的應(yīng)用[J]．中國(guó)農(nóng)村水利水電，2011(2):50-52．[ZHOU Z M，ZHAO M L，LI L．The Application of GM(1，1)Model to the Prediction ofUnderground Water Levelin the Downstream of the Luanhe River[J]．China Rural Water and Hydropower，2011(2):50-52．(in Chinese)]

[5] 楊艷，賈三滿，王海剛．北京平原區(qū)地面沉降現(xiàn)狀及發(fā)展趨勢(shì)分析[J]．上海地質(zhì)，2010，31(4):23-28．[YANG Y，JIA S M，WANG H G．The Status and Development of Land Subsidence in Beijing Plain[J]．Shanghai Geology，2010，31(4):23-28．(in Chinese)]

[6] 賈三滿，王海剛，趙守生，等．北京地面沉降機(jī)理研究初探[J]．城市地質(zhì)，2007，2(1):20-26．[JIA S M，WANG H G，ZHAO S S，et al．A Tentative Study of the Mechanism of land Subsidence in Beijing[J]．City Geology，2007，2(1):20-26．(in Chinese)]

[7] 劉予．北京市地面沉降區(qū)含水層和壓縮層組劃分及地面沉降自動(dòng)監(jiān)測(cè)系統(tǒng)[D]．長(zhǎng)春:吉林大學(xué)，2004．[LIU Y． Divided Water-bearing Zones and Compressible Zones of Beijing Land Subsidence Area and Land Subsidence Automatic Monitoring System[D]．Changchun:Jilin University，2004．(in Chinese)]

[8] 黃景銳，胡安焱，張煥楚，等．基于非等間距序列GM(1，1)模型的地下水溫度預(yù)測(cè)[J]．水文地質(zhì)工程地質(zhì)，2013，40(1):48-52．[HUANG J R，HU A Y，ZHANG H C，etal．Prediction ofshallow groundwater temperature based on GM(1，1)model for non-equidistant sequence[J]．Hydrogeology ＆Engineering Geology，2013，40(1):48-52．(in Chinese)]

[9] 杜釗鋒，宮輝力，王灑，等．短時(shí)空基線PS-InSAR在北京地面沉降監(jiān)測(cè)中的應(yīng)用[J]．水文地質(zhì)工程地質(zhì)，2012，39(5):116-120．[DU Z F，GONG H L，WANG S，et al．Application of small spatio-temporal baseline PS-InSAR to the study of land subsidence in Beijing[J]．Hydrogeology ＆ Engineering Geology，2012，39(5):116-120．(in Chinese)]

[10] 潘網(wǎng)生．基于灰色—馬爾可夫模型的西安市主要城市建設(shè)用地需求量預(yù)測(cè)[D]．西安:長(zhǎng)安大學(xué)，2008．[PAN W S．Based on the Theory of Gray-Markov Forecasting Model for Main Land-Demand of Xi'an City Foundation[D]． Xi'an:Chang'an University，2008．(in Chinese)]