CN107103090B - 栅格数据局部奇异性迭代分析方法及装置 - Google Patents
栅格数据局部奇异性迭代分析方法及装置 Download PDFInfo
- Publication number
- CN107103090B CN107103090B CN201710310280.XA CN201710310280A CN107103090B CN 107103090 B CN107103090 B CN 107103090B CN 201710310280 A CN201710310280 A CN 201710310280A CN 107103090 B CN107103090 B CN 107103090B
- Authority
- CN
- China
- Prior art keywords
- data
- integral
- index
- preset
- augmented
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
- G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
- G06F16/29—Geographical information databases
Landscapes
- Engineering & Computer Science (AREA)
- Databases & Information Systems (AREA)
- Theoretical Computer Science (AREA)
- Remote Sensing (AREA)
- Data Mining & Analysis (AREA)
- Physics & Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Image Generation (AREA)
Abstract
本发明提供一种栅格数据局部奇异性迭代分析方法及装置，涉及地学数据分析领域。该方法包括：根据栅格数据构置对应的积分图；根据积分图、多个尺度各不相同且同中心的预设方形窗口及预设的正规化尺度参数，建立局部奇异性模型；重复对最新获得的局部系数栅格数据进行积分图构置及局部奇异性模型建立的迭代处理，以获得迭代后的局部系数栅格数据及迭代后的奇异性指数余维数栅格数据，直至迭代后的奇异性指数余维数栅格数据满足预设条件或迭代次数达到预设次数；根据奇异性指数余维数栅格数据及每次迭代获得迭代后的奇异性指数余维数栅格数据，输出迭代分析结果。实现从海量地学数据中分析出的局部奇异性更准确，且计算时间更快速提高分析效率。
Description
技术领域
本发明涉及地学数据分析领域，具体而言，涉及栅格数据局部奇异性迭代分析方法及装置。
背景技术
地学数据是数据的一种类型，通常指与地球参考空间(二维或三维)位置有关的、表达与地球科学研究中各种实体和过程状态属性的数据。地学数据来源具有多源性，包括地质、地球物理、地球化学、遥感、地理等数据。从数据存储格式上，地学数据可以分为矢量地学数据和栅格地学数据两大类，矢量数据格式和栅格数据格式在地理信息系统支持下可以相互转换。对地学数据的分析是人类探索地球、认识地球的主要途径。通过对地学数据的分析，可以获得地质资源(例如，矿产资源)的分布情况。在资源告急的现如今，只有获得准确的地质资源分布，才能做到在保护地球的同时对资源进行合理的开发，才能保证可持续发展。
时空分布的不均匀性和多尺度性是地学数据的重要特征，分形/多重分形方法在地学数据分析得到广泛应用，“奇异性-广义自相似性-分形谱系”为核心的多重分形非线性矿产定量预测评价理论和方法体系，在地质找矿中取得了良好的应用效果。
局部奇异性分析方法是目前获得国内外高度肯定的并在实际地质找矿中取得了良好的应用效果的地学数据分析方法。其通过获取地学数据的奇异性指标，从而揭示矿产资源的内在关联性。需要说明的是，奇异性是指在很小的时间-空间范围内具有巨大能量释放或巨量物质形成的现象。多重分形性是用以描述非均匀分布物体分形维数的测度指标体系。具有多重分形分布特征的地学数据通常意味着具有局部奇异性特征。在开展局部奇异性分析时使用的地学数据应符合或经过适当的变换后符合其取值范围是非负的，具有可加性。
但是，发明人发现在大数据时代，地质勘查数据越来越丰富，同时找矿难度也随之越来越大。面对海量的地学数据，通过现有的局部奇异性分析方法获得的奇异性指标不够精准的问题越来越凸显。同时，现有的地学数据分析方法对海量的地学数据的处理速度过低。局部奇异性模型涉及到多个不同尺度窗口滑动平均值的计算，采用对窗口内数据累加求和再除以数据个数获得平均值这种常规计算方法，窗口数量越多，窗口越大，计算速度越慢。
发明内容
本发明的一个目的是提供一种栅格数据局部奇异性迭代分析方法，以解决上述问题。
本发明的另一个目的是提供一种栅格数据局部奇异性迭代分析装置，以解决上述问题。
为了实现上述目的，本发明实施例采用的技术方案如下：
本发明实施例提供一种栅格数据局部奇异性迭代分析方法，所述方法包括：获得栅格数据；根据所述栅格数据构置对应的积分图；根据所述积分图、多个尺度各不相同且同中心的预设方形窗口及预设的正规化尺度参数，建立局部奇异性模型，以获得局部系数栅格数据和奇异性指数余维数栅格数据；重复对最新获得的局部系数栅格数据进行积分图构置及所述局部奇异性模型建立的迭代处理，以获得迭代后的局部系数栅格数据及迭代后的奇异性指数余维数栅格数据，直至所述迭代后的奇异性指数余维数栅格数据满足预设条件或迭代次数达到预设次数；根据所述奇异性指数余维数栅格数据及每次迭代获得的所述迭代后的奇异性指数余维数栅格数据，输出所述迭代分析结果。
本发明实施例还提供一种栅格数据局部奇异性迭代分析装置，所述装置包括：栅格数据输入模块，用于获得栅格数据；积分图构置模块，用于根据所述栅格数据构置对应的积分图；局部奇异性计算模块，用于根据所述积分图、多个尺度各不相同且同中心的预设方形窗口及预设的正规化尺度参数，建立局部奇异性模型，以获得局部系数栅格数据和奇异性指数余维数栅格数据；迭代模块，用于重复对最新获得的局部系数栅格数据进行积分图构置及所述局部奇异性模型建立的迭代处理，以获得迭代后的局部系数栅格数据及迭代后的奇异性指数余维数栅格数据，直至所述迭代后的奇异性指数余维数栅格数据满足预设条件或迭代次数达到预设次数；奇异性分析结果输出模块，用于根据所述奇异性指数余维数栅格数据及每次迭代获得的所述迭代后的奇异性指数余维数栅格数据，输出所述迭代分析结果。
与现有技术相比，本发明提供一种栅格数据局部奇异性迭代分析方法及装置。所述方法包括：根据所述栅格数据构置对应的积分图；根据所述积分图、多个尺度各不相同且同中心的预设方形窗口及预设的正规化尺度参数，建立局部奇异性模型，以获得局部系数栅格数据和奇异性指数余维数栅格数据；重复对最新获得的局部系数栅格数据进行积分图构置及所述局部奇异性模型建立的迭代处理，以获得迭代后的局部系数栅格数据及迭代后的奇异性指数余维数栅格数据，直至所述迭代后的奇异性指数余维数栅格数据满足预设条件或迭代次数达到预设次数。输出所述迭代分析结果。以实现从海量地学数据中分析出的局部奇异性更精准，同时还提高了对海量的地学数据的分析速度。
为使本发明的上述目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附附图，作详细说明如下。
附图说明
为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。
图1为本发明的应用环境示意图。
图2为图1中处理机的方框示意图。
图3为本发明较佳实施例提供的栅格数据局部奇异性迭代分析装置的功能模块示意图。
图4为图3中示出的积分图构置模块的功能子模块示意图。
图5为图3中示出的局部奇异性计算模块的功能子模块示意图。
图6为多个预设方形窗口在栅格数据上的覆盖示意图。
图7为栅格数据与栅格数据的增广积分图的对应关系示意图。
图8为预设方形窗口在栅格数据上出现越界时，预设方形窗口在增广积分图的对应覆盖示意图。
图9为图3中示出的迭代模块的功能子模块示意图。
图10为本发明较佳实施例提供的栅格数据局部奇异性迭代分析方法的流程图。
图11为图10中示出的步骤S103的子步骤流程图。
图12为图10中示出的步骤S104的子步骤流程图。
图13为图10中示出的步骤S105的子步骤流程图。
图标：10-并行处理机系统；11-主处理机；12-从处理机；100-处理机；101-存储器；102-存储控制器；103-处理器；104-外设接口；105-通信单元；200-栅格数据局部奇异性迭代分析装置；201-预设参数模块；202-栅格数据输入模块；203-积分图构置模块；2031-增广积分图创建子模块；2032-增广积分图构置子模块；204-局部奇异性计算模块；2041-窗口平均值计算子模块；2042-正规化处理子模块；2043-奇异性模型拟合子模块；2044-获得子模块；205-迭代模块；2051-积分图构置子模块；2052-迭代子模块；206-奇异性分析结果输出模块。
具体实施方式
下面将结合本发明实施例中附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。
本发明实施例所提供的栅格数据局部奇异性迭代分析方法及装置可应用于如图1所示的并行处理机系统10中。如图1所示，并行处理机系统10包括主处理机11及多个从处理机12。主处理机11与从处理机12通过总线连接，实现数据交互。
于本发明实施例中，主处理机11及从处理机12结构均相同。但在实现本发明提供的栅格数据局部奇异性迭代分析的过程中承担不同的工作。具体地，主处理机11用于接收数据，接收数据包括预设参数及栅格数据，其中预设参数包括多个尺度不同的预设方形窗口的参数、预设正规化尺度参数、迭代精度及预设次数等。需要说明的是，正规化尺度参数要大于栅格的分辨率，是正整数，也可以是更一般的正实数，优选实施例中，正规化尺度参数选用所有预设方形窗口尺度的几何平均值。预设方形窗口的长和宽均以栅格数据中数据单元的边长为其度量单位，且预设方形窗口的长和宽均奇数。数据单元为正方形。可选地，主处理机11根据不同的预设方形窗口分别激活对应数量的从处理机12。主处理机11与从处理机12的总个数与预设方形窗口数量相同。使主处理机11与每个从处理机12分别被分配一个尺度不同的预设方形窗口。需要说明的是，当预设方形窗口中存在尺度为1x1，则尺度为1x1的预设方形窗口分配给主处理机11，其他预设方形窗口分别分配给激活的从处理机12。
在其他实施例中，栅格数据局部奇异性迭代分析方法及装置还可应用于串行处理系统。预设方形窗口同心覆盖于栅格数据上，串行处理系统依次计算每个预设方形窗口对应的平均值，并根据正规化尺度参数完成栅格数据局部奇异性迭代分析。
图2示出本发明较佳实施例提供的处理机100的方框示意图。处理机100可根据用途的不同作为主处理机11及从处理机12。所述处理机100包括栅格数据局部奇异性迭代分析装置200、存储器101、存储控制器102、处理器103、外设接口104、通信单元105。
所述存储器101、存储控制器102、处理器103、外设接口104、通信单元105各元件相互之间直接或间接地电性连接，以实现数据的传输或交互。例如，这些元件相互之间可通过一条或多条通讯总线或信号线实现电性连接。所述栅格数据局部奇异性迭代分析装置200包括至少一个可以软件或固件(firmware)的形式存储于所述存储器101中或固化在所述处理机100的操作系统(operating system，OS)中的软件功能模块。所述处理器103用于执行存储器101中存储的可执行模块，例如所述栅格数据局部奇异性迭代分析装置200包括的软件功能模块或计算机程序。
其中，存储器101可以是，但不限于，随机存取存储器(Random Access Memory，RAM)，只读存储器(Read Only Memory，ROM)，可编程只读存储器(Programmable Read-OnlyMemory，PROM)，可擦除只读存储器(Erasable Programmable Read-Only Memory，EPROM)，电可擦除只读存储器(Electric Erasable Programmable Read-Only Memory，EEPROM)等。其中，存储器101用于存储程序，所述处理器103在接收到执行指令后，执行所述程序，所述处理器103以及其他可能的组件对存储器101的访问可在所述存储控制器102的控制下进行。
处理器103可能是一种集成电路芯片，具有信号的处理能力。上述的处理器103可以是通用处理器，包括中央处理器(Central Processing Unit，简称CPU)、网络处理器(Network Processor，简称NP)等；还可以是数字信号处理器(DSP)、专用集成电路(ASIC)、现场可编程门阵列(FPGA)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件。可以实现或者执行本发明实施例中的公开的各方法、步骤及逻辑框图。通用处理器可以是微处理器或者该处理器103也可以是任何常规的处理器103等。
所述外设接口104将各种输入/输出装置耦合至处理器103以及存储器101。在一些实施例中，外设接口104，处理器103以及存储控制器102可以在单个芯片中实现。在其他一些实例中，他们可以分别由独立的芯片实现。
通信单元105可以通过有线或者无线的方式与其他通讯终端建立通信连接，实现数据的交互。
第一实施例
请参考图3本发明较佳实施例提供的栅格数据局部奇异性迭代分析装置200的功能模块示意图。栅格数据局部奇异性迭代分析装置200包括：预设参数模块201、栅格数据输入模块202、积分图构置模块203、局部奇异性计算模块204、迭代模块205及奇异性分析结果输出模块206。
预设参数模块201用于预设参数。其中预设参数包括预设方形窗口的个数、每个预设方形窗口的尺度参数、预设方形窗口的正规化尺度参数、迭代精度及预设次数等。
栅格数据输入模块202用于获得栅格数据。栅格数据就是将空间划分成大小相等分布均匀、紧密相连的网格，每一个网格称为一个单元，并在各单元上赋予相应的属性值来表示实体的一种数据形式。栅格数据中每一个数据单元均有一个相对位置坐标，以其在栅格数据中所处的行索引和列索引来表示，遵循C语言、Python等流行编程语言约定，所述栅格数据的索引从0开始。如下表1所示：
表1：
需要说明的是，栅格数据是地学数据中最简单、最常用的、最直观的空间数据。对栅格数据可以很方便地进行邻域函数计算。方形窗口是最常用的邻域，方形窗口的宽和长(以单元边长为度量单位)通常设置为奇数，矩形窗口中心所对应的栅格数据单元即为焦元。需要指出的是，局部奇异性分析中，方形窗口中又以正方形窗口最为简便，它不涉及各向异性的问题。
本实施例中，栅格数据输入模块202用于首先获取原始数据。原始数据即为需要被分析的地学数据。再判断所述原始数据是否为栅格格式数据。当所述原始数据不是栅格格式，根据所述原始数据，利用空间插值模型生成所述栅格数据，其中，所述空间插值模型可以但不限于是距离平方反比插值模型(IDW)、克立格空间插值模型(Kriging)等。当原始数据本身就是栅格数据时，需要判断其是否存在缺失数据，并在存在数据缺失的时候，采用距离平方反比插值、克立格等空间插值模型生成对缺失值进行估值，使得栅格数据各个位置都具有有效数值。对于连片较大范围的缺失数据，插值具有较大的不确定性，同样先进行插值，在完成局部奇异性分析后对所得计算结果中的原缺失值位置的结果设置为缺失值，这样便于计算过程简单化而不失一般性。
主处理机11的积分图构置模块203，用于根据所述栅格数据构置对应的积分图。积分图的格式与栅格数据相同。本实施例中，积分图推荐使用增广积分图。如图4所示，积分图构置模块203包括增广积分图创建子模块2031及增广积分图构置子模块2032。
增广积分图创建子模块2031，用于创建初始增广积分图，其中，所述初始增广积分图为一空白矩阵，且所述初始增广积分图的行数比所述栅格数据多一行及所述初始增广积分图的列数比所述栅格数据多一列。
增广积分图构置子模块2032，用于将所述初始增广积分图的首行及首列均预置数值为零的积分图值。其中，所述初始增广积分图的首行的行索引为0，所述初始增广积分图的首列的列索引为0。
增广积分图构置子模块2032，还用于根据所述栅格数据，利用公式：
SAT(i,j)＝ρ(i-1,j-1)+SAT(i-1,j)+SAT(i,j-1)-SAT(i-1,j-1)，
填充所述初始增广积分图除了首行及首列以外的其他单元的积分图值，以构置所述增广积分图，其中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的整数。当i＝1或j＝1时，SAT(i-1,j-1)代表首行或首列中的积分图值，例如，SAT(0,j)代表首行中列索引为j的积分图值，SAT(i,0)代表首列中行索引为i的积分图值。增广积分图构置时，行索引i的变化范围为1到M，列索引j的变化范围为1到N，行索引和列索引的索引增量为1，其中M和N分别为原始栅格数据的总行数和总列数。遍历所有索引，按前述公式进行计算完成增广积分图的构置。需要说明的是，增广积分图构置时，当i＝1或j＝1时，行索引为(i-1)、列索引为j的积分图值SAT(i-1,j)、行索引为i、列索引为(j-1)的积分图值SAT(i,j-1)正好就是积分图创建时预置在矩阵0行、0列的0值。
需要说明的是，公式：
与前述增广积分公式等效，0行及0列的值预置为0，i及j均为大于等于1的整数，j′及i′均为大于等于0的整数。其区别在于前述公式利用前面步骤所得的积分图值直接代入到后继计算中使得增广积分图的构置速度更快，减少运算量，提高整个栅格数据局部奇异性迭代分析的效率，其中SAT(i,j)表示置于所述初始增广积分图行索引为i列索引为j列的积分图值，ρ(i′,j′)代表所述栅格数据中行索引为i′列索引为j′的数据，
代表对所有满足条件i′<i且j′<j的ρ(i′,j′)求和，例如，SAT(1,1)＝ρ(0,0)，SAT(2,2)＝ρ(0,0)+ρ(0,1)+ρ(1,0)+ρ(1,1)，SAT(3,2)＝ρ(0,0)+ρ(0,1)+ρ(1,0)+ρ(1,1)+ρ(2,0)+ρ(2,1)。
以下表2为例，表2为积分图构置模块203根据表1所示的栅格数据构置的增广积分图：
表2
局部奇异性计算模块204，用于根据所述积分图、多个尺度各不相同且同中心的预设方形窗口及预设的正规化尺度参数，建立局部奇异性模型，以获得局部系数栅格数据和奇异性指数余维数栅格数据。本实施例中，如图5所示，局部奇异性计算模块204包括窗口平均值计算子模块2041、正规化处理子模块2042、奇异性模型拟合子模块2043及获得子模块2044。
窗口平均值计算子模块2041，用于根据所述预设方形窗口及所述积分图，计算每个所述预设方形窗口所覆盖的所述栅格数据的平均值，其中，所有的所述预设方形窗口在栅格数据对应的焦元位置相同，所述焦元为所述预设方形窗口中心对应的所述栅格数据的数据单元。其中，主处理机11将构置的增广积分图拷贝给被激活的每个从处理机12。被激活的每个从处理机12分别根据预先分配得到的预设方形窗口尺度及复制所得的增广积分图计算对应的平均值。需要说明的是，主处理机11及被激活的每个从处理机12对应的不同尺度的预设方形窗口在栅格数据对应的焦元位置相同，所述焦元为所述预设方形窗口中心对应的所述栅格数据的数据单元。这就使如果把不同的预设方形窗口放置在同一份栅格数据上时，预设方形窗口呈同心嵌套。需要说明的是，当预设计算方框中存在尺度为1x1这种极端情形时，则尺度为1x1的预设计算方框直接分配给主处理机11，其所需计算的平均值就是当前焦元位置的栅格数据值。以图6为例，图6所示4个不同尺度的预设方形窗口覆盖于表1所示的栅格数据上，4个预设方形窗口同心嵌套且分别覆盖了栅格数据中的1个数据，9个数据、25个数据及49个数据，四个预设方形窗口的中心在栅格数据对应的焦元的行索引为3，列索引为3。在本实施例中，计算每个所述预设方形窗口对应的平均值时，主处理机11与被激活的从处理机12同时计算，提高速度。
在其他实施例中，即当运用于串行处理系统时，则依次计算每个预设方框对应的平均值。
在本实施例中，窗口平均值计算子模块2041执行所述计算每个所述预设方形窗口所覆盖的所述栅格数据的平均值的方式包括：
根据所述预设方形窗口、所述预设方形窗口对应的所述焦元位置及增广积分图，利用公式：
SUM(i,j,w,h)＝SAT(r2,c2)+SAT(r1,c1)-SAT(r1,c2)-SAT(r2,c1)、
计算所述预设方形窗口所覆盖的所述栅格数据的和，其中，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)则为预设方形窗口在增广积分图中对应的四个积分图值。
以图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。
预设方形窗口在栅格数据图中进行滑动计算，预设方形窗口对应的焦元靠近栅格数据边缘位置时，会出现预设方形窗口越出栅格数据坐标范围的情况，预设方形窗口实际上只有其中一部分覆盖有栅格数据。预设方形窗口在增广积分图中对应的四个积分图值中至少一个积分图值的索引(包括行索引和/或列索引)出现越界，进而无法根据索引从增广积分图中获得有效的积分图值。因此当预设方形窗口出现越界情况时，需要对与预设方形窗口对应的四个积分图值的索引进行预处理。预设方形窗口在增广积分图中对应的四个积分图值所涉及四个坐标，即行索引r1、r2和列索引c1、c2，需要按如下四种基本情形对索引进行修正：
①当计算窗口向上越界时，即行索引r1<0，则行索引修改为0，也即r1＝0，此时行索引为r1的积分图值就是增广积分图中第0行中预置的0值；
②当计算窗口向左越界时，即列索引c1<0，则列索引修改为0，也即c1＝0，此时列索引c1的积分图值就是增广积分图中第0列中预置的0值；
③当计算窗口向下越界时，即行索引r2大于栅格数据总行数M，则行索引r2修改为M，也即r2＝M，此时的积分图取值为SAT(M,j)，其中j为[0,N]区间的任意值，N为栅格数据总列数；
④当计算窗口向右越界时，即列索引c2大于栅格数据列数N，则列索引c2修改为N，也即c2＝N，此时的积分图取值为SAT(i,N)，其中i为[0,M]区间的任意值。
需要指出的是，上述四个坐标的预处理需要逐个判别并根据需要进行相应的修正，在全部完成之后，才能根据上述积分图公式计算获得所述预设计算方框实际覆盖数据的和值。
计算所述预设方形窗口所覆盖的所述栅格数据的平均值，计算方法为根据所述预设方形窗口所覆盖的所述栅格数据的和除以所述栅格数据的个数，具体为公式：
其中，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)时，预设方形窗口内的覆盖栅格数据的个数。统计所述预设方形窗口内覆盖栅格数据的个数为两种情形：
(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个。需要指出的是，此种情形下对某个大小的窗口内数据个数是个固定值，其个数统计只需预先计算一次，以后该尺度预设方形窗口的平均值计算都可以直接可以使用该结果，无需每次计算平均值时重复统计个数。
(2)当预设方形窗口越出栅格数据坐标范围时，也即行索引r1、r2和列索引c1、c2在预处理中有修正索引的情形，需要统计预设方形窗口内覆盖的栅格数据的有效个数。覆盖个数不仅与预设方形窗口的尺度有关，还与窗口位置有关。具体为，根据公式：
T(i,j,w,h)＝(r2-r1)×(c2-c1)，
其中，r1和r2为宽为w长为h的预设方形窗口中心对应焦元的索引为(i,j)时，对应的四个积分图值经过预处理修正后的行索引，c1和c2为经过预处理修正后的列索引。需要指出的是，上述公式是通用的，对预设方形窗口完全覆盖于栅格数据情形也适用，当然用此公式相比之下会增加计算量。
以图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。
再例如，图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。
需要指出的是，当预设方框覆盖的栅格单元数等于1，也即窗口的宽为1且长为1时，此时平均值即为原来输入栅格数据对应焦元位置的值，只有当预设方框覆盖的栅格单元数大于1时，计算平均值时才需要借助积分图。
在计算预设方形窗口所覆盖的栅格数据的和时，无论窗口是大是小，利用积分图都归结为四个积分图值的加减法，比起直接对窗口内的数值简单相加求和，大大减少了运算次数，特别是窗口越大，计算效率更明显。
正规化处理子模块2042，用于根据所述正规化尺度参数对所述预设方形窗口的尺度进行正规化处理，以获得每个所述预设方形窗口的正规化尺度大小。具体为，根据所述正规化尺度参数，利用公式：
st＝lt/L，
获得每个所述预设方形窗口的正规化尺度大小。其中，t的取值可以是小于或等于m，且大于0的整数，m为预设方形窗口数量，st代表第t个预设方形窗口的正规化后的尺度大小，lt代表第t个预设方形窗口的尺度，L代表正规化尺度参数。约定lt大小为所述预设方形窗口面积的开方，也即
wt代表第t个预设方形窗口的宽，ht代表第t个预设方形窗口的长。对正方形窗口，lt即为所述预设方框的边长，也即lt＝wt＝ht。需要指出的是，在迭代模型中L应大于1，否则迭代模型不收敛不能得到预期的精度要求；对于常规的非迭代模型，L可以等于1，从而得到更为简化的局部奇异性模型。L的取值大小不影响局部奇异性指数余维数的期望值，但影响局部系数的估计期望值以及局部系数总体个别值预测区间。当正规化尺度取值为所有预设计算方框的尺寸的几何平均值时，这使得局部系数在迭代过程中产生的估计误差最小，可以减少迭代计算中的不确定性的传递。
奇异性模型拟合子模块2043，用于根据所有的所述预设方形窗口对应的平均值及对应的所述正规化尺度大小，采用双对数坐标系下的加权线性回归方法，以获得所述焦元位置所对应的局部奇异性指数余维数及局部系数，其中，所述加权线性回归方法中各个样本点对应的权重赋值为所述样本点对应的所述预设方形窗口的平均值的平方。具体为对任一焦元(i,j)位置，多尺度下的均值及对应的尺度存在如下关系：
MEAN(i,j,wt,ht)＝c(i,j)st -Δα(i,j)，
其中，m为预设方形窗口数量。MEAN(i,j,wt,ht)为所述位于(i,j)第t个宽为wt长为ht的预设方形窗口对应的平均值，st表示第t个预设方形窗口对应的正规化尺度大小。Δa(i,j)和c(i,j)为待定参数，分别称为位于(i,j)的局部奇异性指数余维数和局部系数。
局部奇异性模型实质为非线性的幂律模型，可选的，可以通过幂律模型转换为双对数坐标系下的加权线性回归模型来获得上述公式中所述的Δa(i,j)和c(i,j)。幂律模型转换为双对数坐标系下的加权线性回归模型可以提高Δa(i,j)和c(i,j)的预测精度，同时又节省计算量，具体地，通过采用一种特殊的加权方式进行线性模型拟合来逼近原始幂律模型拟合的结果。回归模型形如：
yt(i，j)＝ln(c(i，j))+(-Δa(i，j))×xt，
其中，t＝1,2,3...,m。xt＝ln(st)，yt＝ln(MEAN(i,j,wt,ht))，
(xt,yt)为第t个预设方形窗口对应的数据点。应用加权最小二乘法，使得二元函数
其中Wt为第t个预设方形窗口对应的数据点(xt,yt)的权重。对焦元在(i,j)位置处的加权回归，为达到逼近原始幂律模型拟合效果。具体地，利用公式：
Wt(i,j)＝[MEAN(i,j,wt,ht)]2，
获取每个预设方形窗口对应的权重值。其中，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。可以采用残差平方和的大小来判定相对优劣，残差平方和越小，说明拟合效果越好。残差平方和的计算公式为：
其中，SSE(i,j)为预设方形窗口对应焦元的索引为(i,j)时对应的残差平方和。
以图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)。
用幂律模型直接进行非线性拟合时，在i＝3，j＝3位置处，局部奇异性指数余维数Δα(3,3)＝1.0699，局部系数c(3,3)＝10.4569，残差平方和SSE(3,3)＝2.3232。
用双对数坐标系下采用常规线性回归方法时，可得，Δα(3,3)＝0.9945，c(3,3)＝10.8024，SSE(3,3)＝5.5322。
用双对数坐标系下采用加权线性回归方法时，图6中四个尺度由小到达预设方形窗口对应权重赋值为分别为1156.0、87.0489、47.3344、23.2324，可得，Δα(3,3)＝1.0548，c(3,3)＝10.6337，SSE(3,3)＝2.8354。
从预设方形窗口对应焦元索引为(3,3)的局部奇异性分析结果表明，该位置处具有明显的正奇异性，随着尺度由大变小，局部邻域中的平均值(例如代表某种物质的密度或某种地质过程的能量)越来越高，而且呈幂率增长。加权线性回归方法所得的SSE比常规线性回归方法所得的SSE更小，更加接近幂率模型非线性拟合所得的SSE。但是非线性拟合法需要预先提供相对合理的初始值，通过迭代来不断改进估值直至获得最优解，计算比较耗时。因此，本发明中提出使用双对数坐标系下加权线性回归方法来高效高精度地获得局部奇异性指数余维数及局部系数的估计值。
获得子模块2044，用于所述预设方形窗口的中心同步遍历每一个所述栅格数据的数据单元，并重复根据所有的所述预设方形窗口对应的平均值，与正规化尺度大小建立局部奇异性模型，以获得所述栅格数据中每一个数据单元对应的局部奇异性指数余维数及局部系数，进而获得所述局部系数栅格数据及所述奇异性指数余维数栅格数据。
迭代模块205，用于重复对最新获得的局部系数栅格数据进行积分图构置及所述局部奇异性模型建立的迭代处理，以获得迭代后的局部系数栅格数据及迭代后的奇异性指数余维数栅格数据，直至所述迭代后的奇异性指数余维数栅格数据满足预设条件或迭代次数达到预设次数。需要说明的是，预设条件为迭代后的奇异性指数余维数栅格数据的所有单元的数据的平方之和小于预设迭代精度。
在本实施例中，如图9所示，迭代模块205包括以下功能子模块：
积分图构置子模块2051，用于根据最新获得的所述局部系数栅格数据构置对应的积分图。原理同积分图构置模块203相同，在此不再赘述。
迭代子模块2052，用于根据所述最新获得的局部系数栅格数据对应的积分图、所述预设方形窗口及预设的正规化尺度参数，建立局部奇异性模型，以获得所述迭代后的局部系数栅格数据及迭代后的奇异性指数余维数栅格数据。
迭代子模块2052，还用于判断迭代后的所述奇异性指数余维数栅格数据是否满足所述预设条件或所述迭代次数是否达到预设次数。
奇异性分析结果输出模块206，用于根据所述奇异性指数余维数栅格数据及每次迭代获得的所述迭代后的奇异性指数余维数栅格数据，输出所述迭代分析结果。
不失一般性，奇异性分析结果可以是以奇异性指数余维数栅格数据的形式输出，具体地，根据所述奇异性指数余维数栅格数据及每次迭代获得的所述迭代后的奇异性指数余维数栅格数据，利用公式：
获得最终输出的奇异性分析结果，其中，Δα*(i,j)为最终输出奇异性指数余维数栅格数据中行索引为i，列索引为j位置的数值，n为迭代总次数，Δαk(i,j)为第k次迭代后的奇异性指数余维数栅格数据中行索引为i，列索引为j位置的数值，当k＝0时，Δαk(i,j)代表局部奇异性计算模块204根据栅格数据获得的所述奇异性指数余维数栅格数据中行索引为i，列索引为j的栅格值。当n为0时，意味着局部奇异性迭代模型退化为常规的非迭代模型。
不失一般性，奇异性分析结果可以是以奇异性指数栅格数据的形式输出。利用一下公式：
获得最终输出的奇异性分析结果，其中，α*(i,j)为最终输出的奇异性分析结果即最终奇异性指数栅格数据中行索引为i，列索引为j位置的数值；E代表栅格数据的维数，在本实施例中采用的是二维的栅格数据，因此E取值为2；Δαk(i,j)为第k次迭代后的奇异性指数余维数栅格数据中行索引为i，列索引为j位置的数值，当k＝0时，Δαk(i,j)代表局部奇异性计算模块204根据栅格数据获得的所述奇异性指数余维数栅格数据中行索引为i，列索引为j的栅格值。
第二实施例
请参照图10，图10为本发明较佳实施例提供的一种栅格数据局部奇异性迭代分析方法的流程图。所述方法包括：
步骤S101，预设参数。其中预设参数包括预设方形窗口的个数、每个预设方形窗口的尺度参数、预设方形窗口的正规化尺度参数、迭代精度及预设次数等。
在本发明实施例中，步骤S101可以由预设参数模块201执行。
步骤S102，获得栅格数据。
在本发明实施例中，步骤S102可以由栅格数据输入模块202执行。具体地，可以是获取原始数据；判断所述原始数据是否为栅格格式数据；当所述原始数据不是栅格格式，根据所述原始数据，利用空间插值模型生成所述栅格数据，其中，所述空间插值模型可以但不限于是距离平方反比插值模型、克立格空间插值模型等；当原始数据本身就是栅格数据时，需要继续判断其是否存在缺失数据，并在存在数据缺失的时候，采用距离平方反比插值、克立格等空间插值模型生成对缺失值进行估值，使得栅格数据各个位置都具有有效数值。对于连片较大范围的缺失数据，插值具有较大的不确定性，同样先进行插值，在完成局部奇异性分析后对所得计算结果中的原缺失值位置的结果设置为缺失值，这样便于计算过程简单化而不失一般性。
步骤S103，根据所述栅格数据构置对应的积分图。
在本发明实施例中，步骤S103可以由积分图构置模块203执行。如图11所示，步骤S103包括以下子步骤：
子步骤S1031，创建初始增广积分图，其中，所述初始增广积分图为一空白矩阵，且所述初始增广积分图的行数比所述栅格数据多一行及所述初始增广积分图的列数比所述栅格数据多一列。
在本发明实施例中，步骤S1031可以由增广积分图创建子模块2031执行。
子步骤S1032，将所述初始增广积分图的首行及首列均预置数值为零的积分图值。
在本发明实施例中，步骤S1032可以由增广积分图构置子模块2032执行。首行的行索引为0，首列的列索引为0。
子步骤S1033，填充所述初始增广积分图。具体地，根据所述栅格数据，利用公式：
SAT(i,j)＝ρ(i-1,j-1)+SAT(i-1,j)+SAT(i,j-1)-SAT(i-1,j-1)，
填充所述初始增广积分图除了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的整数。
在本发明实施例中，步骤S1033可以由增广积分图构置子模块2032执行。
步骤S104，根据所述积分图、多个尺度各不相同且同中心的预设方形窗口及预设的正规化尺度参数，建立局部奇异性模型，以获得局部系数栅格数据和奇异性指数余维数栅格数据。
在本发明实施例中，步骤S104可以由局部奇异性计算模块204执行。如图12所示，步骤S104包括以下子步骤：
子步骤S1041，根据所述预设方形窗口及所述积分图，计算每个所述预设方形窗口所覆盖的所述栅格数据的平均值，其中，所有的所述预设方形窗口在栅格数据对应的焦元位置相同，所述焦元为所述预设方形窗口中心对应的所述栅格数据的数据单元。
本发明实施例中，子步骤S1041可以由窗口平均值计算子模块2041执行。具体为，根据所述预设方形窗口、所述预设方形窗口对应的所述焦元位置及积分图，利用公式：
SUM(i,j,w,h)＝SAT(r2,c2)+SAT(r1,c1)-SAT(r1,c2)-SAT(r2,c1)、
计算所述预设方形窗口所覆盖的所述栅格数据的和，其中，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，按如下四种情形对索引进行修正：
①行索引r1<0，则行索引修改为0，也即r1＝0；
②列索引c1<0，则列索引修改为0，也即c1＝0；
③行索引r2大于栅格数据总行数M，则行索引r2修改为M，也即r2＝M；
④列索引c2大于栅格数据列数N，则列索引c2修改为N，也即c2＝N。
统计所述预设计算方框内覆盖的所述栅格数据的个数，根据公式：
T(i,j,w,h)＝w×h，或
T(i,j,w,h)＝(r2-r1)×(c2-c1)
其中，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)。
然后，根据所述预设方形窗口所覆盖的所述栅格数据的和及所述栅格数据的个数，计算所述预设方形窗口所覆盖的所述栅格数据的平均值。
子步骤S1042，根据所述正规化尺度参数对所述预设方形窗口的尺度进行正规化处理，以获得每个所述预设方形窗口的正规化尺度大小。
本发明实施例中，子步骤S1042可以由正规化处理子模块2042执行。
子步骤S1043，根据所有的所述预设方形窗口对应的平均值及对应的所述正规化尺度大小，建立局部奇异性模型。具体为，采用双对数坐标系下的加权线性回归方法，以获得所述焦元位置所对应的局部奇异性指数余维数及局部系数，其中，所述加权线性回归方法中各个样本点对应的权重赋值为所述样本点对应的所述预设方形窗口的平均值的平方。
本发明实施例中，子步骤S1043可以由奇异性模型拟合子模块2043执行。
子步骤S1044，获得所述局部系数栅格数据及所述奇异性指数余维数栅格数据。具体地，所述预设方形窗口的中心同步遍历每一个所述栅格数据的数据单元，并重复根据所有的所述预设方形窗口对应的平均值及正规化尺度大小建立局部奇异性模型，以获得所述栅格数据中每一个数据单元对应的局部奇异性指数余维数及局部系数，进而获得所述局部系数栅格数据及所述奇异性指数余维数栅格数据。
本发明实施例中，子步骤S1044可以由获得子模块2044执行。
步骤S105，重复对最新获得的局部系数栅格数据进行积分图构置及所述局部奇异性模型建立的迭代处理，以获得迭代后的局部系数栅格数据及迭代后的奇异性指数余维数栅格数据，直至所述迭代后的奇异性指数余维数栅格数据满足预设条件或迭代次数达到预设次数。
在本发明实施例中，步骤S105可以由迭代模块205执行。如图13所示，步骤S105包括以下子步骤：
子步骤S1051，根据最新获得的所述局部系数栅格数据构置对应的积分图。
在本发明实施例中，子步骤S1051可以由积分图构置子模块2051执行。
子步骤S1052，根据所述最新获得的局部系数栅格数据对应的积分图及所述预设方形窗口及预设的正规化尺度参数，建立局部奇异性模型，以获得所述迭代后的局部系数栅格数据及迭代后的奇异性指数余维数栅格数据。
在本发明实施例中，子步骤S1052可以由迭代子模块2052执行。
步骤S1053，判断迭代后的所述奇异性指数余维数栅格数据是否满足所述预设条件或所述迭代次数是否达到预设次数。当迭代后的所述奇异性指数余维数栅格数据满足所述预设条件和/或所述迭代次数达到预设次数，结束迭代，流程进入S106。当迭代后的所述奇异性指数余维数栅格数据不满足所述预设条件且所述迭代次数未达到预设次数，则重复对迭代后最新获得的所述局部系数栅格数据进行迭代，即流程回到S1051，重复对最新获得的局部系数栅格数据进行迭代。
在本发明实施例中，步骤S1053可以由迭代子模块2052执行。
步骤S106，根据所述奇异性指数余维数栅格数据及每次迭代获得的所述迭代后的奇异性指数余维数栅格数据，输出所述迭代分析结果。
在本发明实施例中，步骤S106可以由转换模块206执行。
综上所述，本发明提供本发明提供一种栅格数据局部奇异性迭代分析方法及装置。其中，所述方法包括：获得栅格数据。根据所述栅格数据构置对应的积分图。根据所述积分图、多个尺度各不相同且同中心的预设方形窗口及预设的正规化尺度参数，建立局部奇异性模型，以获得局部系数栅格数据和奇异性指数余维数栅格数据。重复对最新获得的局部系数栅格数据进行积分图构置及所述局部奇异性模型建立的迭代处理，以获得迭代后的局部系数栅格数据及迭代后的奇异性指数余维数栅格数据，直至所述迭代后的奇异性指数余维数栅格数据满足预设条件或迭代次数达到预设次数。根据所述奇异性指数余维数栅格数据及每次迭代获得的所述迭代后的奇异性指数余维数栅格数据，输出所述迭代分析结果。以实现从海量地学数据中获得更准确的局部奇异性指数余维数，同时还提高了对海量的地学数据的分析速度。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。
Claims (7)
1.一种栅格数据局部奇异性迭代分析方法，其特征在于，所述方法包括：
获得栅格数据；
根据所述栅格数据构置对应的积分图；其中，所述积分图为增广积分图；所述增广积分图的首行及首列均为零；
根据所述积分图、多个尺度各不相同且同中心的预设方形窗口及预设的正规化尺度参数，建立局部奇异性模型，以获得局部系数栅格数据和奇异性指数余维数栅格数据；具体包括：根据所述预设方形窗口、所述预设方形窗口对应的焦元位置及积分图，利用公式：
SUM(i,j,w,h)＝SAT(r2,c2)+SAT(r1,c1)-SAT(r1,c2)-SAT(r2,c1)、
计算所述预设方形窗口所覆盖的所述栅格数据的和，其中，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处的积分图值；统计所述预设方形窗口内覆盖的所述栅格数据的个数；根据所述预设方形窗口所覆盖的所述栅格数据的和及所述栅格数据的个数，计算所述预设方形窗口所覆盖的所述栅格数据的平均值；其中，所有的所述预设方形窗口在栅格数据对应的所述焦元位置相同，所述焦元为所述预设方形窗口中心对应的所述栅格数据的数据单元；根据所述正规化尺度参数对所述预设方形窗口的尺度进行正规化处理，以获得每个所述预设方形窗口的正规化尺度大小；根据所有的所述预设方形窗口对应的平均值及对应的所述正规化尺度大小，采用双对数坐标系下的加权线性回归方法，以获得所述焦元位置所对应的局部奇异性指数余维数及局部系数，其中，所述加权线性回归方法中各个样本点对应的权重赋值为所述样本点对应的所述预设方形窗口所覆盖的所述栅格数据的平均值的平方；具体地，利用公式：
Wt(i,j)＝[MEAN(i,j,wt,ht)]2，
获取每个预设方形窗口对应的权重值；其中，Wt(i,j)为第t个预设方形窗口在积分图中对应焦元为(i,j)时对应的数据点(xt,yt)的权重，MEAN(i,j,wt,ht)为第t个预设方形窗口在积分图中对应焦元为(i,j)时对应的平均值；所述预设方形窗口的中心同步遍历每一个所述栅格数据的数据单元，并重复根据所有的所述预设方形窗口对应的平均值，计算所述栅格数据中每一个数据单元对应的局部奇异性指数余维数及局部系数，获得所述局部系数栅格数据及所述奇异性指数余维数栅格数据；在所述预设方形窗口遍历过程中出现越界时，对r1、r2、c1、c2中至少一个进行校正；
重复对最新获得的局部系数栅格数据进行积分图构置及所述局部奇异性模型建立的迭代处理，以获得迭代后的局部系数栅格数据及迭代后的奇异性指数余维数栅格数据，直至所述迭代后的奇异性指数余维数栅格数据满足预设条件或迭代次数达到预设次数；
根据所述奇异性指数余维数栅格数据及每次迭代获得的所述迭代后的奇异性指数余维数栅格数据，输出所述迭代分析结果。
2.如权利要求1所述的栅格数据局部奇异性迭代分析方法，其特征在于，所述根据所述栅格数据构置对应的积分图包括：
创建初始增广积分图，其中，所述初始增广积分图为一空白矩阵，且所述初始增广积分图的行数比所述栅格数据多一行及所述初始增广积分图的列数比所述栅格数据多一列；
将所述初始增广积分图的首行及首列均预置数值为零的积分图值；
根据所述栅格数据，利用公式：
SAT(i,j)＝ρ(i-1,j-1)+SAT(i-1,j)+SAT(i,j-1)-SAT(i-1,j-1)，
填充所述初始增广积分图，以构置所述增广积分图，其中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的整数。
3.如权利要求1所述的栅格数据局部奇异性迭代分析方法，其特征在于，所述重复根据最新获得的局部系数栅格数据进行所述积分图构置及所述局部奇异性模型建立的迭代处理，以获得迭代后的局部系数栅格数据及迭代后的奇异性指数余维数栅格数据，直至所述迭代后的奇异性指数余维数栅格数据满足预设条件或迭代次数达到预设次数的步骤包括：
根据最新获得的所述局部系数栅格数据构置对应的积分图；
根据所述最新获得的局部系数栅格数据对应的积分图、所述预设方形窗口及预设的正规化尺度参数，建立局部奇异性模型，以获得所述迭代后的局部系数栅格数据及迭代后的奇异性指数余维数栅格数据；
判断迭代后的所述奇异性指数余维数栅格数据是否满足所述预设条件或所述迭代次数是否达到预设次数；
当迭代后的所述奇异性指数余维数栅格数据满足所述预设条件或所述迭代次数达到预设次数，结束迭代；
当迭代后的所述奇异性指数余维数栅格数据不满足所述预设条件且所述迭代次数未达到预设次数，则重复对迭代后最新获得的所述局部系数栅格数据进行迭代。
4.如权利要求1所述的栅格数据局部奇异性迭代分析方法，其特征在于，所述获得栅格数据的步骤包括：
获取原始数据；
判断所述原始数据是否为栅格格式数据；
当所述原始数据不是栅格格式，根据所述原始数据，利用空间插值模型生成所述栅格数据，其中，所述空间插值模型包括距离平方反比插值模型或克立格空间插值模型。
5.一种栅格数据局部奇异性迭代分析装置，其特征在于，所述装置包括：
栅格数据输入模块，用于获得栅格数据；
积分图构置模块，用于根据所述栅格数据构置对应的积分图；其中，所述积分图为增广积分图；所述增广积分图的首行及首列均为零；
局部奇异性计算模块，用于根据所述积分图、多个尺度各不相同且同中心的预设方形窗口及预设的正规化尺度参数，建立局部奇异性模型，以获得局部系数栅格数据和奇异性指数余维数栅格数据；其中，所述局部奇异性计算模块包括窗口平均值计算子模块，用于根据所述预设方形窗口、所述预设方形窗口对应的焦元位置及积分图，利用公式：
SUM(i,j,w,h)＝SAT(r2,c2)+SAT(r1,c1)-SAT(r1,c2)-SAT(r2,c1)、
计算所述预设方形窗口所覆盖的所述栅格数据的和，其中，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处的积分图值；统计所述预设方形窗口内覆盖的所述栅格数据的个数；根据所述预设方形窗口所覆盖的所述栅格数据的和及所述栅格数据的个数，计算所述预设方形窗口所覆盖的所述栅格数据的平均值；其中，所有的所述预设方形窗口在栅格数据对应的所述焦元位置相同，所述焦元为所述预设方形窗口中心对应的所述栅格数据的数据单元；正规化处理子模块，用于根据所述正规化尺度参数对所述预设方形窗口的尺度进行正规化处理，以获得每个所述预设方形窗口的正规化尺度大小；奇异性模型拟合子模块，用于根据所有的所述预设方形窗口对应的平均值及对应的所述正规化尺度大小，采用双对数坐标系下的加权线性回归方法，以获得所述焦元位置所对应的局部奇异性指数余维数及局部系数，其中，所述加权线性回归方法中各个样本点对应的权重赋值为所述样本点对应的所述预设方形窗口的平均值的平方；具体地，利用公式：
Wt(i,j)＝[MEAN(i,j,wt,ht)]2，
获取每个预设方形窗口对应的权重值；其中，Wt(i,j)为第t个预设方形窗口在积分图中对应焦元为(i,j)时对应的数据点(xt,yt)的权重，MEAN(i,j,wt,ht)为第t个预设方形窗口在积分图中对应焦元为(i,j)时对应的平均值；
获得子模块，用于所述预设方形窗口的中心同步遍历每一个所述栅格数据的数据单元，并重复根据所有的所述预设方形窗口对应的平均值，计算所述栅格数据中每一个数据单元对应的局部奇异性指数余维数及局部系数，获得所述局部系数栅格数据及所述奇异性指数余维数栅格数据；当所述预设方形窗口遍历过程中出现越界时，对r1、r2、c1、c2中至少一个进行校正；
迭代模块，用于重复对最新获得的局部系数栅格数据进行积分图构置及所述局部奇异性模型建立的迭代处理，以获得迭代后的局部系数栅格数据及迭代后的奇异性指数余维数栅格数据，直至所述迭代后的奇异性指数余维数栅格数据满足预设条件或迭代次数达到预设次数；
奇异性分析结果输出模块，用于根据所述奇异性指数余维数栅格数据及每次迭代获得的所述迭代后的奇异性指数余维数栅格数据，输出所述迭代分析结果。
6.如权利要求5所述的栅格数据局部奇异性迭代分析装置，其特征在于，所述积分图构置模块包括：
增广积分图创建子模块，用于创建初始增广积分图，其中，所述初始增广积分图为一空白矩阵，且所述初始增广积分图的行数比所述栅格数据多一行及所述初始增广积分图的列数比所述栅格数据多一列；
增广积分图构置子模块，用于将所述初始增广积分图的首行及首列均预置数值为零的积分图值；及还用于根据所述栅格数据，利用公式：
SAT(i,j)＝ρ(i-1,j-1)+SAT(i-1,j)+SAT(i,j-1)-SAT(i-1,j-1)，
填充所述初始增广积分图，以构置所述增广积分图，其中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的整数。
7.如权利要求5所述的栅格数据局部奇异性迭代分析装置，其特征在于，所述迭代模块包括：
积分图构置子模块，用于根据最新获得的所述局部系数栅格数据构置对应的积分图；
迭代子模块，用于根据所述最新获得的局部系数栅格数据对应的积分图、所述预设方形窗口及预设的正规化尺度参数，建立局部奇异性模型，以获得所述迭代后的局部系数栅格数据及迭代后的奇异性指数余维数栅格数据；
迭代子模块，还用于判断迭代后的所述奇异性指数余维数栅格数据是否满足所述预设条件或所述迭代次数是否达到预设次数。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710310280.XA CN107103090B (zh) | 2017-05-05 | 2017-05-05 | 栅格数据局部奇异性迭代分析方法及装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710310280.XA CN107103090B (zh) | 2017-05-05 | 2017-05-05 | 栅格数据局部奇异性迭代分析方法及装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107103090A CN107103090A (zh) | 2017-08-29 |
| CN107103090B true CN107103090B (zh) | 2020-04-10 |
Family
ID=59657624
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710310280.XA Active CN107103090B (zh) | 2017-05-05 | 2017-05-05 | 栅格数据局部奇异性迭代分析方法及装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107103090B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113065091B (zh) * | 2021-04-12 | 2021-11-05 | 中国地质科学院地质力学研究所 | 地学信息各向异性分布规律分析方法及装置、电子设备 |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101694671A (zh) * | 2009-10-27 | 2010-04-14 | 中国地质大学（武汉） | 一种基于地学栅格图像的空间加权主成分分析的方法 |
-
2017
- 2017-05-05 CN CN201710310280.XA patent/CN107103090B/zh active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101694671A (zh) * | 2009-10-27 | 2010-04-14 | 中国地质大学（武汉） | 一种基于地学栅格图像的空间加权主成分分析的方法 |
Non-Patent Citations (5)
| Title |
|---|
| A batch sliding window method for local singularity mapping and its application for geochemical anomaly identification;FanXiao 等;《Computers&Geosciences》;20151102;全文 * |
| 基于栅格数据模型的局部奇异性分析迭代算法;陈志军 等;《2011 年第十届全国数学地质与地学信息学术研讨会论文集》;20111125;全文 * |
| 基于栅格数据的局部奇异性分析迭代方法及其软件实现;陈志军 等;《地质学刊》;20141231;第38卷(第4期);第616-619页 * |
| 积分图及其应用;阿玛尼迪迪;《https://www.cnblogs.com/codingmengmeng/p/6567124.html》;20170317;第1、3页 * |
| 陈志军 等.基于栅格数据的局部奇异性分析迭代方法及其软件实现.《地质学刊》.2014,第38卷(第4期), * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107103090A (zh) | 2017-08-29 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| O’Brien et al. | A fast and objective multidimensional kernel density estimation method: fastKDE | |
| CN110232471B (zh) | 一种降水传感网节点布局优化方法及装置 | |
| US20170338802A1 (en) | Actually-measured marine environment data assimilation method based on sequence recursive filtering three-dimensional variation | |
| CN113837450B (zh) | 基于深度学习的河网密集流域水情趋势预测方法及其应用 | |
| KR20210027269A (ko) | 깊이 이미지 보완 방법 및 장치, 컴퓨터 판독 가능한 저장 매체 | |
| CN104123457B (zh) | 一种稳健的卫星遥感影像有理函数模型参数估计方法 | |
| CN102509304A (zh) | 基于智能优化的摄像机标定方法 | |
| CN111611541B (zh) | 基于Copula函数的无资料地区降水数据推求方法及系统 | |
| CN110929862B (zh) | 定点化的神经网络模型量化装置和方法 | |
| Yang et al. | Multiobjective sensitivity analysis and optimization of distributed hydrologic model MOBIDIC | |
| CN107103090B (zh) | 栅格数据局部奇异性迭代分析方法及装置 | |
| CN114822033A (zh) | 基于特征金字塔网络的路网交通流量数据修复方法及系统 | |
| CN117891962B (zh) | 一种城市分布式光伏系统数据的图数据库构建方法及应用 | |
| Palacios-Rodríguez et al. | Generalized Pareto processes for simulating space-time extreme events: an application to precipitation reanalyses | |
| De Hoop et al. | An exact redatuming procedure for the inverse boundary value problem for the wave equation | |
| CN110807428A (zh) | 煤类样品的识别方法、装置、服务器及存储介质 | |
| Tian et al. | Impacts of surface model generation approaches on raytracing-based solar potential estimation in urban areas | |
| CN117218551B (zh) | 一种基于误差分析的估测算法优化方法及装置 | |
| CN114998530B (zh) | 一种基于实景三维地形的水体监测方法和装置 | |
| CN111199092A (zh) | 太阳辐射遥感估算方法、系统及数据处理装置 | |
| Dutta et al. | Variogram calculations for random fields on regular lattices using quadrature methods | |
| Wang et al. | A geostatistical approach to the change-of-support problem and variable-support data fusion in spatial analysis | |
| Li et al. | Extraction of digital terrain model based on regular mesh generation in mountainous areas | |
| CN110794469A (zh) | 基于最小地质特征单元约束的重力反演方法 | |
| CN115951405B (zh) | 一种基于稀疏化卷积和自适应矩估计的速度变化成像方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |