CN102262789B - Random dot data gridding method by using Laplace equation - Google Patents

Random dot data gridding method by using Laplace equation Download PDF

Info

Publication number
CN102262789B
CN102262789B CN 201110184310 CN201110184310A CN102262789B CN 102262789 B CN102262789 B CN 102262789B CN 201110184310 CN201110184310 CN 201110184310 CN 201110184310 A CN201110184310 A CN 201110184310A CN 102262789 B CN102262789 B CN 102262789B
Authority
CN
China
Prior art keywords
gridding
data
value
node
point
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.)
Expired - Fee Related
Application number
CN 201110184310
Other languages
Chinese (zh)
Other versions
CN102262789A (en
Inventor
黄翼坚
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Northwest University
Original Assignee
Northwest University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Northwest University filed Critical Northwest University
Priority to CN 201110184310 priority Critical patent/CN102262789B/en
Publication of CN102262789A publication Critical patent/CN102262789A/en
Application granted granted Critical
Publication of CN102262789B publication Critical patent/CN102262789B/en
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Abstract

The invention discloses a random dot data gridding method by using a Laplace equation. To boundary nodes, the method provided by the invention utilizes an inverse distance to a power to compute a gridding value according to extracted partial random dot data. To adjacent internal codes comprising known random dot data, the method utilizes the inverse distance to a power to compute a gridding value according to the known data which is adjacent to the dot. To adjacent internal codes comprising no known data, the method utilizes a finite difference approximate expression iteration of the Laplace equation to compute a gridding value. To grid nodes which are overlapped with the known random dot data, the method directly utilizes the known data as a gridding value.

Description

A kind of random point data gridding method that is used for geophysics field of utilizing Laplace's equation
Technical field
The invention belongs to random point data gridding technical field, relate to a kind of computer graphics, geology, geophysics and geochemical discrete data gridding technology, especially a kind of random point data gridding method that is used for geophysics field of utilizing Laplace's equation.
Background technology
Before the computer graphics field is carried out isoplethes drawing and in geology, geophysics and geochemistry observation data is further processed, need to carry out gridding to the given data of stochastic distribution, namely calculate data value on the regular grid node according to known random point data value, this is a data Interpolation Process.
Gridding method commonly used has anti-distance weighted (Inverse distance to a power) method, Ke Lijin (Kriging) method of interpolation, minimum curvature (Minimum curvature) method etc.These methods have their own characteristics each.The extrapolability of anti-distance weighted method is not strong, the computing velocity of golden method of interpolation is relatively slow in the gram, the minimum-curvature method computing velocity is relatively very fast and the result is smooth, but in number of data points situation large or repeatedly gridding, its computing velocity still needs further to improve.The present invention is the improvement on the minimum-curvature method basis, and computing velocity is faster and result of calculation is smooth.
Summary of the invention
The object of the invention is to overcome the shortcoming of above-mentioned prior art, a kind of random point data gridding method that is used for geophysics field of utilizing Laplace's equation is provided, the method can solve from stochastic distribution but the given data value of non-overlapping copies calculates the data value problem on the regular grid node, and computing velocity is faster and result of calculation is smooth.
The objective of the invention is to solve by the following technical programs:
The present invention utilizes the random point data gridding method that is used for geophysics field of Laplace's equation, may further comprise the steps:
1) the gridding zone is expanded carried out gridding behind the border again and calculate, and two row and the two row nodes of ragged edge in all grid nodes are called boundary node, all the other nodes are called internal node;
2) to boundary node, when the number of all given datas greater than 100 the time, extract 100 data from all given datas, accordingly the data inverse distance weighted interpolation method computing net value of formatting; When the number of all given datas is less than or equal to 100, adopt the inverse distance weighted interpolation method computing net value of formatting according to all given datas;
The node that peripheral region in the internal node is included given data is called the Constrained point, and Constrained is put the bound data that given data that the peripheral region comprises is called this Constrained point; Above-described peripheral region refers to a rectangular region centered by this node, and this rectangular length and width are respectively half of Gridding length and width;
3) to Constrained point bound data employing inverse distance weighted interpolation method computing net value of formatting according to this point; The node that peripheral region in the internal node is not comprised given data is called without obligatory point; Above-described peripheral region refers to a rectangular region centered by this node, and this rectangular length and width are respectively half of Gridding length and width;
4) to having or not obligatory point give null value, then adopt the finite-difference approximation formula iterative computation gridding value of Laplace's equation;
5) to having or not the obligatory point moving window method of average calculate one by one last gridding value.
In the above step 1), the width that expands the border is 5 sizing grids, and the data in gridding zone before the border are not expanded in last only output.
Above step 2) in, adopt the formula of inverse distance weighted interpolation method computation bound node gridding value to be according to the part given data that extracts or whole given data:
Figure GDA00002840734500031
Wherein, u 1The gridding value of expression boundary node, N represents the number of used given data, k represents the sequence number of used given data, f kThe value that represents k used given data, d kRepresent that k used given data is to the distance of calculation level.
In the above step 3), adopt the computing formula of the inverse distance weighted interpolation method computing net value of formatting to be to Constrained point according to the bound data of this point:
Figure GDA00002840734500032
Wherein, u 2The gridding value of expression Constrained point, N represents total number of bound data, k represents the sequence number of bound data, f kThe value that represents k bound data, d kRepresent that k bound data is to the distance of calculation level.
In the above step 4), the expression formula of described Laplace's equation is:
∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 = 0 - - - ( 3 )
Wherein, u represents the gridding value function, and x and y represent two components of Cartesian coordinates; Utilize the differential in the Using Second-Order Central Difference replacement Laplace's equation (3), can obtain the gridding iterative computation formula without obligatory point
u j i = ( u j + 1 i + u j - 1 i ) dx + ( u j i + 1 + u j i - 1 ) dy 2 ( dx + dy ) - - - ( 4 )
Wherein, dx and dy represent respectively length and the width of grid, and i and j represent respectively the grid node sequence number on x axle and the z direction of principal axis,
Figure GDA00002840734500042
Expression is positioned on the x direction of principal axis gridding value of j Nodes on i node and the z direction of principal axis.
In the above step 5), the described employing moving window method of average calculate one by one have or not the computing formula of obligatory point gridding value to be:
w jj ii = 1 25 [ Σ i = ii - 2 ii + 2 Σ j = jj - 2 jj + 2 u j i ] - - - ( 5 )
Wherein, ii and jj represent respectively the grid node sequence number on x axle and the z direction of principal axis,
Figure GDA00002840734500044
Represent that last result of calculation is positioned on the x direction of principal axis gridding value of j Nodes on i node and the z direction of principal axis, i and j represent respectively the grid node sequence number on x axle and the z direction of principal axis,
Figure GDA00002840734500045
Representing the 4th) step calculates the gridding value that is positioned on the x direction of principal axis j Nodes on i node and the z direction of principal axis of gained.
The present invention has following beneficial effect:
(1) the present invention carries out gridding calculating after the border is expanded in the gridding zone again, and process on such expansion limit can reduce the border obtaining value method to inner gridding result's impact;
(2) the present invention utilizes at most 100 data to come the gridding value of computation bound node by the inverse distance weighted interpolation method, when the number of known random point data is larger, can determine rapidly the boundary node value;
(3) only according to the data inverse distance weighted interpolation method computing net value of formatting on little square region around this point, computing velocity is very fast to Constrained point in the present invention;
(4) the present invention is after obligatory point is given null value to having or not, adopt the finite-difference approximation formula iterative computation gridding value of two-dimentional Laplace's equation, on the one hand, the two dimension Laplace's equation has the Second Order Partial differential term to x and y coordinate, can make the gridding result smooth based on this equation, on the other hand, the two dimension Laplace's equation only has two partial differential items, form is extremely simple, and the finite-difference formula of setting up based on this equation is also very simple, can calculate rapidly the gridding value without obligatory point;
(5) the present invention to having or not obligatory point adopt the moving window method of average to calculate one by one last gridding value, can further increase the smoothness of result of calculation, and not change the data value with the given data coincide point.
Description of drawings
Fig. 1 is the flat distribution map of one group of actual measurement geophysics gravity anomaly random point data, the wherein position of "+" expression eyeball.
Fig. 2 is the flat distribution map of one group of regular grid node, wherein the position of "+" expression node.
The isogram of Fig. 3 for adopting the present invention that the represented actual measurement geophysics gravity anomaly data of Fig. 1 are drawn after by regular node gridding shown in Figure 2.
Embodiment
The present invention utilizes the random point data gridding method that is used for geophysics field of Laplace's equation, specifically may further comprise the steps:
1) the gridding zone is expanded carried out gridding behind the border again and calculate, and two row and the two row nodes of ragged edge in all grid nodes are called boundary node, all the other nodes are called internal node; In preferred embodiment of the present invention, the width that this step expands the border is 5 sizing grids, and the data in gridding zone before the border are not expanded in last only output.
2) to boundary node, when the number of all given datas greater than 100 the time, extract 100 data from all given datas, accordingly the data inverse distance weighted interpolation method computing net value of formatting; When the number of all given datas is less than or equal to 100, adopt the inverse distance weighted interpolation method computing net value of formatting according to all given datas;
The node that peripheral region in the internal node is included given data is called the Constrained point, and Constrained is put the bound data that given data that the peripheral region comprises is called this Constrained point; Above-described peripheral region refers to a rectangular region centered by this node, and this rectangular length and width are respectively half of Gridding length and width;
In this step, adopt the formula of inverse distance weighted interpolation method computation bound node gridding value to be according to the part given data that extracts or whole given data:
Figure GDA00002840734500061
Wherein, u 1The gridding value of expression boundary node, N represents the number of used given data, k represents the sequence number of used given data, f kThe value that represents k used given data, d kRepresent that k used given data is to the distance of calculation level.
3) to Constrained point bound data employing inverse distance weighted interpolation method computing net value of formatting according to this point; The node that peripheral region in the internal node is not comprised given data is called without obligatory point; Above-described peripheral region refers to a rectangular region centered by this node, and this rectangular length and width are respectively half of Gridding length and width;
Wherein, adopt the computing formula of the inverse distance weighted interpolation method computing net value of formatting to be to Constrained point according to the bound data of this point:
Figure GDA00002840734500071
Wherein, u 2The gridding value of expression Constrained point, N represents total number of bound data, k represents the sequence number of bound data, f kThe value that represents k bound data, d kRepresent that k bound data is to the distance of calculation level.
4) to having or not obligatory point give null value, then adopt the finite-difference approximation formula iterative computation gridding value of Laplace's equation; The expression formula of Laplace's equation is suc as formula (3):
∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 = 0 - - - ( 3 )
Wherein, u represents the gridding value function, and x and y represent two components of Cartesian coordinates; Utilize the differential in the Using Second-Order Central Difference replacement Laplace's equation (3), can obtain the gridding iterative computation formula without obligatory point
u j i = ( u j + 1 i + u j - 1 i ) dx + ( u j i + 1 + u j i - 1 ) dy 2 ( dx + dy ) - - - ( 4 )
Wherein, dx and dy represent respectively length and the width of grid, and i and j represent respectively the grid node sequence number on x axle and the z direction of principal axis,
Figure GDA00002840734500074
Expression is positioned on the x direction of principal axis gridding value of j Nodes on i node and the z direction of principal axis.
5) to having or not obligatory point adopt the moving window method of average to calculate one by one last gridding value, wherein adopt the moving window method of average calculate one by one have or not the computing formula of obligatory point gridding value to be:
w jj ii = 1 25 [ Σ i = ii - 2 ii + 2 Σ j = jj - 2 jj + 2 u j i ] - - - ( 5 )
Wherein, ii and jj represent respectively the grid node sequence number on x axle and the z direction of principal axis,
Figure GDA00002840734500082
Represent that last result of calculation is positioned on the x direction of principal axis gridding value of j Nodes on i node and the z direction of principal axis, i and j represent respectively the grid node sequence number on x axle and the z direction of principal axis,
Figure GDA00002840734500083
Representing the 4th) step calculates the gridding value that is positioned on the x direction of principal axis j Nodes on i node and the z direction of principal axis of gained.
Be applied to a specific embodiment of geophysics field below in conjunction with the present invention, the present invention done describing in further detail.
Fig. 1 is the flat distribution map of one group of actual measurement geophysics gravity anomaly data point, and wherein the position of "+" number expression observation data has 10773 data points.The minimum and maximum value of data is respectively 5.1 and-2.5.Data point distributes very close at the x direction of principal axis, and not exclusively point-blank, it is comparatively sparse then to distribute at the y direction of principal axis.Because of the natural conditions restriction, some place does not obtain observation data.
Carry out the processing such as spectrum analysis and wavenumber domain inverting for drawing isoline on computers or to these data, need to survey thus the random point data value and calculate numerical value on the regular grid node.Suppose the plane distribution of this regular grid node as shown in Figure 2.In Fig. 2, the position of "+" number expression grid node.The scope of this rule mesh zone on x axle and y direction of principal axis all is 0 ~ 30km, and its sizing grid on x axle and y direction of principal axis is 0.3km all, has 101 * 101=10201 node.
Gridding performing step of the present invention is as follows:
1. read in all random point data as shown in Figure 1;
2. read in the gridding zone at the axial maximal value of x and minimum value and sizing grid, be designated as respectively xmax=0, xmin=30 and dx=0.3;
3. read in the gridding zone at the axial maximal value of y and minimum value and sizing grid, be designated as respectively ymax=0, ymin=30 and dy=0.3;
4. add up minimum and maximum value poor of known random point data, remember that this difference is f_max=5.1-(2.5)=7.6, add up total number of known random point data, remember this total number N=10773;
5. the 1.5km width is expanded respectively on the axial both sides of x in the gridding zone to reading in, and expands respectively the 1.5km width on the axial both sides of y, obtains new gridding zone;
6. total data is extracted by the N/100=107 interval and obtain 100 data, according to these data the boundary node in new grid type zone is utilized formula (1) computing net value of formatting;
7. line by line all internal nodes in new grid type zone are had or judge without obligatory point, to wherein give null value without obligatory point;
8. to the Constrained point, according to corresponding bound data according to formula (2) computing net value of formatting;
9. to without obligatory point, adopt successively formula (4) to upgrade the gridding value, record simultaneously the absolute value of the difference of this nodal value before and after upgrading, be called gridding value renewal amount, be designated as d;
10. the maximal value of obligatory point renewal amount that statistics has or not is designated as d_max;
11. if d_max is more than or equal to 7.6 * 10 -6, returned for the 9th step and continue to calculate; Otherwise adopt the moving window method of average to calculate one by one the gridding value that institute has or not obligatory point according to formula (5), and according to the gridding value of all grid nodes of gridding zone output that read in, calculate end.
Above step is that the microcomputer of 1.7GHz is tested in dominant frequency, and gridding used computing time is 1.7 seconds, and Fig. 3 is the isogram of drawing according to the gridding result, and this figure isoline is smooth.This explanation computing velocity of the present invention is fast, and the gridding result is smooth, has realized the re-set target of invention.

Claims (6)

1. an actual measurement random point data gridding method that is used for geophysics field of utilizing Laplace's equation is characterized in that, may further comprise the steps:
1) the gridding zone is expanded carried out gridding behind the border again and calculate, and two row and the two row nodes of ragged edge in all grid nodes are called boundary node, all the other nodes are called internal node;
2) to boundary node, when the number of all known actual measurement random point data greater than 100 the time, extract 100 data from all given datas, accordingly the data inverse distance weighted interpolation method computing net value of formatting; When the number of all known actual measurement random point data is less than or equal to 100, adopt the inverse distance weighted interpolation method computing net value of formatting according to all given datas;
The node that peripheral region in the internal node is included known actual measurement random point data is called the Constrained point, and Constrained is put the bound data that given data that the peripheral region comprises is called this Constrained point; Above-described peripheral region refers to a rectangular region centered by this node, and this rectangular length and width are respectively half of Gridding length and width;
3) to Constrained point bound data employing inverse distance weighted interpolation method computing net value of formatting according to this point; The node that peripheral region in the internal node is not comprised known actual measurement random point data is called without obligatory point; Above-described peripheral region refers to a rectangular region centered by this node, and this rectangular length and width are respectively half of Gridding length and width;
4) to having or not obligatory point give null value, then adopt the finite-difference approximation formula iterative computation gridding value of Laplace's equation;
5) the having or not obligatory point employing moving window method of average is calculated last gridding value one by one.
2. the actual measurement random point data gridding method that is used for geophysics field of utilizing Laplace's equation according to claim 1, it is characterized in that, in the step 1), the width that expands the border is 5 sizing grids, and the data in gridding zone before the border are not expanded in last only output.
3. the actual measurement random point data gridding method that is used for geophysics field of utilizing Laplace's equation according to claim 1, it is characterized in that, step 2) in, adopt the formula of inverse distance weighted interpolation method computation bound node gridding value to be according to the part given data that extracts or whole given data:
Figure FDA00002865349600021
Wherein, u 1The gridding value of expression boundary node, N represents the number of used given data, k represents the sequence number of used given data, f kThe value that represents k used given data, d kRepresent that k used given data is to the distance of calculation level.
4. the actual measurement random point data gridding method that is used for geophysics field of utilizing Laplace's equation according to claim 1, it is characterized in that, in the step 3), adopt the computing formula of the inverse distance weighted interpolation method computing net value of formatting to be to Constrained point according to the bound data of this point:
Figure FDA00002865349600022
Wherein, u 2The gridding value of expression Constrained point, N represents total number of bound data, k represents the sequence number of bound data, f kThe value that represents k bound data, d kRepresent that k bound data is to the distance of calculation level.
5. the actual measurement random point data gridding method that is used for geophysics field of utilizing Laplace's equation according to claim 1 is characterized in that in the step 4), the expression formula of described Laplace's equation is:
∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 = 0 - - - ( 3 )
Wherein, u represents the gridding value function, and x and y represent two components of Cartesian coordinates; Utilize the differential in the Using Second-Order Central Difference replacement Laplace's equation (3), can obtain the gridding iterative computation formula without obligatory point
u j i = ( u j + 1 i + u j - 1 i ) dx + ( u j i + 1 + u j i - 1 ) dy 2 ( dx + dy ) - - - ( 4 )
Wherein, dx and dy represent respectively length and the width of grid, and i and j represent respectively the grid node sequence number on x axle and the z direction of principal axis,
Figure FDA00002865349600032
Expression is positioned on the x direction of principal axis gridding value of j Nodes on i node and the z direction of principal axis.
6. the actual measurement random point data gridding method that is used for geophysics field of utilizing Laplace's equation according to claim 1, it is characterized in that, in the step 5), the described employing moving window method of average calculate one by one have or not the computing formula of obligatory point gridding value to be:
w jj ii = 1 25 Σ i = ii - 2 ii + 2 Σ j = jj - 2 jj + 2 u j i
Wherein, ii and jj represent respectively the grid node sequence number on x axle and the z direction of principal axis, Represent that last result of calculation is positioned on the x direction of principal axis gridding value of j Nodes on i node and the z direction of principal axis, i and j represent respectively the grid node sequence number on x axle and the z direction of principal axis,
Figure FDA00002865349600035
The expression step 4) is calculated the gridding value that is positioned on the x direction of principal axis j Nodes on i node and the z direction of principal axis of gained.
CN 201110184310 2011-07-04 2011-07-04 Random dot data gridding method by using Laplace equation Expired - Fee Related CN102262789B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN 201110184310 CN102262789B (en) 2011-07-04 2011-07-04 Random dot data gridding method by using Laplace equation

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN 201110184310 CN102262789B (en) 2011-07-04 2011-07-04 Random dot data gridding method by using Laplace equation

