CN113537329B - 一种逐位置快速估算各类地物概率分布的方法 - Google Patents

Abstract

本发明属于空间信息技术领域，具体涉及一种逐位置快速估算各类地物概率分布的方法。本发明经过步骤1.给定估算的所需参数；步骤2.根据目标区域大小，设定目标规则格网TG，并根据目标地物的空间结构特征选择和目标规则格网具有相同空间分辨率的训练图像TI；步骤3.将TI缩放到从1到G的不同粒度大小，步骤4.对每个TI_g中的网格单元p_g，TI中与其对应的像元组成一个集合；步骤5.将TG缩放到从1到G的不同粒度大小等12个步骤。本发明解决现有多点模拟方法难以根据少数样本数据在研究区域内快速逐位置计算各类地物概率分布的问题。

Description

一种逐位置快速估算各类地物概率分布的方法

技术领域

本发明属于空间信息技术领域，具体涉及一种逐位置快速估算各类地物概率分布的方法。

背景技术

模拟一直是研究地理现象的一种常用技术和方法。例如对在气象、地质等很多方面，很多研究都是关于如何从有限样本点模拟出研究对象在研究区域的实际分布。20世纪90年代Guardiano，Srivastava和Joumel等提出了使用多点统计学来对进行地理现象建模的这一理论方法，即多点模拟。这一方法可以弥补传统地质统计学变差函数的不足：仅考虑两点之间的统计关系。基于这一思想，专家们提出了很多相关算法并进行了验证。目前这些模拟算法已经被广泛的用于石油和矿物探测、超分辨率制图以及遥感图像分类等各个方面。尤其是在矿物储层模拟方面，通常只能通过有限的几十个样本去估计整个区域不同位置的石油分布情况。

目前，多点模拟方法可以分为三个类别。第一个类别是基于正则方程的方法，例如ENESIM，SNESIM以及IMPALA。这些方法通过训练图像来学习不同空间近邻条件下各类别的空间分布。这些方法主要针对的是类别数据。第二个类别是基于模式匹配的方法。这些方法从训练图像中直接学习结构模式，然后将这些模式直接粘贴到模拟格网中。代表性的方法有FILTERSIM，SIMPAT 和RPAFSIM。第三种类型的方法是基于抽样或者信息论的方法，例如DS算法。这类算法通常需要大量的运算才能获得比较好的模拟结果。

在实际应用中，通常单次的模拟结果是没有任何参考价值的，必须要通过多次模拟来获得研究区域不同位置各个类别出现的概率。通常需要成百上千次的模拟才能较为准确的估算不同位置各类别出现的概率。然而，如此多的模拟需要大量的计算时间。然而在很多情况下，时间效率尤为关键，因此需要新方法能够快速直接的推算出不同位置的各类别概率。

本发明提出了快速估估算不同类别地物复杂空间概率分布的方法。这一方法首先将训练图像和目标目标规则格网分别转换到从粗到细多个不同分辨率规则格网上。不同分辨率规则格网上的每个像元的值被设置为一个类别矢量。矢量的每个元素都是这一元素所对应类别在这一像元空间范围内出现的数目。

其次，该方法分别在不同分辨率规则格网上，根据训练图像，在给定的距离范围内，分别统计每个近邻位置出现不同类别组合时，中心像元出现各类别的数目。然后，在每个分辨率规则格网上，将这一统计结果存储下来。

最后，从最粗的规则格网开始，根据样本数据构建当前像元的条件数据，并以此为依据估计当前格网的每个像元上不同位置近邻条件下各个类别的出现的概率。然后，把各类别的边缘概率作为各类别先验概率，再融合所有位置近邻的推断结果，并得出当前像元各类别出现的概率。最后将这一概率作为下一层较细规则格网对应像元的先验概率，使用同样的方法推导出这一层规则格网的不同位置出现不同类别的概率。依次类推，得到最细规则格网上每个像元不同类别的概率分布。

发明内容

本发明的目的是面向对具有复杂几何特征的多类别地物，解决现有多点模拟方法难以根据少数样本数据在研究区域内快速逐位置计算各类地物概率分布的问题，提供一种逐位置快速估算各类地物概率分布的方法。

