基于线性GSI二维半变异函数各向异性结构建模及估计研究——以DEM数据为例

删除或更新信息，请邮件至freekaoyan#163.com(#换成@)

本站小编 Free考研考试/2021-12-29

高歆许昌学院城市与环境学院,许昌 461000

Anisotropic modeling and estimation for a two-dimensional semi-variogram based on the linear GSI Model: Taking DEM data as an example

GAO XinCollege of Urban and Environmental Sciences, Xuchang University, Xuchang 461000, Henan, China

收稿日期:2019-09-12修回日期:2020-03-11网络出版日期:2020-11-20

基金资助:

河南省高校人文社会科学研究一般项目.2021-ZDJH-0340
河南省高校人文社会科学研究一般项目.2020-ZZJH-426
辽宁省自然科学基金.2019-MS-342
河南省科技攻关项目.182102310924
河南省高等学校重点科研项目.21A170019

Received:2019-09-12Revised:2020-03-11Online:2020-11-20

摘要
鉴于传统各向异性二维半变异函数各向同性化方法未充分考虑或无力精确描述其内部结构信息的缺陷,本研究通过引入线性广义尺度不变（GSI）模型,以DEM数据作为验证对象,对二维半变异函数各向异性结构信息进行多尺度建模,并采用旋转椭圆法、两步搜索作图法等方法对系统参数进行估计,最后以球状模型为例对理论半变异函数的估计精度,及其在空间数据插值中的应用效果进行对比研究。结果表明：各向异性普遍存在于地形数据的空间变异中,有证据表明,这种各向异性结构中处处显现出不同的变形特征,但是也存在着某种规则性的成分,如各向同性圆形或近圆形等值线,因此,在对坐标进行各向同性化处理时不适合采用“一刀切”的方式去处理;GSI系统参数皆能得到较高精度的估计,如决定系数R²普遍达到了0.99以上,间接证明了GSI模型对地形数据各向异性结构处理的有效性和适用性;通过理论模型估计和插值结果对比,线性GSI坐标转换法比传统坐标转换法有了明显的精度提升,并且展现出了较高的边缘信息恢复能力,但也表现出了一定的局限性和不稳定性。
关键词： 广义尺度不变性;半变异函数;各向异性;数字高程模型;坐标转换

Abstract

Anisotropy has been found to widely exist in the nature, and also regarded as one of the several essential attributes of geographical phenomena and processes. Therefore, it needs some complex models and methods to analyze and explain these phenomena, and deal with such problems as the optimization and interpolation for discrete monitoring points, or the uncertainty analysis by stochastic simulation over the space for some regional variables. The traditional treatment on an anisotropic modeling in kriging interpolation based on coordinate transformation does not fully consider or cannot accurately describe the internal structures of a two-dimensional anisotropic semi-variogram. Therefore, this study introduced a model named as linear generalized scale invariance (GSI) to simulate the anisotropic information of a 2D semi-variogram by using DEM as input data, while the system parameters were estimated by using the rotating ellipse and two-step search mapping method, and the comparisons of the two methods including GSI and traditional coordinate transformation applied in the fittings to the spherical model and the corresponding kriging interpolation were also made. The results firstly showed that anisotropy is common and ubiquitous in the spatial variability for topographic data, and the complexity is characterized by a change for an anisotropic ration when the corresponding scale changes, i.e. the different deformation behaviors over the whole semi-variogram maintained. However, some evidences showed that there are some regular features like an isotropic component existing in the anisotropic mechanism, such as a circular scale or circular contour. When facing such a complex structure, simple and rough treatment is obviously not enough. Secondly, the related parameters in GSI model can be estimated with high accuracy, such as the values of R², which are all almost over 0.99 for the six regions. This indirectly proved that the validity and applicability of GSI model in the treatment to the anisotropic structure. In addition, for the fittings of the theoretical spherical model, the GSI model showed huge advantage over the traditional transformation and isotropic methods. Finally, as the enhanced effect originated from application of the GSI model in the interpolation processes, the coordinate transformation based on the linear GSI model had better improvement in accuracy than the traditional coordinate transformation, as well as a high ability of edge information recovery, although it exhibited some limitations and instability due to its complex covariance structure.

Keywords：generalized scale invariance;semi-variogram;anisotropy;DEM;coordinate transformation

PDF (6772KB)元数据多维度评价相关文章导出EndNote|Ris|Bibtex 收藏本文
本文引用格式
高歆. 基于线性GSI二维半变异函数各向异性结构建模及估计研究——以DEM数据为例. 地理研究[J], 2020, 39(11): 2607-2625 doi:10.11821/dlyj020190794
GAO Xin. Anisotropic modeling and estimation for a two-dimensional semi-variogram based on the linear GSI Model: Taking DEM data as an example. Geographical Research[J], 2020, 39(11): 2607-2625 doi:10.11821/dlyj020190794

1 引言

在自然界中,地理或者环境监测变量取值常常布满整个二维或者三维空间,但是受制于技术和成本等因素,人们总习惯于在一些重点区域布设点位进行观测,显然这些离散不连续的观测值无法反映地理环境因子在空间上的分布和变异特征。如果想要获取其他非监测点位的属性值,只有通过数学或者物理模型的方法去估计。其中,插值是一种将点源数据转化为面状数据常见的数学手段,常用的插值方法有反距离权重、克里金和径向基函数等方法^[1,2,3],其中克里金插值因其无偏估计和误差度量成为最受欢迎的方法之一。对于传统的克里金插值,常见的做法是首先在所有方向上搜索不同距离的坐标对以绘制半变异函数曲线,然后采用理论函数进行拟合,确定出权重后进行估计,显然这种做法忽略了方向信息,它将所有已知点的权重仅仅视为距离的函数。然而,自然界中的地理现象在空间分布上除了存在一定的连续性以外,还呈现着明显的差异性和不规则性,例如,在各种内外营力作用下,地球表面在不同方向上表现出截然不同的地形地貌特征^[4,5,6],每一点的高程梯度也因位置和方向的不同而不同。因此,将方向信息融入到半变异函数模型中将对克里金插值的精度提升具有重要意义^[7]。

各向异性克里金插值的关键是对各向异性半变异函数的处理,过去经验表明,自然界中的各向异性变量或者过程几乎都可以划分为几何各向异性或者带状各向异性,前者代表变程在不同方向上变化（各个方向拥有相同的基台值）,而后者意味着在不同方向上存在着不同的基台值^[8]。Zimmerman D L对此提出了不同的看法,他从更一般的角度进行了研究,认为将各向异性半变异函数类型划分为变程各向异性、基台各向异性和块金各向异性更为合理^[9]。以几何各向异性为例,对各向异性处理的方式主要有两种,第一种是传统方法坐标转换法^[10,11],基本思想是通过变程椭圆的长短轴比和倾角等参数建立各向同性转换模型,然后使用转换后的已知点坐标值和属性值进行参数估计获取理论模型参数,这种方法适用于多种各向异性类型。例如,Eriksson M等从椭圆几何视角以一般的方式基于坐标转换对变程、块金和基台等多种各向异性类型的处理方式进行了较为全面的总结和研究,尤其是从公式对比方面指出了已有程序在处理幂律各向异性时可能存在的问题^[12]。由于半变异函数是以两点之间的欧氏距离（向量的模）作为自变量,显然对于任意两点之间的坐标差△x和△y,都需要将其处理成距离,对于各向同性来说,处理还算简单,但是对于各向异性来说,某种程度上限制了半变异函数的扩展,从而也对变程等值线的形状和各向异性模型的使用构成了一定的束缚。第二种方法是可分离函数乘积法^[13],其简化了各向异性处理操作,基本思想是采用可化为x轴和y轴一维函数乘积的二维协方差函数直接估计,无需将各个方向的变程化为相等,但需要将x轴或y轴旋转至主变形方向。近年来相关理论也得到了一定程度的发展,例如,Ecker M D,Shen Y等提出了一种基于贝叶斯理论的各向异性参数估计模型,这种模型在参数估计精度提升实施策略方面显示出了一定的优势^[14,15];Hristopulos D T等以超椭圆函数形式对协方差函数进行参数化,根据阶数n的不同,其形状也各不相同,并且通过各向异性随机场协方差二阶导数张量和协方差函数中各向异参数之间的恒等关系提出了一种非参数估计方法^[16]。虽然可分离函数乘积法定义较为简单,但是与其对应的随机场却较为稀少,相反,实际生活中碰到的绝大多数随机场,其协方差函数都是不可分离的^[17]。为了摆脱传统各向异性模型的束缚,Allard D等革新了其构造思想,他们从半变异函数或者协方差函数谱函数表达式出发,通过构建方向积分函数,并将方向测度表达成一组核函数的和,从而根据核函数的类型使一组新各向异性模型的建立成为可能^[18]。经过近几十年的发展,各向异性协方差函数或者半变异函数相关理论中融入了许多新思想和新理论,它的一些基本概念和内涵需要重新认识和思考,De Iaco S等综述了半变异函数理论中一些基本概念如可分离性、各向同性、对称性和严格正定性等特征之间的逻辑关系,并得到一些有意义的结论,例如,定义在部分重叠定义域上的协方差函数可以用于对非几何各向异性建模,这类函数是严格正定的,但是不可分离的^[19]。

综上所述,不管是坐标转换法还是可分离函数法,均以如何构造变程等值线的形状作为其基本出发点,然后将这一套变异机制用于所有的已知点,并未顾及该变异机制与内部点的变形特征相符与否。如图1所示,假设有一模拟半变异函数等值线（图1a,c = -0.1, e = -0.2, f = 0.1）,显然传统处理方式仅考虑了变程等值线（最外层等值线）的变形特征（图1b）,即处理过后仅有最外层等值线为圆形,完全没有顾及变程内部等值线的方向信息,那么这样拟合后的结果其内部仍是各向异性的,并没有真实反映实际半变异函数（图1a各向同性化后的曲线应为图1c）。显然,这种“一刀切”的方式是不合理的,如果这种变形特征是一种多尺度的行为,即变形特征随尺度的变化而变化,那又应当如何处理呢?鉴于这一认识,本研究将从多尺度坐标转换的视角来探讨这个问题,以坐标转换处理方式和常见变形类型——几何各向异性为例来进行说明。

图1

新窗口打开|下载原图ZIP|生成PPT
图1传统坐标转换法和GSI坐标转换法对各向异性半变异函数等值线处理对比

Fig. 1A comparison between the traditional method and the GSI method for treating isolines of an anisotropic semi-variogram before the performation of an estimation

2 研究数据和方法

2.1 研究区域和数据

从全国主要山地中选择六个区域进行验证,其位置、所属山系及地形地貌特征总结见表1。研究数据采用由美国航天局（NASA）与日本经济产业省（METI）共同推出的ASTER GDEM（Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model,先进星载热发射和反射辐射仪全球数字高程模型）,其分辨率为30 m,验证区域大小统一设置为512×512像元,面积约236 km²（图2）。

Tab. 1
表1
表1研究区域的位置和地形地貌特征
Tab. 1The locations and terrain features of the six study regions

研究区域	位置	所属山系	经纬度
区域1	陕西省北部	黄土高原丘陵沟壑区	37°N 109°E
区域2	陕西省和湖北省交界处	秦岭山地	33°N 110°E
区域3	陕西省和重庆市交界处	西北-东南走向的大巴山地	32°N 108°E
区域4	湖北省和湖南省交界处	东北-西南延伸的武陵山地	29°N 108°E
区域5	云南省北部	青藏高原横断山地	28°N 98°E
区域6	福建省西南部	武夷山地	26°N 117°E

新窗口打开|下载CSV

图2

新窗口打开|下载原图ZIP|生成PPT
图2研究区域DEM数据

Fig. 2The DEM data of the six study regions

2.2 GSI模型

GSI模型（Generalized Scale Invariance）,也叫广义尺度不变模型,是一种幂律模型,它是Schertzer D等在使用分形理论研究湍流现象中的尺度不变特征提出来的,他们认为在自然界中除了一些简单的自相似现象以外（如两个物体经过放大和缩小变得相同）,还存在一些复杂的自相似现象（经过旋转、拉伸和缩放后变得相同）^[20,21]。如前所述,坐标转换的核心是旋转和尺度缩放操作,同样这两种操作也是GSI的核心思想,它能够通过这两种操作兼顾全部等值线的变形信息,从而使半变异函数结构从整体上得到控制。GSI模型中比较难理解的是尺度概念,同其他幂律模型相比,GSI中尺度定义的方式有所不同,一般幂律模型中尺度是定义在某个支撑区间上,如污染物浓度值采样的某个时间段、矿石品位估计的矿块体积等物理量上,而GSI将位于同一条等值线（如图1中的圆和椭圆,也可以是很复杂的曲线）上的点视为同一尺度,即一条等值线代表一个尺度,相同尺度上点位的函数值都是相等的。那么这种尺度应如何度量呢?对于分布在圆上的点,可以利用该点到原点的距离作为其尺度的度量,即半径（圆上每个点都有一个尺度值,且都相等）。按照这种定义,对于一些复杂几何曲线如椭圆上的点来说,它们到原点的距离是不等的,显然这种定义在GSI模型中是行不通的,因此需要提出一种新的尺度度量函数或范式。通常在操作欧氏距离中的点时,关注的重点是坐标值和其坐标构成的向径距离（尺度）,不同的是,在这里操作这些点时,既要习惯于其在欧氏坐标系中的坐标,同时还要关注其位于某条等值线上,因为后者决定了其尺度的大小（非欧氏距离）。显然,可以很轻易地将这种尺度定义与整个区域中的单元格一一对应起来,对于一个研究区域来说,用无数条等值线去覆盖整个区域,每条曲线代表了一个尺度,也就意味着对于研究区域中的任意一点,都能找到划过该点的等值线,同时也就确定了它的尺度。如果该系统是自相似的,那么这些不同尺度上的函数值就可以通过幂律准则建立起相互之间的关联,意味着一旦已知其中某个尺度上的函数值,就能计算出所有尺度上的值。

下面给出GSI模型的符号描述,如上所述,要建立一个GSI模型,则需要确定三个要素：① 尺度改变算子T_λ;② 一个基准圆或者球E₁和一组由其生成的等值线E_λ;③ 等值线尺度的度量函数φ。幂律关系是湍流等地球物理领域中的一个基本准则,例如在湍流域中,结构函数满足S(λ^-1x)=λ^-2^HS(x),S为结构函数,λ是尺度比（注意与Λ区分）,有λ=φ(E₁)/φ(E_λ),H代表幂阶数,x为位置或者距离向量。上述模型可视为各向同性的,Schertzer D等在此基础上,将其调整为S(T_λx)=λ^-2^HS(x),其中T_λ=λ^-^G是尺度算子,G是一矩阵,λ^-^G为一矩阵函数,需要将其化为指数形式然后通过幂级数展开进行计算^[20,21]。可以看出,一般结构函数中的各向同性尺度被一个更一般的尺度算子λ^-^G所取代,意味着一组同心圆等值线被一组可以任意复杂的曲线所代替。GSI模型本质上就是将一个尺度转换算子λ^-^G不断地作用在一个圆或者其他曲线上,随着λ的变化从而生成一组同心等值线的过程,这组等值线仅仅用G中少数的几个参数就关联了起来。下面看一下φ的定义,对于一个圆或者球尺度φ(E₁)（基准尺度）来说,很容易找到一个正数来定义,如圆或者球的半径,而对于一个不规则形状的曲线,应如何定义呢?从圆或者球半径的定义受到启发,已知圆或者球半径与圆面积或者球体积之间满足幂函数关系,如其可以表示为φ^D(E₁),D为该半径φ(E₁)的幂。因此,这里可以采用对不规则曲线的面积或者体的体积进行幂运算来获取其尺度,且结果肯定是一正数,如果该不规则曲线为一椭圆,其面积为A,则其尺度可以定义为A^1/2。

