一种狭长图斑分块融解方法

李成名, 殷勇, 武鹏达, 吴伟

李成名, 殷勇, 武鹏达, 吴伟. 一种狭长图斑分块融解方法[J]. 武汉大学学报 ( 信息科学版), 2018, 43(12): 2259-2268. DOI: 10.13203/j.whugis20180274
引用本文: 李成名, 殷勇, 武鹏达, 吴伟. 一种狭长图斑分块融解方法[J]. 武汉大学学报 ( 信息科学版), 2018, 43(12): 2259-2268. DOI: 10.13203/j.whugis20180274
LI Chengming, YIN Yong, WU Pengda, WU Wei. A Partitioned Dissolution Method for Long and Narrow Patches[J]. Geomatics and Information Science of Wuhan University, 2018, 43(12): 2259-2268. DOI: 10.13203/j.whugis20180274
Citation: LI Chengming, YIN Yong, WU Pengda, WU Wei. A Partitioned Dissolution Method for Long and Narrow Patches[J]. Geomatics and Information Science of Wuhan University, 2018, 43(12): 2259-2268. DOI: 10.13203/j.whugis20180274

一种狭长图斑分块融解方法

基金项目: 

国家基础测绘项目 A1705

国家自然科学基金 41871375

详细信息
    作者简介:

    李成名, 博士, 研究员, 主要从事全息时空、时空大数据、序码地图、智慧城市研究。cmli@casm.ac.cn

    通讯作者:

    殷勇, 博士。yinyong.@casm.ac.cn

  • 中图分类号: P283

A Partitioned Dissolution Method for Long and Narrow Patches

Funds: 

The National Basic Surveying and Mapping Project A1705

the National Natural Science Foundation of China 41871375

