[0054] 与前述增广积分公式等效，0行及0列的值预置为0，i及j均为大于等于1的整数，j′及i′均为大于等于0的整数。其区别在于前述公式利用前面步骤所得的积分图值直接代入到后继计算中使得增广积分图的构置速度更快，减少运算量，提高整个栅格数据局部奇异性迭代分析的效率，其中SAT(i,j)表示置于所述初始增广积分图行索引为i列索引为j列的积分图值，ρ(i′,j′)代表所述栅格数据中行索引为i′列索引为j′的数据， 代表对所有满足条件i′
[0055] 以下表2为例，表2为积分图构置模块203根据表1所示的栅格数据构置的增广积分图：
[0056] 表2
[0057]
[0058]
[0059] 局部奇异性计算模块204，用于根据所述积分图、多个尺度各不相同且同中心的预设方形窗口及预设的正规化尺度参数，建立局部奇异性模型，以获得局部系数栅格数据和奇异性指数余维数栅格数据。本实施例中，如图5所示，局部奇异性计算模块204包括窗口平均值计算子模块2041、正规化处理子模块2042、奇异性模型拟合子模块2043及获得子模块2044。
[0060] 窗口平均值计算子模块2041，用于根据所述预设方形窗口及所述积分图，计算每个所述预设方形窗口所覆盖的所述栅格数据的平均值，其中，所有的所述预设方形窗口在栅格数据对应的焦元位置相同，所述焦元为所述预设方形窗口中心对应的所述栅格数据的数据单元。其中，主处理机11将构置的增广积分图拷贝给被激活的每个从处理机12。被激活的每个从处理机12分别根据预先分配得到的预设方形窗口尺度及复制所得的增广积分图计算对应的平均值。需要说明的是，主处理机11及被激活的每个从处理机12对应的不同尺度的预设方形窗口在栅格数据对应的焦元位置相同，所述焦元为所述预设方形窗口中心对应的所述栅格数据的数据单元。这就使如果把不同的预设方形窗口放置在同一份栅格数据上时，预设方形窗口呈同心嵌套。需要说明的是，当预设计算方框中存在尺度为1x1这种极端情形时，则尺度为1x1的预设计算方框直接分配给主处理机11，其所需计算的平均值就是当前焦元位置的栅格数据值。以图6为例，图6所示4个不同尺度的预设方形窗口覆盖于表1所示的栅格数据上，4个预设方形窗口同心嵌套且分别覆盖了栅格数据中的1个数据，9个数据、25个数据及49个数据，四个预设方形窗口的中心在栅格数据对应的焦元的行索引为3，列索引为3。在本实施例中，计算每个所述预设方形窗口对应的平均值时，主处理机11与被激活的从处理机12同时计算，提高速度。
[0061] 在其他实施例中，即当运用于串行处理系统时，则依次计算每个预设方框对应的平均值。
[0062] 在本实施例中，窗口平均值计算子模块2041执行所述计算每个所述预设方形窗口所覆盖的所述栅格数据的平均值的方式包括：
[0063] 根据所述预设方形窗口、所述预设方形窗口对应的所述焦元位置及增广积分图，利用公式：
[0064] SUM(i,j,w,h)＝SAT(r2,c2)+SAT(r1,c1)-SAT(r1,c2)-SAT(r2,c1)、
[0065]
[0066]
[0067] 及
[0068]
[0069] 计算所述预设方形窗口所覆盖的所述栅格数据的和，其中，SUM(i,j,w,h)代表宽为w且长为h的所述预设方形窗口在所对应的焦元为所述栅格数据中行索引为i且列索引为j的数据单元时所覆盖的所述栅格数据的和，r1及r2为在所述增广积分图获取积分图值的行索引，c1及c2为在所述增广积分图中获取积分图值的列索引，SAT(r2,c2)为所述增广积分图在行索引为r2列索引为c2处的积分图值，SAT(r1,c1)为所述增广积分图在行索引为r1列索引为c1处的积分图值，SAT(r1,c2)为所述增广积分图在行索引为r1列索引为c2处的积分图值，SAT(r2,c1)为所述增广积分图在行索引为r2列索引为c1处的积分图值。SAT(r2,c2)、SAT(r1,c1)、SAT(r1,c2)及SAT(r2,c1)则为预设方形窗口在增广积分图中对应的四个积分图值。
[0070] 以图7为例，图7中左侧为栅格数据外侧为索引，横向的为列索引，纵向为行索引。右侧为与左侧栅格数据对应的增广积分图，增广积分图外侧为索引，横向的为列索引，纵向为行索引。左侧为栅格数据有10行10列(即总行数M＝10，总列数N＝10)，右侧增广积分图总行数为(M+1)＝11、总列数为(N+1)＝11。数据单元的边长为1，预设方形窗口的宽w＝5，长h＝5，它的中心在栅格数据中对应的焦元行列索引为3列索引为3，即索引为(3,3)，i＝3、j＝
3。预设方形窗口在增广积分图中对应四个积分图值的行索引分别为r1＝i-(h-1)/2＝3-2＝1、r2＝i+(h+1)/2＝3+3＝6、列索引分别为c1＝j-(w-1)/2＝3-2＝1及c2＝j+(w+1)/2＝3+
3＝6。根据增广积分图可获得SAT(r1,c1)＝SAT(1,1)＝7，SAT(r1,c2)＝SAT(1,6)＝22，SAT(r2,c1)＝SAT(6,1)＝24，SAT(r2,c2)＝SAT(6,6)＝211，则此时SUM(3,3,5,5)＝SAT(6,6)+SAT(1,1)-SAT(1,6)-SAT(6,1)＝211+7-22-24＝172。
[0071] 预设方形窗口在栅格数据图中进行滑动计算，预设方形窗口对应的焦元靠近栅格数据边缘位置时，会出现预设方形窗口越出栅格数据坐标范围的情况，预设方形窗口实际上只有其中一部分覆盖有栅格数据。预设方形窗口在增广积分图中对应的四个积分图值中至少一个积分图值的索引(包括行索引和/或列索引)出现越界，进而无法根据索引从增广积分图中获得有效的积分图值。因此当预设方形窗口出现越界情况时，需要对与预设方形窗口对应的四个积分图值的索引进行预处理。预设方形窗口在增广积分图中对应的四个积分图值所涉及四个坐标，即行索引r1、r2和列索引c1、c2，需要按如下四种基本情形对索引进行修正：
[0072] ①当计算窗口向上越界时，即行索引r1<0，则行索引
修改为0，也即r1＝0，此时行索引为r1的积分图值就是增广积分图中第0行中预置的0值；
[0073] ②当计算窗口向左越界时，即列索引c1<0，则列索引修改为0，也即c1＝0，此时列索引c1的积分图值就是增广积分图中第0列中预置的0值；
[0074] ③当计算窗口向下越界时，即行索引r2大于栅格数据总行数M，则行索引r2修改为M，也即r2＝M，此时的积分图取值为SAT(M,j)，其中j为[0,N]区间的任意值，N为栅格数据总列数；
[0075] ④当计算窗口向右越界时，即列索引c2大于栅格数据列数N，则列索引c2修改为N，也即c2＝N，此时的积分图取值为SAT(i,N)，其中i为[0,M]区间的任意值。
[0076] 需要指出的是，上述四个坐标的预处理需要逐个判别并根据需要进行相应的修正，在全部完成之后，才能根据上述积分图公式计算获得所述预设计算方框实际覆盖数据的和值。
[0077] 计算所述预设方形窗口所覆盖的所述栅格数据的平均值，计算方法为根据所述预设方形窗口所覆盖的所述栅格数据的和除以所述栅格数据的个数，具体为公式：
[0078]
[0079] 其中，MEAN(i,j,w,h)表示宽为w，长为h的预设方形窗口在栅格数据中对应的焦元索引为(i,j)时，预设方形窗口所覆盖的栅格数据的平均值，SUM(i,j,w,h)代表宽为w且长为h的所述预设方形窗口在栅格数据中对应的焦元索引为(i,j)时，所覆盖的所述栅格数据的和；T(i,j,w,h)代表宽为w且长为h的所述预设方形窗口在栅格数据中对应的焦元索引为(i,j)时，预设方形窗口内的覆盖栅格数据的个数。统计所述预设方形窗口内覆盖栅格数据的个数为两种情形：
[0080] (1)当预设方形窗口完全覆盖于栅格数据时，也即r1≥0且r2≤M且c1≥0且c1≥0，行索引r1、r2和列索引c1、c2在预处理中无需修正索引值，每个预设方形窗口内覆盖的栅格数据的个数是固定的，覆盖个数仅与预设方形窗口的尺度有关，T(i,j,w,h)＝w×h。例如，预设方形窗口的宽w＝5，长h＝5(宽和长以数据单元边长为度量单位)，则预设方形窗口覆盖的栅格数据的个数为w×h＝25个。需要指出的是，此种情形下对某个大小的窗口内数据个数是个固定值，其个数统计只需预先计算一次，以后该尺度预设方形窗口的平均值计算都可以直接可以使用该结果，无需每次计算平均值时重复统计个数。
[0081] (2)当预设方形窗口越出栅格数据坐标范围时，也即行索引r1、r2和列索引c1、c2在预处理中有修正索引的情形，需要统计预设方形窗口内覆盖的栅格数据的有效个数。覆盖个数不仅与预设方形窗口的尺度有关，还与窗口位置有关。具体为，根据公式：
[0082] T(i,j,w,h)＝(r2-r1)×(c2-c1)，
[0083] 其中，r1和r2为宽为w长为h的预设方形窗口中心对应焦元的索引为(i,j)时，对应的四个积分图值经过预处理修正后的行索引，c1和c2为经过预处理修正后的列索引。需要指出的是，上述公式是通用的，对预设方形窗口完全覆盖于栅格数据情形也适用，当然用此公式相比之下会增加计算量。
[0084] 以图7为例，该预设方形窗口中心在栅格数据中对应焦元的行索引为i＝3列索引为j＝3，预设方形窗口的宽w＝5，长h＝5，该计算方框完全覆盖在栅格数据范围内，在i＝3，j＝3位置计算所得的积分图值SUM(3,3,5,5)＝172，数据个数T(3,3,5,5)＝5×5＝25，该预设方形窗口对应的平均值MEAN(3,3,5,5)＝172/25＝6.88。
[0085] 再例如，图8所示四个预设方形窗口在增广积分图对应覆盖示意图，所述增广积分图外侧为标注的索引，横向的为列索引，纵向的为行索引。其中由于第0行第0列为增广积分图
构建时添加的，其中设置的0值在统计有效数据时应不计算在内。图8分别示出四个预设方形窗口出现越出栅格数据的不同情形时，预设方形窗口对应在增广积分图的覆盖情形。图9中左上
角的预设方形窗口对应的三个积分图值的索引出现越界，该预设方形窗口中心在栅格数据中对应焦元的行索引为i＝0列索引为j＝0，预设计算方框的宽w＝5，长h＝5，根据前述计算公式，预设计算方框在增广积分图中对应四个积分图值的行索引分别为r1＝i-(h-1)/2＝0-2＝-2、r2＝i+(h+1)/2＝0+3＝3、列索引分别为c1＝j-(w-1)/2＝0-2＝-2及c2＝j+(w+1)/2＝0+3＝3，有3个积分图值SAT(-2,-2)，SAT(-2,3)，SAT(3,-2)是无效的，行索引r1和列c1需要进行修正，预处理后的索引坐标为r1＝0，r2＝3，c1＝0，c2＝3，按顺
时针方向修正后的四个有效的积分图值为SAT(0,0)，SAT(0,3)，SAT(3,3)和SAT(0,3)，于是，在i＝0，j＝0位置计算所得的积分图值SUM(0,0,5,5)＝(29+0-0-0)＝29，数据个数T(0,0,5,5)＝(3-0)×(3-0)＝9，该预设方形窗口对应的平均值MEAN(0,0,5,5)＝29/9＝3.22。再如，图9中其它三个预设方形窗口的情形。图9中左下角的预设方形窗口对应的三个积分图值的索引出现越界，SAT(7,-2)、SAT(12,-2)、SAT(12,3)积分图值无效，在增广积分图中预处理后的索引坐标为r1＝7、r2＝10、c1＝0、c2＝3，平均值为(92+0-63-0)/(3×3)＝29/9＝3.22。图
9中右下角的预设方形窗口对应的三个积分图值的索引出现越界，SAT(7,12)、SAT(12,12)、SAT(12,7)积分图值无效，在增广积分图中预处理后的索引坐标为r1＝0、r2＝3、c1＝7、c2＝
10，平均值为(134+0-105-0)/(3×3)＝29/9＝3.22。图9中右上角的预设方形窗口对应的三个积分图值的索引出现越界，SAT(-2,7)、SAT(-2,12)、SAT(3,12)积分图值无效，在增广积分图中预处理后的索引坐标为r1＝7、r2＝10、c1＝7、c2＝10，平均值为(455+236-325-327)/
9＝39/9＝4.33。
[0086] 需要指出的是，当预设方框覆盖的栅格单元数等于1，也即窗口的宽为1且长为1时，此时平均值即为原来输入栅格数据对应焦元位置的值，只有当预设方框覆盖的栅格单元数大于1时，计算平均值时才需要借助积分图。
[0087] 在计算预设方形窗口所覆盖的栅格数据的和时，无论窗口是大是小，利用积分图都归结为四个积分图值的加减法，比起直接对窗口内的数值简单相加求和，大大减少了运算次数，特别是窗口越大，计算效率更明显。
[0088] 正规化处理子模块2042，用于根据所述正规化尺度参数对所述预设方形窗口的尺度进行正规化处理，以获得每个所述预设方形窗口的正规化尺度大小。具体为，根据所述正规化尺度参数，利用公式：
[0089] st＝lt/L，
[0090] 获得每个所述预设方形窗口的正规化尺度大小。其中，t的取值可以是小于或等于m，且大于0的整数，m为预设方形窗口数量，st代表第t个预设方形窗口的正规化后的尺度大小，lt代表第t个预设方形窗口的尺度，L代表正规化尺度参数。约定lt大小为所述预设方形窗口面积的开方，也即 wt代表第t个预设方形窗口的宽，ht代表第t个预设方形窗口的长。对正方形窗口，lt即为所述预设方框的边长，也即lt＝wt＝ht。需要指出的是，在迭代模型中L应大于1，否则迭代模型不收敛不能得到预期的精度要求；对于常规的非迭代模型，L可以等于1，从而得到更为简化的局部奇异性模型。L的取值大小不影响局部奇异性指数余维数的期望值，但影响局部系数的估计期望值以及局部系数总体个别值预测区间。当正规化尺度取值为所有预设计算方框的尺寸的几何平均值时，这使得局部系数在迭代过程中产生的估计误差最小，可以减少迭代计算中的不确定性的传递。
[0091] 奇异性模型拟合子模块2043，用于根据所有的所述预设方形窗口对应的平均值及对应的所述正规化尺度大小，采用双对数
坐标系下的加权线性回归方法，以获得所述焦元位置所对应的局部奇异性指数余维数及局部系数，其中，所述加权线性回归方法中各个样本点对应的权重赋值为所述样本点对应的所述预设方形窗口的平均值的平方。具体为对任一焦元(i,j)位置，多尺度下的均值及对应的尺度存在如下关系：
[0092] MEAN(i,j,wt,ht)＝c(i,j)st-Δα(i,j)，
[0093] 其中，m为预设方形窗口数量。MEAN(i,j,wt,ht)为所述位于(i,j)第t个宽为wt长为ht的预设方形窗口对应的平均值，st表示第t个预设方形窗口对应的正规化尺度大小。Δa(i,j)和c(i,j)为待定参数，分别称为位于(i,j)的局部奇异性指数余维数和局部系数。
[0094] 局部奇异性模型实质为非线性的幂律模型，可选的，可以通过幂律模型转换为双对数坐标系下的加权线性回归模型来获得上述公式中所述的Δa(i,j)和c(i,j)。幂律模型转换为双对数坐标系下的加权线性回归模型可以提高Δa(i,j)和c(i,j)的预测精度，同时又节省计算量，具体地，通过采用一种特殊的加权方式进行线性模型拟合来逼近原始幂律模型拟合的结果。回归模型形如：
[0095] yt(i，j)＝ln(c(i，j))+(-Δa(i，j))×xt，
[0096] 其中，t＝1,2,3...,m。xt＝ln(st)，yt＝ln(MEAN(i,j,wt,ht))，
[0097] (xt,yt)为第t个预设方形窗口对应的数据点。应用加权最小二乘法，使得二元函数其中Wt为第t个预设方形窗口对应的数据点(xt,yt)的权重。对焦元在(i,j)位置处的加权回归，为达到逼近原始幂律模型拟合效果。具体地，利用公式：
[0098] Wt(i,j)＝[MEAN(i,j,wt,ht)]2，
[0099] 获取每个预设方形窗口对应的权重值。其中，Wt(i,j)为第t个预设方形窗口在积分图中对应焦元为(i,j)时对应的数据点(xt,yt)的权重，MEAN(i,j,wt,ht)为第t个预设方形窗口在积分图中对应焦元为(i,j)时对应的平均值。在此加权方式下获得两个回归参数，也即斜率p和截距q，于是可得Δa(i,j)＝-p，c(i,j)＝eq。可以采用残差平方和的大小来判定相对优劣，残差平方和越小，说明拟合效果越好。残差平方和的计算公式为：
[0100]
[0101] 其中，SSE(i,j)为预设方形窗口对应焦元的索引为(i,j)时对应的残差平方和。
[0102] 以图6为例，图6所示4个不同尺度的预设方形窗口覆盖于表1所示的栅格数据上，4个同心嵌套的预设计算方框的焦元行索引i＝3，列索引j＝3，对应的预设方形窗口尺度分别为1、3、5、7，尺度正规化参数L的推荐参考值为窗口尺度序列的几何平均值，即本例中我们设置为近似值3，预设方形窗口对应的正规化尺度大小分别是0.3333、1.0000、1.6667及2.3333。正规化尺度大小的值在用户进行参数设置后，可以预先计算好，在任意焦元位置处实施奇异性分析可以直接复用。根据增广积分图在行索引为3且列索引为3位置处可获得各个窗口的平均值分别为34.00、9.33、6.88、4.82。图6中的四个由小到大的预设方形窗口在双坐标系下对应的数据点分别为(-1.0986,3.5264)，(0.0000,2.2332)，(0.5108,1.9286)及(0.8473,1.5728)。
[0103] 用幂律模型直接进行非线性拟合时，在i＝3，j＝3位置处，局部奇异性指数余维数Δα(3,3)＝1.0699，局部系数c(3,3)＝10.4569，残差平方和SSE(3,3)＝2.3232。
[0104] 用双对数坐标系下采用常规线性回归方法时，可得，Δα(3,3)＝0.9945，c(3,3)＝10.8024，SSE(3,3)＝5.5322。
[0105] 用双对数坐标系下采用加权线性回归方法时，图6中四个尺度由小到达预设方形窗口对应权重赋值为分别为1156.0、87.0489、47.3344、23.2324，可得，Δα(3,3)＝1.0548，c(3,3)＝10.6337，SSE(3,3)＝2.8354。
[0106] 从预设方形窗口对应焦元索引为(3,3)的局部奇异性分析结果表明，该位置处具有明显的正奇异性，随着尺度由大变小，局部邻域中的平均值(例如代表某种物质的
密度或某种地质过程的能量)越来越高，而且呈幂率增长。加权线性回归方法所得的SSE比常规线性回归方法所得的SSE更小，更加接近幂率模型非线性拟合所得的SSE。但是非线性拟合法需要预先提供相对合理的初始值，通过迭代来不断改进估值直至获得最优解，计算比较耗时。因此，本发明中提出使用双对数坐标系下加权线性回归方法来高效高精度地获得局部奇异性指数余维数及局部系数的估计值。
[0107] 获得子模块2044，用于所述预设方形窗口的中心同步遍历每一个所述栅格数据的数据单元，并重复根据所有的所述预设方形窗口对应的平均值，与正规化尺度大小建立局部奇异性模型，以获得所述栅格数据中每一个数据单元对应的局部奇异性指数余维数及局部系数，进而获得所述局部系数栅格数据及所述奇异性指数余维数栅格数据。
[0108] 迭代模块205，用于重复对最新获得的局部系数栅格数据进行积分图构置及所述局部奇异性模型建立的迭代处理，以获得迭代后的局部系数栅格数据及迭代后的奇异性指数余维数栅格数据，直至所述迭代后的奇异性指数余维数栅格数据满足预设条件或迭代次数达到预设次数。需要说明的是，预设条件为迭代后的奇异性指数余维数栅格数据的所有单元的数据的平方之和小于预设迭代精度。
[0109] 在本实施例中，如图9所示，迭代模块205包括以下功能子模块：
[0110] 积分图构置子模块2051，用于根据最新获得的所述局部系数栅格数据构置对应的积分图。原理同积分图构置模块203相同，在此不再赘述。
[0111] 迭代子模块2052，用于根据所述最新获得的局部系数栅格数据对应的积分图、所述预设方形窗口及预设的正规化尺度参数，建立局部奇异性模型，以获得所述迭代后的局部系数栅格数据及迭代后的奇异性指数余维数栅格数据。
[0112] 迭代子模块2052，还用于判断迭代后的所述奇异性指数余维数栅格数据是否满足所述预设条件或所述迭代次数是否达到预设次数。
[0113] 奇异性分析结果输出模块206，用于根据所述奇异性指数余维数栅格数据及每次迭代获得的所述迭代后的奇异性指数余维数栅格数据，输出所述迭代分析结果。
[0114] 不失一般性，奇异性分析结果可以是以奇异性指数余维数栅格数据的形式输出，具体地，根据所述奇异性指数余维数栅格数据及每次迭代获得的所述迭代后的奇异性指数余维数栅格数据，利用公式：
[0115]
[0116] 获得最终输出的奇异性分析结果，其中，Δα*(i,j)为最终输出奇异性指数余维数栅格数据中行索引为i，列索引为j位置的数值，n为迭代总次数，Δαk(i,j)为第k次迭代后的奇异性指数余维数栅格数据中行索引为i，列索引为j位置的数值，当k＝0时，Δαk(i,j)代表局部奇异性计算模块204根据栅格数据获得的所述奇异性指数余维数栅格数据中行索引为i，列索引为j的栅格值。当n为0时，意味着局部奇异性迭代模型退化为常规的非迭代模型。
[0117] 不失一般性，奇异性分析结果可以是以奇异性指数栅格数据的形式输出。利用一下公式：
[0118]
[0119] 获得最终输出的奇异性分析结果，其中，α*(i,j)为最终输出的奇异性分析结果即最终奇异性指数栅格数据中行索引为i，列索引为j位置的数值；E代表栅格数据的维数，在本实施例中采用的是二维的栅格数据，因此E取值为2；Δαk(i,j)为第k次迭代后的奇异性指数余维数栅格数据中行索引为i，列索引为j位置的数值，当k＝0时，Δαk(i,j)代表局部奇异性计算模块204根据栅格数据获得的所述奇异性指数余维数栅格数据中行索引为i，列索引为j的栅格值。
[0120] 第二实施例
[0121] 请参照图10，图10为本发明较佳实施例提供的一种栅格数据局部奇异性迭代分析方法的流程图。所述方法包括：
[0122] 步骤S101，预设参数。其中预设参数包括预设方形窗口的个数、每个预设方形窗口的尺度参数、预设方形窗口的正规化尺度参数、迭代精度及预设次数等。
[0123] 在本发明实施例中，步骤S101可以由预设参数模块201执行。
[0124] 步骤S102，获得栅格数据。
[0125] 在本发明实施例中，步骤S102可以由栅格数据输入模块202执行。具体地，可以是获取原始数据；判断所述原始数据是否为栅格格式数据；当所述原始数据不是栅格格式，根据所述原始数据，利用空间插值模型生成所述栅格数据，其中，所述空间插值模型可以但不限于是距离平方反比插值模型、克立格空间插值模型等；当原始数据本身就是栅格数据时，需要继续判断其是否存在缺失数据，并在存在数据缺失的时候，采用距离平方反比插值、克立格等空间插值模型生成对缺失值进行估值，使得栅格数据各个位置都具有有效数值。对于连片较大范围的缺失数据，插值具有较大的不确定性，同样先进行插值，在完成局部奇异性分析后对所得计算结果中的原缺失值位置的结果设置为缺失值，这样便于计算过程简单化而不失一般性。
[0126] 步骤S103，根据所述栅格数据构置对应的积分图。
[0127] 在本发明实施例中，步骤S103可以由积分图构置模块203执行。如图11所示，步骤S103包括以下子步骤：
[0128] 子步骤S1031，创建初始增广积分图，其中，所述初始增广积分图为一空白矩阵，且所述初始增广积分图的行数比所述栅格数据多一行及所述初始增广积分图的列数比所述栅格数据多一列。
[0129] 在本发明实施例中，步骤S1031可以由增广积分图创建子模块2031执行。
[0130] 子步骤S1032，将所述初始增广积分图的首行及首列均预置数值为零的积分图值。
[0131] 在本发明实施例中，步骤S1032可以由增广积分图构置子模块2032执行。首行的行索引为0，首列的列索引为0。
[0132] 子步骤S1033，填充所述初始增广积分图。具体地，根据所述栅格数据，利用公式：
[0133] SAT(i,j)＝ρ(i-1,j-1)+SAT(i-1,j)+SAT(i,j-1)-SAT(i-1,j-1)，
[0134] 填充所述初始增广积分图除了0行及0列以外的其他单元的积分图值，以构置所述增广积分图，其中SAT(i,j)表示所述增广积分图行索引为i且列索引为j的积分图值，ρ(i-1,j-1)代表所述栅格数据中行索引为(i-1)且列索引为(j-1)的栅格值，SAT(i-1,j)代表所述增广积分图行索引为(i-1)且列索引为j的积分图值，SAT(i,j-1)代表所述增广积分图行索引为i且列索引为(j-1)的积分图值，SAT(i-1,j-1)代表所述增广积分图行索引为(i-1)且列索引为(j-1)的积分图值，i和j均为大于或等于1的整数。
[0135] 在本发明实施例中，步骤S1033可以由增广积分图构置子模块2032执行。
[0136] 步骤S104，根据所述积分图、多个尺度各不相同且同中心的预设方形窗口及预设的正规化尺度参数，建立局部奇异性模型，以获得局部系数栅格数据和奇异性指数余维数栅格数据。
[0137] 在本发明实施例中，步骤S104可以由局部奇异性计算模块204执行。如图12所示，步骤S104包括以下子步骤：
[0138] 子步骤S1041，根据所述预设方形窗口及所述积分图，计算每个所述预设方形窗口所覆盖的所述栅格数据的平均值，其中，所有的所述预设方形窗口在栅格数据对应的焦元位置相同，所述焦元为所述预设方形窗口中心对应的所述栅格数据的数据单元。
[0139] 本发明实施例中，子步骤S1041可以由窗口平均值计算子模块2041执行。具体为，根据所述预设方形窗口、所述预设方形窗口对应的所述焦元位置及积分图，利用公式：
[0140] SUM(i,j,w,h)＝SAT(r2,c2)+SAT(r1,c1)-SAT(r1,c2)-SAT(r2,c1)、
[0141]
[0142]
[0143] 及
[0144]
[0145] 计算所述预设方形窗口所覆盖的所述栅格数据的和，其中，SUM(i,j,w,h)代表宽为w且长为h的所述预设方形窗口在所对应的焦元为所述栅格数据中行索引为i且列索引为j的数据单元时所覆盖的所述栅格数据的和，r1及r2为在所述增广积分图获取积分图值的行索引，c1及c2为在所述增广积分图中获取积分图值的列索引，SAT(r2,c2)为所述增广积分图在行索引为r2列索引为c2处的积分图值，SAT(r1,c1)为所述增广积分图在行索引为r1列索引为c1处的积分图值，SAT(r1,c2)为所述增广积分图在行索引为r1列索引为c2处的积分图值，SAT(r2,c1)为所述增广积分图在行索引为r2列索引为c1处的积分图值；对预设方形窗口在增广积分图中对应的四个积分图值所涉及四个坐标，即行索引r1，r2和列索引c1，c2，按如下四种情形对索引进行修正：
[0146] ①行索引r1<0，则行索引修改为0，也即r1＝0；
[0147] ②列索引c1<0，则列索引修改为0，也即c1＝0；
[0148] ③行索引r2大于栅格数据总行数M，则行索引r2修改为M，也即r2＝M；
[0149] ④列索引c2大于栅格数据列数N，则列索引c2修改为N，也即c2＝N。
[0150] 统计所述预设计算方框内覆盖的所述栅格数据的个数，根据公式：
[0151] T(i,j,w,h)＝w×h，或
[0152] T(i,j,w,h)＝(r2-r1)×(c2-c1)
[0153] 其中，MEAN(i,j,w,h)代表宽为w且长为h的所述预设计算方框在所对应的焦元为所述栅格数据中行索引为i且列索引为j的数据单元时所覆盖的所述栅格数据的个数。当r1，r2∈[0,M]且c1，c2∈[0,N]时，使用公式T(i,j,w,h)＝w×h；当r1、r2、c1、c2中至少有一个被预处理修正过，则使用公式T(i,j,w,h)＝(r2-r1)×(c2-c1)。
[0154] 然后，根据所述预设方形窗口所覆盖的所述栅格数据的和及所述栅格数据的个数，计算所述预设方形窗口所覆盖的所述栅格数据的平均值。
[0155] 子步骤S1042，根据所述正规化尺度参数对所述预设方形窗口的尺度进行正规化处理，以获得每个所述预设方形窗口的正规化尺度大小。
[0156] 本发明实施例中，子步骤S1042可以由正规化处理子模块2042执行。
[0157] 子步骤S1043，根据所有的所述预设方形窗口对应的平均值及对应的所述正规化尺度大小，建立局部奇异性模型。具体为，采用双对数坐标系下的加权线性回归方法，以获得所述焦元位置所对应的局部奇异性指数余维数及局部系数，其中，所述加权线性回归方法中各个样本点对应的权重赋值为所述样本点对应的所述预设方形窗口的平均值的平方。
[0158] 本发明实施例中，子步骤S1043可以由奇异性模型拟合子模块2043执行。
[0159] 子步骤S1044，获得所述局部系数栅格数据及所述奇异性指数余维数栅格数据。具体地，所述预设方形窗口的中心同步遍历每一个所述栅格数据的数据单元，并重复根据所有的所述预设方形窗口对应的平均值及正规化尺度大小建立局部奇异性模型，以获得所述栅格数据中每一个数据单元对应的局部奇异性指数余维数及局部系数，进而获得所述局部系数栅格数据及所述奇异性指数余维数栅格数据。
[0160] 本发明实施例中，子步骤S1044可以由获得子模块2044执行。
[0161] 步骤S105，重复对最新获得的局部系数栅格数据进行积分图构置及所述局部奇异性模型建立的迭代处理，以获得迭代后的局部系数栅格数据及迭代后的奇异性指数余维数栅格数据，直至所述迭代后的奇异性指数余维数栅格数据满足预设条件或迭代次数达到预设次数。
[0162] 在本发明实施例中，步骤S105可以由迭代模块205执行。如图13所示，步骤S105包括以下子步骤：
[0163] 子步骤S1051，根据最新获得的所述局部系数栅格数据构置对应的积分图。
[0164] 在本发明实施例中，子步骤S1051可以由积分图构置子模块2051执行。
[0165] 子步骤S1052，根据所述最新获得的局部系数栅格数据对应的积分图及所述预设方形窗口及预设的正规化尺度参数，建立局部奇异性模型，以获得所述迭代后的局部系数栅格数据及迭代后的奇异性指数余维数栅格数据。
[0166] 在本发明实施例中，子步骤S1052可以由迭代子模块2052执行。
[0167] 步骤S1053，判断迭代后的所述奇异性指数余维数栅格数据是否满足所述预设条件或所述迭代次数是否达到预设次数。当迭代后的所述奇异性指数余维数栅格数据满足所述预设条件和/或所述迭代次数达到预设次数，结束迭代，流程进入S106。当迭代后的所述奇异性指数余维数栅格数据不满足所述预设条件且所述迭代次数未达到预设次数，则重复对迭代后最新获得的所述局部系数栅格数据进行迭代，即流程回到S1051，重复对最新获得的局部系数栅格数据进行迭代。
[0168] 在本发明实施例中，步骤S1053可以由迭代子模块2052执行。
[0169] 步骤S106，根据所述奇异性指数余维数栅格数据及每次迭代获得的所述迭代后的奇异性指数余维数栅格数据，输出所述迭代分析结果。
[0170] 在本发明实施例中，步骤S106可以由转换模块206执行。
[0171] 综上所述，本发明提供本发明提供一种栅格数据局部奇异性迭代分析方法及装置。其中，所述方法包括：获得栅格数据。根据所述栅格数据构置对应的积分图。根据所述积分图、多个尺度各不相同且同中心的预设方形窗口及预设的正规化尺度参数，建立局部奇异性模型，以获得局部系数栅格数据和奇异性指数余维数栅格数据。重复对最新获得的局部系数栅格数据进行积分图构置及所述局部奇异性模型建立的迭代处理，以获得迭代后的局部系数栅格数据及迭代后的奇异性指数余维数栅格数据，直至所述迭代后的奇异性指数余维数栅格数据满足预设条件或迭代次数达到预设次数。根据所述奇异性指数余维数栅格数据及每次迭代获得的所述迭代后的奇异性指数余维数栅格数据，输出所述迭代分析结果。以实现从海量地学数据中获得更准确的局部奇异性指数余维数，同时还提高了对海量的地学数据的分析速度。
[0172] 以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉
本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以
权利要求的保护范围为准。