Publications (2)

Publication Number Publication Date
CN102262789A CN102262789A (en) 2011-11-30
CN102262789B true CN102262789B (en) 2013-05-29

Family

ID=45009404

Family Applications (1)

Application Number Title Priority Date Filing Date
CN 201110184310 Expired - Fee Related CN102262789B (en) 2011-07-04 2011-07-04 Random dot data gridding method by using Laplace equation

Country Status (1)

Country Link
CN (1) CN102262789B (en)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104376135A (en) * 2013-08-14 2015-02-25 复旦大学 Plane boundary surface charge density extraction method combined with boundary integral equation method and random method
CN103901475B (en) * 2014-03-31 2017-03-08 中国石油天然气股份有限公司 One attribute isogram method for drafting and device
CN108269268B (en) * 2018-01-31 2020-04-21 杭州师范大学 Method for automatically extracting typhoon high-wind-speed cloud system area based on microwave scatterometer data

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR100351674B1 (en) * 1999-10-29 2002-09-10 한국과학기술원 Surface element layer formation method for improving hexahedral mesh shape in finite element method
CN101383047B (en) * 2007-09-03 2011-05-04 鸿富锦精密工业（深圳）有限公司 Curved surface meshing method

Also Published As

Publication number Publication date
CN102262789A (en) 2011-11-30

