留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

球面三角格网单元面积的分形计算

王政 赵学胜 罗富丽 李亚路 孙中昶

王政, 赵学胜, 罗富丽, 李亚路, 孙中昶. 球面三角格网单元面积的分形计算[J]. 武汉大学学报 ● 信息科学版, 2020, 45(10): 1541-1546. doi: 10.13203/j.whugis20180478
引用本文: 王政, 赵学胜, 罗富丽, 李亚路, 孙中昶. 球面三角格网单元面积的分形计算[J]. 武汉大学学报 ● 信息科学版, 2020, 45(10): 1541-1546. doi: 10.13203/j.whugis20180478
WANG Zheng, ZHAO Xuesheng, LUO Fuli, LI Yalu, SUN Zhongchang. A Fractal Method for Area Calculation of the Spherical Triangular Grid[J]. Geomatics and Information Science of Wuhan University, 2020, 45(10): 1541-1546. doi: 10.13203/j.whugis20180478
Citation: WANG Zheng, ZHAO Xuesheng, LUO Fuli, LI Yalu, SUN Zhongchang. A Fractal Method for Area Calculation of the Spherical Triangular Grid[J]. Geomatics and Information Science of Wuhan University, 2020, 45(10): 1541-1546. doi: 10.13203/j.whugis20180478

球面三角格网单元面积的分形计算

doi: 10.13203/j.whugis20180478
基金项目: 

国家自然科学基金 41671394

国家重点研发计划 2018YFB0505301

详细信息
    作者简介:

    王政,博士,主要从事全球离散格网的性质与应用研究。wangzheng@ncst.edu.cn

    通讯作者: 赵学胜,博士,教授,博士生导师。zxs@cumtb.edu.cn
  • 中图分类号: P228

A Fractal Method for Area Calculation of the Spherical Triangular Grid

Funds: 

The National Natural Science Foundation of China 41671394

the National Key Research and Development Program of China 2018YFB0505301

More Information
    Author Bio:

    WANG Zheng, PhD candidate, specializes in the theories and applications of discrete global grid. E-mail: wangzheng@ncst.edu.cn

    Corresponding author: ZHAO Xuesheng, PhD, professor. E-mail:zxs@cumtb.edu.cn
图(5) / 表(2)
计量
  • 文章访问数:  838
  • HTML全文浏览量:  167
  • PDF下载量:  60
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-12-27
  • 刊出日期:  2020-10-05

球面三角格网单元面积的分形计算

doi: 10.13203/j.whugis20180478
    基金项目:

    国家自然科学基金 41671394

    国家重点研发计划 2018YFB0505301

    作者简介:

    王政,博士,主要从事全球离散格网的性质与应用研究。wangzheng@ncst.edu.cn

    通讯作者: 赵学胜,博士,教授,博士生导师。zxs@cumtb.edu.cn
  • 中图分类号: P228

摘要: 球面三角格网是一种多用途的全球离散格网,是数字地球领域空间信息表达的一种有效方案。针对球面三角格网单元面积计算的一般方法准确度不高的问题,提出一种面积计算的分形方法。在同一剖分层中,把球面三角格网重新编码为0或1,计算相同编码的格网单元的面积,拟合多重分形的面积计算函数,以此作为球面三角格网的面积计算模型。用相对中误差计算该模型的准确度,并用算例进行验证。结果显示, 随着多重分形的次数增加, 格网单元的面积计算准确度变高,面积计算准确度是可控的。

English Abstract

