基于AR模型的上海地区地面沉降预测分析

书书书第２８卷第５期２００９年１０月大地测量与地球动力学ＪＯＵＲＮＡＬＯＦＧＥＯＤＥＳＹＡＮＤＧＥＯＤＹＮＡＭＩＣＳＶｏｌ．２９Ｎｏ．５　Ｏｃｔ．，２００９　　文章编号：１６７１５９４２（２００９）０５０１２１０５基于ＡＲ模型的上海地区地面沉降预测分析焉建国１）　陈正松２）　罗志才３）　李　琼３）１）武汉大学测绘遥感信息工程国家重点实验室，武汉　４３００７９２）中国地震局地震研究所，武汉　４３００７１３）地球空间环境与大地测量教育部重点实验室，武汉　４３００７９摘　要　利用ＡＲ模型，对上海地区地面沉降作了拟合评估，依据ＡＩＣ推出了ＡＲ系列预报模型。计算结果表明，ＡＲ（４）模型的结果最符合真实值。用ＡＲ（４）模型对上海地区未来１０年的地面沉降值进行了预测。关键词　ＡＲ模型；ＡＩＣ准则；预测；地面沉降；自回归中图分类号：Ｐ２０７　　　　文献标识码：ＡＡＮＡＬＹＳＩＳＡＮＤＰＲＥＤＩＣＴＩＯＮＯＦＬＡＮＤＳＵＢＳＩＤＥＮＣＥＩＮＳＨＡＮＧＨＡＩＢＡＳＥＤＯＮＡＲＭＯＤＥＬＹａｎＪｉａｎｇｕｏ１），ＣｈｅｎＺｈｅｎｇｓｏｎｇ２），ＬｕｏＺｈｉｃａｉ３）ａｎｄＬｉＱｉｏｎｇ３）１）ＳｔａｔｅＫｅｙＬａｂｏｒａｔｏｒｙｏｆＩｎｆｏｒｍａｔｉｏｎＥｎｇｉｎｅｅｒｉｎｇｉｎＳｕｒｖｅｙｉｎｇ，ＭａｐｐｉｎｇａｎｄＲｅｍｏｔｅＳｅｎｓｉｎｇ，ＷｕｈａｎＵｎｉｖｅｒｓｉｔｙ，Ｗｕｈａｎ４３００７９２）ＩｎｓｔｉｔｕｔｅｏｆＳｅｉｓｍｏｌｏｇｙ，ＣｈｉｎａＥａｒｔｈｑｕａｋｅＡｄｍｉｎｉｓｔｒａｔｉｏｎ，ＷｕＨａｎ　４３００７１３）ＫｅｙＬａｂｏｒａｔｏｒｙｏｆＧｅｏｓｐａｃｅＥｎｖｉｏｎｍｅｎｔａｎｄＧｅｏｄｅｓｙ，ＭｉｎｉｓｔｒｙｏｆＥｄｕｃａｔｉｏｎ，Ｗｕｈａｎ　４３００７９Ａｂｓｔｒａｃｔ　ＴｈｅＡＲｍｏｄｅｌｉｓｕｓｅｄｔｏｍａｋｅｆｉｔｔｉｎｇｅｖａｌｕａｔｉｏｎｏｆｌａｎｄｓｕｂｓｉｄｅｎｃｅｉｎＳｈａｎｇｈａｉａｒｅａ．ＡＲｓｅｒｉａｌｐｒｅｄｉｃｔｉｏｎｍｏｄｅｌｉｓｃｅｒｉｖｅｄａｃｃｏｒｄｉｎｇｔｏＡＩＣｌａｗ，ａｎｄｔｈｅｃｏｍｐｕｔｅｄｒｅｓｕｌｔｉｎｄｉｃａｔｅｓｔｈａｔｔｈｅｒｅｓｕｌｔｓｆｒｏｍＡＲ（４）ｍｏｄｅｌａｒｅｅｘａｃｔｌｙｉｎａｃｃｏｒｄａｎｃｅｗｉｔｈａｃｔｕａｌｒｅｓｕｌｔｓａｎｄｗｉｔｈｔｈｉｓｍｏｄｅｌｔｈｅｐｒｅｄｉｃｔｉｏｎｏｆｌａｎｄｓｕｂｓｉｄｅｎｃｅｉｎｆｕｔｕｒｅ１０ｙｅａｒｓｏｆＳｈａｎｇｈａｉａｒｅａｗａｓｂｅｅｎｇｉｖｅｎ．Ｋｅｙｗｏｒｄｓ：ＡＲｍｏｄｅｌ；ＡＩＣｌａｗ；ｐｒｅｄｉｃｔｉｏｎ；ｌａｎｄｓｕｂｓｉｄｅｎｃｅ；ａｕｔｏｒｅｇｒｅｓｓｉｏｎ１　地面沉降概述地面沉降是指区域性地面高程下降的现象，一般是通过测量地面水准点的高程变化而发现的。根据地面沉降的成因，可以将地面沉降分为两大类，即自然的和人为的［１］。通常所说的与研究的地面沉降均指人为因素产生的地面沉降。目前，国内外对地面沉降的理论研究主要集中在对地下水的开采引起地面沉降的机理及依据数理统计建立预测模型［２］，以地下水开采为主要原因的地面沉降模拟和预测始终是地面沉降研究的热点。地面沉降过程包含了影响其变化的各种确定性因素和随机因素的信息，可分为两类，即考虑沉降机理的确定性模型研究和随机统计模型研究。地面沉降模型的研究始终受到高度重视，许多学者致力于这方面的研究，提出了许多计算模型［３－７］。目前，地面沉收稿日期：２００９０１０８基金项目：国家９７３计划（２００７ＣＢ７１４４０５）；国家教育部新世纪优秀人才计划（ＮＣＥＴ－０７－０６３５）作者简介：焉建国，男，１９７９年生，讲师，主要研究方向为卫星定轨与月球重力场．Ｅ－ｍａｉｌ：[e-mail]大地测量与地球动力学２９卷降预测计算模型可分为两类：一类为土水模型，另一类为沉降量与时间关系模型。其中，土水模型包括两个方面：一方面是确定含水层水位与开采（回灌）量之间的关系，即水位模型；另一方面是计算由于含水层水位变化而引起黏性土及含水层本身的变形，即土力学模型。水位模型和土力学模型按照不同的作用关系便形成了地面沉降预测计算模型。沉降量与时间关系模型［８，９］：如生命旋回模型，主要从地面沉降的整个发展过程来考虑，直接由沉降量与时间之间的关系构成；又如泊松旋回模型［１０］、Ｖｅｒｈｕｌｓｔ生物模型和灰色预测模型等。除此之外，神经网络模型也在地面沉降预测中得到应用，用于预测地面沉降的趋势，为地面沉降的综合治理提供理论依据。本文介绍一种新的预测方法———自回归预测模型［４］。２　自回归模型原理［５］２．１　多元线性自回归假设系统有已知的ｍ个自变量（即预测因子）ｕ１，ｕ２，…，ｕｍ，预测对象记为ｙ，则多元线性回归方程为：＾ｙ（ｋ）＝ｂ０＋ｂ１ｕ１（ｋ）＋ｂ２ｕ２（ｋ）＋…＋ｂｍｕｍ（ｋ）（１）其中＾ｙ（ｋ）为预测对象ｙ（ｋ）的回归估计值，ｂ０，ｂ１，…，ｂｍ称为回归系数，ｕ１（ｋ）为前一项值，ｕ２（ｋ）为前第二项值，ｕｍ（ｋ）为前第ｍ项值。回归系数ｂｉ（ｉ＝０，１，…，ｍ）的最优估计值可由最小二乘法原理得到，即ｂｉ的确定满足残差平方和Ｑ达到最小。残差平方和为：Ｑ＝∑ｎｋ＝１［ｙ（ｋ）－＾ｙ（ｋ）］２＝∑ｎｋ＝１｛ｙ（ｋ）－［ｂ０＋ｂ１ｕ１（ｋ）＋ｂ２ｕ２（ｋ）＋…＋ｂｍｕｍ（ｋ）］｝２（２）由微分学中的极值原理知，式（２）满足如下方程组：Ｑ（ｂ０，ｂ１，…，ｂｍ）（ｂ０）＝０Ｑ（ｂ０，ｂ１，…，ｂｍ）（ｂ１）＝０……Ｑ（ｂ０，ｂ１，…，ｂｍ）（ｂｍ）＝０（３）即Ｑ对各回归系数ｂｉ的偏微商为零，由Ｑｂ０得：∑ｎｋ＝１ｙ（ｋ）－［ｂ０＋ｂ１ｕ１（ｋ）＋ｂ２ｕ２（ｋ）＋…＋ｂｍｕｍ（ｋ）］＝０即可得：ｂ０＝ｙ－［ｂ１珔ｕ１＋ｂ２珔ｕ２＋…＋ｂｍ珔ｕｍ（４）其中：珋ｙ＝１ｎ∑ｎｋ＝１ｙ（ｋ），珔ｕｉ＝１ｎ∑ｎｋ＝１ｕｉ（ｋ）　ｉ＝１，２，…，ｍ将式（４）代入式（３）可得下列解ｍ元回归系数的正规线性方程组：ｓ１１ｂ１＋ｓ１２ｂ２＋…＋ｓ１ｍｂｍ＝ｓ１ｙｓ２１ｂ１＋ｓ２２ｂ２＋…＋ｓ２ｍｂｍ＝ｓ２ｙ……ｓｍ１ｂ１＋ｓｍ２ｂ２＋…＋ｓｍｍｂｍ＝ｓｍｙ（５）其中，ｓｉｊ＝ｓｊｉ＝∑ｎｉ＝１［ｕｉ（ｋ）－珔ｕｉ］［ｕｊ（ｋ）－珔ｕｊ］（６）ｓｉｙ＝∑ｎｉ＝１［ｕｉ（ｋ）－珔ｕｉ］［ｙ（ｋ）－珋ｙ］（７）ｓｉｊ称为ｕｉ与ｕｊ之间的相关矩，ｓｉｙ称为ｕｉ与ｙ之间的相关矩。在样本数据给定的前提下，各ｓｉｊ和ｓｉｙ的数值均可算出。因此，正规线性方程组式（３）是以回归系数ｂ０，ｂ１，…，ｂｍ为未知数的ｍ元联立方程组，可用消去法求解。解出ｂ０，ｂ１，…，ｂｍ之后，由式（２）解出ｂ０，从而就可列出多元线性回归方程。如果引进无量纲的单相关系数ｒｉｊ来代替相关矩ｓｉｊ求解正规线性方程，由于ｒｉｊ的量值差异要比ｓｉｊ小，使得舍入误差要小，在计算机上计算更是如此。由于相关矩与相关系数有如下关系：ｒｉｊ＝ｓｉｊｓｉｉｓ槡ｊｊ，ｒｉｙ＝ｓｉｙｓｉｉｓ槡ｙｙ　ｉ＝１，２，…，ｍ因此正规线性方程组也可写成如下形式：ｒ１１ｂ′１＋ｒ１２ｂ′２＋…＋ｒ１ｍｂ′ｍ＝ｒ１ｙｒ２１ｂ′１＋ｒ２２ｂ′２＋…＋ｒ２ｍｂ′ｍ＝ｒ２ｙ　　　　　　　　ｒｍ１ｂ′１＋ｒｍ２ｂ′２＋…＋ｒｍｍｂ′ｍ＝ｒｍｙ（８）其中，ｂ′ｉ＝ｓ槡ｉｉｓ槡ｙｙｂｉ　　ｉ＝１，２，…，ｍ（９）２．２　模型阶数的确定对于＾ｙ（ｋ）＝ｂ０＋ｂ１ｕ１（ｋ）＋ｂ２ｕ２（ｋ）＋…＋ｂｍｕｍ（ｋ），当数据越多时，回归系数ｂｋ越多，得到的数据越精确。然而，在实际中ｂｋ越多有时误差反而越大，所以ｂｋ的选择需要跟数据的多少联系。因此必须找到一个最优的方法。ＡＩＣ准则就是用于识别模型的系数个数［５］：ＡＩＣ（ｋ）＝ｎｌｎ（σ２ｅ）＋２Ｋ（１０）Ｋ为回归系数的个数，也称为阶数，一般为Ｋ＝１，２，…，ｎ／４，σ２ｅ是阶数为Ｋ时残差的方差，例如：ＡＩＣ（１）＝３，ＡＩＣ（２）＝３．１，…，ＡＩＣ（１０）＝３．２２２１　第５期焉建国等：基于ＡＲ模型的上海地区地面沉降预测分析根据式（１０）在计算得出的值中挑出ＡＩＣ（Ｋ）的最小值，即为模型的阶数。计算过程主要分以下几个步骤：１）计算样本序列的自相关系数ｒｋ＝１ｎ－ｋ∑ｎ－ｋｊ＝１ｕｊｕｊ＋ｋ１ｎ（１１）其中ｒｋ为自相关系数，ｋ为阶数，一般数据少时用１ｎ以下。２）利用递推法计算自回归系数ｂｋ＋１ｋ＋１＝ｒｋ＋１－ｂｋ１ｒｋ－Ｂｋ２ｒｋ－１…－ｂｋｋｒ１１－ｂｋ１ｒ１－ｂｋ２ｒ２－…－ｂｋｋｒｋｂｋ＋１ｊ＝ｂｋｊ－ｂｋ＋１ｋ＋１ｂｋｋ＋１－}ｊ（１２）３）计算残差的方差σ２ｅ＝１－ｂｋ１ｒ１－ｂｋ２ｒ２－…－ｂｋｋｒｋ（１３）４）根据不同阶数Ｋ，计算ＡＩＣＡＩＣ（ｋ）＝ｎｌｎ（σ２ｅ）＋２Ｋ（１４）５）选择ＫＡＩＣ（Ｋ）＝ｍｉｎＡＩＣ（Ｋ）（１５）３　实例计算与分析采用的数据来源于中国地质调查局和上海城市地质调查报告［６，７］。依据ＡＩＣ准则，通过计算可以得到ｋ阶的ＡＩＣ（ｋ）值。如ｋ＝２时，模型为ＡＲ（１），即ｋ阶对应的是ＡＲ（ｋ－１）模型，通过计算，表１为ｋ取不同值时对应的ＡＩＣ（ｋ）值。表１　ＡＩＣ（Ｋ）值Ｔａｂ．１　ＡＩＣ（Ｋ）ｖａｌｕｅｓｋ２３４５６７８９１０１１ＡＩＣ（ｋ）１７２．６６３２１６５．４６６０１５３．２０７１４０．１４４２１４２．７８９８１４１．６８２４１４３．９９４４１４４．５６０１１４６．３０６７１４８．７０３６　　在显著水平为０．０５的条件下，根据ＡＩＣ准则计算得到ＡＲ（ｋ）各个阶次的参数如表２所示，由于后面几种模型比较复杂，而且ＡＩＣ值比较大，实用价值不是很大，故没有计算其系数值。根据ＡＲ理论得到的预测模型，由于阶次越高，计算也越复杂，并不是阶次越高，预测模型的精度越高，因此我们计算了前面几种阶次模型的系数值（表２）。表２　ＡＲ（ｋ）模型系数Ｔａｂ．２　ＰａｒａｍｅｔｅｒｓｏｆＡＲ（ｋ）ｍｏｄｅｌＡＲ（１）ＡＲ（２）ＡＲ（３）ＡＲ（４）ＡＲ（５）ＡＲ（６）参数ｂ０１１．６２３９．９５２９．１７３９．５９７１０．３５８９．３９４ｂ１０．７０４１．０８００．９０２０．９２２０．９３１０．９７７ｂ２－０．４９００．４００．０７７０．６０－０．０１１ｂ３－０．５１０－０．６３７－０．５０２－０．４９１ｂ４０．１３９－０．１５４－０．２５４ｂ５０．２６９０．６０９ｂ６－０．３４２依据表２得到的各ＡＲ（ｋ）模型的系数值以及原始的沉降值（本文所涉及到的数据全部为上海城区年平均地面沉降量），在显著水平为０．０５的条件下，利用１９７１—２００１年原始年平均沉降量的数据依据得到一系列的预测值ＡＲ（１）、ＡＲ（２）、ＡＲ（３）、ＡＲ（４）、ＡＲ（５）、ＡＲ（６），利用上述预测值分别与原始值作差进行比较分析，得到系列残差图如图１所示。从图１可以看出，每一种模型的预测值与真实值之间的残差在０～１ｍｍ之间变化，其中ＡＲ（４）模型得到的残差值最平稳，残差的最大值最小；同时对年沉降量的估计值与真实值之间的残差值进一步统计分析可以得到各个模型的残差平方和、残差最大值（表３）。表３　模型残差统计表（单位：ｍｍ）Ｔａｂ．３　Ｓｔａｔｉｓｔｉｃｓｏｆｔｈｅｒｅｓｉｄｕａｌｓｏｆｔｈｅｍｏｄｅｌ模型残差平方和残差最大值ＡＲ（１）６．８１２２０．８２３４ＡＲ（２）５．６４９３０．５６３５ＡＲ（３）４．３５６４０．４２６５ＡＲ（４）２．５６３５０．２５３５ＡＲ（５）３．２５６１０．３０３６ＡＲ（６）３．５６８６０．３８２３从表３中可看出，ＡＲ（４）模型的残差平方和最小，结合表１的ＡＩＣ值以及图１的残差值可知，对于上海城区地面沉降，在上面的几个模型中，ＡＲ（４）模型的预测结果最可靠；由于选择的试验区同为上海，地面沉降量的自然因素相同，故可以利用已有数据检验ＡＲ（４）模型的准确性。利用１９９８—２００１年的数据来预测２００２—２００６年的沉降量（表４），并与真实值进行比对。从表４可知ＡＲ（４）预测值与真实值相差不大，残差以