Similar Documents

Publication Publication Date Title
CN103713315B (en) A kind of seismic anisotropy parameter full waveform inversion method and device
Camera et al. Evaluation of interpolation techniques for the creation of gridded daily precipitation (1× 1 km2); Cyprus, 1980–2010
CN102937721B (en) Limited frequency tomography method for utilizing preliminary wave travel time
CN105549080B (en) A kind of relief surface waveform inversion method based on auxiliary coordinates
CN103499835A (en) Method for inverting near-surface velocity model by utilizing preliminary waveforms
CN107015274A (en) One kind missing seismic exploration data recovery and rebuilding method
CN107765308B (en) Reconstruct low-frequency data frequency domain full waveform inversion method based on convolution thought Yu accurate focus
CN110058302A (en) A kind of full waveform inversion method based on pre-conditional conjugate gradient accelerating algorithm
CN105319581A (en) Efficient time domain full waveform inversion method
CN106483559A (en) A kind of construction method of subsurface velocity model
CN102778699A (en) Electromagnetic data terrain correction method
CN104459784A (en) Two-dimensional Lg wave Q value tomographic imaging method based on single station data, double station data and double event data
CN111580163B (en) Full waveform inversion method and system based on non-monotonic search technology
CN102262789B (en) Random dot data gridding method by using Laplace equation
CN105761310A (en) Simulated analysis and image display method of digital map of sky visible range
CN108318921A (en) A kind of quick earthquake stochastic inversion methods based on lateral confinement
CN104573333A (en) Method for optimizing of model selection based on clustering analysis
CN105607122A (en) Seismic texture extraction and enhancement method based on total variation seismic data decomposition model
CN105184297A (en) Polarized SAR image classification method based on tensor and sparse self-coder
CN101793977A (en) Estimation method of hydrogeological parameters
CN103308945B (en) Simulating generating and forecasting method for first arriving former noise for land exploration
CN102269823A (en) Wave field reconstruction method based on model segmentation
CN103077274A (en) High-precision curve modeling intelligent method and device
CN105676292A (en) 3D earthquake data de-noising method based on 2D curvelet transform
CN103631990A (en) Simulated scene model establishment method and system for SAR irradiation region

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20130529

Termination date: 20140704

EXPY Termination of patent right or utility model