根据矩阵G的形式,GSI模型可以分为线性GSI和非线性GSI,有研究表明^[22]¹⁴,线性GSI基本上能够解释自然界中大多数现象的各向异性变形信息,这里只讨论线性GSI,有关非线性GSI可以参考其他文献^[23,24,25]。这里的“线性”跟线性代数中的概念一致,即一种可以将转换过程近似为线性转换的特殊形式,每一次操作意味着对向量空间中的变量进行线性转换。在一个二维线性GSI模型中,G可以表达成四个简单矩阵的线性组合：

（1）

G = d E + c K + e I + f J

式中：四个矩阵分别是E=

[\begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix}]

;K=

[\begin{matrix} 1 & 0 \\ 0 & - 1 \end{matrix}]

;I=

[\begin{matrix} 0 & 1 \\ 1 & 0 \end{matrix}]

;J=

[\begin{matrix} 0 & 1 \\ - 1 & 0 \end{matrix}]

;c决定了椭圆长短轴的比例;e决定了旋转;f起着拉伸和压缩的作用;d为空间域维数有关的参数,d=1。关于c、e和f三个参数的拉伸和旋转效应,可参考图3（以图1a为基准）。要找出GSI的显式形式,需要对幂矩阵函数λ^-^G的指数形式exp(-Glnλ)进行幂级数展开,展开结果如下：

（2）

T_{λ} = \overset{\infty}{\sum_{n = 0}} (\ln λ^{- 1})^{n} G^{n} / n!

式中：λ为尺度比,即圆尺度和各级等值线尺度的比。

图3

新窗口打开|下载原图ZIP|生成PPT
图3GSI中参数的变化特征

Fig. 3The spatial variations of the elliptical contours with the parameters changing in the GSI model

下面看一下将上述尺度转换机理引入到半变异函数中的幂函数形式,对于半变异函数,其估计公式为：

（3）

γ^{*} (h) = \frac{1}{2 N (h)} \overset{N (h)}{\sum_{i = 1}} [Z (x_{i}) - Z (x_{i} {+ h)]}^{2}

式中：Z(x_i)为某点的函数值,这里指的是区域某点的高程值（m）;N为某确定间隔的样点对数量;x_i为某点的向量坐标（x_i, y_i）;h为两点之间的增量,有h=(△x_i,△y_i)。如果半变异函数具有幂律特征,那么变换后的公式为γ(T_λh)=λ^-^sγ(h),s为幂阶数。

总之,GSI模型使用少量参数来关联各级等值线,将它们纳入到一个统一的框架中,使从全局上控制半变异函数成为可能。这里需要强调的是,GSI模型不仅仅是一种幂律模型,更重要的是,它提供了一种坐标转换的思路和方法,一种全局各向异性结构各向同性化的方法。如前所述,GSI系统结构是尺度算子λ^-^G不断地作用在一个圆形曲线上从而生成的一组各向异性复杂曲线,那么它的逆运算就是将一组各向异性的复杂曲线转换成了同一个圆,然后再分别除以各级等值线的尺度比,从而就实现了各向同性化,后面将会把这种坐标转换方法应用于理论半变异函数模型的拟合和空间数据插值中。

3 结果与讨论

3.1 二维半变异函数计算

与一维半变异函数相比,二维函数不仅需要计算两点之间的相关性,同时还要考虑方向上的变化特征。对于一些稀疏采样数据集,存在某个方向或某个位置上因样本过少而估计不稳定的问题,为了提高其估计的稳定性,在计算之前往往需要按照一定的步长对角度和距离进行分组合并。与之相比,本研究数据为固定间隔且连续分布的栅格数据,不存在采样点稀疏的问题,为了充分发挥计算机矩阵计算的优点,这里采用滑动矩阵的方法进行计算。具体计算过程如下：以某4×4矩阵的试验数据M为例,首先分别准备3个矩阵,一个是待填充的半变异函数空矩阵A,边长为2n-1（图4a,n为M的行列栅格个数,此处行列数目相等,A中的任意一点代表一个滞后距）,另外两个是行滑动矩阵B和列滑动矩阵C,B中每个栅格值为其行号,C中每个栅格值为其列号,边长分别为n,如图4b和图4d所示;其次,筛选矩阵A中某给定点满足条件的坐标对及计算该点函数值,以P点为例,计算过程如下：第一步,计算P点同原点之间的行列偏移量x_O-x_P、y_O-y_P;第二步,计算行偏移量x_O-x_P和行滑动矩阵的和,及列偏移量y_O-y_P和列滑动矩阵的和,从而得到两个新矩阵（图4c和图4e）,即行滑动增量矩阵B'和列滑动增量矩阵C';第三步,在B'和C'中找出既落在行定义域范围（图4c中的红色框）又落在列定义域范围内（图4e中的红色框）的交集（图4c和图4e中的绿色框）,落在交集中的栅格值即为以P点为滞后距的终点坐标,而与此对应的行列滑动矩阵B和C中的栅格值则为起点坐标（图4b和图4d中的绿色框）;最后,从M中取出与坐标对相对应的值代入到半变异函数公式中计算即可;以此类推,可以得出全部滞后距函数值的估计。

图4

新窗口打开|下载原图ZIP|生成PPT
图4二维半变异函数计算

Fig. 4A diagram to depict how to calculate a 2D semi-variogram by using a sliding matrix method

为了凸显研究结果,研究对象的工作范围需要有所限制,这里采用变程作为确定外边界的参考依据。以区域1为例,具体实施如下：为了提升各方向上半变异函数在原点附近的拟合精度,这里采用二级套合球状模型进行拟合（图5）,在拟合曲线上根据最小基台值计算每个方向上的变程,用所有方向变程中的最大值作为参考边界的下限。经计算,区域1的最大变程约为678.61 m,位于π/2方向上。另外,由于在球状模型拟合过程中,会产生平滑效应,尤其会造成大尺度上变异信息的损失,因此,在具体实施的时候,研究范围可以适当放大。关于区域非平稳性对半变异函数的影响,这里采用一阶多项式去趋势,所有的估计结果都是建立在对数据去趋势的基础上。图6显示了6个区域的估计结果（显示范围不小于上述确定的最大变程值）,从图上可以看出半变异函数呈现出明显的各向异性,并且这种特征随着尺度的变化而变化,证明了地形空间异向分布的普遍性和复杂性,因此,在使用理论模型估计和插值的时候,需要重视这些因素。

图5

新窗口打开|下载原图ZIP|生成PPT
图5区域10~π间8个方向上的半变异函数

Fig. 5The semi-variograms along the eight directions within 0 and π for the region 1 based on nested structures

图6

新窗口打开|下载原图ZIP|生成PPT
图6六个区域半变异函数估计

Fig. 6The estimates of semi-variograms for the six regions

3.2 系统参数估计

由于系统参数较多,并且涉及到矩阵函数运算,因此GSI参数估计不能通过解析模型直接进行,而是需要依靠数值优化的方法分阶段逐步完成^[22,26]^45,40。在参考前人成果的基础上,这里对估计过程进行了简化,实施过程分为两步：第一步,先采用旋转椭圆法对半变异函数进行拟合,计算不同尺度上点集最优拟合椭圆的几何参数,得到的这组椭圆称为样本椭圆,从样本椭圆的几何参数中就能得出尺度等参数的估计;第二步,在第一步参数估计的基础上,继续采用旋转椭圆法或者加权旋转椭圆法对样本椭圆进行拟合,拟合后的结果称为GSI椭圆,根据样本椭圆和GSI椭圆的一致性建立目标函数,就可以得出其余参数的估计。