王政, 赵学胜, 罗富丽, 李亚路, 孙中昶. 球面三角格网单元面积的分形计算[J]. 武汉大学学报 ● 信息科学版, 2020, 45(10): 1541-1546. doi: 10.13203/j.whugis20180478
引用本文: 王政, 赵学胜, 罗富丽, 李亚路, 孙中昶. 球面三角格网单元面积的分形计算[J]. 武汉大学学报 ● 信息科学版, 2020, 45(10): 1541-1546. doi: 10.13203/j.whugis20180478
WANG Zheng, ZHAO Xuesheng, LUO Fuli, LI Yalu, SUN Zhongchang. A Fractal Method for Area Calculation of the Spherical Triangular Grid[J]. Geomatics and Information Science of Wuhan University, 2020, 45(10): 1541-1546. doi: 10.13203/j.whugis20180478
Citation: WANG Zheng, ZHAO Xuesheng, LUO Fuli, LI Yalu, SUN Zhongchang. A Fractal Method for Area Calculation of the Spherical Triangular Grid[J]. Geomatics and Information Science of Wuhan University, 2020, 45(10): 1541-1546. doi: 10.13203/j.whugis20180478
  • 球面三角格网灵活、可扩展,适用于多种地球系统模型[1-3],是数字地球或虚拟地球空间信息表达的一种有效方案[4-5]。由于(椭)球面是流形空间,在数学上与欧氏空间不同胚,直接在(椭)球面上无法递归剖分出像平面栅格那样完全相同的格网单元[6-7]。国内外学者研究了球面三角格网单元的面积计算方法。鉴于均匀格网的面积容易计算,很多学者致力于改进格网面积均匀性的研究。这种改进的格网主要有两种,即等积投影格网和面积优化格网[8-9]。等积投影格网指在平面上构造均匀格网,然后按照某种保持均匀性的地图投影方法把该平面格网逆投影到地球表面,形成面积相等的格网[10-11]。逆等积投影计算复杂,形成的格网在一定程度上改变了单元的性质。面积格网指把(椭)球面格网单元的大圆弧线即大地线用一些非大圆弧线取代,对格网单元的面积进行优化,形成面积近似相等的格网[7, 12]。不过,这种近似等积的格网只能满足气象预警、洋流模拟等面积精确度要求不高的应用需求,对格网计算精确度较高的应用需求就略显不足。

    为了解决上述问题,本文提出一种基于分形编码的方法计算球面三角格网单元的面积。格网单元的编码是格网系统处理问题的核心,格网单元编码运算以现代矩阵理论为基础,具有快速索引和高效计算的功能[13-14]。因此球面三角格网单元的面积计算必然以格网单元编码为基础。

    • 地理编码隐含地理实体的空间位置信息,在全球离散格网统一坐标系中,它是格网单元的唯一标识。地理编码即地理位置坐标X与几何属性Y之间具有直接的函数关系,即Y=fX)。本文构建球面三角格网单元面积Y与地理编码X之间的函数关系f

      本文选择正八面体作为基础多面体。对球面三角形的剖分采用直接球面剖分法,即分别取初始球面三角形的边弧中点,若依次用大圆弧连接相应中点,则构成4个新的子球面三角形,它具有递归嵌套的本质特征[15]。然后按照Goodchild等[16]提出的地理编码方法对格网单元进行初始地理编码,如图 1所示。

      图  1  球面三角形直接剖分

      Figure 1.  Direct Subdivision Method

    • 球面是一个不可展曲面,它具有与平面不同的几何性质。球面三角形的面积与其内角有关,面积表达式为$ \left( {\mathop \sum \limits_{i = 1}^3 {\alpha _i} - {\rm{\pi }}} \right) \times {R^2}$。直接剖分球面等边三角形得到的是大小和形状不完全相等的4个子球面三角形,经面积公式计算,其中1个子球面三角形面积较大,另外3个子球面三角形面积较小,如图 1所示。因此,直接剖分球面三角形的方法是不等面积的格网剖分法,剖分所得球面子三角形存在大小和形状的差异。

      然而,在应用中若按照不等面积格网计算其单元的面积,因为存在数量巨大的格网单元,那么格网单元的面积计算必然难以进行。因此,为了降低应用复杂度,在应用中把格网单元按照均匀格网计算[7, 17-18]

      $$ f\left( n \right) = {S_i}/{S_{ij}}, i = 0, 1, 2 \ldots {4^n} - 1;j = 0, 1, 2, 3$$ (1)

      根据式(1)分析格网Si与子格网Sij的面积数量关系,计算结果如图 2所示,图 2(a)表示第nfn)的最大值、最小值和平均值,图 2(b)表示第nfn)的平均值及变化率。

      图  2  球面三角形与剖分所得球面三角形面积比

      Figure 2.  Area Ratio of Spherical Triangle to That of Spherical Sub-triangles

      图 2可知,随着剖分层次的增加,fn)的平均值越来越趋于常数4,格网单元面积变化收敛。从格网单元面积的极值和平均值的变化来看,对一个球面三角形的剖分即第一次剖分,所得到的格网单元面积的差异最大,随着剖分次数增加,所得到的网格单元面积的差异逐渐减小。故格网单元面积的差异性与其父格网单元面积的大小具有正相关性。面积初始值相同的球面三角形可以归为一个子集,分别计算子集中的三角形面积变化。

      分别把球面三角形面积较大的子集记为BG、面积较小的记为SL,如图 3所示,BG对应初始地理编码X中以数字0开头的三角形,即X∈[0,4n-1-1],SL对应初始地理编码X中以数字1、2和3开头的三角形,即X∈[4n-1,4n-1]。

      图  3  球面三角格网单元的子集特征

      Figure 3.  Cell Subsets of Spherical Triangular Grids

    • 根据上述特征,分别计算第1层球面三角格网单元组成的两个子集。由于子集内的三角格网单元的面积Li差异变小,因此把子集中的格网单元视为等面积单元,把面积的平均值L或数学期望E作为该组格网单元面积的真值$\widetilde L $。

      取CGCS2000(China geodetic coordinate system 2000)的长轴作为均匀球体半径R,则球面三角形第一次剖分的面积分别为SBG=0.551R2SSL=30.339 8R2,第n层格网中两个子集的球面三角形个数分别为NBG=4n-1NSL=34n-1,则第n层球面三角形面积的数学期望分别如下:

      $$ \begin{array}{l} {{\tilde L}_{{\rm{BG}}}} = {E_{{\rm{BG}}}}\left( n \right) = {S_{{\rm{BG}}}}/{N_{{\rm{BG}}}} = 2.205 \times {R^2} \times {4^{ - n}} = \\ \;\;\;\;\;\;\;8.971 \times {10^{13}}{{\rm{e}}^{ - 1.386n}}, n = 1, 2, 3 \cdots \end{array}$$ (2)
      $$ \begin{array}{l} {{\tilde L}_{{\rm{SL}}}} = {E_{{\rm{SL}}}}\left( n \right) = {S_{{\rm{SL}}}}/{N_{{\rm{SL}}}} = 1.359 \times {R^2} \times {4^{ - n}} = \\ \;\;\;\;\;\;\;5.529{\rm{}}9 \times {10^{13}}{{\rm{e}}^{ - 1.386n}}, n = 1, 2, 3 \cdots \end{array}$$ (3)

      式(2)和式(3)的计算准确度用中误差()和相对中误差()表示。对于第n层格网,球面三角形的面积观测值为Li,由式(2)和式(3)计算的面积值为$\widetilde L $,则观测值Li与计算值$\widetilde L $的差值为,则中误差公式为$\sigma = \sqrt {\left[ {\Delta \Delta } \right]/m} $(m为该层三角形总数),相对中误差公式$\gamma = \sigma /\tilde L \times 100{\rm{\% }}$,计算结果见表 1σAvgγAvg分别表示均匀格网法计算的面积中误差和相对中误差。

      表 1  子集模型与均匀格网模型的准确度

      Table 1.  Accuracy of the Subset Model and the Uniform Grid Model

      层次 σBG/m σSL/m σAvg/m γBG/% γSL/% γAvg/%
      1 1 431 543.379 3 341 614.090 4 300 933 108 340.890 0 0 26.922
      2 408 424 739 196.440 186 702 379 955.430 998 634 254 578.553 7.285 5.402 25.005
      3 105 219 774 196.952 48 056 323 978.280 244 182 169 430.843 7.507 5.562 24.456
      4 26 352 081 047.338 12 035 043 972.553 60 693 648 203.963 7.520 5.571 24.315
      5 6 588 752 511.296 3 009 087 071.344 15 151 275 607.137 7.521 5.572 24.280
      6 1 647 199 554.119 752 276 858.055 3 786 433 169.730 7.521 5.572 24.271
      7 411 800 067.058 188 069 294.077 946 521 649.861 7.521 5.572 24.268
      8 102 950 019.569 47 017 324.776 236 624 996.769 7.521 5.572 24.268
      9 25 737 504.940 11 754 331.217 59 155 910.703 7.521 5.572 24.268
      10 6 434 376.237 2 938 582.806 14 788 956.520 7.521 5.572 24.268

      表 1可见,两种计算方法计算的相对中误差都是收敛的,依次近似收敛于7.521%和5.572%以及24.268%。基于子集法计算的格网单元面积的中误差和相对中误差均小于均匀格网法对其计算的结果。

    • 准备阶段:在地球上均匀地选择4个区域(非行政区域,分别为GL、US、CH、AS)作为试验对象,如图 4所示。图 4中的绿色区域表示在UTM投影下的试验对象及其位置。试验对象的经纬度范围见表 2

      图  4  UTM投影下的试验对象试验对象示意图

      Figure 4.  Test Objects in UTM Projection

      表 2  试验对象范围及面积预测相对误差

      Table 2.  Test Objects Ranges and Their Relative Error of Areas

      实体 试验对象位置范围/(°) 面积预测相对误差/%
      东经 北纬 东经 北纬 平均值法 子集法 文献[7]
      GL -73.042 83.627 -11.312 59.977 18.249 2.331 4.222
      US -124.733 49.389 -66.955 24.947 9.239 0.603 1.667
      CH 73.499 53.558 135.087 18.161 5.279 0.485 3.921
      AS 113.147 -10.693 153.639 -39.139 29.767 2.926 1.818

      地理实体格网化:在本试验中,选择一种均匀性较好的等积格网以便与本文方法作比较,这里选择文献[7]中等积格网。分别用第7层球面三角格网单元和第7层文献[7]中等积格网单元各自表达试验对象。

      面积计算:统计覆盖试验对象的格网数量,分别计算试验对象的面积。计算结果见图 5,图中横坐标表示不同的计算方法;纵坐标表示试验对象的面积(单位:m2),右侧标记了试验对象名称。

      图  5  基于不同方法计算的试验对象的面积

      Figure 5.  Areas Calculated by Different Methods

      准确度分析:计算试验对象的统计面积与其真实面积的相对误差,以此描述结果的准确度,计算结果见表 2

      试验结果显示本文方法计算的试验对象的面积普遍接近其真实值,其次是采用文献[7]等积格网的计算方法,球面三角格网的直接平均值法计算的试验对象面积与其真实面积的差异最大。

    • 为了更加清楚地表达分形,把格网单元的地理编码转换为用0和1表达的分形编码。对于第1层格网单元,上文中已阐述分形的方法,即分为BG和SL两个子集,此处把BG称为分形编码为0的子集,把SL称为分形编码为1的子集。分形编码为0的格网单元的进一步剖分得到分形编码为00和01的两个子集;分形编码为1的格网单元进一步剖分得到格网单元分形编码为10和11的两个子集。以此类推,对其他层的分形编码的方法相似。

      对于编码为A的任意格网单元,若进一步提高计算准确度,则必然产生分形编码A0和A1的格网单元。若给出AA0的格网单元的分形维数,则可以在A已知时计算任意A0的格网单元面积;同理,可计算任意A1的格网单元面积。以分形编码0为例,它由N位0编码的球面三角形构成,其面积组成的分形为$ 0, 00 \cdots \underbrace {0 \cdots 0}_N$。把面积作为纵坐标(m2)、格网单元所在的层作为横坐标,则可以拟合出格网单元的面积函数。

    • 分形编码面积计算模型的准确度用中误差和相对中误差来表示,对于像地理编码为0的这种区域内的准确度计算方法与上文基于单元子集的计算方法相同;对于像地理编码为1、2、3的这种跨区域的准确度,采用误差传播定律计算。

    • 根据分形编码的计算方法,可以计算任意分形次数和任意编码的格网单元的面积。首先根据分形编码计算初始剖分面积。初始剖分面积越小,剖分所得格网单元面积预测越准确。以分形编码为0开头的格网单元的面积计算为例予以说明。

      算例:按照第3次和第5次分形计算球面三角形0203012的表面积。

      地理编码转分形编码。本算例中把0203012转为0101011。

      根据分形次数计算初始面积。本算例中初始分形编码分别为010和01010,根据初次分形编码格网单元的面积模型计算得到。

      第3次分形编码法初始面积为:

      S010=1 381 791 553 080.890 m2

      第5次分形编码法初始面积为:

      S01010=8 632 105 890.887 m2

      根据格网单元子集是均匀单元的判据,计算最终球面三角形的面积。本算例中,球面三角形编码为7位即第7层格网,第3次分形即把第3层相应编码的格网单元作为初始剖分单元,第7层球面三角形个数为44;同理把第5层格网单元作为初始剖分单元,第7层球面三角形个数为42。按照子集内的单元是均匀的假设分别计算0101011的面积。

      第3次分形后:

      S0101011=S010/44=5 397 623 254.222 m2

      第5次分形后:

      S0101011=S01010/42=5 395 066 181.430 m2

      换算为原来地理编码的面积。球面三角形0203012的估算面积如下:

      第3次分形编码法:

      S0203012=5 397 623 254.222 m2

      第5次分形编码法:

      S0203012=5 395 066 181.430 m2

      评估相对误差。地理编码为0203012的真实面积S=5 386 470 314.976 m2,则估算面积与真实面积的两种分形编码法相对误差分别如下。

      第3次分形编码法:

      $$ \begin{array}{l} {\delta _2} = \left( {\left| {{S_{0203012}} - S} \right|/S} \right) \times 100{\rm{\% }} = 0.207{\rm{\% }} < \\ \;\;\;\;1.805{\rm{\% }} = {\sigma _{01}} \end{array} $$

      第5次分形编码法:

      $$ \begin{array}{l} {\delta _4} = \left( {\left| {{S_{0203012}} - S} \right|/S} \right) \times 100{\rm{\% }} = 0.160{\rm{\% }} < \\ \;\;\;\;0.441{\rm{\% }} = {\sigma _{0101}} \end{array} $$

      上述算例中,采用了两种分形公式计算了球面三角形的面积,相对中误差高于模型误差。

    • 基于子集法计算的试验对象面积的相对误差值与均匀格网法以及文献[7]法计算的结果作比较。采用子集法计算的结果误差比均匀格网法小,前者大约是后者的1/10,即本文子集法比均匀格网法计算的面积准确度提高了大约10倍。子集法在初始球面三角形的中心位置的计算准确度更高,而越接近初始边角位置面积的准确度越低。经比较试验对象的面积准确度,该方法比文献[7]方法普遍较优,只有试验对象AS的面积误差比文献[7]方法计算的大。

      子集法计算面积的相对中误差分别收敛于7.521%和5.572%,这个数值不一定能满足高精确计算的需要。本文进一步提出分形编码的方法,以便提高面积计算的准确度满足更高精度计算的需要。然而基于子集法的面积计算是分形编码法的基础。根据分形编码模型的准确度变化情况,发现分形次数越多,格网单元计算的面积越接近其真实面积,因此分形的次数和面积计算的准确度具有正相关性。

    • 采用子集法计算球面三角网格单元面积能提高其面积计算准确度。在应用中,该方法把球面三角形分成两个子集,在每个子集内认为格网单元是均匀的,按照均匀格网来计算单元面积。该方法比直接均匀格网法计算球面三角格网单元的准确度提高了大约10倍,提高效果显著。

      该模型对球面三角形的中心区域适用性较好。球面三角格网模型是一种表面模型,剖分该表面所得格网单元的面积差异也与被剖分单元的内角的大小和差异有关。因为初始球面三角形的顶角固定,随着剖分的进行,包含该角的格网单元的内角差异越来越大,单元形状变化增大,容易造成计算误差的增加。

      本文研究了基于球面的格网单元的面积计算模型。鉴于椭球面是大地测量和计算的基准面,且中国现已应用和推广的CGCS2000坐标系也是以参考椭球面为基准面的坐标系,故基于椭球面的面积计算模型是下一步要研究的一个重要问题。

参考文献 (18)

目录

    /

    返回文章
    返回