为解决上述技术问题，本发明采用的技术方案包含以下步骤：

步骤1.给定估算的所需参数，包括最大粒度G，近邻搜索半径τ，以及包含m个样本的样本数据集合H＝{h₁，h₂，…，h_m}，给出研究区域的所有地物类别，并将C个类别组成的集合记为{s₁，s₂，…，s_C}，G为一个正整数参数，τ为一个正实数参数；

步骤2.根据目标区域大小，设定目标规则格网TG，并根据目标地物的空间结构特征选择和目标规则格网具有相同空间分辨率的训练图像TI，其中空间分辨率为每个格网单元所表示的地面实际尺寸的大小；

步骤3.将TI缩放到从1到G的不同粒度大小，并分别记为TI₁，TI₂，…， TI_G；在将TI缩放到粒度1≤g≤G过程中，选定TI的某个顶点格网单元o为坐标原点，其坐标为(0，0，0)；对于每个TI中的像元p，将p在三个坐标方向上在TI 规则格网上和o相差的格网单元数记为p的坐标(p_x，p_y，p_z)；同时，选定TI_g中对应于o的顶点格网单元o_g为坐标原点，其坐标为(0，0，0)；对于每个TI_g中的像元 p_g，将p_g在三个坐标方向上在TI_g规则格网上同o_g相差的格网单元数记为p_g的坐标(p_gx，p_gy，p_gz)；

步骤4.对每个TI_g中的网格单元p_g，TI中与其对应的像元组成一个集合 F(p_g)＝{p∈TI|2^g-1×p_gx≤p_x＜2^g-1×(p_gx+1)Λ2^g-1×p_gy≤p_y＜ 2^g-1×(p_gy+1)Λ2^g-1×p_gz≤p_z＜2^g-1×(p_gz+1)}；当计算出F(p_g)后，将网格单元p_g对应的值设定为(s₁(p_g)，s₂(p_g)，…，s_C(p_g))，其中s_i(p_g)表示 F(p_g)中被标识为s_i类别的网格单元个数有多少；

步骤5.将TG缩放到从1到G的不同粒度大小，并分别记为TG₁，TG₂，…， TG_G；在将TG缩放到粒度1≤g≤G过程中，选定TG的某个顶点格网单元o′为坐标原点，其坐标为(0，0，0)；对于每个TG中的像元q，将q在三个坐标方向上在TG规则格网上和o′相差的格网单元数记为q的坐标(q_x，q_y，q_z)；同时，选定TG_g中对应于o′的顶点格网单元o′_g为坐标原点，其坐标为(0，0，0)；对于每个TG_g中的像元q_g，将q_g在三个坐标方向上在TG_g规则格网上同o′_g相差的格网单元数记为q_g的坐标(q_gx，q_gy，q_gz)；

步骤6.将样本数据集合H中的每个样本h指定到距离其最近的TG网格单元中；对每个TG中的网格单元q，样本数据中被指定到网格单元q中的样本集合记为H_q，q所对应的值为一个向量(s₁(q)，s₂(q)，…，s_C(q))，其中s_i(q)表示H_q中被标识为s_i类别的样本数有多少；

步骤7.对每个TG_g中的网格单元q_g，TG中与其对应的像元组成一个集合 F(q_g)＝{q∈TG|2^g-1×q_gx≤q_x＜2^g-1×(q_gx+1)Λ2^g-1×q_gy≤q_y＜ 2^g-1×(q_gy+1)Λ2^g-1×q_gz≤q_z＜2^g-1×(q_gz+1)}；当计算出F(q_g)后，将网格单元q_g对应的值设定为(s₁(q_g)，s₂(q_g)，…，s_C(q_g))，其中s_i(q_g)表示 F(q_g)中具有s_i类别标识的网格单元个数有多少；

步骤8.根据给定的近邻搜索半径τ，初始化近邻模板T；T由距离坐标原点不超过τ个网格单元的整数位置矢量构成；对于T中的每个元素t，初始化类别矩阵

为0；

每一行表示一个类别组合；例如有两个类别s₁和 s₂，则

有三行，分别表示t中只出现类别s₁，只出现类别s₂以及同时出现类别s₁和s₂；

第i列代表各种组合下s_i类别在中心像元出现的个数；

步骤9.设定g＝G，并将TI中各类别的比例作为各类别在训练图像中的概率质量函数，记为{P_ti(s₁)，P_ti(s₂)，…，P_ti(s_C)}，同时设置TG_g中每个网格单元q_g各类别的先验概率质量函数为

步骤10.初始化所有类别矩阵

为0，t∈T；对于TI_g上的每个网格单元p_g，在TI_g中寻找同p_g距离不超过τ个网格单元集合P_g；对于P_g中的每个网格单元δ，其坐标为(δ_x，δ_y，δ_z)，在近邻模板T中找到δ对应的元素t^δ＝(δ_x- p_gx，δ_y-p_gy，δ_z-p_gz)；在

中找到网格单元δ所有出现类别所对应的类别组合所在行i，并将p_g所对应的值(s₁(p_g)，s₂(p_g)，…，s_C(p_g))累加到

的i行；

步骤11.对于TG_g中的每个网格单元q_g计算各类别的概率质量函数，具体方法如下：

(1)将(P(s₁(q_g))，P(s₂(q_g))，…，P(s_C(q_g)))初始化为0；

(2)在TG_g中寻找同q_g距离不超过τ个网格单元集合TQ_g；

(3)对于TQ_g中的每个网格单元γ，其坐标为(γ_x，γ_y，γ_z)；如果 (s₁(γ)，s₂(γ)，…，s_C(γ))为0向量，则从TQ_g中去除γ；如果(s₁(γ)，s₂(γ)，…，s_C(γ)) 不为0向量，则在近邻模板T中找到γ对应的元素t^γ＝(γ_x-p_gx，γ_y-p_gy，γ_z- p_gz)；在

中找到网格单元γ所有出现类别所对应的类别组合所在行j，然后设置

为

第j行的值。

(4)对每个类别s_i使用公式

融合TQ_g中不同网格单元所推导出的类别概率，其中

计算出所有的P(s_i(q_g))后，使用

对P(s_i(q_g))进行归一化得到每个类别的概率。

(5)如果g≠1，在TG_g-1中找到和q_g对应的网格单元集合Q_g-1＝ {q_g-1∈TG_g-1|2×q_gx≤q_(g-1)x＜2×(q_gx+1)Λ2×q_gy≤q_(g-1)y＜ 2×(q_gy+1)Λ2×q_gz≤q_(g-1)z＜2×(q_gz+1)}，并将q_g-1各类别的先验概率质量函数

设置为 {P(s₁(q_g))，P(s₂(q_g))，…，P(s_C(q_g))}，设定g＝g-1，跳转到步骤10；

步骤12.将TG₁，TG₂，…，TG_G上每个网格单元的各类别概率质量函数结果输出。

附图说明

图1为本发明的流程图；

图2为训练图像；

图3为粒度为2的训练图像(TI₂)构建方法示意图；

图4为粒度为2的目标规则格网构建方法示意图；

图5为粒度为2的目标规则格网在每个网格单元上的各类别概率质量估计结果。

具体实施方式

本实施实例中使用了一个多类别训练图像进行说明，为了方便演示起见，采用了二维情况进行说明，但本方法同样适用于三维情况。同时在二维情况估算过程中，设置所有像元的z轴坐标为0，并不影响其计算结果。结合图 1，具体包括以下步骤：

步骤1.给定估算所需的最大粒度G＝2，近邻搜索半径τ＝20，研究区域的地物类别数目C＝3，记这些类别分别为{s₁，s₂，s₃}。研究区域包含m＝20个样本的样本数据集合，样本数据如下所示：

为了表示和计算方便，样本数据的最后一列使用了s₁，s₂，s₃的下标来表示这三个类别。

步骤2.根据目标区域大小，我们设定了目标规则格网TG为大小100×100 像元，并根据目标地物的空间结构特征选择和目标规则格网相同空间分辨率的训练图像TI，其大小为250×250像元。训练图像和目标规则格网每个像元所表示的地面实际尺寸的大小相同。训练图像如图2所示：

步骤3.将TI缩放到从1到2的不同粒度大小，并分别记为TI₁和TI₂。将TI 缩放到粒度1≤g≤2过程中，选定TI的某个顶点格网单元o为坐标原点，其坐标为(0，0，0)；对于每个TI中的像元p，将p在三个坐标方向上在TI规则格网上和θ相差的格网单元数记为p的坐标(p_x，p_y，p_z)；同时，选定TI_g中对应于o的顶点格网单元o_g为坐标原点，其坐标为(0，0，0)；对于每个TI_g中的像元p_g，将p_g在三个坐标方向上在TI_g规则格网上和o_g相差的格网单元数记为p_g的坐标 (p_gx，p_gy，p_gz)；

步骤4.对每个TI_g中的网格单元p_g，TI中与其对应的像元组成一个集合 F(p_g)＝{p∈TI|2^g-1×p_gx≤p_x＜2^g-1×(p_gx+1)Λ2^g-1×P_gy≤p_y＜ 2^g-1×(p_gy+1)Λ2^g-1×p_gz≤p_z＜2^g-1×(p_gz+1)}；当计算出F(p_g)后，将网格单元p_g对应的值设定为(s₁(p_g)，s₂(p_g)，s₃(p_g))，其中s_i(p_g)表示F(p_g)中被标识为s_i类别的网格单元个数有多少，如图3例子所示。

图3左面子图是一个简单的训练图像，包含B(黑色)，W(白色)两个类别，每一个黑色或白色方框代表一个像元。上图右面子图是左面子图缩放到粒度为2的训练图像，每个方框代表该粒度下的一个像元。都以最左上角顶点作为原点，那么两个子图中左上角顶点的坐标都为(0，0)，并且右面子图中左上角顶点对应的F((0，0))集合为{(0，0)，(0，1)，(1，0)，(1，1)}这四个左面子图中的像元。因此粒度为2的训练图像中左上角顶点的所对应的值为矢量(1，3)。这一矢量的含义是出现1个黑色类别和3个白色类别。

步骤5.将TG缩放到从1到G的不同粒度大小，并分别记为TG₁，TG₂，…， TG_G；在将TG缩放到粒度1≤g≤G过程中，选定TG的某个顶点格网单元o′为坐标原点，其坐标为(0，0，0)；对于每个TG中的像元q，将q在三个坐标方向上在TG规则格网上和o′相差的格网单元数记为q的坐标(q_x，q_y，q_z)；同时，选定TG_g中对应于o′的顶点格网单元o′_g为坐标原点，其坐标为(0，0，0)；对于每个TG_g中的像元q_g，将q_g在三个坐标方向上在TG_g规则格网上和o′_g相差的格网单元数记为q_g的坐标(q_gx，q_gy，q_gz)；

步骤6.将样本数据集合H中的每个样本h指定到距离其最近的TG网格单元中；对每个TG中的网格单元q，样本数据中被指定到网格单元q中的样本集合记为H_q，q所对应的值为一个向量(s₁(q)，s₂(q)，s₃(q))，其中s_i(q)表示H_q中被标识为s_i类别的样本数有多少。

本实例的这一步骤中20个样本分别被指定到其对应的20个不同的像元中。例如第一个样本其x和y坐标为15和7，类别为s₂，这一样本被赋值到 TG坐标为(15，7)的像元处，该像元对应的值设置为(0，1，0)。

步骤7.对每个TG_g中的网格单元q_g，TG中与其对应的像元组成一个集合 F(q_g)＝{q∈TG|2^g-1×q_gx≤q_x＜2^g-1×(q_gx+1)Λ2^g-1×q_gy≤q_y＜ 2^g-1×(q_gy+1)Λ2^g-1×q_gz≤q_z＜2^g-1×(q_gz+1)}；当计算出F(q_g)后，将网格单元q_g对应的值设定为(s₁(q_g)，s₂(q_g)，s₃(q_g))，其中s_i(q_g)表示F(q_g)中具有s_i类别标识的网格单元个数有多少。

以图4为例，该图左面子图为已经将样本指定之后的目标规则格网，其中虚线框表示的像元不包含任何样本，白色实线框表示的像元表示包含只出现一次W类地物，黑色实线框表示的像元表示只出现一次B类地物。图4右面子图为粒度为2的目标规则格网。每个方框代表该粒度下的一个像元。都以最左上角顶点作为原点，那么两个子图中左上角顶点的坐标都为(0，0)，并且图4右面子图中左上角顶点对应的F((0，0))集合为{(0，0)，(0，1)，(1，0)，(1，1)}这四个左面子图中的像元。以左上角为例，其值为一个向量(0，1)，表示只出现W类地物。以右面子图中间位置在坐标为的像元为例，其在粒度为2的目标规则格网中的坐标为(1，2)，F((1，2))集合为{(2，4)，(3，4)，(2，5)，(3，5)}。左面子图中这四个像元中B 类地物和W类地物分别出现一次，因此在粒度为2的目标规则格网中的坐标为(1，2)的像元对应的值为一个向量(1，1)，表示W类地物和B类地物各出现一次。

步骤8.根据给定的近邻搜索半径，初始化近邻模板T；T由距离坐标原点不超过τ个网格单元的整数位置矢量构成；对于T中的每个元素t，初始化类别矩阵

为0；

每一行表示一个类别组合；例如有两个类别s₁和 s₂，则

有三行，分别表示t中只出现类别s₁，只出现类别s₂以及同时出现类别s₁和s₂；

第i列代表各种组合下s_i类别在中心像元出现的个数；

二维情况下，当近邻搜索半径为20时，该模板共有1257个元素。对于每个元素，需要初始化一个类别矩阵，由于类别数为3，其大小为7×3。从第一行到第七行分别表示类别组合{s₁}，{s₂}，{s₃}，{s₁，s₂}，{s₁，s₃}，{s₂，s₃}和 {s₁，s₂，s₃}。以第一行为例，其表示{s₁}类别组合出现时，中心像元出现s₁，s₂，s₃三个类别的数目分别是多少。

步骤9.设定g＝2，并将TI中各类别的比例作为各类别在训练图像中的概率质量函数，其值为{P_ti(s₁)，P_ti(s₂)，…，P_ti(s_C)}＝{0.3122，0.4333，0.2545}，同时设置TG_g中每个网格单元q_g各类别的先验概率质量函数

步骤10.初始化所有类别矩阵

为0，t∈T；对于TI_g上的每个网格单元p_g，在TI_g中寻找同p_g距离不超过τ个网格单元集合P_g；对于P_g中的每个网格单元δ，其坐标为(δ_x，δ_y，δ_z)，在近邻模板T中找到δ对应的元素t^δ＝(δ_x- p_gx，δ_y-p_gy，δ_z-p_gz)；在

中找到网格单元δ所有出现类别所对应的类别组合所在行i，并将p_g所对应的值(s₁(p_g)，s₂(p_g)，s₃(p_g))累加到

的i行。

例如，对于训练图像的每个像元如果发现当其x轴右侧侧紧邻(模板中对应的矢量为t^δ＝(1，0))的像元类别组合为{s₁}时，则将该像元的三个类别出现次数的值累加到这一模板元素对应的

的第1行的三个类别数目的值上。

步骤11.对于TG_g中的每个网格单元q_g计算各类别的概率质量函数，具体方法如下：

(6)将(P(s₁(q_g))，P(s₂(q_g))，P(s₃(q_g)))初始化为0；

(7)在TG_g中寻找同q_g距离不超过τ个网格单元集合TQ_g；

(8)对于TQ_g中的每个网格单元γ，其坐标为(γ_x，γ_y，γ_z)；如果 (s₁(γ)，s₂(γ)，s₃(γ))为0向量，则从TQ_g中去除γ；如果(s₁(γ)，s₂(γ)，s₃(γ))不为 0向量，则在近邻模板T中找到γ对应的元素t^γ＝(γ_x-p_gx，γ_y-p_gy，γ_z-p_gz)；在

中找到网格单元γ所有出现类别所对应的类别组合所在行j，然后设置

为

第j行的值。

(9)对每个类别s_i使用公式

融合TQ_g中不同网格单元所推导出的类别概率，其中

计算出所有的P(s_i(q_g))后，使用

(10)对P(s_i(q_g))进行归一化得到每个类别的概率。如果g≠1，在TG_g-1中找到和q_g对应的网格单元集合Q_g-1＝{q_g-1∈TG_g-1|2×q_gx≤q_(g-1)x＜ 2×(q_gx+1)Λ2×q_gy≤q_(g-1)y＜2×(q_gy+1)Λ2×q_gz≤q_(g-1)z＜ 2×(q_gz+1)}，并将q_g-1各类别的先验概率质量函数

设置为{P(s₁(q_g))，P(s₂(q_g))，P(s₃(q_g))}，设定 g＝g-1，跳转到步骤10；

在这一步中，每个粒度为g的目标规则格网中像元上所估算出的概率质量函数都作为粒度为g-1的目标规则格网中相应像元在下一次估算中的先验。

步骤12.将TG₁，TG₂上每个网格单元的各类别概率质量函数结果输出。图5为TG₂上每个网格单元的各类别概率质量函数结果。该图从左到右依次为 s₁，s₂，s₃三个类别的概率估算结果。最后一个图例为概率值。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的权利要求范围之内。

Claims (1)

1.一种逐位置快速估算各类地物概率分布的方法，其特征在于：包括以下步骤：

步骤1.给定估算的所需参数，包括最大粒度G，近邻搜索半径τ，以及包含m个样本的样本数据集合H＝{h₁，h₂，…，h_m}，给出研究区域的所有地物类别，并将C个类别组成的集合记为{s₁，s₂，…，s_C}，G为一个正整数参数，τ为一个正实数参数；

步骤2.根据目标区域大小，设定目标规则格网TG，并根据目标地物的空间结构特征选择和目标规则格网具有相同空间分辨率的训练图像TI，其中空间分辨率为每个格网单元所表示的地面实际尺寸的大小；

步骤3.将TI缩放到从1到G的不同粒度大小，并分别记为TI₁，TI₂，…，TI_G；在将TI缩放到粒度1≤g≤G过程中，选定TI的某个顶点格网单元o为坐标原点，其坐标为(0，0，0)；对于每个TI中的像元p，将p在三个坐标方向上在TI规则格网上和o相差的格网单元数记为p的坐标(p_x，p_y，p_z)；同时，选定TI_g中对应于o的顶点格网单元o_g为坐标原点，其坐标为(0，0，0)；对于每个TI_g中的像元p_g，将p_g在三个坐标方向上在TI_g规则格网上同o_g相差的格网单元数记为p_g的坐标(p_gx，p_gy，p_gz)；

步骤4.对每个TI_g中的网格单元p_g，TI中与其对应的像元组成一个集合F(p_g)＝{p∈TI|2^g-1×p_gx≤p_x＜2^g-1×(p_gx+1)∧2^g-1×p_gy≤p_y＜2^g-1×(p_gy+1)∧2^g-1×p_gz≤p_z＜2^g-1×(p_gz+1)}；当计算出F(p_g)后，将网格单元p_g对应的值设定为(s₁(p_g)，s₂(p_g)，…，s_C(p_g))，其中s_i(p_g)表示F(p_g)中被标识为s_i类别的网格单元个数有多少；

步骤5.将TG缩放到从1到G的不同粒度大小，并分别记为TG₁，TG₂，…，TG_G；在将TG缩放到粒度1≤g≤G过程中，选定TG的某个顶点格网单元o′为坐标原点，其坐标为(0，0，0)；对于每个TG中的像元q，将q在三个坐标方向上在TG规则格网上和o′相差的格网单元数记为q的坐标(q_x，q_y，q_z)；同时，选定TG_g中对应于o′的顶点格网单元o′_g为坐标原点，其坐标为(0，0，0)；对于每个TG_g中的像元q_g，将q_g在三个坐标方向上在TG_g规则格网上同o′_g相差的格网单元数记为q_g的坐标(q_gx，q_gy，q_gz)；

步骤6.将样本数据集合H中的每个样本h指定到距离其最近的TG网格单元中；对每个TG中的网格单元q，样本数据中被指定到网格单元q中的样本集合记为H_q，q所对应的值为一个向量(s₁(q)，s₂(q)，…，s_C(q))，其中s_i(q)表示H_q中被标识为s_i类别的样本数有多少；

步骤7.对每个TG_g中的网格单元q_g，TG中与其对应的像元组成一个集合F(q_g)＝{q∈TG|2^g-1×q_gx≤q_x＜2^g-1×(q_gx+1)∧2^g-1×q_gy≤q_y＜2^g-1×(q_gy+1)∧2^g-1×q_gz≤q_z＜2^g-1×(q_gz+1)}；当计算出F(q_g)后，将网格单元q_g对应的值设定为(s₁(q_g)，s₂(q_g)，…，s_C(q_g))，其中s_i(q_g)表示F(q_g)中具有s_i类别标识的网格单元个数有多少；

步骤8.根据给定的近邻搜索半径τ，初始化近邻模板T；T由距离坐标原点不超过τ个网格单元的整数位置矢量构成；对于T中的每个元素t，初始化类别矩阵

为0；

每一行表示一个类别组合；例如有两个类别s₁和s₂，则

有三行，分别表示t中只出现类别s₁，只出现类别s₂以及同时出现类别s₁和s₂；

第i列代表各种组合下s_i类别在中心像元出现的个数；

步骤9.设定g＝G，并将TI中各类别的比例作为各类别在训练图像中的概率质量函数，记为{P_ti(s₁)，P_ti(s₂)，…，P_ti(s_C)}，同时设置TG_g中每个网格单元q_g各类别的先验概率质量函数为

步骤10.初始化所有类别矩阵

为0，t∈T；对于TI_g上的每个网格单元p_g，在TI_g中寻找同p_g距离不超过τ个网格单元集合P_g；对于P_g中的每个网格单元δ，其坐标为(δ_x，δ_y，δ_z)，在近邻模板T中找到δ对应的元素t^δ＝(δ_x-p_gx，δ_y-p_gy，δ_z-p_gz)；在

中找到网格单元δ所有出现类别所对应的类别组合所在行i，并将p_g所对应的值(s₁(p_g)，s₂(p_g)，…，s_C(p_g))累加到

的i行；

步骤11.对于TG_g中的每个网格单元q_g计算各类别的概率质量函数，具体方法如下：

(1)将(P(s₁(q_g))，P(s₂(q_g))，…，P(s_C(q_g)))初始化为0；

(2)在TG_g中寻找同q_g距离不超过τ个网格单元集合TQ_g；

(3)对于TQ_g中的每个网格单元γ，其坐标为(γ_x，γ_y，γ_z)；如果(s₁(γ)，s₂(γ)，…，s_C(γ))为0向量，则从TQ_g中去除γ；如果(s₁(γ)，s₂(γ)，…，s_C(γ))不为0向量，则在近邻模板T中找到γ对应的元素t^γ＝(γ_x-p_gx，γ_y-p_gy，γ_z-p_gz)；在

中找到网格单元γ所有出现类别所对应的类别组合所在行j，然后设置

为

第j行的值；

(4)对每个类别s_i使用公式

融合TQ_g中不同网格单元所推导出的类别概率，其中

计算出所有的P(s_i(q_g))后，使用

对P(s_i(q_g))进行归一化得到每个类别的概率；

(5)如果g≠1，在TG_g-1中找到和q_g对应的网格单元集合Q_g-1＝{q_g-1∈TG_g-1|2×q_gx≤q_(g-1)x＜2×(q_gx+1)∧2×q_gy≤q_(g-1)y＜2×(q_gy+1)∧2×q_gz≤q_(g-1)z＜2×(q_gz+1)}，并将q_g-1各类别的先验概率质量函数

设置为{P(s₁(q_g))，P(s₂(q_g))，…，P(s_C(q_g))}，设定g＝g-1，跳转到步骤10；

步骤12.将TG₁，TG₂，…，TG_G上每个网格单元的各类别概率质量函数结果输出。

Images