3.2.1 尺度参数估计首先介绍一下旋转椭圆法的基本思想,如图7所示,假设有一个以圆点O为中心的椭圆对某样本集进行拟合,该椭圆方程的一般形式为Cx²+2Dxy+Fy²=1,C、D和F是椭圆方程系数,也是待估参数。对于一样本点P(x_P, y_P),有OP与x轴之间的夹角为φ_P,以圆点为中心沿该角度绘制的直线与椭圆相交于P',OP'为极径ρ,P'点坐标为x_P'=ρcosφ_P、y_P'=ρsinφ_P,该点的极径ρ用参数的形式可表示为ρ(φ_P)=(Ccos²φ_P+2Dsinφ_Pcosφ_P+Fsin²φ_P)^-1/2。P到原点之间的距离OP等于

\sqrt[]{x_{P}^{2} + y_{P}^{2}}

,记为k_P,OP和OP'之间的偏差可以表示为

d_{P} = k_{P}^{2} - ρ^{2} (φ_{P})

,如果采样点数目为n,那么本样本集的平均误差可以表示为

E^{2} = \frac{1}{n} \overset{n}{\sum_{P = 1}} d_{P}^{2}

。

图7

新窗口打开|下载原图ZIP|生成PPT
图7样本椭圆拟合

Fig. 7A diagram to depict how to fit the samples on a scale with an ellipse

根据以上误差方程和最小二乘法基本原理,分别对三个参数进行求导,从而可以得出以下参数C、D和F的估计方程^[22]⁸⁰：

（4）

\begin{array}{l} (\begin{matrix} \sum_{P} {k_{P}}^{4} co s^{4} φ_{P} & \sum_{P} {k_{P}}^{4} \sin φ_{P} co s^{3} φ_{P} & \sum_{P} {k_{P}}^{4} si n^{2} φ_{P} co s^{2} φ_{P} \\ \sum_{P} {k_{P}}^{4} \sin φ_{P} co s^{3} φ_{P} & \sum_{P} {k_{P}}^{4} si n^{2} φ_{P} co s^{2} φ_{P} & \sum_{P} {k_{P}}^{4} si n^{3} φ_{P} \cos φ_{P} \\ \sum_{P} {k_{P}}^{4} si n^{2} φ_{P} co s^{2} φ_{P} & \sum_{P} {k_{P}}^{4} si n^{3} φ_{P} \cos φ_{P} & \sum_{P} {k_{P}}^{4} si n^{4} φ_{P} \end{matrix}) \times (\begin{array}{l} C \\ 2 D \\ F \end{array}) = \\ (\begin{array}{l} \sum_{P} {k_{P}}^{2} co s^{2} φ_{P} \\ \sum_{P} {k_{P}}^{2} \sin φ_{P} \cos φ_{P} \\ \sum_{P} {k_{P}}^{2} si n^{2} φ_{P} \end{array}) \end{array}

式中：k_P为样点到原点的距离;φ_P为样点原点连线与x轴的夹角。根据公式（4）就可以计算出椭圆方程中的基本参数,从而就可以得出椭圆的长短半径和倾角,公式如下^[22]⁸⁰：

（5）

\{\begin{array}{l} θ = 1 / 2 \arctan [2 D / (C - F)] \\ K_{A} = (C co s^{2} θ + 2 D \sin θ \cos θ + F si n^{2} {θ)}^{- \frac{1}{2}} \\ K_{B} = (C si n^{2} θ - 2 D \sin θ \cos θ + F co s^{2} {θ)}^{- \frac{1}{2}} \end{array}

式中：

θ

为该椭圆的倾角;K_A和K_B为该椭圆长短轴上的半径（长短半径视具体情况而定）。之后,就可以得到该椭圆的尺度,记为Λ,有φ(E_λ)=Λ=

\sqrt[]{K_{A} K_{B}}

,以及该椭圆上的各向异性比,记为K_B/K_A。

下面以区域1为例看一下尺度等参数的估计过程,详细步骤如下：首先将该区域上样本区间按照等差数列分为若干级,并对每一级样本点按照公式（4）计算C、D和F三个参数,然后根据公式（5）计算每一级椭圆的倾角θ、长短半径K_A和K_B、尺度Λ和各向异性比K_B/K_A,拟合结果见图8。另外,通过各向异性比K_B/K_A与尺度之间的关系可以确定圆尺度,从图9可以看出,与K_B/K_A=1相交的点代表其长短半径相等,因此对应的尺度值就是圆尺度,可以得出六个区域都存在圆尺度或者近圆尺度。为了精确得出圆尺度数值,这里采用二阶多项式对各向异性比随尺度的变化曲线进行拟合,如区域1圆尺度值约为410.4 m,其他区域的圆尺度估计结果参见表2。计算出圆尺度以后,通过圆尺度与每一级尺度Λ的比就可以得出各级尺度比λ。

图8

新窗口打开|下载原图ZIP|生成PPT
图8区域1样本椭圆拟合结果

Fig. 8A fitting to the semi-variogram of region 1by using the rotating ellipse method

图9

新窗口打开|下载原图ZIP|生成PPT
图9六个区域各向异性比随尺度的变化

Fig. 9Anisotropic ratios corresponding to the scales of DEM data of the six regions

Tab. 2
表2
表2六个区域GSI模型相关参数估计结果
Tab. 2The estimates of the related parameters in the GSI model for the six regions

研究数据	最小基台(m²)	最大变程(m)	圆尺度(m)	c	e	f	R_c²	R_e²	R_f²
区域1	1556	678.61	410.37	-0.1972	-0.1696	0.1924	0.9987	0.9999	0.9968
区域2	6234	1093.35	572.31	0.1011	-0.4123	-0.2226	0.9997	1.0000	0.9996
区域2	6234	1093.35	923.18	0.2017	-0.8014	-0.4199	0.9991	0.9997	0.9990
区域3	13678	2122.37	232.70	0.1172	-0.2011	-0.2056	0.9956	0.9987	0.9885
区域4	4197	1261.67	199.14	-0.2703	-0.1836	0.1915	0.9907	0.9986	0.9927
区域5	79592	4213.85	29.48	-0.0029	-0.3980	-0.1945	0.9995	0.9993	1.0000
区域6	15001	4423.97	49.58	-0.1761	-0.4026	-0.0979	1.0000	0.9991	0.9981

新窗口打开|下载CSV

3.2.2 c、e和f估计虽然找到一组椭圆简化了半变异函数,但是所需参数仍太多,每一个椭圆方程都需要专门的参数去描述,显然这不是最终目的。经过第一次建模以后,还需要对问题域进行第二次建模,即把这些椭圆用少数的参数进行关联,纳入到统一的框架中,下面看一下具体的建模过程。所谓第二次建模本质上是对尺度转换算子中其余参数的估计,除了尺度参数以外,待估参数还剩下c、e和f,这三个参数是关联整个椭圆系统的关键参数,控制着椭圆系统的旋转和拉伸。通过分析可以得出,上述问题本质上是一个在三维问题域中寻找最优解的问题,只要知道三个参数的大致变化范围,然后设置一定的步长,通过遍历整个定义域就可以得到最优解。为了获得最优解,首先需要找到目标函数的解析形式,这里仍然采用上述旋转椭圆法的思想进行获取,基本原理就是采用GSI模型对样本椭圆进行拟合,本次拟合后得到的椭圆就是GSI椭圆。具体获取过程如下：如图10所示,假设在圆尺度上有一点P,其坐标为(x_P, y_P),经过坐标转换λ^-G作用以后,变为GSI椭圆上的P'点,有

[\begin{array}{l} x_{P'} \\ y_{P'} \end{array}] = λ^{- G} [\begin{array}{l} x_{P} \\ y_{P} \end{array}]

。以原点O为起点,作过P'的射线,有

OP' = \sqrt[]{x_{P'}^{2} + y_{P'}^{2}}

,OP'与x轴之间的夹角为φ_P_',另有该射线与样本椭圆交于点P'',有OP''=(Ccos²φ_P_'+2Dsinφ_P_'cosφ_P'+Fsin²φ_P')^-1/2,因此,OP'-OP''就构成这个样本点的拟合误差,如果参与计算的样点覆盖所有椭圆,假设样点数目是n,那么就能得出本问题的目标函数,即

E_{GSI}^{2} = \frac{1}{n} \overset{n}{\sum_{P = 1}} (OP'^{2} - OP''^{2})^{2}

。为了加速收敛和得到较高精度的估计,这里采用两步搜索算法进行最优值的获取,即首先采用大步长得到初步估计,其次在初步估计附近再次采用小步长计算,最后将通过多项式对误差参数曲线拟合得到最优值。根据Pflug和Lewis等的研究结果,三个参数c、e和f的变化一般介于-1~1、-1.5~1.5以及-1~1之间^[22,26]^43,65。

以区域1为例,下面看一下c、e和f的估计过程,首先对圆尺度圆周进行均匀分割以获取参与计算的样点,这里将采样间隔设为π/20,并利用极坐标计算出每一个点在圆尺度上的坐标值;其次,按照公式

[\begin{array}{l} x_{P'} \\ y_{P'} \end{array}] = λ^{- G} [\begin{array}{l} x_{P} \\ y_{P} \end{array}]

,将按照步长0.1序列化的c、e和f值带入到坐标转换算子λ^-G中,计算样本点经过尺度转换后在各个尺度上对应点的坐标值,并根据转换后坐标与x轴之间的夹角寻找同方向上样本椭圆对应点的坐标;最后,计算目标函数,通过对比找出最小误差对应的参数值,经计算,初步估计结果为c=-0.2、e=-0.2和f=0.2。在第二轮搜索中,步长统一调整为0.01,优化过程按照顺序逐次进行,具体实施如下：首先优化c,保持e、f第一轮估计值不变,将c的区间设置为-0.3~-0.1,根据上述模型计算目标函数值,计算完以后绘制函数值关于c的二阶多项式拟合图（图11a）,然后通过一阶求导计算最优值;其次,将得到的c的最优值带入模型,保持f上轮估计不变,将e的搜索范围设置为-0.3~-0.1,与c的优化步骤相同,从而得到e的估计值（图11b）;最后将c、e的最优值代入到模型中,f区间设置为0.1~0.3,优化过程与其他两个参数类似（图11c）。通过两轮搜索以后,三个参数的估计值分别是c=-0.1972、e=-0.1696和f=0.1924。六个区域的估计结果汇总于表2和图12中（最小基台和最大变程计算结果,请参考3.1中确定工作范围的方法）,从图中可以看出,除了区域2采用了两次GSI建模（分割点约为1.5 km,见图9）,对于大多数区域仅仅通过一次GSI建模就能较好地刻画地形数据半变异函数中的各向异性特征。如果当前步长得到的R²较小,可以通过10的n次幂缩小步长重复上述过程,直到决定系数取值满意为止,例如在区域5和区域6的计算过程中,n=1的情况下,上述逐步搜索过程重复了三次。

图10

新窗口打开|下载原图ZIP|生成PPT
图10GSI模型对样本椭圆拟合

Fig. 10A diagram showing the fitting of the GSI model to the previous fitted ellipse

图11

新窗口打开|下载原图ZIP|生成PPT
图11基于二阶多项式对区域1半变异函数c、e和f的估计

Fig. 11The estimates of c, e and f of region 1 based on second-order polynomials

图12

新窗口打开|下载原图ZIP|生成PPT
图12六个区域样本椭圆和GSI椭圆叠合显示

Fig. 12The overlay show of the fitted ellipses to the estimated semi-variogram and the ellipses generated by the GSI model for the six regions

3.3 顾及各向异性的理论半变异函数模型拟合

首先来看处理几何各向异性的传统做法,一般计算过程如下,根据绘制的二维半变异函数图确定变程椭圆,通过椭圆方程计算该椭圆的长短半径以及各向异性比K_B/K_A,同时获取该椭圆的倾角θ,然后根据公式（6）对所有已知点坐标进行转换,将转换后的坐标值代入到加权多项式回归拟合模型中进行估计：

（6）

[\begin{array}{l} x_{u}^{'} \\ x_{v}^{'} \end{array}] = [\begin{matrix} \cos θ & - \sin θ \\ \sin θ & \cos θ \end{matrix}] [\begin{matrix} 1 & 0 \\ 0 & \frac{K_{B}}{K_{A}} \end{matrix}] [\begin{matrix} \cos θ & \sin θ \\ - \sin θ & \cos θ \end{matrix}] [\begin{array}{l} x_{u} \\ x_{v} \end{array}]

式中：x_u、x_v是转换前的已知点坐标;x'_u、x'_v是转换后已知点的坐标值。这里以球状模型和区域4为例,按照上述步骤得到的估计结果为,块金常数C₀=143.1752 m²,变程a=5932.374 m和拱高C=8978.4 m²。获取到理论模型估计以后,就能计算出二维平面中任何一点的各向异性半变异函数值。具体操作步骤如下：首先生成一个301×301像元的空矩阵（这里将区域4的工作范围限制在301×301像元内）,对每个单元格的坐标值按照公式（6）进行坐标转换,将转换后的坐标代入到估计模型中就能得到二维半变异函数图,结果见图13a。另外,区域4的各向同性模型估计参数分别为,C₀=988.2847m²,a=2802.525 m,C=6777.4 m²,半变异函数模拟结果见图13b。

图13

新窗口打开|下载原图ZIP|生成PPT
图13区域4三种方法的球状半变异函数模型估计

Fig. 13The estimates of the spherical model by using the three methods including traditional coordinate transformation, isotropic method and the GSI coordinate transformation for region 4

下面再看一下基于线性GSI的坐标转换处理过程,处理过程跟上面传统处理方法类似,首先需将已知点坐标进行各向同性化,向同性化的公式为：

（7）

[\begin{array}{l} x_{u, λ}^{'} \\ x_{v, λ}^{'} \end{array}] = λ^{G} [\begin{array}{l} x_{u, λ} \\ x_{v, λ} \end{array}] / λ

式中：x_u,λ、x_v,λ是在λ尺度上转换前的已知点坐标;x'_u,λ、x'_v,λ是在λ尺度上转换后已知点的坐标值,坐标转换后,同样采用加权多项式回归拟合模型进行估计,可得估计结果为C₀= 115.6784 m²,a=2169.795 m,C=7654.1 m²。在生成二维半变异函数图的过程中,难点在于每个单元格尺度比λ的估计,想直接得到λ,需要求解每个像元关于λ包含三个方程的方程组,由于存在矩阵函数,该方程组中含有超越方程,直接求解会比较困难。为了避免直接求解超越方程,这里采用正向逐步搜索法进行估计,具体实施如下：首先按照固定间隔生成一个λ序列,将λ序列代入到尺度转换模型中生成一组等值线去覆盖事先准备的空矩阵,然后从原点做一条过待填像元的射线,分别交该像元最近邻内外两条等值线于两点,比较该像元至这两点之间的距离,如果全部大于1个像元,就计算这两条等值线对应λ值的平均值,并生成一条位于两者中间的等值线,将该像元到原点的距离与新等值线继续比较,如此反复,直到离最近的一条等值线的距离小于1个像元为止,那么最近等值线λ值就是该像元的估计值,λ估计结果见图14。估计出每个像元的λ值后,那么就得到了每个像元的各向同性坐标,将其代入到含有估计值的理论函数模型中,可以得到二维半变异函数分布图（图13c）。以原点为中心,以正方形为外部边界,不断增加正方形的边长,统计位于该正方形边界上像元估计误差的均值,这样就绘制出一条平均估计偏差关于正方形边长的曲线（图15）,从结果上可以看出,几乎在整个变程区间内,基于线性GSI坐标转换法比原坐标转换法和各向同性方法在精度方面有了明显提升,尤其是在中距离上提升了近一倍以上。

图14

新窗口打开|下载原图ZIP|生成PPT
图14区域4尺度比 λ 矩阵

Fig. 14A matrix of the scale ratio λ for region 4

图15

新窗口打开|下载原图ZIP|生成PPT
图15区域4三种方法半变异函数估计误差随尺度的变化

Fig. 15Variations of the estimation biases on the spherical model with scales by using the three methods for region 4

3.4 各向异性半变异函数增强估计后的插值精度影响分析

下面以六个区域为例,对使用GSI模型增强估计后的半变异函数进行克里金插值的精度影响进行分析,并同时与各向同性和传统各向异性估计作对比。这里采用两种不同样本容量的数据集进行验证,一为数目较大随机采样的样本集,用于半变异函数估计（目的是确保半变异函数短距离处的估计精度及降低计算成本）,样本容量设为10000;另一为分布规则采样间隔固定的稀疏样本（目的是确保估计的稳定性和平滑性）,用于克里金插值,样本容量设为2500,样本数量为100。之所以这样做,通过试验发现,对于拥有复杂结构的半变异函数,由于受影响的因素增多,其克里金估计的不稳定性和不确定性显著增强。下面以区域4为例给出一次GSI增强型克里金估计的完整过程,也可参考图16：

图16

新窗口打开|下载原图ZIP|生成PPT
图16GSI坐标转换克里金插值流程

Fig. 16The kriging interpolation process by using the GSI coordinate transformation

（1）准备数据,通过生成随机数的方法从区域4中产生一个大样本,并且按照10×10像元间隔获得一个均匀采样的小样本。

（2）按照文中3.1的方法,对大样本去趋势后计算二维半变异函数。

（3）方向半变异函数提取,以π/16作为间隔,以某方向上相邻的10个单元格为一组进行合并,在0~π之间提取16个方向半变异函数曲线（图17）。

图17

新窗口打开|下载原图ZIP|生成PPT
图17区域4不同方向上的半变异函数曲线

Fig. 17The semi-variogram curves along the 16 directions ranging from 0 to π for region 4

（4）按照一定的间隔获取不同尺度上的等值线样点集,根据对二维半变异函数的观察,这里将间隔设置为200 m²,最大半变异函数值设置为3000 m²（与3.3章节中的估计结果相差较大,目的是消除长距离估计的不稳定性）,与等值线对应的样点坐标采用分段线性插值方法确定（图18）。

图18

新窗口打开|下载原图ZIP|生成PPT
图18区域4等值线样点和样本椭圆估计

Fig. 18Sample points on the contours at the scales and the corresponding estimates for the experimental semi-variogram for region 4

（5）与旋转椭圆法相比,加权旋转椭圆法以距离平方的倒数作为权重构建目标函数,这样做更好地保证了短距离半变异函数值的精度（在克里金插值中,短距离权重占比更高）。根据加权旋转椭圆法对上述等值线样点进行样本椭圆拟合,并绘制各向异性比随尺度变化的曲线,以及计算圆尺度。经计算,圆尺度约为65.556 m。

（6）对（5）中的结果进行GSI拟合,计算参数c、e和f,经估计,c=0.0101、e=0.2901和f=0.1943。然后根据此3参数的估计结果进行每个单元格尺度比λ的生成,该参数是计算GSI球状模板矩阵每个单元格半变异函数值的必备条件,这里采取逐步搜索法进行计算,具体计算方法参考3.3章节。

（7）将GSI转换机制用于球状模型参数的估计中,得出球状模型估计参数为：C₀=-11.4688 m²,a=1060.662 m,C=2941.1 m²;而传统坐标转换法估计结果为：C₀=146.2938 m²,a=2740.578 m,C=3172.7 m²;各向同性模型的估计结果为：C₀=653.8444 m²,a=1729.371 m,C=2693 m²。

（8）将参数c、e、f、C₀_、a、C和λ矩阵等代入半变异函数公式中生成GSI半变异函数模板矩阵。

（9）对小样本数据进行边缘化处理,然后滑动模板矩阵分别覆盖离散矩阵中的每个栅格,获取落入邻域范围内的样点及权重并进行插值（图19a~图19c）。

图19

新窗口打开|下载原图ZIP|生成PPT
图19区域1和区域4三种方法克里金插值

Fig. 19The kriging interpolation by using the three methods including traditional coordinate transformation, isotropic method and the GSI coordinate transformation for regions 1 and 4

按照上述流程对六个区域的所有样本依次进行插值,然后将得到的所有估计与真值比较,并计算均方根误差RMSE（Root Mean Square Error）的均值和标准差,所有结果汇总于表3中。从表3可以看出,对于区域4,与传统坐标转换方法相比较,经过GSI估计后的插值精度有了较大程度的改善,这与半变异函数的估计优势是相吻合的,同样的结论也在其余五个区域上得到了验证,证明了GSI模型在拟合半变异函数结构上的优势,因此在使用传统坐标转换的时候,一定不能忽视其内部结构的影响。然而,除了区域4以外,包括GSI和传统方法的各向异性插值在各向同性方法面前并没有显现出特有方法上的优势,甚至在有些区域计算结果上呈现出较大的劣势。究其原因,首先,从半变异函数的定义来看,它是一个利用空间其余位置的样点值来弥补重复观测得到的全局统计量,某点的函数值都是所有观测值综合后的结果,这种统计平均不足以准确描述某点的局部相关性,可想而知,每一点的插值都受到了其余点位的影响,尤其是受到数据主体结构的影响;其次,克里金插值还与半变异函数的估计精度、结构特征、滑动窗口大小以及样本观测值的数量和分布等因素也有着密切的关系。另外,对于各向异性插值,其影响机制则更为复杂,最大的不同在于非欧几何度量彻底改变了协方差插值矩阵的结构,降低了平滑性,尤其对于簇聚程度高的样点,极易导致异常和极端权系数的出现,给插值结果带来极大的不稳定性和不确定性,这也是上述验证采用稀疏规则分布样本的原因。例如,位于某个待插点滑动窗口边缘的样点,协方差矩阵中该点与其余点构成的协方差向量容易出现一值独大的情形（由边界对于样点组合的限制引起）,在非欧距离中体现的则更为明显,从而导致插值波动幅度过大,形成条带效应。虽然从均方根误差这个指标上来看GSI插值的稳定度较低,但是发现GSI模型在几乎所有的图像上都表现出了较好的边缘信息恢复和细节描述能力,这也是各向同性方法所不具备的。例如,以区域1为例（图19d~图19f）,虽然三种方法计算出的均方根误差相差不大,但是很明显能够看出GSI模型的边缘信息恢复是最好的。综上所述,在不考虑上述因素影响的前提下,就均方根误差而言各向同性插值体现出了较好的稳定性和精度,主要是由其受到的限制条件较少和简单,而各向异性插值表现出较强的不稳定性和不确定性,原因是其受到的影响因素更为复杂,因此,若想使用各向异性克里金插值获得理想的结果,除了选择合适的各向异性模型及准确的半变异函数理论模型估计以外,自然也不能忽视其他因素的影响如新的距离度量导致的异常和极端权值或者来自数据自身斑块结构特征的影响。

Tab. 3
表3
表3六区域传统各向异性、各向同性和GSI克里金估计均方根误差的均值和标准差
Tab. 3The means and standard deviations of RMSE concerned with kriging estimates versus the real values by using the traditional anisotropic, isotropic and GSI methods for the six regions

	均值（m）			标准差（m）
	传统各向异性	各向同性	GSI	传统各向异性	各向同性	GSI
区域1	21.6778	22.1308	21.1273	0.1592	0.1546	0.1852
区域2	30.5789	28.7977	28.1073	0.2777	0.2901	0.3217
区域3	42.3000	28.5649	35.9589	0.2304	0.3804	0.2535
区域4	59.4448	28.4447	23.2694	0.2821	0.1654	0.3421
区域5	104.3453	42.1847	56.1315	1.0085	0.6435	0.5204
区域6	31.9274	20.0796	23.7915	0.2454	0.2706	0.1671

新窗口打开|下载CSV

3.5 讨论

作为一个理论模型,旨在采用少数的已知点对较少的参数进行估计,从而推断出覆盖在定义域上的全部信息。当然,在保证精度的前提下,需要的已知点和待估计参数越少越好,模型和参数越简单、越容易估计越好。与传统观念矛盾的是,GSI模型原理和估计过程却较为复杂,参数和需要的已知点也较多,但是随着遥感技术的快速发展和广泛应用,借助遥感数据中的结构或者属性信息来辅助插值变得越来越普遍^[27,28],因此样本点稀少带来的影响将逐渐变小。纵然是连续分布的面状遥感数据,一般也不会直接采取其半变异函数进行插值,也需要通过理论模型估计或者正定化处理后才能使用^[27]。另外,传统的处理方式仅仅采用了最外层椭圆的各向异性信息得出了各向同性转换机制,显然其并未考虑内部结构信息,例如圆形等值线在这种转换机制下就变成了各向异性的（图1a和图1b）,因此,它并不是一种真正的各向同性转换机制,而GSI模型不仅考虑了各级椭圆的方向信息,也对整体结构进行了控制（图1c）,虽然并不能完全反映全部信息,但是比原有的转换机制有了较大程度的提升。从尺度方面来看,传统方式转换后的椭圆是以长轴为半径的圆,而GSI则从面积方面定义尺度来进行约束,从而保证了对真实半变异函数邻域信息的合理覆盖。为了避免非平稳性带来的影响,在估计之前对原数据进行了一阶多项式去趋势处理,显然这种操作会对使用半变异函数分析原始图像的尺度特征产生困难。如果非平稳性很弱或者可以忽略,那么就能直接采用原数据进行估计。从三种方法的半变异函数估计对克里金插值精度的影响分析来看,顾及内部结构的GSI建模其插值效果较传统各向异性的提升是明显的,然而,相对于各向同性方法来说,各向异性插值并没有体现出多少方法上的优势,对数据本身似乎也有着某种特定的要求和限制,尽管本次所有的验证区域都是各向异性的。可见选择合适的方法进行克里金插值至关重要,如果不去分析总结这些方法的影响因素和适用条件,一味盲目地使用,效果只会适得其反。当然上述结论是在采用稀疏规则样本插值和半变异函数准确估计的前提下得出的,更多一般的、普适性的结论还需要进一步分析研究。

4 结论

从六个研究区域二维半变异函数估计的空间分布上来看,变异特征随着尺度的变化而变化,呈现出比较复杂的状况,传统的坐标转换处理方法必然顾此失彼,这里需要一种多尺度的坐标转换模型。当然这种变化也全然不是无规可循的,例如在六个区域都发现了圆形或者近圆形等值线,正是这些特殊尺度构成了本次建模的起点。由于GSI模型的原理较为复杂,直接获取参数估计比较困难,这里将估计过程分为了两个阶段,重复使用了旋转椭圆法,即首先估计出尺度相关参数,然后在此基础上再次通过两步正向搜索作图法估计其余参数,从估计结果来看,本方法可以获得最优解。不同于传统坐标转换模型,GSI坐标转换模型是一种全局多尺度控制模型,它为各向异性半变异函数各向同性化提供了一种新的思路和方法,结果表明,它对半变异函数理论模型整体估计精度提升具有明显作用,对于某些区域克里金插值也表现出一定的优势,且显现出了良好的边缘信息恢复能力,为进一步丰富和完善地统计学中各向异性结构分析和应用的理论和方法提供了依据。

致谢：

真诚感谢二位匿名评审专家在论文评审中所付出的时间和精力,评审专家针对本文内容和结构提出了许多建设性的意见,使本文获益匪浅。

参考文献原文顺序
文献年度倒序
文中引用次数倒序
被引期刊影响因子

[1]

冯波, 陈明涛, 岳冬冬, 等. 基于两种插值算法的三维地质建模对比
吉林大学学报(地球科学版), 2019,49(4):1200-1208.

[本文引用: 1]

[ Feng

, Chen

Mingtao

, Yue

Dongdong

, et al. Comparison of 3D geological modeling based on two different interpolation methods
Journal of Jilin University (Earth Science Edition), 2019,49(4):1200-1208.]

[本文引用: 1]

[2]

吕海洋, 盛业华, 李佳, 等. 基于RASM的紧支撑径向基函数自适应并行地形插值方法
武汉大学学报(信息科学版), 2017,42(9):1316-1322.

[本文引用: 1]

[ Lv

Haiyang

, Sheng

Yehua

, Li

Jia

, et al. An adaptive parallel CSRBF terrain interpolation method based on RASM
Geomatics and Information Science of Wuhan University, 2017,42(9):1316-1322.]

[本文引用: 1]

[3]

聂磊, 舒红, 刘艳. 复杂地形地区月平均气温(混合)地理加权回归克里格插值
武汉大学学报(信息科学版), 2018,43(10):1553-1559.

[本文引用: 1]

[ Nie

Lei

, Shu

Hong

, Liu

Yan

. Interpolation of monthly average temperature by using (mixed) geographically weighted regression kriging in the complex terrain region
Geomatics and Information Science of Wuhan University, 2018,43(10):1553-1559.]

[本文引用: 1]

[4]

夏伟杰, 周建江, 姚楠. 各向异性的分形地形生成方法研究
中国图象图形学报, 2009,14(11):2356-2361.

DOI:10.11834/jig.20091126 URL [本文引用: 1]

传统的分形地形生成方法得到的地形是各向同性的，为了使生成的分形地形具有各向异性特征，提出了一种新的分形地形生成方法，该方法利用组合分形布朗曲面模型，将具有不同特征的两种分形布朗曲面相融合，使得生成的地形具有各向异性特征。对组合分形布朗曲面算法进行了仿真实验，生成了最终的分形地形。对仿真结果的分析表明，生成的分形地形的特征具有各向异性特征，和实际的自然地形特征相符合，从而证明了该算法的有效性。

[ Xia

Weijie

, Zhou

Jianjiang

, Yao

Nan

. Anisotropic fractal terrain generation
Journal of Image and Graphics, 2009,14(11):2356-2361.]

[本文引用: 1]

[5]

Lovejoy

, Schertzer

. Scaling and multifractal fields in the solid earth and topography
Nonlinear Processes in Geophysics, 2007,14(4):465-502.

DOI:10.5194/npg-14-465-2007 URL [本文引用: 1]

[6]

Merwade V

, Maidment D

, Goff J

. Anisotropic considerations while interpolating river channel bathymetry
Journal of Hydrology, 2006,331(3-4):731-741.

DOI:10.1016/j.jhydrol.2006.06.018 URL [本文引用: 1]

[7]

Gyasi-Agyei

. Assessment of radar-based locally varying anisotropy on daily rainfall interpolation
Hydrological Sciences Journal, 2016,61(10):1890-1902.

[本文引用: 1]

[8]

Budrikaite,

, Ducinskas,

. Modelling of geometric anisotropic spatial variation
In: Mathematical Modeling and Analysis, Proceedings of the 10th International Conference MMA2005&CMAM2. Trakai, 2005: 361-366.

[本文引用: 1]

[9]

Zimmerman D

. Another look at anisotropy in geostatistics
Mathematical Geology, 1993,25(4):453-470.

DOI:10.1007/BF00894779 URL [本文引用: 1]

[10]

Isaaks E

, Srivastava R

. An Introduction to Applied Geostatistics. New York: Oxford University Press, 1989: 149-154.

[本文引用: 1]

[11]

Goovaerts

. Geostatistics for Natural Resources Evaluation. New York: Oxford University Press, 1997: 480-485.

[本文引用: 1]

[12]

Eriksson

, Siska P

. Understanding anisotropy computations
Mathematical Geology, 2000,32(6):683-700.

DOI:10.1023/A:1007590322263 URL [本文引用: 1]

Most descriptions of anisotropy make reference to reduced distances and conversion of anisotropic models to isotropic counterparts and equations are presented for a certain class of range-anisotropic models. Many descriptions state that sill anisotropy is modelled using a range-anisotropic structure having a very elongated ellipse. The presentations typically have few or no intervening steps. Students and applied researchers rarely follow these presentations and subsequently regard the programs that compute anisotropic variograms as black-boxes, the contents of which are too complex to try to understand. We provide the geometry necessary to clarify those computations. In so doing, we provide a general way to model any type of anisotropy (range, sill, power, slope, nugget) on an ellipse. We note cases in the literature in which the printed descriptions of anisotropy on an ellipse do not match the stated or coded models. An example is provided in which both range- and sill-anisotropic structures are fitted to the experimental variogram values from an elevation data set using the provided equations and weighted nonlinear regression. The original variogram values are plotted with the fitted surfaces to view the fit and anisotropic structure in many directions at once.]]>