More Information
    Author Bio:

    LI Chengming, PhD, professor, specialige in digital city, smart city and map generalization automatically. E-mail:cmli@casm.ac.cn

    Corresponding author:

    YIN Yong, PhD. E-mail:yinyong@casm.ac.cn

  • 摘要: 融解是图斑综合过程中的一种常见操作,涉及大量计算。面对传统方法受机器计算能力所限,难以对大范围海量图斑进行处理的问题,引入分块策略开展图斑融解,并就块与块之间边界处狭长图斑分裂线拓扑变化问题,提出了一种狭长图斑分块融解方法。首先归纳了狭长图斑分块时在格网边界处出现的4种拓扑变化模式,并针对每种模式提出了相应的分裂线拓扑变化恢复方法;然后利用中国贵州省赤水市地理国情普查实际数据进行了验证。试验结果表明,该方法不仅可以处理海量图斑数据,可以极大地提高融解效率,而且融解结果与整体处理结果保持高度一致,具有良好的可行性。
    Abstract: The dissolution operation is a common operation in patch generalization and it involves a large number of computations. Due to limitations in computational power, it is very difficult to process a large number of patches over a large area using traditional dissolution methods. To overcome this limitation, this paper proposes a partitioned dissolution method for long and narrow patches (LN patches) by introducing the block strategy and focusing on the topological changes that occur in the partition lines of LN patches at partition-cell boundaries. Firstly, the topological changes that appear around partition-cell boundaries during the dissolution of LN patches are summarized into four patterns. For each of these patterns, a corresponding method is formulated for reconciling the topological changes of the partition line. Our approach is then validated by using the national geographical conditions data of Chishui City, Guizhou Province. It is experimentally demonstrated that our method is capable of processing large-scale patch data and greatly improves the efficiency of patch dissolution opera-tions. The results of our partitioned dissolution method are also found to be strongly consistent with the unpartitioned method of dissolution. Our method is therefore viable for practical applications.
  • 在直线、平面拟合、空间直角坐标变换、自回归模型求解中,系数矩阵和右端观测量由随机和非随机元素组成,且同一随机元素会在不同的位置出现[1-4]。常用的变量含误差(errors-in-variables,EIV)模型及相应的整体最小二乘(total least squares,TLS)方法假设系数阵中所有元素含误差。因此,顾及增广误差矩阵元素的随机特性及其相互关系的结构EIV(structured EIV,SEIV)模型和结构加权整体最小二乘(structured weighted TLS,SWTLS)方法应运而生[5-10]。文献[5]采用极大似然估计方法从含有误差的时间序列数据中确定动态线性系统的参数,被认为是SWTLS问题研究的开端。文献[6]定义了结构整体最小二乘(structured TLS,STLS)这一术语,并将其转换为非线性广义奇异值分解问题求解。文献[7]将结构增广误差矩阵用其中的独立随机元素表示,并命名为约束TLS方法。对于某些列为固定元素而另外一些列为随机元素的混合LS-TLS问题,文献[8]引入正交三角分解计算模型参数。对于误差为非等权相关的情形,文献[9]发展了WLS-WTLS的迭代算法。文献[10]将SEIV模型用部分EIV(partial EIV,PEIV)模型表达,在加权最小二乘准则下推导了参数估计的Gauss-Newton型迭代算法,并进一步给出了参数的一阶近似方差、非线性信赖域区间及偏差。文献[11]给出了PEIV模型的两种迭代算法,其计算效率与系数矩阵中随机元素的数量有关。文献[12]提出了顾及系数阵元素和右端向量相关的广义PEIV模型。文献[13]将系数误差矩阵表示成已知矩阵和独立误差向量的乘积,推导了附有线性和二次约束的STLS问题的迭代算法。文献[14]和文献[15]采用变量投影法将增广系数矩阵表示成仿射结构矩阵与独立随机变量的乘积,然后将STLS问题转换为非线性等式约束优化问题求解。文献[16]研究了含多个右端观测向量的STLS问题。文献[17]提出了结构整体最小范数(structured total least norm,STLN)方法,这一方法可以最小化误差向量的范数,是STLS问题的重要拓展[17-18]。文献[19-20]研究了STLN问题的快速算法。文献[21]将STLN拓展到系数阵和观测向量具有共同元素的情形。针对自回归(auto-regression, AR)模型这一特定STLS问题,文献[3-4]提出了两种新解法。

    在STLS算法推导中,不同文献使用了不同的平差准则,其差异在于是否考虑独立误差的重复次数。如二维坐标转换中,源系统的纵横坐标值在系数矩阵中出现了2次,而目标系统的纵横坐标值在右端向量中只出现了1次。文献[1,7,10,13-15]给出的平差准则未考虑重复次数,而文献[3-4,17,21]均考虑了重复次数。不同的平差准则必定会得到不同的平差结果,其参数估值在统计意义上孰优孰劣,目前尚未给出明确结果。本文从函数模型和数值模拟两方面入手,证明了不考虑独立误差元素重复次数的平差准则能够得到统计意义下更优的解。

    EIV模型的函数表达式为[1]

    y+ey=A+EAx (1a)

    式中,yey分别表示n维观测向量及其误差;AEA分别表示n×m维系数矩阵及其误差矩阵;xm维参数向量。误差向量的随机模型为:

    e=eAey~00,σ02Qe=00,σ02QAAQAyQyAQyy  (1b)

    式中,eA=vecEA,vec表示矩阵向量化算子,即将n×m维矩阵的每一列从左至右叠加成一个nm维列向量;σ02是单位权方差;QAAQyy分别是eAey的对称正定协因数矩阵;QyA=QAyT表示两者的协因数矩阵。当向量e的协因数矩阵Qe可逆时,采用WTLS准则求参数的最优估值:

    Φe=eTPee (2)

    式中,Pe=Qe-1e的权矩阵。当系数矩阵具有某种结构时,Qe为秩亏矩阵无凯利逆。许多学者将EIV函数模型进行改化,一般是从结构误差矩阵或增广误差矩阵中提取独立误差向量,进一步构造目标函数求解。常用的几种EIV模型修正方法如下:

    1) PEIV模型及其平差准则。PEIV模型选取系数阵中独立随机量的真值a¯作为待求量,其函数模型为[10]

    y=xTInh+Ma¯+ey (3a)
    a=a¯+ea (3b)

    式中,a是系数矩阵中随机元素所构成的t维列向量;a¯ea分别是相应的真值和误差向量;In表示n维单位矩阵;h是已知的nm维常数向量,其元素包含系数阵中的非随机元素;Mnm×t维常数矩阵,其形式由系数矩阵中非随机元素的个数及元素间的相关性确定。若ay相互独立[10],即cov(ea,ey)=0,且两者的方差矩阵为D(ea)=σ02ω-1,D(ey)=σ02W-1,其中ωW分别为eaey的权矩阵。相应的平差准则为:

    Φa¯,x=eaTωea+eyTWey (4)

    可见PEIV的平差准则中没有考虑独立随机误差的重复次数。如果将eaey组合成独立随机误差向量γ=eaTeyTT,权矩阵Pγ=diag(ω,W)=(diag(ω-1,W-1))-1,则平差准则式(4)可以写为:

    Φ(a¯,x)=γTPγγ (5)

    因此PEIV模型的SWTLS解是令独立随机误差的加权平方和最小的参数估值。文献[12]中的GPEIV模型、文献[13]中约束结构SWTLS方法和文献[14-15]中的变量投影法,其本质上都是采用式(5)所示的平差准则。

    2) STLN模型及其平差准则。STLN方法是定义一个与参数x相关的矩阵X,提取系数矩阵中的独立随机误差ea,使下式成立[17]

    Xea=EAx (6)

    联立式(1a)和式(6),将观测值残差r表达成系数阵独立误差ea和参数x的函数:

    r(ea,x)=-ey=y-Ax-Xea (7)

    将式(7)在eax的近似值处线性化,舍去二次项后得到线性模型,然后采用下列平差准则[17]

    Φea,x=0.5eyTey+0.5eaTDa2 (8)

    式中,Da=diagd1,d2,,dt,dii=1,2,,t)表示ea中的第i个元素在误差矩阵EA中重复出现的次数。尽管STLN 方法没有考虑误差的权值,但可以很方便地将观测误差的权纳入平差准则。仍令γ=eaTeyTT,定义Ds=diagDa,In,由于式(8)中常数项对求极值无影响,则平差准则式(8)等价于:

    Φx=γTDs2γ (9)

    可见,STLN方法考虑了独立误差ea的重复次数,且以重复次数的平方将独立误差纳入平差模型。

    3) 虚拟误差模型及其平差准则。文献[3]以AR模型参数估计为背景,将式(1a)所示的EIV模型在观测值的近似值处线性化。设系数阵真值A˜的近似值为A0,改正数为A,参数x的近似值为x0,改正数为x,则线性化方程为:

    v=A0x0+A0x+Ax0-y (10)

    式中,v表示y对应的改正数。然后通过矩阵等价变换得到Ax0=A10A20xBv,A10A20表示由x0按照一定规则构造的近似矩阵;xB表示设计矩阵中独立观测值对应的改正数向量。将式(10)化为:

    v=In-A20-1A0x+In-A20-1A10xB-In-A20-1y-A0x0 (11)

    组成虚拟观测值误差方程为:

    vB=xB (12)

    联立式(11)和式(12)得到:

    vg=Bgxg-lg (13)

    式中,vg为所有随机观测值改正数向量;Bg为等效设计矩阵;xg为所有待估参数和虚拟参数的改正数;lg表示y-A0x0的线性变换[3]。采用如下平差准则:

    Φx=vgTPgvg (14)

    式中,vg实际上是EA|ey中的独立误差向量γPg=DgP由两部分组成,其中Dg=diagd1,d2,,dt+n,dii=1,2,,t+n,且t+n为独立误差的个数)表示vg中的第i个元素在(EA|ey)中的重复次数,P为观测值向量的权矩阵。可见,虚拟误差模型解AR模型这一典型SWTLS问题考虑了重复次数。与STLN方法的差异在于,它采用的是重复次数本身,而不是STLN中重复次数的平方。

    上述3种典型的STLS方法采用了不同的平差准则,不同的平差准则一定会得到不同的平差结果。为了比较不同准则下的平差结果,下面将上述3种典型的平差准则下的解纳入到一个统一的平差模型中。

    EIV模型(式(1a))又可以表示为[1]

    y-Ax-Be=0 (15)

    式中,B=xTIn-Inn×nm+n矩阵,其中为Kronecker积符号,定义为GH=gijH,其中G=gij,且H为任意矩阵。设γEA|ey中的t个独立随机误差向量,则eAey可以分别表示为:

    eA=H1γ (16a)
    ey=H2γ (16b)

    式中,H1H2分别为nm×tn×t矩阵,其元素均为常数。联立式(16a)和式(16b):

    e=eAey=H1H2γ=Hγ (17)

    将式(17)代入式(15),并令G=BH,则EIV模型可以写为:

    y-Ax-Gγ=0 (18)

    独立随机误差的随机模型为:

    Eγ=0, Dγ=σ02Qγ=σ02Pγ-1 (19)

    式中,Qγ表示γ的协因数矩阵;Pγ表示其权矩阵。令γ=γ1 γ2  γtT,且γi出现的次数为dii=1,2,,t)。定义D=diagk1,k2,,kt,组成综合权矩阵:

    Psγ=DPγ (20)

    相应的平差准则为:

    ψγ=γTPsγγ=min (21)

    易知,当ki=1时,式(21)等价于以PEIV模型为代表的准则式(5);当ki=di2时,式(21)等价于以STLN模型为代表的准则式(9);当ki=di时,式(21)等价于以虚拟误差模型为代表的准则式(14)。以结构化函数模型(18)为基础,基于平差准则式(21),构造如下Lagrange目标函数:

    φx,γ,λ=γTPsγγ-2λTy-Ax-Gγ (22)

    式中,λn维Lagrange乘子向量。令目标函数式(22)对各待定量xγλ的偏导数为0,可得:

    12φx,γ,λxx^,γ^,λ^=A+E^ATλ^=0 (23a)
    12φx,γ,λγx^,γ^,λ^=Psγγ^+GTλ^=0 (23b)
    12φx,γ,λλx^,γ^,λ^=y-Ax^-Gγ^=0 (23c)

    由式(23b)得到独立误差向量γ的估值为:

    γ^=-Psγ-1GTλ^=-QsγGTλ^ (24)

    式中,Qsγ=Psγ-1,表示综合权矩阵的逆;x^γ^λ^分别为xγλ的估值。将式(24)代入式(23c)并移项可得:

    λ^=-GQsγGT-1y-Ax^ (25)

    将式(25)代入式(23a)可以得到法方程:

    (A+E^A)T(GQsγGT)-1Ax^=(A+E^A)T(GQsγGT)-1y (26)

    则参数的估值为:

    x^=[(A+E^A)T(GQsγGT)-1A]-1(A+E^A)T(GQsγGT)-1y (27)

    由于式(26)法方程矩阵不对称,若在式(26)两端加上A+E^ATGQsγGT-1E^Ax^,那么x^可由下式估计:

    x^=[(A+E^A)T(GQsγGT)-1(A+E^A)]-1(A+E^A)T(GQsγGT)-1(y+E^Ax^) (28)

    根据上述推导过程,可以得出通用模型计算SWTLS问题的迭代步骤为:

    1) 给定AyQγ,根据γ的重复次数和不同平差准则的计入方式确定对角阵D,根据(EA|ey)的结构确定H1H2。采用式(20)计算综合权阵Psγ和对应的协因数矩阵Qsγ=Psγ-1。计算Qyy=H2QγH2T和参数的初值x^0=x^LS=ATQyy-1A-1ATQyy-1y

    2) 根据初值x^0计算矩阵BG。采用式(25)和式(24)分别计算λ^γ^,然后由式(16a)计算e^A的值,从而有E^A=ivece^A,ivecvec的逆运算,表示将nm向量恢复成n×m维矩阵;

    3) 将GQsγE^Ax^0分别代入式(27)或式(28),计算参数估值x^

    4) 如果最后两次的估值x^足够接近于给定的阈值,终止迭代。否则,转向步骤2)。

    尽管由式(18)和式(19)组成的SWTLS一般模型能通过选取不同的D矩阵获得不同准则下的迭代最优解,但无法获得参数的方差或均方误差(mean square error,MSE)等精度评定指标,只能在模拟实验中检验参数精度。由文献[22-23]可知,WTLS并不是一种新的平差方法,仅仅是LS框架下的另外一种平差模型,且经典平差理论中的Gauss-Helmert(GH)模型(附有参数的条件平差模型)是EIV模型的一个特例。因此WTLS问题的解可由非线性GH模型导出。若给定x的近似值x0γ的近似值γ0,式(18)可线性化为:

    fx,γ=y-Ax-Gγ=y-Ax0-G0γ0+fx,γxx0,γ0x+fx,γγx0,γ0γ (29)

    式中,Δx=x-x0Δγ=γ-γ0G0EA0为近似值x0γ0处的GEA的值。且有:

    fx,γxx0,γ0=-(A+EA0)=-(A+ivec(H1γ0)) (30)
    fx,γγx0,γ0=-G0=-x0TIn-InH (31)

    将式(30)和式(31)代入式(29)并作适当变换可以得到:

    A+EA0x+G0γ-y+EA0x0=0 (32)

    可见SEIV模型的线性化形式(式(32))是一个标准GH模型,该模型不再体现γ的重复次数,实际上重复次数已经由投影矩阵H表达了。重复次数体现了EIV模型的结构性,而投影矩阵H正是描述上述结构性的量。根据经典LS理论,其平差准则应为:

    Φγ,x=γTPγγ (33)

    因此,从SWTLS问题的线性化形式证明了其平差准则应满足式(33)。由线性化模型式(32)得到的最小二乘解是在给定x0γ0情形下的一个近似解,为了得到更严密的解,将式(32)进一步表达为:

    l=A+EA0x+G0γ+G0γ0 (34)

    式中,l=y-Ax0。为了与§2中符号一致,记Psγ=ItPγ=PγItt维单位矩阵)。根据式(34)及平差准则式(33),构造如下Lagrange乘子函数:

    min:φx,γ,λ=γ0+γTPsγγ0+γ+2λTl-A+EA0x-G0γ-G0γ0 (35)

    分别求φx,γxγλ的偏导数,并令其值为0,可得:

    12φx,γx=-A+EA0Tλ=0 (36)
    12φx,γγ=Psγγ0+γ-G0Tλ=0 (37)
    12φx,γλ=l-A+EA0x-G0γ-G0γ0=0 (38)

    由式(37)可得:

    γ=QsγG0Tλ-γ0  (39)

    将式(39)代入式(38)得:

    l-A+EA0x-G0QsγG0Tλ=0 (40)

    由式(40)可得Lagrange乘子向量为:

    λ=G0QsγG0T-1l-A+EA0x (41)

    式(40)两边同乘以A+EA0T并顾及式(36),则有:

    x=[(A+EA0)T(G0QsγG0T)-1(A+EA0)]-1(A+EA0)T(G0QsγG0T)-1(y-Ax0) (42)
    x=x0+Δx=[(A+EA0)T(G0QsγG0T)-1(A+EA0)]-1(A+EA0)T(G0QsγG0T)-1(y+EA0x0) (43)

    由GH模型推导得到的参数估计(式(43))与§2中通用模型得到的解(式(28))具有完全相同的形式。GH模型的算法流程如下:

    1) 给定AyQγ,根据(EA|ey)的结构确定H1H2。计算Qyy=H2QγH2T和参数的初值x^0=x^LS=ATQyy-1A-1ATQyy-1y,使用上标(i(i=0,1,2,,N,N为迭代次)对不同迭代次数下的各参数进行区分,则γEA的初值分别为γ0=0EA0=ivecH1γ0=0

    2) 将x^i代入式(31)计算BiGi,由式(41)计算λi,由式(39)计算γi并更新γi=γi-1+γi,由式(30)计算EAi

    3) 由式(42)计算Δxi-1

    4) 如果Δx小于给定的阈值,则终止迭代;否则,更新xi=Δxi-1+xi-1,转向步骤2)。

    注意,若不限定式(35)中Psγ=Pγ,那么采用线性化的GH模型,在准则式(21)下,即可选取不同类型D矩阵的条件下,可以得到和§2中通用模型完全一致的解。采用GH模型进行算法推导有3个方面的意义:(1)证明了SEIV模型的平差准则中不应该再考虑重复矩阵。(2)证明了由通用模型(式(18))和线性化GH模型作为条件,在相同的平差准则下得到的参数估值是一致的。(3)将SEIV模型线性化后,便于采用协方差传播定律计算参数估值的近似精度。

    根据§2、§3的算法得到SWTLS解以后,可计算独立误差向量的估值γ^,单位权方差计算式为:

    σ^02=γ^TPsγγ^n-m (44)

    SWTLS解的近似协因数矩阵可以根据线性化GH模型的解(式(42))求解。虚拟观测值l可以写成:

    l=y-Ax0=-x0TInInvecAy=-B0Hη=-G0η (45)

    式中,η(A|y)中的t个独立随机观测向量。根据协因数传播律,向量l的协因数矩阵为:

    Qll=G0QsγG0T (46)

    实际计算中近似值x0可以用SWTLS估值x^代替,相应的残差矩阵近似值EA0用对应的估值E˜A代替。如果忽略式(42)中A˜=A+EA0这一项中A的随机性,那么可由协因数传播律得到SWTLS解的近似协因数矩阵为:

    Qx^x^=QΔxΔx=A˜TQll-1A˜-1A˜TQll-1QllQll-1A˜A˜TQll-1A˜-1=A˜TQll-1A˜-1 (47)

    从而得到SWTLS解的近似方差为:

    Dx^x^=σ^02A˜TQll-1A˜-1 (48)

    首先采用二维仿射变换实例验证算法的正确性和可行性。设第i个公共点Pi在源坐标系和目标坐标系中的坐标分别为xi,yiXi,Yi,仿射变换模型为Xi=a0+a1xi+a2yi,Yi=b0+b1xi+b2yi,其中a0b0为两个坐标系的平移量,a1a2b1b2分别表示两个坐标系之间的旋转和尺度变换参数。假设有3个或以上的公共点,仿射变换模型可以用EIV模型(式(1))表示,且有:

    A=1x1y10000001x1y11x2y20000001x2y21xnyn0000001xnyn,x=a0a1a2b0b1b2 (49)

    给定12个公共点在源坐标系中的坐标分别为(-1,1)、(1.2,-3.0)、(-2.6,3.0)、(3.0,1.5)、(-4.8,-1.0)、(5.2,0.2)、(6.0,5.5)、(-7.2,2.2)、(7.8,-2.0)、(8.5,2.2)、(-9.5,-5.0)、(10.0,-0.8),参数的真值为x˜=[10 4 -2 -10 1 3]T,根据仿射变换模型计算出12个点在目标系中的真实坐标分别为(4,-8)、(20.8,-17.8)、(-6.4,-3.6)、(19.0,-2.5)、(-7.2,-17.8)、(30.4,-4.2)、(23.0,12.5)、(-23.2,-10.6)、(45.2,-8.2)、(39.6,5.1)、(-18.0,-34.5)、(51.6,-2.4)。对所有坐标值添加方差为D48=σ02I48=I48I4848×48的单位矩阵)的随机误差。首先,不考虑系数矩阵中源坐标重复2次而目标坐标重复1次的差异,取重复矩阵D=D1=I48时,分别采用本文方法、PEIV模型[10]和变量投影法[14]进行求解,得到的参数估计结果如表1所示。

    表  1  不同平差方法的SWTLS解(D=D1
    Table  1.  SWTLS Solutions with Different Adjustment Methods (D=D1
    平差方法â0â1â2b̂0b̂1b̂2
    PEIV模型[10]10.412 0554.065 602-2.136 739-9.641 7471.116 1702.912 190
    变量投影法[14]10.412 0554.065 602-2.136 739-9.641 7471.116 1702.912 190
    本文方法10.412 0554.065 602-2.136 739-9.641 7471.116 1702.912 190
    下载: 导出CSV 
    | 显示表格

    当取D=D2=diag(4,,4,1,,1)(4和1的个数都为24)时,采用本文方法和STLN方法[17]进行求解,得到的参数估计结果如表2所示。

    表  2  不同平差方法的SWTLS解(D=D2
    Table  2.  SWTLS Solutions with Different Adjustment Methods (D=D2
    平差方法â0â1â2b̂0b̂1b̂2
    STLN方法[17]10.418 6814.056 913-2.101 082-9.639 4401.118 3362.855 364
    本文方法10.418 6814.056 913-2.101 082-9.639 4401.118 3362.855 364
    下载: 导出CSV 
    | 显示表格

    表1可以看出,由于PEIV模型[10]、变量投影法[14]均未考虑源系统中坐标的重复次数,两者的平差准则等价于本文提出的通用模型中重复矩阵取单位矩阵时的准则(式(21)),尽管三者对结构误差的处理采用了不同的函数模型,但在相同的准则下都得到了完全一致的平差结果,说明本文提出的通用平差算法是可行有效的。由表2可以看出,STLN方法[17]顾及了系数阵中随机观测值的重复次数,并且是将观测重复数的平方纳入平差模型,和本文通用模型在同样的平差准则下也得到了一致的平差结果,进一步证明了本文算法的可行性和有效性。

    为了验证不同平差准则下SWTLS解的统计性质,模拟系数阵元素和右端项含有公共元素的SEIV模型,系数阵和右端项的真值及误差的结构分别为:

    A˜|y˜=3.62-1.030.754.050.75-4.59-4.161.252.752.50-38.881.55-18.89-4.505.455.288.25-10.854.55-16.701.553.203.177.956.001.252.500.20-1.72-3.210.850.210.550.48-3.008.420.8517.501.55-6.182.19-8.014.55-2.210.1201.320.284.85-8.726.75-5.362.502.906.751.421.985.50-1.8015.003.077.20-7.932.5015.004.0617.500.456.2524.353.202.490.65-10.6532.6023.370.756.751.25-1.802.500.8515.001.5517.509.3620.2052.758.254.5524.353.2012.806.0032.607.9514.785.8710.9027.60 (50)
    EA|ey=0γ4γ30γ60γ20γ6000γ20γ30γ70000γ40000γ5γ60γ8000γ9γ700γ1000γ9γ600γ8000γ140γ100γ150γ150000γ16γ90γ17γ1700000γ2100γ190γ20γ1γ2γ3γ4γ5γ6γ7γ8γ9γ10γ11γ12γ13γ14γ15γ16γ17γ18γ19γ20γ21γ22γ23γ24γ25 (51)

    待估参数的真值为x˜=[152]T,独立误差向量γ=[γ1 γ2  γ25]T是零均值独立同分布的随机变量,其方差为Dγ=σ02I25,其中I25表示25×25阶单位矩阵,γii=1,2,,25)的重复次数di可以EA|ey的结构获得。分别采用3种平差准则计算参数的估值,准则式(21)中对应的D矩阵分别取D1=I25,D2=diag(d1,d2,,d25)=diag(1,3,1,3,,1,3,1),D3=diag(d12,d22,,d252)=diag(1,9,1,9,,1,9,1)。根据式(44)和式(48)计算单位权方差和参数估值的方差。由于参数真值已知,可以计算均方误差MSE(x^)=x^-x˜Tx^-x˜,以此来检核3种平差准则下估值的精度,并与近似方差对照。分别取误差水平σ02=0.25σ02=1进行计算。为了体现结果的统计性,将不同误差水平的实验各进行10 000次,将每次运算的结果取平均值,结果分别见表3表4

    表  3  不同平差准则下的SWTLS结果(σ02=0.25
    Table  3.  SWTLS Results Under Different Adjustment Principles (σ02=0.25
    统计项统计量D=D1D=D2D=D3
    估计参数值x̂10.999 9981.000 2810.999 975
    x̂25.002 3895.007 3595.009 953
    x̂31.999 5961.999 3201.998 190
    单位权方差σ̂020.239 7370.537 3491.463 591
    参数方差估值σ̂x129.43×10-61.72×10-55.76×10-5
    σ̂x225.23×10-49.14×10-43.19×10-3
    σ̂x328.28×10-61.60×10-55.05×10-5
    σ̂Σ25.40×10-49.47×10-43.30×10-3
    MSEm̂x129.76×10-61.10×10-51.63×10-5
    m̂x225.44×10-46.21×10-49.40×10-4
    m̂x328.50×10-69.39×10-61.42×10-5
    m̂Σ25.63×10-46.41×10-49.71×10-4
    下载: 导出CSV 
    | 显示表格
    表  4  不同平差准则下的SWTLS结果(σ02=1
    Table  4.  SWTLS Results Under Different Adjustment Principles (σ02=1
    统计项统计量D=D1D=D2D=D3
    估计参数值x̂11.000 2201.003 7701.007 659
    x̂24.997 8385.022 9635.049 989
    x̂32.000 7252.004 0002.007 714
    单位权方差σ̂020.960 3412.146 4735.839 121
    参数方差估值σ̂x123.78×10-56.87×10-52.30×10-4
    σ̂x222.09×10-33.65×10-31.27×10-2
    σ̂x323.32×10-56.41×10-52.02×10-4
    σ̂Σ22.17×10-33.78×10-31.32×10-2
    MSEm̂x124.02×10-54.50×10-56.67×10-5
    m̂x222.21×10-32.53×10-33.85×10-3
    m̂x323.53×10-53.88×10-55.82×10-5
    m̂Σ22.29×10-32.62×10-33.97×10-3
    下载: 导出CSV 
    | 显示表格

    表3可以看出,不考虑重复次数的平差准则(采用D1)获得的参数估值更接近于真值。参数各分量的MSE最小,从数值上验证了不考虑系数阵误差重复系数的平差准则是最优的。单位权方差的估值略小于真值,是由于式(43)没有考虑非线性模型线性化引入的偏差项。同理,各参数分量的方差均小于对应的MSE,原因是SEIV模型的非线性特性会引入估计偏差,方差加上偏差的平方和等于MSE。但是,对比方差和MSE的差异可知,参数3个分量的方差估值和MSE之间的差异((σ^xi2-m^xi2)/m^xi2)分别为3.38%、3.86%和2.59%,说明当误差较小时,参数方差估值和MSE的差异很小,可以作为精度评定的指标。因此,本文不再采用非线性最小二乘偏差修正或Monte Carlo模拟的思想求解偏差值[10,24]

    从上述理论分析可知,将SEIV模型线性化后,其平差准则不应顾及误差重复次数。式(44)和式(48)只适用于D为单位阵的情形。从后面两列可以看出,当平差准则中计入误差重复的次数时,单位权方差估值与真值以及参数分量的方差估值与MSE均有较大的差异。当D=D3时,单位权方差估值约为真值的6倍,方差估值约为对应MSE的3.5倍,这是由于单位权方差公式中重复计算残差平方和所致,由此进一步证明不能采用这两种平差准则用线性近似方法求单位权中误差及方差。

    表4可知,增大观测误差的方差,能够得到与表3一致的结论,进一步验证了最优平差准则应该选择重复矩阵D为单位阵。由式(48)给出的方差估值与MSE在3个分量上的偏差分别为5.97%、5.43%、5.95%,进一步证明了本文的近似精度评定方法在误差较小的情况下是可行有效的。

    结构EIV模型系数矩阵中的随机元素重复出现的次数是否应计入平差准则以及如何计入平差准则,目前尚未形成定论。本文从模型分析和数值验证两方面入手,证明重复次数不应计入平差准则。主要贡献如下:

    1)总结了已有的3种处理SEIV模型的平差准则,指出不同的平差准则会得到不同的平差结果。提出了一种通用的SWTLS平差模型,通过选取不同的综合权矩阵Psγ得到的目标函数等价于上述3种不同准则下的目标函数。采用Lagrange乘子法推导了通用模型的解并给出了计算步骤。

    2)迭代方法无法给出参数的统计性质,本文将通用模型线性化得到GH模型,从理论上分析了误差重复次数不应计入平差准则的原因。推导了GH模型的算法,证明其与通用模型得到的结果是一致的。根据误差传播律得到了参数的近似方差估值。

    3)通过实例验证了本文提出的算法与已有方法结果一致,证明了本文方法可行有效。通过模拟计算证明了不考虑误差重复次数的平差准则得到的解在MSE意义下最优,且参数的近似方差是MSE的良好近似,可以作为精度评定的指标。

  • 图  1   狭长图斑融解方法

    Figure  1.   Dissolution Method for Long and Narrow Patches

    图  2   块与块边界处狭长图斑分裂线拓扑变化模式

    Figure  2.   Topological Changes Around Partition-cell Boundaries

    图  3   狭长图斑分裂线拓扑变化恢复

    Figure  3.   Reconciling the Topological Changes of the Partition Line

    图  4   各个狭长图斑向邻近不同地类图斑和邻近图斑转换的面积对比图

    Figure  4.   Comparison of the Transformation Areas for Adjacent Different Classes of Patches and Adjacent Each Patch

    图  5   转换面积对比分析

    Figure  5.   Comparison of the Transformation Areas

    图  6   2号狭长图斑转换面积差异值折线图

    Figure  6.   Line Chart of Transformation Areas Difference for Adjacent Each Patch of No.2 Patch

    图  7   局部区域融解结果对比

    Figure  7.   Comparison of Local Regional Dissolution Results

    图  8   分块格网划分实验结果

    Figure  8.   Experimental Results of Partitioned Grid Division

    图  9   典型区域的融解结果

    Figure  9.   Typical Area Dissolution Results

    表  1   向邻近不同地类图斑转换的面积对比

    Table  1   Comparison of Transformation Areas for Adjacent Different Classes of Patches

    狭长图斑序号 转换的面积及其面积差/m2
    y1/x1/
    x1y1
    y2/x2/
    x2y2
    y3/x3/
    x3y3
    y4/x4/
    x4y4
    y5/x5/
    x5y5
    y6/x6/
    x6y6
    y7/x7/
    x7y7
    y8/x8/
    x8y8
    y9/x9/
    x9y9
    y10/x10/
    x10y10
    1 2 028.26/
    2 028.26/
    0
    0/0/0 62 198.50/
    62 198.50/
    0
    125.54/
    125.54/
    0
    79.71/
    79.71/
    0
    0/0/0 1 057.30/
    1 057.30/
    0
    0/0/0 88.57/
    88.57/
    0
    205.10/
    205.10/
    0
    2 111 727.05/
    111 706.63/
    -20.43
    21 910.79/
    21 911.94/
    1.14
    679 283.60/
    679 158.49/
    -125.12
    227 419.37/
    227 425.03/
    5.67
    16 539.35/
    16 539.35/
    0
    6.29/
    23.60/
    17.31
    18 003.06/
    18 110.01/
    106.95
    16 703.98/
    16 704.21/
    0.23
    12 187.03/
    12 187.03/
    0
    30.52/
    30.52/
    0
    3 0/0/0 0/0/0 2 897.04/
    2 897.04/
    0
    0/0/0 0/0/0 0/0/0 0/0/0 0/0/0 0/0/0 0/0/0
    34 20 910.38/
    20 910.47/
    0.09
    0/0/0 127 935.87/
    127 935.77/
    -0.09
    6 977.15/
    6 977.15/
    0
    1 171.53/
    1 171.54/
    0.01
    0/0/0 677.88/
    677.88/
    0
    0/0/0 396.97/
    396.97/
    0
    0/0/0
    35 43 206.76/
    43 210.07/
    3.31
    3 468.64/
    3 468.64/
    0
    243 393.11/
    243 421.22/
    28.11
    3 371.70/
    3 371.70/
    0
    21 383.76/
    21 352.33/
    -31.43
    27.12/
    27.12/
    0
    3 732.69/
    3 732.69/
    0
    375.92/
    375.93/
    0
    964.26/
    964.26/
    0
    0/0/0
    36 1 093.49/
    1 093.49/
    0
    0/0/0 4 859.06/
    4 859.06/0
    103.70/
    103.70/
    0
    0/0/0 7.08/
    7.08/0
    0/0/0 0/0/0 0/0/0 0/0/0
    37 16 937.85/
    16 937.85/
    0
    0/0/0 31 604.67/
    31 604.67/
    0
    64 869.86/
    64 869.86/
    0
    972.07/
    972.07/0
    166.49/
    166.49/
    0
    3 768.17/
    3 768.17/
    0
    13 166.14/
    13 166.14/
    0
    0/0/0 0/0/0
    下载: 导出CSV

    表  2   向邻近各个图斑转换的面积对比

    Table  2   Comparison of the Transformation Areas for Adjacent Each Patch

    狭长图斑编号 邻近图斑数量n s/m2 t/m2
    1 43 15.68 -18.36
    2 1 650 14.52 103.69
    3 4 0 0
    4 14 0 0
    5 7 0 0
    6 2 0 0
    7 4 0 0
    8 7 0 0
    9 1 0 0
    10 56 0.11 1.78
    11 16 0.97 -2.76
    12 2 0.01 0.11
    13 2 0 0
    14 8 0.16 0.80
    15 8 0 0
    16 12 0 0
    17 5 0 0
    18 61 0 0
    19 7 0 0
    20 136 0 0
    21 4 0 0
    22 5 0 0
    23 1 0 0
    24 59 0 0
    25 15 0 0
    26 277 0.35 7.30
    27 46 0 0
    28 51 0 0
    29 36 0 0
    30 6 0 0
    31 196 0.14 3.75
    32 5 0 0
    33 3 0 0
    34 298 0.79 8.02
    35 740 10.70 -51.22
    36 16 0 0
    37 106 0.33 4.46
    下载: 导出CSV

    表  3   11号狭长图斑向邻近图斑转换的面积对比

    Table  3   Comparison of the Transformation Areas for Adjacent Each Patch of No.11 Patch

    邻近图斑编号 xi/m2 yi/m2 (xiyi)/m2
    1 48.22 48.22 0
    2 103.65 103.65 0
    3 47.51 47.51 0
    4 32.62 32.62 0
    5 184.83 184.83 0
    6 141.62 141.62 0
    7 108.01 108.01 0
    8 76.36 76.36 0
    9 50.13 50.13 0
    10 823.80 823.80 0
    11 784.99 787.75 -2.76
    12 258.27 258.27 0
    13 634.70 634.70 0
    14 332.37 329.61 2.76
    15 153.67 154.07 -0.40
    16 59.27 58.87 0.40
    下载: 导出CSV

    表  4   2号狭长图斑向邻近图斑转换的面积对比

    Table  4   Comparison of the Transformation Areas for Adjacent Each Patch of No.2 Patch

    邻近图斑编号 xi/m2 yi/m2 (xiyi)/m2
    1 508 313.42 315.41 -1.99
    1 513 196.05 198.00 -1.96
    1 514 842.10 830.60 11.50
    1 515 864.15 869.26 -5.11
    1 516 456.82 457.93 -1.11
    1 517 1 176.72 1 170.95 5.77
    1 518 698.78 725.01 -26.23
    1 521 293.25 294.21 -0.97
    1 523 49.92 40.52 9.39
    1 525 7 957.00 7 949.82 7.18
    1 526 7 277.28 7 287.94 -10.65
    1 530 278.70 284.47 -5.77
    1 534 1 821.34 1 833.81 -12.47
    1 535 831.02 727.33 103.69
    1 536 2 678.10 2 674.33 3.77
    1 539 8 748.64 8 835.53 -86.89
    1 540 3 069.17 3 067.09 2.07
    1 542 295.75 296.07 -0.33
    1 545 556.18 553.35 2.83
    1 547 808.82 808.62 0.20
    1 550 16 954.12 16 954.96 -0.84
    1 551 1 725.15 1 725.41 -0.26
    1 552 87.23 83.64 3.59
    1 553 158.90 159.10 -0.20
    1 555 531.19 529.39 1.80
    1 556 596.62 594.32 2.30
    1 562 149.26 158.50 -9.24
    1 563 3 441.52 3 439.16 2.36
    1 564 639.18 638.85 0.33
    1 567 1 885.65 1 889.17 -3.53
    1 572 4 473.22 4 477.37 -4.15
    1 573 6 119.59 6 108.94 10.65
    1 574 1 231.00 1 230.70 0.30
    1 575 9 051.01 9 051.32 -0.30
    1 577 4 373.62 4 385.39 -11.76
    1 578 15 622.03 15 631.76 -9.74
    1 579 19 679.41 19 680.78 -1.37
    1 582 401.74 405.17 -3.43
    1 585 244.87 245.58 -0.71
    1 587 2 388.91 2 387.67 1.24
    1 589 191.30 191.05 0.26
    1 590 173.16 172.20 0.96
    1 592 2 009.89 2 015.18 -5.29
    1 595 616.13 618.89 -2.77
    1 600 4 560.03 4 559.41 0.61
    1 601 12 528.97 12 527.43 1.54
    1 607 92.67 95.42 -2.74
    1 608 713.40 712.25 1.14
    1 611 3 118.38 3 115.64 2.74
    1 612 986.80 963.11 23.69
    1 614 3 139.87 3 140.07 -0.20
    1 616 94.91 93.51 1.40
    1 617 378.69 373.93 4.76
    1 618 624.93 629.69 -4.76
    1 621 3 834.37 3 830.62 3.76
    1 626 9 327.31 9 321.17 6.15
    1 628 12 940.37 12 961.99 -21.62
    1 630 31.96 32.13 -0.17
    1 631 214.39 214.21 0.17
    1 634 15 700.92 15 707.63 -6.70
    1 636 433.29 430.04 3.26
    1 637 789.65 784.88 4.78
    1 638 771.48 768.72 2.77
    1 639 1 020.26 1 029.07 -8.81
    1 640 9 217.77 9 216.60 1.16
    1 641 391.98 393.52 -1.54
    1 642 6 292.59 6 302.44 -9.85
    1 643 16 364.30 16 360.85 3.45
    1 645 696.46 698.53 -2.06
    1 648 3 230.55 3 213.54 17.01
    1 649 4 078.84 4 061.89 16.95
    总计 243 533.03 243 533.03 0
    下载: 导出CSV
  • [1] 艾廷华, 杨帆, 李精忠.第二次土地资源调查数据建库中的土地利用图综合缩编[J].武汉大学学报·信息科学版, 2010, 35 (8): 887-891 http://ch.whu.edu.cn/CN/abstract/abstract1024.shtml

    Ai Tinghua, Yang Fan, Li Jingzhong. Land-use Data Generalization for the Database Construction of the Second Land Resource Survey[J]. Geomatics and Information Science of Wuhan University, 2010, 35(8): 887-891 http://ch.whu.edu.cn/CN/abstract/abstract1024.shtml

    [2]

    Aichholzer O, Aurenhammer F, Alberts D, et al. A Novel Type of Skeleton for Polygons[M]. Berlin, Heidelberg: Springer, 1996

    [3]

    Haunert J H, Sester M. Area Collapse and Road Centerlines Based on Straight Skeletons[J]. GeoInformatica, 2008, 12(2): 169-191 doi: 10.1007/s10707-007-0028-x

    [4]

    Lee D T. Medial Axis Transformation of a Planar Shape[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1982, PAMI-4 (4): 363-369 doi: 10.1109/TPAMI.1982.4767267

    [5]

    Cloppet F, Oliva J M, Stamon G. Angular Bisector Network, a Simplified Generalized Voronoi Diagram: Application to Processing Complex Intersections in Biomedical Images[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2000, 22(1): 120-128 doi: 10.1109/34.824824

    [6]

    Ruas A. Multiple Paradigms for Automating Map Generalization: Geometry, Topology, Hierarchical Partitioning and Local Triangulation[C]. ACSM/ASPRS Annual Convention and Exposition, Charlotte, USA, 1995

    [7]

    Ware J M, Jones C B, Bundy G L. A Triangulated Spatial Model for Cartographic Generalization of Areal Objects[C]// Kraak M J, Molenaar M. Advance in GIS Research Ⅱ (the 7th Int Symposium on Spatial Data Handling). London: Taylor & Francis, 1997a: 173-192

    [8]

    Delucia A A, Black R T. A Comprehensice Approach to Automatic Feature Generalization[C]. The 13th International Cartographic Conference, Morelia, Mexico, 1987

    [9] 艾廷华, 郭仁忠.支持地图综合的面状目标约束Delaunay三角网剖分[J].武汉大学学报·信息科学版, 2000, 25(1): 35-41 http://ch.whu.edu.cn/CN/abstract/abstract5070.shtml

    Ai Tinghua, Guo Renzhong. A Constrained Delaunay Partitioning of Areal Objects to Support Map Generalization[J]. Geomatics and Information Science of Wuhan University, 2000, 25(1): 35-41 http://ch.whu.edu.cn/CN/abstract/abstract5070.shtml

    [10] 艾廷华, 刘耀林.土地利用数据综合中的聚合与融合[J].武汉大学学报·信息科学版, 2002, 27(5): 486-492 http://ch.whu.edu.cn/CN/abstract/abstract4997.shtml

    Ai Tinghua, Liu Yaolin. Aggregation and Amalgamation in Land-use Data Generalization[J]. Geoma-tics and Information Science of Wuhan University, 2002, 27(5): 486-492 http://ch.whu.edu.cn/CN/abstract/abstract4997.shtml

    [11]

    Penninga F, Verbree E, Quak W, et al. Construction of the Planar Partition Postal Code Map Based on Cadastral Registration[J]. GeoInformatica, 2005, 9(2): 181-204 doi: 10.1007/s10707-005-6430-3

    [12]

    Jones C B, Bundy G L, Ware M J. Map Generalization with a Triangulated Data Structure[J]. Cartography and Geographic Information Systems, 1995, 22(4): 317-331 doi: 10.1559/152304095782540221

    [13]

    Uitermark H, Vogels A, van Oosterom P. Semantic and Geometric Aspects of Integrating Road Networks[M]. Berlin, Heidelberg :Springer, 1999

    [14]

    Touya G. Relevant Space Partitioning for Collaborative Generalization[C]. The 13th Workshop of the ICA Commission on Generalisation and Multiple Representation, Zürich, Switzertand, 2010

    [15]

    Touya G, Berli J, Lokhat I, et al. Experiments to Distribute and Parallelize Map Generalization Processes[J]. The Cartographic Journal, 2017, 54(4): 322-332 doi: 10.1080/00087041.2017.1413787

    [16]

    Thiemann F, Werder S, Globig T, et al. Investigations into Partitioning of Generalization Processes in a Distributed Processing Framework[C]. The 26th International Cartographic Conference, Dresden, Germany, 2013

图(9)  /  表(4)
计量
  • 文章访问数:  1563
  • HTML全文浏览量:  340
  • PDF下载量:  229
  • 被引次数: 0
出版历程
  • 收稿日期:  2018-07-14
  • 发布日期:  2018-12-04

目录

/

返回文章
返回