CN103268423A - 一种地理现象多点模拟空间尺度选择的方法 - Google Patents

Abstract

一种地理现象多点模拟空间尺度选择的方法：（1）计算模板中心和边界的最短距离；（2）计算训练图像空间相关性小于指定阈值的临界距离；（3）通过模板中心和边界的最短距离以及空间相关性临界距离确定所需考虑的空间尺度大小即多重网格数的大小。本发明能够在进行多点模拟前根据空间相关性确定所需空间尺度的大小，避免了由专家经验选择空间尺度所带来的主观性和不确定性，为生产生活实践中解决相关问题，例如石油储层模拟，提供多重网格数参数设定的方法。本发明可以应用到地质矿产、气象、遥感、海洋研究、军事侦察及环境监测等领域。

Description

一种地理现象多点模拟空间尺度选择的方法

技术领域

本发明属于空间信息技术领域，具体地说，涉及一种对多点模拟中使用多重网格进行模拟时多重网格数参数的确定方法。

技术背景

模拟一直是研究地理现象的一种常用技术和方法。例如对在气象、地质等很多方面，很多研究都是关于如何从有限样本点模拟出研究对象在研究区域的实际分布。20世纪90年代Guardiano,Srivastava和Journel等提出了使用多点统计学来对进行地理现象建模的这一理论方法，即多点模拟。这一方法可以弥补传统地质统计学变差函数的不足：仅考虑两点之间的统计关系。基于这一思想，专家们提出了很多相关算法并进行了验证。这些算法的基本方法都是从训练图像中提取空间结构模式并将其应用到目标图像上。早期提出的多点模拟算法，不论进行估计还是模拟，其时间复杂度较高。Strebelle在其博士论文中提出的SNESIM算法首次成功的解决了这一问题。SNESIM算法引入搜索树这一存储结构。在模拟之前，SNESIM算法将训练图像中的空间结构信息通过一次扫描存储到搜索树中。在进行模拟时，SNESIM算法对每个像元都根据周边条件数据从搜索树中提取条件概率，从而降低了时间复杂度。SNESIM算法目前已经得到了广泛应用，并且有了进一步的发展。

虽然目前SNESIM算法已经有了比较深入的研究。但是这一方法需要大量的参数，对于这些参数的设置通常都需要使用者的经验来确定，例如多重网格数、模板、训练图像，最少条件数据数等等。对于这些参数，现有的对这方面的研究基本上是进行灵敏度分析，例如：Strebelle在其博士论文中对不同参数就进行了比较细致的对比分析。在此基础上Liu发表相关论文对SNESIM算法进行了灵敏度分析。但是这方面的研究仍然没有一种方法可以确定这些参数。多重网格数是一个关系到模拟质量的至关重要的参数，多重网格数代表了模拟过程中所需考虑最大空间尺度的大小，该参数直接关系到模拟的质量和速度，如果使用过小的网格数，那么就无法有效模拟出目标地物的结构特征，但是如果使用比较大的网格数，由于需要多次建立搜索树，那么算法的效率就会降低，并且Liu讨论过，在大于某一特定网格数之后，模拟的质量提升就非常有限。如果能够在模拟前根据地物空间分布特征确定这一参数，就可以在保证模拟质量的前提下尽可能的提高模拟效率。

空间相关性是空间分析中的最为重要的技术之一。针对名义变量，目前被广泛使用和接受的有效的空间相关性探测方法是连接数统计方法，这一方法主要是将实际观测到的不同类别邻接对的数目同理论上随机分布的这一数目进行比较，如果实测值大于理论值，也就是说不同地物更容易邻接，那么地物分布在空间上就是负相关的。反之，如果实测值小于理论值，也就是说同类地物更容易邻接，那么地物分布在空间上就是正相关的。但是Goodchild指出这一方法并没有像Moran’s I等连续值变量那样给出一个特定的指标来反映空间相关程度，只是给出判断是否相关。本发明中对这一方法进行了归一化，给出了类似Moran’s I但是针对分类变量的空间相关性系数。

空间相关性为探测目标地物是否具有大尺度的空间结构特征提供了一个有效的工具。如果大尺度上的空间相关性是随机的，不同地物都是随机分布的，那么就可以得出这一尺度上就不具有空间相关性。也就没有必要再多点模拟的过程中对这一尺度的信息进行建模，而根据地理学第一定律：“任何事物都相关，只是相近的事物关联更紧密”，大于这一尺度的空间相关性也一定较弱，考虑大于这一尺度的结构特征也并不会显著提高多点模拟质量。

基于上述原理，本发明提出了一种地理现象多点模拟空间尺度选择的方法。这一方法主要是通过计算目标地物训练图像的空间相关性来确定在何种尺度上进行多点模拟。所确定的这一尺度同时结合模板结构就能够在多点模拟中自动确定多重网格数，使用这一多重网格数可以确保模拟结果的质量。

发明内容

本发明的技术解决问题：在进行石油储层模拟等应用领域中进行多点模拟时根据地物特征自动确定所需要的空间尺度大小。该方法使用名义变量的空间相关性作为依据，使用了规范化的连接数统计来计算空间相关性系数，并且根据这一系数结合模板结构来确定所需的多重网格数。本发明提出的多重网格数确定方法可以在保证模拟质量的前提下，尽可能提高模拟效率，避免专家经验所具有的不确定性，同时又无需比较设置不同参数的模拟结果就能得出质量满足要求的模拟结果，便于多点模拟的实际应用。

本发明的技术解决方案：一种地理现象多点模拟空间尺度选择的方法，包括以下基本步骤：

步骤1、确定k阶邻接矩阵W_k；

k阶邻接矩阵W_k为一个n×n大小的矩阵，n为训练图像的像元数，邻接矩阵中第i个像元和第j个像元如果k阶相邻那么其对应W_k中的元素w_k(i,j)的值为1，否则w_k(i,j)的值为0，在不同应用中可以从以下三种邻接矩阵选择一个：

（1）当综合考虑所有方向上的空间尺度时，使用

其中如果|i_x-j_x|+|i_y-j_y|=k，那么

否则

（2）当仅考虑东西方向上的空间尺度时，使用

其中如果|i_x-j_x|=0且|i_y-j_y|=k，那么 $w_k^{\mathrm{EW}}(i,j)=1,$ 否则 $w_k^{\mathrm{EW}}(i,j)-0;$

（3）当仅考虑南北方向上的空间尺度时，使用其中如果|i_x-j_x|=k且|i_y-j_y|=0，那么 $w_k^{\mathrm{NS}}(i,j)=1,$ 否则 $w_k^{\mathrm{NS}}(i,j)=0;$

在上述三种邻接矩阵的定义中i_x和i_y分别为i号像元的横纵坐标，j_x和j_y分别为j号像元的横纵坐标；

步骤2、计算模板τ中心和边界的最短距离d；

（1）设置d'初始值为1；

（2）计算向量集合

其中D(x,y)是向量(x,y)与模板中心的邻接的阶数，如果w_m((0,0),(x,y))=1，那么D(x,y)=m，V_d'的具体计算过程为对于每一个满足|x'|+|y'|小于等于d'的向量计算其D(x',y')，如果D(x',y')小于等于d'，就将其放入到集合V_d'中；

（3）如果V_d'包含于τ，那么d'=d'+1并重复步骤（2）和（3），否则d=d'-1；

步骤3、计算训练图像Ti中空间相关性小于指定阈值a的临界距离s；

（1）设置s'的初值为1；

（2）计算s'阶训练图像空间相关性SAC_s'(Ti)，其中k阶空间相关性SAC_k(Ti)计算方法如下：

（I）根据邻接矩阵W_k，选择出所有相邻的像元对集合J_all={(i,j)|w_k(i,j)=1}，这里W_k(i,j)表示邻接矩阵第i行第j列元素的值；

（II）在J_all中选择出所属类别不相同的邻接像元对组成集合

这里C_i和C_j分别表示像元i和j所属的类别；

（III）计算J_rs中元素个数同J_all中元素个数的比例

（IV）根据不同类别C={C₁,C₂,...,C_m}像元所占比例R(C₁),R(C₂),...,R(C_m)，计算不同类别像元对在空间随机分布情况下J_all中所属类别不相同的邻接像元对组成集合

所占比例 $R_{\mathrm{ran}}\left(\mathrm{rs}\right)=\underset{C_i\≠C_j}{\mathrm{\Σ}}R\left(C_i\right)R\left(C_j\right);$

（V）按照公式（3）计算空间相关性系数：

${\mathrm{SAC}}_k\left(\mathrm{Ti}\right)=\left\{\begin{array}{cc}\frac{R_{\mathrm{ran}}\left(\mathrm{rs}\right)-R\left(\mathrm{rs}\right)}{R_{\mathrm{ran}}\left(\mathrm{rs}\right)},& R_{\mathrm{ran}}\left(\mathrm{rs}\right)-R\left(\mathrm{rs}\right)>0\\ \frac{R_{\mathrm{rna}}\left(\mathrm{rs}\right)-R\left(\mathrm{rs}\right)}{1-R_{\mathrm{ran}}\left(\mathrm{rs}\right)},& R_{\mathrm{ran}}\left(\mathrm{rs}\right)-R\left(\mathrm{rs}\right)\≤0\end{array}\right.---\left(1\right)$

（3）如果SAC_s'(Ti)大于指定阈值a，那么s'=s'+1，并重复执行（2）和（3），否则s=s'-1。

步骤4、根据计算得到的s和d按照公式（2）确定多点模拟的空间尺度即多重网格数g：

本发明与现有技术相比的优点在于：

（1）使用空间相关性来确定多点模拟参数；

（2）提出并使用了一种规范化的名义变量空间相关性计算方法；

（3）不依赖于专家经验自动确定多点模拟的多重网格数。

（4）根据需要确定阈值a，从空间相关性的意义上去控制模拟结果的质量。

附图说明

图1    本发明的实现流程图

图2    本发明实施实例1所使用的训练图像

图3    本发明实施实例1和实施实例2所使用的模板

图4    本发明实施实例1所生成的模拟结果

图5    本发明实施实例2所使用的训练图像

图6    本发明实施实例2所生成的模拟结果

具体实施方式

图1为本发明确定空间尺度大小方法的流程图，下面将结合实例给出两种具体实施例子。

实施例1

第一个实例对典型的各项同性地物进行模拟，这里使用的目标地物假定为圆形，这里简单的对两种地物分别根据其颜色命名为白色地物和黑色地物。图2是进行多点模拟时需要的训练图像，大小为400×300，图3是进行多点模拟时需要的模板。

整个实施过程可以分为四个步骤:

（1）选择邻接矩阵。

这一实例中的目标地物是圆形，是典型的各向同性地物，因此需要综合考虑各个方向的空间尺度。因此采用了邻接矩阵

其中如果|i_x-j_x|+|i_y-j_y|=k，那么 $w_k^{\mathrm{dis}}(i,j)=1,$ 否则 $w_k^{\mathrm{dis}}(i,j)=0,$

（2）计算模板边界元素中距离中心最近的距离。

这一步需要通过前面定义的

来确定d。首先设置d'=1，在d'阶邻接范围内所有的矢量为{(0,1),(1,0),(-1,0),(0,-1)}，这些矢量都是在模板范围内，然后将d'加1，所有的矢量仍然在模板范围内。直到d'=4时，至少有矢量(0,4)落在模板范围外，所以这里的d值为3。

（3）根据阈值计算空间相关性的临界距离。

通过所选择的邻接矩阵可以计算各阶空间相关性，这里以一阶空间相关性系数为例，首先统计出总的互相连接的像元对数目为|J_all|=239300，其中不同类别连接对的数目为|J_rs|=2392。训练图像中一共有两种地物，白色地物比例为0.715，黑色地物比例为0.285，因此R_ran(rs)=0.204，这样通过公式(1)可以计算得到SAC₁(Ti)=0.951。这里设定阈值为a=0.1。当k=18时，相关性系数才小于0.1，因此最大尺度s为17。

（4）最终确定所需多重网格数。

直接根据计算得到的s和d，通过公式(2)可以计算得到所需的多重网格数

本实例最后在图4中不仅给出了g=4时的多点模拟结果，模拟图像的大小为200×200。图4同时给出了多重网格数为从1到6的多点模拟结果。从图4可以看出，当多重网格数大于等于4时，模拟结果的质量变化不大，也说明了本发明的有效性。

实施例2

第一个实例对典型的各项异性地物进行模拟，这里使用的目标地物假定为具有条带化特征，这里简单的对两种地物分别根据其颜色命名为白色地物和黑色地物。图5是进行多点模拟时需要的训练图像，大小为250×250。图3是进行多点模拟时需要的模板。

整个实施过程同样分为四个步骤：

（1）选择邻接矩阵。

这一实例中的目标地物具有典型的东西方向条带化特征，因此仅需考虑东西方向的空间尺度即可。因此采用了邻接矩阵

其中如果|i_x-j_x|=0且|i_y-j_y|=k，那么 $w_k^{\mathrm{ES}}(i,j)=1,$ 否则 $w_k^{\mathrm{ES}}(i,j)=0.$

（2）计算模板边界元素中距离中心最近的距离。

这一步需要通过前面选择的

来确定d。首先设置d'=1，在d'阶邻接范围内所有的矢量为{(1,0)，(-1,0)}，这些矢量都是在模板范围内，然后将d'加1，所有的矢量仍然在模板范围内。直到d'=4时，至少有矢量(0,4)落在模板范围外，所以这里的d值为3。

（3）根据阈值计算空间相关性的临界距离。

通过所选择的邻接矩阵可以计算各阶空间相关性，这里以一阶空间相关性系数为例，首先统计出总的互相连接的像元对数目为|J_all|=62250（由邻接矩阵决定只统计东西方向），其中不同类别连接对的数目为|J_rs|=830。训练图像中一共有两种地物，白色地物比例为0.722，黑色地物比例为0.278，因此R_ran(rs)=0.201，这样通过公式(1)可以计算得到SAC₁(Ti)=0.935。这里设定阈值为a=0.1。当k=20时，相关性系数小于0.1。因此最大尺度s为19。

（4）最终确定所需多重网格数。

直接根据计算得到的s和d，通过公式(2)可以计算得到所需的多重网格数

本实例最后在图6中不仅给出了g=4时的多点模拟结果，模拟图像的大小为200×200。图6同时给出了多重网格数为从1到6的多点模拟结果。从图6可以看出，当多重网格数大于等于4时，模拟结果的质量变化不大，也说明了本发明的有效性。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，本发明也同样适用于各种不同类型地物以及3维地物的模拟。此外为了避免样本点的影响，本发明实例并没有使用样本点，但使用样本点也同样适用于本发明。本发明的实施实例使用了SNESIM算法，本发明同样适用于其他使用模板和训练图像的多点模拟算法来设定多重网格数。对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的权利要求范围之内。

Claims (1)

1.一种地理现象多点模拟空间尺度选择的方法，其特征在于包括以下步骤：

步骤1、确定k阶邻接矩阵W_k；

步骤2、计算模板τ中心和边界的最短距离d；

步骤3、计算训练图像Ti中空间相关性小于指定阈值α的临界距离s；

步骤4、根据计算得到的s和d按照下面公式确定多点模拟的空间尺度即多重网格数g：

所述步骤1中k阶邻接矩阵W_k为一个n×n大小的矩阵，n为训练图像的像元数，邻接矩阵中第i个像元和第j个像元如果k阶相邻那么其对应W_k中的元素w_k(i,j)的值为1，否则w_k(i,j)的值为0，在不同应用中可以从以下三种邻接矩阵选择一个：

（1）当综合考虑所有方向上的空间尺度时，使用其中如果|i_x-j_x|+|i_y-j_y|=k，那么 $w_k^{\mathrm{dis}}(i,j)=1,$ 否则 $w_k^{\mathrm{dis}}(i,j)=0;$

（2）当仅考虑东西方向上的空间尺度时，使用

其中如果|i_x-j_x|=0且|i_y-j_y|=k，那么 $w_k^{\mathrm{EW}}(i,j)=1,$ 否则 $w_k^{\mathrm{EW}}(i,j)=0;$

（3）当仅考虑南北方向上的空间尺度时，使用

其中如果|i_x-j_x|=k且|i_y-j_y|=0，那么 $w_k^{\mathrm{NS}}(i,j)=1,$ 否则 $w_k^{\mathrm{NS}}(i,j)=0;$

在上述三种邻接矩阵的定义中i_x和i_y分别为i号像元的横纵坐标，j_x和j_y分别为j号像元的横纵坐标；

所述步骤2中按照下面步骤计算模板τ中心和边界的最短距离d：

（1）设置d'初始值为1；

（2）计算向量集合V_d'={(x,y)|D(x,y)≤d',x∈Z,y∈Z}，其中D(x,y)是向量(x,y)与模板中心的邻接的阶数，如果w_m((0,0),(x,y))=1，那么D(x,y)=m，V_d'的具体计算过程为对于每一个满足|x'|+|y'|小于等于d'的向量计算其D(x',y')，如果D(x',y')小于等于d'，就将其放入到集合V_d'中；

（3）如果V_d'包含于τ，那么d'=d'+1并重复步骤（2）和（3），否则d=d'-1；

所述步骤3中按照下面步骤计算训练图像Ti的k阶空间相关性SAC_k(Ti)：

（1）根据邻接矩阵W_k，选择出所有相邻的像元对集合J_all={(i,j)|w_k(i,j)=1}；

（2）在J_all中选择出所属类别不相同的邻接像元对组成集合J_rs={(i,j)∈J_all|C_i≠C_j}，这里C_i和C_j分别表示像元i和j所属的类别；

（3）计算J_rs中元素个数同J_all中元素个数的比例

（4）根据不同类别C={C₁,C₂,...,C_m}像元所占比例R(C₁),R(C₂),...,R(C_m)，计算不同类别像元对在空间随机分布情况下J_all中所属类别不相同的邻接像元对组成集合

所占比例 $R_{\mathrm{ran}}\left(\mathrm{rs}\right)=\underset{C_i\≠C_j}{\mathrm{\Σ}}R\left(C_i\right)R\left(C_j\right);$

（5）按照下面公式计算空间相关性系数：

${\mathrm{SAC}}_k\left(\mathrm{Ti}\right)=\left\{\begin{array}{cc}\frac{R_{\mathrm{ran}}\left(\mathrm{rs}\right)-R\left(\mathrm{rs}\right)}{R_{\mathrm{ran}}\left(\mathrm{rs}\right)},& R_{\mathrm{ran}}\left(\mathrm{rs}\right)-R\left(\mathrm{rs}\right)>0\\ \frac{R_{\mathrm{rna}}\left(\mathrm{rs}\right)-R\left(\mathrm{rs}\right)}{1-R_{\mathrm{ran}}\left(\mathrm{rs}\right)},& R_{\mathrm{ran}}\left(\mathrm{rs}\right)-R\left(\mathrm{rs}\right)\≤0\end{array}\right.$

所述步骤3中按照下面步骤计算训练图像Ti中空间相关性小于指定阈值α的临界距离s：

（1）设置s'的初值为1；

（2）计算s'阶训练图像空间相关性SAC_s'(Ti)；

（3）如果SAC_s'(Ti)大于指定阈值α，那么s'=s'+1，并重复执行（2）和（3），否则s=s'-1。

Images