[13]

Cressie N

. Statistics for Spatial Data
New York: John Wiley, 1993: 83-90.

[本文引用: 1]

[14]

Ecker M

, Gelfand A

. Bayesian modeling and inference for geometrically anisotropic spatial data
Mathematical Geology, 1999,31(1):67-83.

[本文引用: 1]

[15]

Shen

, Gelfand A

. Exploring geometric anisotropy for point-referenced spatial data
Spatial Statistics, 2019,32(4):100370.

DOI:10.1016/j.spasta.2019.100370 URL [本文引用: 1]

[16]

Hristopulos D

. New anisotropic covariance models and estimation of anisotropic parameters based on the covariance tensor identity
Stochastic Environmental Research and Risk Assessment, 2002,16(1):43-62.

DOI:10.1007/s00477-001-0084-y URL [本文引用: 1]

Many heterogeneous media and environmental processes are statistically anisotropic. In this paper we focus on range anisotropy, that is, stochastic processes with variograms that have direction dependent correlation lengths and direction independent sill. We distinguish between two classes of anisotropic covariance models: Class (A) models are reducible to isotropic after rotation and rescaling operations. Class (B) models can be separated into a product of one-dimensional functions oriented along the principal axes. We propose a new Class (A) model with multiscale properties that has applications in subsurface hydrology. We also present a family of Class (B) models based on non-Euclidean distance metrics that are generated by superellipsoidal functions. Next, we propose a new method for determining the orientation of the principal axes and the degree of anisotropy, i.e., the ratio(s) of the correlation lengths. This information reduces the degrees of freedom of anisotropic variograms and thus simplifies the estimation procedure. In particular, Class (A) models are reduced to isotropic and Class (B) models to one-dimensional functions. Our method is based on an explicit relation between the second-rank slope tensor (SRST), which can be estimated from the data, and the covariance tensor. The procedure is conceptually simple and numerically efficient. It is more accurate for regular (on-grid) data distributions, but it can also be used for sparse (off-grid) spatial distributions. In the case of non-differentiable random fields the method can be extended using generalized derivatives. We illustrate its implementation with numerical simulations.]]>

[17]

Cressie

, Huang H

. Classes of nonseparable, spatio-temporal stationary covariance functions
Journal of the American Statistical Association, 1999,94(448):1330-1339.

DOI:10.1080/01621459.1999.10473885 URL [本文引用: 1]

[18]

Allard

, Senoussi

, Porcu

. Anisotropy models for spatial data
Mathematical Geosciences, 2016,48(3):305-328.

DOI:10.1007/s11004-015-9594-x URL [本文引用: 1]

[19]

De Iaco

, Posa

, Cappello

, et al. Isotropy, symmetry, separability and strict positive definiteness for covariance functions: A critical review
Spatial Statistics, 2019,29(1):89-108.

DOI:10.1016/j.spasta.2018.09.003 URL [本文引用: 1]

[20]

Schertzer

, Lovejoy

. Generalized scale invariance in turbulent phenomena
Physicochemical Hydrodynamics, 1985,6:623-635.

[本文引用: 2]

[21]

Lovejoy

, Schertzer

. Generalized scale invariance in the atmosphere and fractal models of rain
Water Resources Research, 1985,21(8):1233-1250.

DOI:10.1029/WR021i008p01233 URL [本文引用: 2]

[22]

Pflug

. Generalized scale invariance, differential rotation and cloud texture
Montreal: Master Dissertation of McGill University, 1993.

[本文引用: 3]

[23]

Lovejoy

, Schertzer

, Pflug

. Generalized scale invariance and differential rotation in cloud radiances
Physica A: Statistical Mechanics and its Applications, 1992,185(1-4):121-128.

DOI:10.1016/0378-4371(92)90446-W URL [本文引用: 1]

[24]

Gagnon J

, Lovejoy

, Schertzer

. Multifractal earth topography
Nonlinear Processes in Geophysics, 2006,13(5):541-570.

DOI:10.5194/npg-13-541-2006 URL [本文引用: 1]

[25]

Schertzer

, Lovejoy

. Multifractals, generalized scale invariance and complexity in geophysics
International Journal of Bifurcation and Chaos, 2011,21(12):3417-3456.

DOI:10.1142/S0218127411030647 URL [本文引用: 1]

The complexity of geophysics has been extremely stimulating for developing concepts and techniques to analyze, understand and simulate it. This is particularly true for multifractals and Generalized Scale Invariance. We review the fundamentals, introduced with the help of pedagogical examples, then their abstract generalization is considered. This includes the characterization of multifractals, cascade models, their universality classes, extremes, as well as the necessity to broadly generalize the notion of scale to deal with anisotropy, which is rather ubiquitous in geophysics.

[26]

Lewis

. The scale invariant generator technique and scaling anisotropy in geophysics
Montreal: Master Dissertation of McGill University, 1993.

[27]

Velasco-Forero C

, Sempere-Torres

, Cassiraga E

, et al. A non-parametric automatic blending methodology to estimate rainfall fields from rain gauge and radar data
Advances in Water Resources, 2009,32(7):986-1002.

DOI:10.1016/j.advwatres.2008.10.004 URL [本文引用: 2]

AbstractQuantitative estimation of rainfall fields has been a crucial objective from early studies of the hydrological applications of weather radar. Previous studies have suggested that flow estimations are improved when radar and rain gauge data are combined to estimate input rainfall fields. This paper reports new research carried out in this field. Classical approaches for the selection and fitting of a theoretical correlogram (or semivariogram) model (needed to apply geostatistical estimators) are avoided in this study. Instead, a non-parametric technique based on FFT is used to obtain two-dimensional positive-definite correlograms directly from radar observations, dealing with both the natural anisotropy and the temporal variation of the spatial structure of the rainfall in the estimated fields. Because these correlation maps can be automatically obtained at each time step of a given rainfall event, this technique might easily be used in operational (real-time) applications. This paper describes the development of the non-parametric estimator exploiting the advantages of FFT for the automatic computation of correlograms and provides examples of its application on a case study using six rainfall events. This methodology is applied to three different alternatives to incorporate the radar information (as a secondary variable), and a comparison of performances is provided. In particular, their ability to reproduce in estimated rainfall fields (i) the rain gauge observations (in a cross-validation analysis) and (ii) the spatial patterns of radar fields are analyzed. Results seem to indicate that the methodology of kriging with external drift [KED], in combination with the technique of automatically computing 2-D spatial correlograms, provides merged rainfall fields with good agreement with rain gauges and with the most accurate approach to the spatial tendencies observed in the radar rainfall fields, when compared with other alternatives analyzed.]]>

[28]

Tandeo

, Autret

, Chapron

, et al. SST spatial anisotropic covariances from METOP-AVHRR data
Remote Sensing of Environment, 2014,141(4):144-148.

DOI:10.1016/j.rse.2013.10.024 URL [本文引用: 1]