CN106340004B - A kind of parallel Cloud-motion wind inversion method that cloud system is pre-processed based on fuzzy clustering - Google Patents
A kind of parallel Cloud-motion wind inversion method that cloud system is pre-processed based on fuzzy clustering Download PDFInfo
- Publication number
- CN106340004B CN106340004B CN201610642651.XA CN201610642651A CN106340004B CN 106340004 B CN106340004 B CN 106340004B CN 201610642651 A CN201610642651 A CN 201610642651A CN 106340004 B CN106340004 B CN 106340004B
- Authority
- CN
- China
- Prior art keywords
- cloud
- wind
- block
- arrow
- image
- 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
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T1/00—General purpose image data processing
- G06T1/20—Processor architectures; Processor configuration, e.g. pipelining
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10032—Satellite or aerial image; Remote sensing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20024—Filtering details
- G06T2207/20032—Median filtering
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Image Processing (AREA)
Abstract
The invention discloses a kind of parallel Cloud-motion wind inversion method that cloud system is pre-processed based on fuzzy clustering, including：Choose three width cloud atlas of adjacent constant duration, data extraction and pretreatment are carried out to cloud atlas image, every width cloud atlas is divided into polylith and follows the trail of block, and tracking block is assigned to the multiple parallel computation processing units specified, Inversion Calculation is carried out to the wind arrow of the tracking block by maximum correlation coefficient in each independent calculation processing unit, the wind speed and wind direction of the tracking block wind arrow is obtained.The present invention effectively improves the actual execution efficiency of Cloud-motion wind inverting, has important application value to the prediction and early warning of basic, normal, high empty global wind field.
Description
Technical field
The present invention relates to satellite derived wind inverting field, and in particular to a kind of to pre-process the parallel of cloud system based on fuzzy clustering
Cloud-motion wind inversion method.
Background technology
, it is necessary to which substantial amounts of Wind Data, the observation of global atmosphere wind field is served in the research of weather and hazard forecasting
Very positive effect.In recent years, with the constantly improve of land wind field observation network, researcher obtains land Wind Data
What is become is more and more easier.However, for the Wind Data of land, the observational network in the area such as ocean, plateau, desert just seems
It is very not enough, obtain these regional Wind Datas more difficult.By satellite remote dada inverting atmospheric wind, undoubtedly to obtaining
Global atmosphere Wind Data serves the supplementary function being highly profitable.
Meanwhile, some land or water body (being referred to as earth's surface) are because overhead is cloudless or partly cloudy and direct " exposed " is in satellite cloud picture
In.Researcher is mainly studied the cloud in cloud atlas, such as cloud classification, cloud identification, Cloud motion wind, therefore realizes cloud sector
The image segmentation of domain and earth surface area is essential.The research work split for satellite cloud picture image deploys already, mesh
The preceding method used generally has following defect：(1) cloud atlas image document over the years is overly dependent upon, it is necessary to spend manpower and materials
Collect and preserve；(2) split to realize using the image based on threshold value, belong to rigid segmentation, accuracy rate is difficult to ensure that；(3) cloud
Figure is a sufficiently complex image, not only comprising cloud and earth's surface, there are some also between cloud and earth's surface and is in " transition rank
The cloud (cloud for for example growing and dissipating) of section "；So, in conventional methods where, easily forced to be divided into a certain class, but this
It is irrational, is also one of direction of research to the rationally segmentation of cloud atlas image.
In general, satellite cloud picture time interval it is shorter, for inverting target cloud is more, wind arrow calculating density more
Height, the Wind Data quality that Cloud motion wind is finally inversed by also can be higher, but this efficiency all to Cloud motion wind inversion algorithm proposes challenge.
Although there is scholar to propose the simple algorithm of Cloud motion wind inversion algorithm, also it is difficult to meet the wind field calculating for requiring real-time increasingly
Efficiency requirements.Therefore, except the inversion algorithm for studying relatively low calculation cost, the architecture of modern computer is made full use of, is ground
Study carefully parallel Cloud motion wind inversion algorithm, be also to lift one of effective way of inversion algorithm actual operating efficiency.
GPU (Graphic Processing Unit), graphics processor is video card " heart ", equivalent to CPU in microcomputer
On effect.GPU has very high transmission bandwidth and numerous computing units, and computing capability is very powerful.It is continuous with GPU
Upgrading is updated, GPU computing capability more comes outstanding, this provides a kind of parallel acceleration thinking for the inverting of Cloud-motion wind；GPU is in height
Possess the advantage of three highly significants in terms of performance computing：First, the parallel processing capability of data is powerful, and NVIDIA companies are current
Newest GPU peak value floating-point operation abilities are more than 3TFLOPS, and this is almost suitable with a Small-sized C PU cluster；Second, GPU possess
Very high data throughput capabilities, K20 band block is more than 200GB/s；3rd, GPU have friendly DLL, NVIDIA companies
The CUDA frameworks of exploitation support a variety of high-level programming languages, such as C, C++ and Fortran language.
To sum up, although in the prior art, article illustrates the satellite Cloud motion wind parallel refutation algorithm based on multi-core CPU,
But due to CPU be good at processing with complex logic calculation procedure task, such as data compression, artificial intelligence, physical simulation, and
GPU is good at processing computation-intensive and the better simply repetitive task of logic, such as image domains relevant treatment operations, and both are in thing
Manage and there is very big difference in structure and data processing, applied satellite Cloud-motion wind inversion method is to need to overcome many on GPU
Although difficult and technology barrier, GPU possesses hundreds and thousands of computing units, but the computing capability of each computing unit is weak, difficult
To perform complicated branch prediction.GPU is dramatically different with CPU to the access mode of video memory, lacks perfect shared storage logic,
Cache capacity at different levels are smaller, and it is real that this brings theory and technology to the complicated algorithm on tradition CPU to the transplanting of GPU parallel algorithms
Difficulty on now.How using computing capability powerful GPU, the satellite derived wind inverting to high density, high real-time is carried out simultaneously
Row processing is a challenge.
The content of the invention
The present invention has designed and developed a kind of parallel Cloud-motion wind inversion method that cloud system is pre-processed based on fuzzy clustering, the present invention
The first purpose be solve in serial Cloud-motion wind inversion method due to iterations it is more caused by Cloud-motion wind inverse time length, reality
The problem of when property is poor, by parallel refutation method by single iteration parallelization, reduces iterations, and then significantly improve Cloud-motion wind
The real-time of inverting.
The present invention goal of the invention two be solve it is long to target cloud mass search time in defined region of search
The problem of caused Cloud-motion wind inverting poor real, by optimizing the method for tracing of spike cloud mass, significantly simplify to the field of search
The search time of domain target cloud mass, and then significantly improve the real-time of Cloud-motion wind inverting.
The third object of the present invention is solved in the prior art to needing the cloud atlas image preprocessing of progress inverting only
Be use the hard dividing method of image, and then exist to be difficult to during image preprocessing obtain a suitable threshold value, and
By such environmental effects, the problem of accuracy rate is difficult to ensure that.
The technical scheme that the present invention is provided is：
A kind of parallel Cloud-motion wind inversion method that cloud system is pre-processed based on fuzzy clustering, including：
Cloud atlas image is divided into polylith and follows the trail of block, and tracking block is assigned to the multiple GPU parallel computations specified
Array, is matched in each independent computing array by maximum correlation coefficient to the target cloud mass in the tracking block
After carry out Inversion Calculation, obtain it is described tracking block wind arrow wind speed and wind direction；
The matching process of target cloud mass is included by the maximum correlation coefficient：Chosen in the tracking block
Spike cloud mass, in the region of search comprising the spike cloud mass, calculate step-size in search more than 1 all with the target cloud
The coefficient correlation of the equivalent block of block pixel, determines maximum correlation coefficient in these coefficient correlations, is and the spike cloud mass
In similarity highest block, the region for increasing pixel centered on the block, step-length is progressively halved and searches again for out phase
Like degree highest block, until searching out the block that step-length is 1 maximum correlation coefficient searched, as target cloud mass；And
Wherein, original cloud atlas image is entered by fuzzy clustering method after separating with racking and determines to carry out Cloud-motion wind inverting
Cloud atlas image, it includes：
The gray feature data of infrared channel in cloud atlas image are gathered, the gray feature data are handled, are obtained
Sample gray feature data, the sample gray feature data are classified, and make each class sample gray feature data point
Not Dui Ying a kind of cloud atlas type, and then determine the cloud atlas image.
Preferably, in addition to, choose three width cloud atlas of adjacent constant duration, the wind of every width cloud atlas wind arrow drawn respectively
Speed and wind direction, by the wind speed of the wind arrow, wind direction and according to sequential by the wind arrow speed difference of adjacent cloud atlas, side in the cloud atlas
It is compared to difference with threshold value, performs following principle：
If wind arrow speed is less than the minimum speed threshold value specified, this wind arrow is rejected；Or
If the difference of two wind arrow speed is more than the maximum speed discrepancy threshold value specified, this wind arrow is rejected；Or
If the difference in two wind arrow directions is more than the poor threshold value in the maximum direction specified, this wind arrow is rejected；Or
If there is no the poor wind arrow for being less than the poor threshold value in the maximum direction specified in direction therewith in the peripheral region of two wind arrows, pick
Except this wind arrow.
Preferably, multi-channel data extraction is carried out to the cloud atlas view data, retains infrared channel data.
Preferably, in addition to：The cloud atlas image is pre-processed, it is included to the longitude and latitude projection such as image progress
Convert and original pie chart is transformed to histogram.
Preferably, the cloud atlas image is pre-processed also includes carrying out grey level enhancement to described image.
Preferably, gray feature data progress processing is included：
The gray feature data input of infrared channel in cloud atlas image is obtained into pixel grey scale value matrix, passed throughLeading diagonal pixel grey scale average in a matrix is calculated, is passed throughCalculate auxiliary in a matrix
Diagonal pixels gray average, and then pass through M=(M1+M2Obtain pixel grey scale average, i.e. sample gray feature data in)/2.
Preferably, the first LONG WAVE INFRARED, the second LONG WAVE INFRARED and short-wave infrared are chosen in the infrared channel
The sample gray feature data are set up 5 gray feature data statistics amounts, Nephogram are entered by sample gray feature data
The segmentation of row image obtains the cloud atlas image.
Preferably, calculated according to the width height and wind arrow density of every width cloud atlas pixel and determine the tracking block number,
It is (Width/ ρ) × (Height/ ρ), in formula, and Width is the width of every width cloud atlas pixel, and Height is every width
The height of cloud atlas pixel, ρ is wind arrow density.
Preferably, the matching process of target cloud mass is included by the maximum correlation coefficient：In the maximum
In correlation coefficient process, spike cloud mass A is chosen in the tracking block, it is N × N pixel image blocks, is including the spike
In cloud mass A M × M pixel image blocks region of search, step-size in search being all equal with the spike cloud mass A pixels for m is calculated
The coefficient correlation of block, maximum correlation coefficient is determined in these coefficient correlations, is and the spike cloud mass A similarity highests
N × N pixel image block A1, in block A1 increasesStep-length, is then kept to by individual pixelIndividual pixel, in Yi Qu
Centered on block A1The coefficient correlation of each N × N pixel images block is calculated in the region of the length of side, in these coefficient correlations
It is middle to determine maximum correlation coefficient and then determine the next and spike cloud mass A similarities highest N × N pixel image block A2,
Step-length is progressively halved until searching out the block that step-length is 1 maximum correlation coefficient searched, as target cloud mass passes through institute
The positional information calculation for stating spike cloud mass and the target cloud mass obtains the wind speed and wind direction of wind arrow.
The present invention is had the advantage that compared with prior art：
1st, the cascading algorithm in satellite derived wind inversion algorithm is changed to parallel algorithm, iterations greatly reduces,
By the computational methods of this single iteration parallelization, so that the wind arrow in whole Zhang Yun's figure can concurrently be performed in GPU, enter
The wind arrow computing of the extensive quantity of row, in the case of the precision and performance for the Cloud motion wind product for not influenceing to be finally inversed by, is saved significantly
The time needed for inverting and wind arrow verification has been saved, and then wind arrow computational efficiency is greatly improved, Cloud-motion wind is significantly improved
The real-time of inverting；
2nd, in the present invention, by the way of variable step is iterated to calculate, it is greatly lowered in region of search to target cloud mass
Search time, and then increase substantially follow the trail of block wind arrow computational efficiency, significantly improve the real-time of Cloud-motion wind inverting.
3rd, in the present invention, more preferable cloud can be obtained by applying to the segmentation pretreatment of cloud atlas image by fuzzy clustering method
Ground separating effect, this method can effectively carry out the segmentation of cloud atlas uncertainty and lifting segmentation accuracy rate, after being separated in cloud
Cloud-motion wind inverting is carried out on image so that efficiency of inverse process is more accurate, while before fuzzy clustering pretreatment is carried out, to gray scale
Characteristic carries out pretreatment by Fast Median Filtering method or diagonal filtering method and eliminates salt-pepper noise, heightens
Picture quality, makes target clear, eliminates the noise jamming produced in the calculating to Cloud-motion wind and identification work, further carries
High efficiency of inverse process.
Brief description of the drawings
Fig. 1 is that tasks of the CPU and GPU in the satellite derived wind inversion method based on GPU platform divides schematic diagram.
Fig. 2 is the FCM result figures of cloud atlas image.
Fig. 3 is the cloud atlas image segmentation flow chart based on FCM.
Basic flow sheets of the Fig. 4 for Cloud-motion wind inversion algorithm after parallel.
Fig. 5 is influence figure of the stream handle quantity to efficiency.
Fig. 6 is the schematic diagram based on the serial and concurrent inversion algorithm accuracy comparison 1 of GPU Cloud motion winds.
Fig. 7 is the schematic diagram based on the serial and concurrent inversion algorithm accuracy comparison 2 of GPU Cloud motion winds.
Fig. 8 is the Cloud motion wind product figure being finally inversed by based on the parallel Cloud motion wind algorithms of GPU.
Personal attendant's thread block increased tendency chart when Fig. 9 is performs.
Figure 10 is speed-up ratio with the increased tendency chart of thread block.
Figure 11 is that spike cloud mass follows the trail of schematic diagram.
Figure 12 is that wind arrow speed derives schematic diagram.
Figure 13 is that wind arrow direction derives schematic diagram.
Figure 14 is the iterative search schematic diagram that step-length is 8.
Figure 15 is the iterative search schematic diagram that step-length is 4.
Embodiment
The present invention is described in further detail below in conjunction with the accompanying drawings, to make those skilled in the art with reference to specification text
Word can be implemented according to this.
As shown in Figure 1 to 4, the invention provides a kind of parallel Cloud-motion wind inverting that cloud system is pre-processed based on fuzzy clustering
Method, comprises the following steps：
Step 1: using three width infrared channel cloud atlas of adjacent constant duration, setting inverted parameters, including satellite remote sensing
Read path, the density of Cloud-motion wind inverting, Cloud-motion wind inversion algorithm relevant parameter, GPU parallel parameter, cloud of the data on disk
Mark wind verification dynamics parameter, the storing path of Cloud-motion wind retrieval products on disk etc.；
Step 2: being decompressed using lossless decompression algorithm to the earth remote sensing data issued from meteorological satellite load
Contracting；
Step 3: to satellite infrared image data extraction, i.e., the satellite that flow generation is decompressed to satellite remote sensing date is more
Channel data carries out data extraction, retains infrared channel data, gives up other channel datas, and opens up that the cloud atlas is wide high to be deposited
Space is stored up, by the cloud atlas data duplication of the infrared channel of extraction to the memory space；
Step 4: pre-processed to the cloud atlas data after extraction, including the enhancing of image projection transformation, gradation of image and
Enter determination cloud atlas image after separating with racking to the cloud atlas image of extraction by fuzzy clustering method；
Step 5: initialized to performing environment, including to the initial of the performing environment of inverted parameters and GPU equipment
Change, it includes：Being calculated according to the width height and wind arrow density of every width cloud atlas pixel needs the wind arrow number of inverting, in GPU
Video memory in open up Array for structural body for preserving Cloud-motion wind, opened up in GPU video memory for preserving each wind arrow using most
Big coefficient correlation and the array of correlation coefficient value produced, are opened up in GPU video memory for storing 3 width infrared cloud image data
Array, CPU copies to three width cloud atlas data in main memory in three arrays that GPU video memorys are opened up successively.
Step 6: carrying out Cloud-motion wind parallel refutation calculating, CPU carries series of parameters and multi-task planning and scheduling rule
Enter GPU equipment, and trigger the parallel refutation of GPU one wind arrow of each thread correspondence to calculate, every width cloud atlas is divided into polylith follows the trail of
Block, and tracking block is assigned to the multiple parallel computation processing units specified, in each independent calculation processing unit
Inversion Calculation is carried out after being matched by maximum correlation coefficient to the target cloud mass in the tracking block, described chase after is obtained
The wind speed and wind direction of track block wind arrow, and be saved in GPU video card；
Wherein, the width of every width cloud atlas pixel, height and wind arrow density, which are calculated, determines the tracking block number, and it is
In (Width/ ρ) * (Height/ ρ), formula, Width is the width of every width cloud atlas pixel, and Height is every width cloud atlas
The height of pixel, ρ is wind arrow density.
Step 7: carrying out quality indicator to Cloud-motion wind, the content of verification is wind arrow speed and wind arrow direction, the object of verification
Be 3 width adjacent time sequences cloud atlas in the result that is finally inversed by of preceding two Zhang Yuns figure and the result that is finally inversed by of rear two Zhang Yuns figure, will
The wind speed of wind arrow, wind direction and the wind arrow speed difference of adjacent cloud atlas, direction difference in three width cloud atlas are carried out with threshold value according to sequential
Compare, performed by following principle：
If wind arrow speed is less than the minimum speed threshold value specified, this wind arrow is rejected；Or
If the difference of two wind arrow speed is more than the maximum speed discrepancy threshold value specified, this wind arrow is rejected；Or
If the difference in two wind arrow directions is more than the poor threshold value in the maximum direction specified, this wind arrow is rejected；Or
If there is no the poor wind arrow for being less than the poor threshold value in the maximum direction specified in direction therewith in the peripheral region of two wind arrows, pick
Except this wind arrow.
Step 8: carrying out data preservation to above-mentioned Cloud-motion wind product, final retrieval products are obtained.
In another embodiment, in step one, remotely-sensed data read path and product in inverted parameters preserve road
Footpath is absolute path；The unit of inverting density is the inverting wind arrow number that every 1 degree of longitude and latitude grid block is included；Inverted parameters
Include the number of pixels of spike cloud mass and trace regions ranks；GPU parallel parameters are the net of the execution unit of GPU kernel functions
Lattice division rule；Cloud-motion wind verification dynamics parameter is Cloud-motion wind credit rating, and credit rating is higher, it is desirable to stricter, is finally inversed by
Wind arrow it is with a high credibility, but number is fewer, otherwise quality is low, and number is more.
In another embodiment, in step 2, satellite remote sensing date decompression algorithm can use and be based on Run- Length Coding
Decompression algorithm, the LZW decompression algorithms based on dictionary encoding technology, the decompression algorithm based on Huffman encoding principle or
Decompression algorithm of the person based on arithmetic coding.
In another embodiment, in step 4, to longitude and latitude projective transformations such as image progress, Mercator's throwing is transformed into
Shadow, to carry out wind arrow inverting, histogram is transformed to by the original pie chart of satellite, facilitates the division and execution of GPU task, to carry
Rise inverting quality；
Grey level enhancement processing is carried out to image can be using linear gradation conversion enhancing algorithm, piecewise linear gray transformation enhancing
Algorithm, histogram equalization enhancing algorithm or local histogram equalization enhancing algorithm strengthen the contrast of image；
In another embodiment, in step 4, original cloud atlas image is entered by fuzzy C-means clustering (FCM)
The cloud atlas image is determined after separating with racking, it comprises the following steps：
Input：Satellite cloud picture image
Output：Eliminate the cloud atlas image of earth surface area
Step a, cloud atlas is pre-processed first, the method for processing is so as to removing cloud atlas image using median filtering method
In noise and interference.
Step b, the target data set for building FCM.Choose infrared cloud image and visible cloud image builds IR-VIS two-dimension spectrums
Feature space, each pixel is stored in infrared, visible images using a two-dimensional array data [256,2] respectively by row
Gray feature value.
Step c, calculating sample general features value.If not having to record 8 kinds of earth's surfaces, gray features of the varieties of clouds in database,
Cloud atlas Sample Storehouse is analyzed, 8 kinds of earth's surfaces, the gray feature of the varieties of clouds is counted and is recorded in database.If so, then ignoring
This step, is continued executing with.
Step d, initial cluster center is set to sample characteristics carries out FCM.8 class earth's surfaces, cloud are obtained from database
The gray feature of class and is clustered as FCM initial cluster center in Infrared-Visible two-dimension spectrum feature space.
Step e, to cloud atlas carry out image segmentation.All pixels are clicked through using the cluster block plan obtained by step (4)
Row judges, if its two-dimentional gray value is in earth's surface gray feature region, for earth's surface pixel, labeled as black；If not existing,
For cloud pixel, retain.
Finally obtain rectangle, grey level enhancement and inactive area propose after cloud atlas.
In another embodiment, in step 6, series of parameters and multi-task planning and scheduling rule are loaded into by CPU
GPU equipment, and the parallel refutation calculating of GPU one wind arrow of each thread correspondence is triggered, GPU performs the anti-of all wind arrows parallel
Drill, the inverting of each wind arrow is performed in independent thread, the inverting of wind arrow calculates target cloud using maximum correlation coefficient
Block and the position relationship for following the trail of cloud atlas, so that it is determined that wind arrow speed and direction, are finally converted to the location of pixels where wind arrow
Actual longitude and latitude position.
In another embodiment, in step 6, Cloud motion wind inverting parallel algorithm uses the Fork- of process-thread
Join models, host process reads cloud atlas data and parameter and generated after network subdivision relationships between nodes, is referred to by OpenMP compilings
Lead sentence, the computational load of n wind arrow is dynamically assigning on m thread, using maximum correlation coefficient matching target cloud mass,
Calculate the steps such as wind speed and direction, wind arrow fixed height, the wind arrow quality indicator of wind arrow independently to be calculated by each thread, reuse maximum phase
Close in Y-factor method Y matching target cloud mass, the global wind arrow array step of write-in, each thread shares cloud atlas datarams region and wind arrow
All wind arrows of array region of memory are calculated finish after, sentence end lines journey is instructed by compiling, host process is by the wind arrow calculated
Data write-in disk is preserved.
In another embodiment, in step 6, maximum correlation coefficient refers to seek in satellite cloud picture trace regions
The similarity between cloud mass similarity highest cloud mass, cloud mass is looked for and followed the trail of then by coefficient correlation as index, i.e., is chased after in cloud atlas
The cloud mass for the coefficient correlation maximum that cloud mass is calculated is found and followed the trail of in track region, specifically includes following steps：
(1) cloud atlas 1 is divided into multiple tracking blocks, cloud is traveled through successively and each follows the trail of block；
(2) in the wind arrow inverting block distributed in current GPU threads in the pixel region of centre of location position, record
Heart point position；
(3) gray variance and standard deviation in region in (2) are calculated；
(4) begun stepping through from first pixel region in the upper left corner of the same block in cloud atlas 2, until traversal knot
Beam, return (1), start on to 0, under from left to right travel through, until traversal terminate, calculate the variance and standard of the same block
Difference；
(5) coefficient correlation of (3) and (4) is calculated, and writes the result into storage matrix；
(6) maximum correlation coefficient in the storage matrix of (5) is found out, its correspondence position is recorded；
(7) according to the position coordinates of (2) and (6), the wind arrow distance and angle of the region unit is calculated, is removed with the distance of wind arrow
With time interval, then wind arrow speed can be obtained；The form that the angle for the wind arrow being finally inversed by is converted into northern drift angle is wind arrow
Direction.
Specifically, in step 6, the matching process of target cloud mass is included by the maximum correlation coefficient：
In the maximum correlation coefficient, spike cloud mass A is chosen in the tracking block, it is N × N pixel image blocks, in bag
In M × M pixel image blocks region of search containing the spike cloud mass A, calculate step-size in search for m all with the spike cloud mass
The coefficient correlation of the equivalent block of A pixels, determines maximum correlation coefficient in these coefficient correlations, is and the spike cloud mass A
Similarity highest N × N pixel image block A1, in block A1 increasesStep-length, is then kept to by individual pixelIndividual picture
Vegetarian refreshments, centered on block A1The coefficient correlation of each N × N pixel images block is calculated in the region of the length of side,
Maximum correlation coefficient is determined in these coefficient correlations and then the next and spike cloud mass A similarities highest N × N pixels are determined
Image block A2, step-length is progressively halved until searching out the block that step-length is 1 maximum correlation coefficient searched, as target
Cloud mass, the wind speed and wind direction of wind arrow are obtained by the positional information calculation of the spike cloud mass and the target cloud mass.
In another embodiment, in the step 8, Cloud-motion wind product storage method is CPU by GPU video memory
Wind arrow data duplication into computer hosting, then export an external memory disk in.
In another embodiment, the data object of inverting of the present invention is the cloud atlas data of the passage of FY-2E satellite infrareds one,
The size of the cloud atlas is that 2581*1399 pixels, i.e. Width are 2581 pixels, and Height is 1399 pixels, chooses constant duration
The cloud atlas of three adjacent passages of FY-2E satellite infrareds one as the present embodiment data processing object.
In another embodiment, inverted parameters in the present invention receive flow, and input satellite remote sensing date reads road
Footpath Path1 and satellite derived wind product storing path Path2；The density p of input reverse；The side of spike cloud mass in input reverse
Long, i.e., target cloud block size is P*P；Input Cloud-motion wind inverting grade D.
In another embodiment, satellite remote sensing date decompression flow in the present invention, according to cloud provided above
Figure read path Path1, reads the remotely-sensed data of FY-2E satellites, using lossless decompression algorithm to FY-2E satellite remote sensing numbers
According to being decompressed, decompression obtains the cloud atlas pixel data of each passage, Temperature Scaling data, longitude and latitude calibration data etc..
In another embodiment, the image data extraction flow of satellite infrared one in the present invention is decompressed from above-mentioned
To each passage cloud atlas data in, the cloud atlas data of satellite infrared one of three Zhang Yun's figures are extracted respectively, copy to microcomputer host open
In array Data1, Data2 and Data3 for warding off.
In another embodiment, gradation of image enhancing processing is contrasted using Nogata equalization algorithm to cloud atlas data
Degree enhancing.
In another embodiment, performing environment initialization flow in the present invention initializes CPU and GPU execution
Environment and internal storage data.Calculated according to the width of density and cloud atlas is high, it is necessary to the Cloud-motion wind number n of inverting；CPU calls CUDA interfaces
The n Array for structural body windRecords for being used to preserve Cloud-motion wind are opened up in GPU shared drive；CPU calls CUDA interfaces
N size is opened up in GPU shared drive to produce using maximum correlation coefficient for preserving each wind arrow for (P+1) * (P+1)
The array of raw (P+1) * (P+1) individual correlation coefficient value；CPU calls CUDA interfaces to open up three use in GPU shared drive
In the array for storing three infrared cloud atlas data, the size of array is the size of cloud atlas according to claim 4
Width*Height；CPU calls CUDA interfaces that three width cloud atlas data in main memory are copied into GPU equipment in step 1 successively
In three arrays that shared drive is opened up.
In another embodiment, Cloud-motion wind inverting flow process in the present invention, CPU is by the inverting of each wind arrow on cloud atlas
It is divided into a mission thread, all thread dividings is performed in grid to a GPU, GPU is according to the task of thread lattice point
Match somebody with somebody and dispatching principle, the inverting to all wind arrows is concurrently performed, and inverting determines the displacement of wind arrow using maximum correlation coefficient, enters
And determine velocity magnitude and the direction of wind arrow.The location of pixels where wind arrow is finally converted into actual longitude and latitude position, will be anti-
Drill the WindRecords arrays opened up in the index position write-in GPU shared drives that obtained wind arrow is distributed by individual execution thread
In.
In another embodiment, flow is preserved in the Cloud-motion wind product of the present invention, CPU falls into a trap GPU shared drive
Cloud-motion wind data duplication after the quality indicator calculated is into the main memory of host's microcomputer, and then CPU is by the Cloud-motion wind number in main memory
According to being preserved, storing path is that above-mentioned inverted parameters receive the Cloud-motion wind product storing path Path2 that flow is received, preservation
Form is：Wind arrow numbering, wind arrow longitude, wind arrow latitude, wind arrow speed and wind arrow direction.
Specific description is done to parallel refutation method and maximum correlation coefficient in conjunction with specific embodiment below.
In another embodiment, as shown in Figure 2 and Figure 3, divide while cloud atlas is entered with density clustering method to rack
From；Pixel is divided according to " similitude " of gray feature between cloud atlas spectral digital image slices vegetarian refreshments come reach cloud sector domain with
The purpose of earth surface area segmentation.Traditional threshold method belongs to the hard dividing method of image, is not only difficult to obtain a suitable threshold
Value, and by such environmental effects, accuracy rate is difficult to ensure that.This method can effectively carry out the segmentation of cloud atlas uncertainty and lifting point
Accuracy rate is cut, Cloud-motion wind inverting is carried out on the image after being separated in cloud so that efficiency of inverse process is more preferably accurate；In satellite digital figure
During generation, processing and projective transformation of picture etc., cloud atlas generally all unavoidably produces some salt-pepper noises, and these noises are not
But picture quality is reduced, also makes objective fuzzy, calculating to Cloud-motion wind and identification work generate very big interference, therefore
It is necessary that noise and reduction interference are removed in satellite cloud picture image.
Gray feature data in satellite cloud picture image are removed at noise and reduction interference by median filter method
Reason obtains gray scale sample characteristics data, specific as follows：
Box filter masterplate is used, following two dimension median filter formula is defined：
Zxy=Med { g (i, j) }=Med { g (x+r, y+s), (r, s) ∈ A (x, y) ∈ I2}
Wherein, Med represents to take median operation；P (x, y) gray value is set to Z during processingxy, as shown in table 1.
The medium filtering module of table 13 × 3
| 0th row | 1st row | 2nd row | |
| 0th row | P0 | P1 | P2 |
| 1st row | P3 | P4 | P5 |
| 2nd row | P6 | P7 | P8 |
Often calculate an intermediate value to be ranked up the gray value of each point in template, according to row or column the pixel in template point
Median calculation is carried out into three groups (every group of 3 pixels), obtained intermediate value is sample gradation data feature, is filtered using Quick Median
Wave method can reduce the number of comparisons between pixel, and by Parallel Implementation, can be greatly enhanced the place of image denoising
Speed is managed, Fast Median Filtering method specifically includes：
Input：Grey scale pixel value in 3 × 3 windows
Output：The gray scale intermediate value of 9 pixels
(1) to every a line in template, maximum, minimum value and intermediate value are asked for respectively, obtain data group A=[Max0,
Max1, Max2], B=[Min0, Min1, Min2], C=[Med0, Med1, Med2], it is as follows.
A：Max0=max [P0, P1, P2], Max1=max [P3, P4, P5], Max2=max [P6, P7, P8]
B：Min0=min [P0, P1, P2], Min1=min [P3, P4, P5], Min2=min [P6, P7, P8]
C：Med0=med [P0, P1, P2], Med1=med [P3, P4, P5], Med2=med [P6, P7, P8]
(2) minimum value in A is sought：Maxmin=min [Max0, Max1, Max2]；
(3) intermediate value in C is sought：Medmed=med [Med0, Med1, Med2]；
(4) maximum in B is sought：Minmax=max [Min0, Min1, Min2]；
(5) med [Maxmin, Medmed, Minmax] is sought, is required intermediate value.
In theory, we are only needed to be partitioned into cloud sector domain and earth surface area from cloud atlas, that is to say, that clusters number is set
It is set to 2.But the cloud atlas segmentation effect being obtained by is extremely undesirable, there is the mistake that cloud sector domain is divided into earth's surface by mistake, therefore
It is necessarily increased clusters number to improve the accuracy rate of cloud atlas image segmentation, in the present embodiment cloud atlas objects in images type
8 classes, i.e. land, water body, cumulonimbus, cumulus congestus, cirrus, medium cloud, cumulus humilis and low clouds are divided into, therefore FCM is clustered herein
Division numbers be set to 8, by such selection, the professional knowledge of meteorological field can be made full use of, constraint and guides poly-
Class process, improves the accuracy rate of cluster.The number of times of circulation is greatly reduced, so that the effect of algorithm is more notable.
In the present embodiment, the sample gradation data feature of infrared (IR) and two passage cloud atlas of steam (WV) is chosen,
The gray feature statistic (averages of some samples) of this 8 class sample is set up, as shown in table 2.
The gray feature statistics of the cloud atlas sample of table 2
It is 8 to determine clusters number, and the initial center that sample gray feature statistic is clustered as FCM, to constrain and
Continuous adjusting and optimizing cluster centre and subordinated-degree matrix in cluster process, cluster iterative process are instructed, two dimensional character is finally obtained
The cluster result in space.
Specific FCM algorithms are as follows：
Input：Clusters number c=8, data bulk n；Fixed value λ, 1≤λ≤∞；Threshold epsilon；
Output：Subordinated-degree matrix U'=[uij]n×n
(1) average of earth's surface and typical varieties of clouds characteristic quality of sample is assigned to V0=[vij]2×8；
(2) cycle-index b=0 is made, is subordinated-degree matrix U using formula in step (5)(0)Assign initial value；
(3)b++；
(4) c cluster centre V is calculatedi (b)(i=1,2 ..., c), formula such as (1-1).
(5) for j=1~n, U is recalculated(b), and make U(b+1)=U(b)；
1. I is calculatedjWith
Ij=i | 1≤i≤c；dij=| | Xj-Vi| |=0 } (1-2)
2. for data Xj, update its degree of membership.IfThen to allMake uij=0,j
=j+1；Otherwise
(6) U is compared(b)And U(b+1).IfPerform step (7)；Otherwise (3) are returned to；
(7) U '=U is made(b+1)；Terminate.
FY-2E meteorological satellites have visible ray, infrared and three spectrum cloud atlas of steam, therefore each pixel in cloud atlas
With visible ray, infrared and three kinds of gray feature attributes of steam；FCM existing more intuitively effects in two dimensional surface space have again
Higher execution efficiency, it is in the present embodiment, infrared and visible ray is made a concrete analysis of from three kind Sexual behavior modes, using infrared logical
Road and visible channel build orthogonal two-dimension spectrum feature space, make on image pixel with the point in two-dimensional feature space one by one
Mapping, i.e. cloud atlas pixel are RI(m,n),RV(m, n), two-dimension spectrum feature space point is SIV(RI,RV), by having described above
Body FCM algorithms carry out FCM clusters to two dimensional gray feature space.
The division of transition cloud system and not distinguishable cloud system is carried out to cloud atlas image：Assuming that the earth's surface, the varieties of clouds in cloud atlas can not be thin
It is divided into c classes：(X1,X2,...,Xc), SIV(AI,AV) in the degree of membership each put be U (u1,u2,...,uc), then can be according to following
Condition determines the classification belonging to each pixel in cloud atlas：If ui=Max (U) ＞ 0.5, then decision-pointCloud
Pixel R in figureI(m,n)、RV(m, n) belongs to classification i；IfThen the pixel is subordinate between two kinds of classifications
Category degree can be classified as transition cloud system higher than bottom threshold and less than upper threshold；If ui≤ 0.3, then the pixel feature abnormalities mould
Paste, it is difficult to distinguish, not distinguishable cloud system can be classified as.Threshold value bound in above-mentioned partiting step can be adjusted according to actual conditions.
The concept of transition cloud system and not distinguishable cloud system is introduced, returning for obscured in image, uncertain pixel is solved
Category problem, carries out FCM clusters to cloud atlas again, obtains Clustering Effect as shown in Figure 2.
It is seen that degree of membership higher than upper threshold (0.5) pixel formed in two-dimensional feature space 8 it is bright
Aobvious cluster areas (A, B, C, D, E, F, G, H), the scope in these regions reflects the gray feature of 8 quasi-representative samples in table 2
Change；X, Y region represent transition cloud system and not distinguishable cloud system respectively；With reference to the spectroscopic data in table 2, according to pixel in two dimension
The distribution of feature space, it may be determined that A is land, B is water body, as shown in Figure 2.
If land, the cluster areas of water body are respectively S in gray scale two-dimensional feature spaceIV(land), SIV(water body), cloud atlas
Middle pixel (x, y) gray value in infrared and visible cloud image is respectively PI(x,y)、PV(x, y), point (x, y) exists after image segmentation
Gray value in infrared and visible cloud image is G respectivelyI(x,y)、GV(x,y).The image segmentation thought based on threshold value is then combined,
The criterion split using the image of FCM methods can be set up, it is as follows：
In another embodiment, gray feature data in satellite cloud picture image are gone by diagonal filtering method
It is specific as follows except noise and reduction interference processing obtain gray scale sample characteristics data：
Cloud atlas represents the radiation that object is sent using gray-scale level, and the gray value that pixel has represents cloud or earth's surface surface
Bright temperature value.The more low then object height above sea level of bright temperature value is higher, such as altostratus；The more high then object height above sea level of bright temperature value is lower, such as land,
Ocean etc..During generation, digital processing and projective transformation etc., satellite cloud picture is logical unavoidably to produce some noises, noise
Picture quality is not only reduced, also makes objective fuzzy, is even more ground to follow-up image segmentation, feature extraction, Cloud-motion wind inverting etc.
Study carefully work and cause inconvenience, therefore it is necessary to remove noise and reduction interference in satellite cloud picture image.
For the feature of Cloud-motion wind inverting main body, the globality especially cloud border of cloud Cloud-motion wind inversion result is influenceed compared with
Greatly, therefore, it is proposed that not only can overcome the disadvantages that cloud cluster noise but also diagonal (major-minor) the mean filter method on clear cloud border can be retained,
Algorithm is as follows：
Diagonal filtering method
Input：Grey scale pixel value matrix A=[a in n × n windowsij]n×n
Output：The gray average M of diagonal pixels
Step1：Seek leading diagonal pixel grey scale average
Step2：Seek auxiliary diagonal pixel grey scale average
Step3：Seek M=(M1+M2)/2；
By the way that the purpose of cloud atlas image classification is had into two, one is cloud and land and water body are separated, it is ensured that Yun Ji
The input of wind algorithm is pure " cloud ".The second is cloud is divided into " high cloud ", " medium cloud ", " low clouds " three major types, it is anti-with Cloud-motion wind
" upper-level winds ", " hollow wind " and " Low level wind " correspondence drilled, obtains more targeted inversion result.
Choose LONG WAVE INFRARED 1 (11.5-12.5 μm), LONG WAVE INFRARED 2 (10.3-11.3 μm) and short-wave infrared (3.5-4.0 μ
M) gray feature of three passage cloud atlas, sets up the gray feature statistics of this 5 class sample of water body, earth's surface, low clouds, medium cloud and high cloud
Amount, as shown in table 3.
The gray feature statistics of the cloud atlas sample of table 3
Determine earth's surface, water body and cloud system total number be 5, using three-dimensional samples gray feature statistic as cluster it is initial in
Continuous adjusting and optimizing cluster centre and subordinated-degree matrix, finally obtain the cluster in three-dimensional feature space in the heart, cluster iterative process
As a result, cloud system clustering algorithm is as follows：
Cloud system clustering algorithm：
Input：Clusters number c=5, Three-Dimensional Gray eigenmatrix value V0=[vij]3×5, data volume n；Algorithm terminates threshold epsilon
Output：Subordinated-degree matrix U'=[uij]n×n
Calculating process：
Step1：V is assigned to 5 class three-dimensional samples characteristic quantity initial values0, to subordinated-degree matrix U(0)Assign as follows
Initial value, b=0；
Step2：b++；Calculate new cluster centre Vi (b)(i=1,2 ..., c)：
Step3：For j=1~n, U is calculated(b), and make U(b+1)=U(b)；
1. I is calculatedjWith
Ij=i | 1≤i≤c；dij=| | Xj-Vi| |=0 }
2. for data Xj, update its degree of membership.
IfThen to allMake uij=0；
Otherwise：
Step4：If | | U(b)-U(b+1)| | ＜ ε, perform step step5；Otherwise Step3 is returned；
Step5：U '=U(b+1)；Terminate
As shown in figure 4, because the calculating of each wind arrow in cloud atlas is all completely self-contained, wind arrow in cloud atlas can be selected to carry out
Coarse grain parallelism, the basic procedure of parallel Cloud-motion wind inversion algorithm is as follows：
(1) decompression satellite issues data；
(2) the satellite infrared passage cloud atlas that three interval times are Dt is taken, temporally cloud atlas 1, cloud is sequentially named as sooner or later
Fig. 2, cloud atlas 3；
(3) gradation of image is strengthened into algorithm respectively and three pictures is subjected to grey level enhancement respectively, more accurately to track
Target cloud mass；
(4) cloud atlas is divided into (Width/ ρ) * (Height/ ρ) block to follow the trail of in block, formula, Width is every width cloud atlas
The width of pixel, Height is the height of every width cloud atlas pixel, and ρ is wind arrow density；
(5) calculation processing unit specified to all blocks distribution on cloud, is used necessarily n parallel processing element
Regular carry out task distribution.Parallel processing element herein includes each stream of the thread and GPU in MPI process, openMP
Processor, and the complete wind arrow inverting task of Charge-de-Mission one in each parallel processing element；
(6) start n parallel processing transaction and parallel computation is carried out to the processing task distributed；
(7) data summarization.The wind arrow of all parallel processing element invertings is collected together, and by effective wind arrow data
Disk is write to preserve；
(8) computing terminates.
In the present embodiment, as shown in figure 5, using the tall and handsome CUDA frameworks up to company, thread is stream handle in GPU
The base unit of execution, and the block comprising multiple thread is the elementary cell of processor scheduling, multiple block compositions one
Individual grid.Multiple thread in same block can share same memory space.The present invention is with a grid as one
Secondary GPU parallel refutations task, performs CUDA kernel function kernel_function<<<M, N>>>(arguments),
Wherein M is that the grid of grid size, i.e., one includes how many block；N is block size, i.e., have in each block
How many thread, each thread distributes the calculating task of a wind arrow.
Cloud motion wind inversion algorithm is accelerated using GPU concurrent techniques, MN wind arrow distribution of computation tasks is arrived GPU's
Performed on thread, the application of thread lattice is one-dimensional variable blocksPerGrid, each for the sake of simplicity by experiment of the invention
Thread Count on threadsPerBlock thread, each block of experiment of the invention is opened up in thread lattice
ThreadsPerBlock is 256, block quantity in thread lattice
BlocksPerGrid=(iMax+threadsPerBlock-1)/threadsPerBlock
The index that each wind arrow inverting task obtains the wind arrow for wanting inverting is
WindIndex=blockDim.x*blockIdx.x+threadIdx.x
Finally call
dev_makeCloud<<<blocksPerGrid,threadsPerBlock>>>(arguments ...) is to each wind arrow
Inverting task is allocated and performed, and the thread in each block is performed serially in a stream handle, and will be finally inversed by
Into the corresponding position windIndex of global wind arrow array opened up.
Accuracy comparison is tested
(1) when inverted parameters are that target cloud block size is 22 × 22 pixels, trace regions size is 44 × 44 pixels, wind arrow
Density is 2/when spending, it have chosen on equator that continuous 10 wind arrows are as comparison other, parallel with the wind of serial algorithm inverting
Swear that accuracy comparison is as shown in table 4.(wherein S represents serial, and T-64 represents 64 thread of distribution, T-128 generations in a thread block
In one thread block block of table distribute 128 thread thread, then the block numbers in thread lattice be general assignment quantity WN divided by
Thread quantity in each block).
Serial and concurrent Cloud motion wind inversion algorithm accuracy comparison 1 of the table 4 based on GPU
(2) inverted parameters are that target cloud block size is 32 × 32 pixels, and trace regions size is 64 × 64, and wind arrow density is
4/degree, continuous 10 wind arrows of 32 degree of north latitude are have chosen as comparison other, parallel with the wind arrow precision of serial algorithm inverting
(wherein S represents serial, and T-64 represents distribution 64 thread, T-128 in a thread block and represents one as shown in table 5 for contrast
128 thread thread are distributed in thread block block, then the block numbers in thread lattice are general assignment quantity WN divided by each
Thread quantity in block).
Serial and concurrent Cloud motion wind inversion algorithm accuracy comparison 2 of the table 5 based on GPU
As can be seen that the parallel Cloud motion wind inverting based on GPU, serial inverting and parallel refutation from figure 6 above and Fig. 7
Obtain the similar numerical result of most of identical, small part, because GPU floating-point operation also has trueness error, therefore take
This result is obtained ideal, it is possible thereby to prove that the parallel method based on GPU is numerically credible and effective.
The Cloud motion wind product figure that parallel refutation algorithm based on GPU is finally inversed by is as depicted in figure 8.
Performance comparison is tested
Experiment is the thread CPU of 4 core 8 using host machine, and the GPU used is the tall and handsome NvidiaGeForce up to company
GT650M, core frequency is 925MHz, and video memory is 1G, and video memory frequency is 5400MHz, a width of 84.6GB/s of video memory band, stream process
Device quantity 768, theoretical calculation ability is 1.42TFLOPs.Respectively to defending in the case of different inverted parameters and thread block
Nebula figure carries out inverting, and algorithm is with respect to speed-up ratio and parallel efficiency change such as table 6 and table 7 with inverted parameters and thread number of blocks
It is shown.
Parallel Cloud motion wind inversion algorithm speed-up ratio and parallel efficiency 1 of the table 6 based on GPU
Parallel Cloud motion wind inversion algorithm speed-up ratio and parallel efficiency 2 of the table 7 based on GPU
By carrying out analysis to Fig. 9 and Figure 10 as can be seen that in the case of inverted parameters identical, thread block is followed substantially
Number is more, and speed-up ratio lifts bigger rule.Work as thread block<Parallel speedup ratio is basic into 2 times of increase trend when 64, when more than
Parallel speedup ratio is relevant with the division of thread block when 64, and two experiment highest speed-up ratios have reached 112 and 103, and acceleration effect is very
Significantly.
As shown in figure 11, maximum correlation coefficient (Maximum Cross-correlation Coefficient, MCC)
Principle can be expressed as choosing a spike cloud mass A in cloud atlas, then cloud mass A (N × N pixel image blocks) is in the mesh comprising A
Mark in field of search S (M × M pixel image blocks) and slide, calculate the coefficient correlation of all tracking cloud masses and spike cloud mass A in S, most
The big corresponding tracking cloud mass B of coefficient correlation is target cloud mass.
(Δ x, Δ y) calculation formula are such as shown in (2-1) by maximum correlation coefficient R：
Wherein,Respectively follow the trail of block and deviate spike cloud mass
Line number and columns, (i, j) is the index value of spike cloud mass, and A (i, j) and B (i, j) are respectively picture of the spike cloud mass with following the trail of block
Plain gray value,Respectively spike cloud mass and follow the trail of block average gray.
Obtain after maximum correlation coefficient, the tracking cloud mass B in target search area is relative to spike cloud mass A on x, y directions
Mobile distance is respectively Δ x and Δ y, if the longitude and latitude of spike cloud mass A central point is (x1,y1), follow the trail of cloud mass B center
The longitude and latitude of point is (x2,y2), target cloud mass B can be calculated relative to spike cloud mass A by setting up earth model (as shown in figure 12)
Apart from d, the speed of wind arrow is just can obtain with d divided by two images observation interval t；Derivation is as follows：
The speed of wind arrow
The wherein distance of A, B point-to-point transmission (RE is earth radius)
D=RE × γ (2-3)
2 points of A, B and the centre of sphere angle of centre of sphere formation can be calculated with space vector scalar product
If vectorialFor (a1,b1,c1), vectorFor (a2,b2,c2), then
a1=RE × cos (x1)×cos(y1) (2-5)
b1=RE × sin (x1)×cos(y1) (2-6)
c1=RE × sin (y1) (2-7)
Similarly
a2=RE × cos (x2)×cos(y2) (2-8)
b2=RE × sin (x2)×cos(y2) (2-9)
c2=RE × sin (y2) (2-10)
Formula (2-5)~(2-10) enters (2-4) can be obtained into centre of sphere angle γ size, the γ obtained is then substituted into (2-3)
Can obtain A, B point-to-point transmission apart from d, d finally is substituted into (2-2) just can obtain the speed υ of wind arrow.
The direction of wind arrow need to calculate derivation by the spherical triangle cosine law.As shown in figure 13, plane AED and the earth
Sphere ABC is tangential on the beginning and end that point A, A, B are respectively wind arrow.
Have to triangle ODE with the triangle cosine law
Because triangle OAE and triangle OAD are right angled triangles, therefore have
OE2=OA2+AE2 (2-11)
OD2=OA2+AD2 (2-12)
(2-12) (2-13) is substituted into (2-11) and had with the triangle cosine law
Cos (α)=cos (β) × cos (γ)+sin (β) × sin (γ) × cos (θ) (2-13)
Wherein
α=x2-x1 (2-14)
β=y2-y1 (2-15)
(2-15) (2-16) and the γ obtained above are substituted into (2-14), θ size can be obtained using arcsin function.
If A points northern area under the line, wind arrow deflection is γ, if x2＜ x1, then wind arrow deflection is 2 π-θ；If A points under the line with
Southern area, then wind arrow deflection is π-θ, if x2＜ x1, then wind arrow deflection is π+θ.
The height of wind arrow can be specified by infrared the one of wind arrow position and the bright temperature value of vapor channel.
Wind arrow calculating needs to carry out wind speed and direction verification after finishing, to reject wind speed and direction in primary Calculation result
The larger wind arrow of deviation.The cloud atlas for defining three continuous time observations is C1、C2、C3.If utilizing C1And C2The wind arrow CV calculated1's
Wind speed and direction is respectively υ1And θ1, utilize C2And C3Calculate wind arrow CV2Wind speed and direction be respectively υ2And θ2, then two wind arrows
The poor Δ θ of wind direction=| θ1-θ2|, the relative mistake Δ υ of wind speed=| υ1-υ2|, when the wind direction difference Δ θ or wind speed difference Δ υ of wind arrow are more than
During given threshold value, then this wind arrow is rejected.
In the present embodiment, the scene that maximum correlation coefficient is used is found out in the wind arrow place block of pixels for wanting inverting
The relation of the maximum 32*32 regions of heart 32*32 pixel regions similarity, similarity and coefficient correlation is monotonic increase, so
The position of target block can be found out using this algorithm, so as to calculate the correlation values of wind arrow, schematic diagram is as shown in figure 11.
In the algorithm, the coefficient correlation time calculated between two blocks is less, and the calculating of coefficient correlation array takes
Longer, if the size of block is N*N, the trace regions where block are M*M, and such as each coefficient correlation is calculated, then needed
(M-N+1) * (M-N+1) secondary coefficient correlation is calculated, the time is calculate a coefficient correlation time more than 1,000 times, therefore, this hair
It is bright this amount of calculation to be optimized.
In the present embodiment, " variable step iterative calculation " is used, variable step-size search method refers in trace regions S, spike
The moving step length of cloud mass by it is original be always that 1 strategy is changed to the strategy of variable step movement, the initial step length of search use compared with
Big value, progressively reduces step-length and reduces hunting zone simultaneously until step-length is 1, every kind of step-length is all found out in the case of current step afterwards
Optimum target matching cloud mass, step-length just have found approximate global optimal matching cloud mass B when being 1.
In this embodiment it is assumed that initial step length is 8, then there is following derivation：
First pass moving step length is 8, calculates (M/8) * (M/8) individual correlation coefficient value, finds out in these correlation coefficient value
Maximum value, so as to find similarity highest N*N block A1, method is as shown in figure 14.
Find after A1, to avoid error, the A1 blocks length of side is increased into step pixel, step-length is then reduced to 8/2
I.e. 4 pixels, i.e. step is 4, to the N*N regions of coefficient correlation maximum in the region of the N+step length of sides centered on A1
Scan for, the coefficient correlation of each N*N blocks and B blocks is calculated respectively, the coefficient matrix size calculated is (N+step)/4*
(N+step)/4, i.e., (N+step)2/ 16, it is the maximum in the case of 4 that this, which just represents these step numbers of iterative search with regard to that can find step,
Coefficient correlation block, searching algorithm is as shown in figure 15.
According to this algorithm, the mode that step progressively halves is successively decreased, the maximum finally searched when step step-lengths are 1
Coefficient correlation block Block be with center block B similarity highest blocks so that according to the positional information calculation of two blocks
Go out the speed and angle of wind arrow.
The former method that target cloud is followed the trail of, which finds optimal similar block, needs our calculating of the number of times of iterative search to understand,
StepCount1=(M-N+1)2
Optimal similar block is found after algorithm optimization needs the number of times of iterative search
StepCount2=(M/step)2+step/2+step/4+…+1
It is that 64, N is that in the case that 32, step is 8, former method needs the search of 1089 steps, after optimization to be revised as M in step-length
Method needs the search of 71 steps, greatly simplifies the time that search needs, can accelerate more than 15 times, and the precision in later stage in theory
Optimization can be directed to the scope searched again for after variable step and do the adjustment related to M sizes.
Although embodiment of the present invention is disclosed as above, it is not restricted in specification and embodiment listed
With it can be applied to various suitable the field of the invention completely, can be easily for those skilled in the art
Other modification is realized, therefore under the universal limited without departing substantially from claim and equivalency range, the present invention is not limited
In specific details and shown here as the legend with description.
Claims (8)
1. a kind of parallel Cloud-motion wind inversion method that cloud system is pre-processed based on fuzzy clustering, it is characterised in that including：
Cloud atlas image is divided into polylith and follows the trail of block, and multiple GPU parallel computations arrays for specifying are assigned to by block is followed the trail of,
It is laggard to the target cloud mass progress matching in the tracking block by maximum correlation coefficient in each independent computing array
Row Inversion Calculation, obtains the wind speed and wind direction of the tracking block wind arrow；
The matching process of target cloud mass is included by the maximum correlation coefficient：Target cloud is chosen in the tracking block
Block, in the region of search comprising the spike cloud mass, calculate step-size in search more than 1 all with the spike cloud mass pixel
The coefficient correlation of equivalent block, determines maximum correlation coefficient, the block with maximum correlation coefficient is in these coefficient correlations
For with the spike cloud mass similarity highest block, the length of side increase of the block is halved into step-length picture centered on the block
In the region of vegetarian refreshments, step-length is progressively halved and searches again for out similarity highest block, searched until searching out step-length for 1
Maximum correlation coefficient block, as target cloud mass；And
Wherein, original cloud atlas image is entered by fuzzy clustering method after separating with racking and determines the cloud for carrying out Cloud-motion wind inverting
Figure image, it includes：
The gray feature data of infrared channel in cloud atlas image are gathered, the gray feature data are handled, sample is obtained
Gray feature data, the sample gray feature data are classified, and make each class sample gray feature data right respectively
A kind of cloud atlas type is answered, and then determines the cloud atlas image.
2. the parallel Cloud-motion wind inversion method as claimed in claim 1 that cloud system is pre-processed based on fuzzy clustering, it is characterised in that
Also include, choose three width cloud atlas of adjacent constant duration, the wind speed and wind direction of every width cloud atlas wind arrow are drawn respectively, by the wind
The wind speed of arrow, wind direction and the wind arrow speed difference of adjacent cloud atlas, direction difference in the cloud atlas are compared with threshold value according to sequential
Compared with the following principle of execution：
If wind arrow speed is less than the minimum speed threshold value specified, this wind arrow is rejected；Or
If the difference of two wind arrow speed is more than the maximum speed discrepancy threshold value specified, this wind arrow is rejected；Or
If the difference in two wind arrow directions is more than the poor threshold value in the maximum direction specified, this wind arrow is rejected；Or
If there is no the poor wind arrow for being less than the poor threshold value in the maximum direction specified in direction therewith in the peripheral region of two wind arrows, this is rejected
Wind arrow.
3. the parallel Cloud-motion wind inversion method as claimed in claim 1 or 2 that cloud system is pre-processed based on fuzzy clustering, its feature is existed
In to cloud atlas view data progress multi-channel data extraction, reservation infrared channel data.
4. the parallel Cloud-motion wind inversion method as claimed in claim 3 that cloud system is pre-processed based on fuzzy clustering, it is characterised in that
Also include：The cloud atlas image is pre-processed, the longitude and latitude projective transformation such as it includes carrying out image and by original disk
Figure is transformed to histogram.
5. the parallel Cloud-motion wind inversion method as claimed in claim 4 that cloud system is pre-processed based on fuzzy clustering, it is characterised in that
The cloud atlas image, which is pre-processed, also to be included carrying out grey level enhancement to described image.
6. the parallel Cloud-motion wind inversion method as claimed in claim 1 that cloud system is pre-processed based on fuzzy clustering, it is characterised in that
Gray feature data progress processing is included：
The gray feature data input of infrared channel in cloud atlas image is obtained into pixel grey scale value matrix, passed through
Leading diagonal pixel grey scale average in a matrix is calculated, is passed throughCalculate auxiliary diagonal pixel ash in a matrix
Average is spent, and then passes through M=(M1+M2Obtain pixel grey scale average, i.e. sample gray feature data in)/2.
7. the parallel Cloud-motion wind inversion method as claimed in claim 6 that cloud system is pre-processed based on fuzzy clustering, it is characterised in that
The sample gray feature data of the first LONG WAVE INFRARED, the second LONG WAVE INFRARED and short-wave infrared are chosen in the infrared channel,
The sample gray feature data are set up into 5 gray feature data statistics amounts, image segmentation is carried out to Nephogram and obtains institute
State cloud atlas image.
8. the parallel Cloud-motion wind inverting side that cloud system is pre-processed based on fuzzy clustering as any one of claim 1,2,4-7
Method, it is characterised in that calculated according to the width height and wind arrow density of every width cloud atlas pixel and determine the tracking block number, it is
In (Width/ ρ) × (Height/ ρ), formula, Width is the width of every width cloud atlas pixel, and Height is every width cloud atlas
The height of pixel, ρ is wind arrow density.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610642651.XA CN106340004B (en) | 2016-08-08 | 2016-08-08 | A kind of parallel Cloud-motion wind inversion method that cloud system is pre-processed based on fuzzy clustering |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610642651.XA CN106340004B (en) | 2016-08-08 | 2016-08-08 | A kind of parallel Cloud-motion wind inversion method that cloud system is pre-processed based on fuzzy clustering |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN106340004A CN106340004A (en) | 2017-01-18 |
| CN106340004B true CN106340004B (en) | 2017-09-01 |
Family
ID=57825212
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610642651.XA Expired - Fee Related CN106340004B (en) | 2016-08-08 | 2016-08-08 | A kind of parallel Cloud-motion wind inversion method that cloud system is pre-processed based on fuzzy clustering |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN106340004B (en) |
Families Citing this family (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108647692B (en) * | 2018-04-11 | 2021-09-28 | 中国海洋大学 | Ocean layer knot extraction method based on LH histogram |
| CN108764563B (en) * | 2018-05-25 | 2021-04-02 | 国网湖南省电力有限公司 | Squall line wind early warning method for power transmission line |
| CN110428453B (en) * | 2019-07-30 | 2020-12-15 | 深圳云天励飞技术有限公司 | Data processing method, data processing device, data processing equipment and storage medium |
| CN111175720B (en) * | 2020-01-15 | 2022-03-08 | 中国科学院国家空间科学中心 | Method and system for quickly inverting on-board sea surface wind field |
| CN111598802B (en) * | 2020-05-12 | 2023-04-25 | 中国科学院合肥物质科学研究院 | Foundation all-sky cloud parameter inversion system and method |
| CN112347956B (en) * | 2020-11-12 | 2022-05-06 | 上海交通大学 | Cloud observation system and method based on multiple unmanned aerial vehicles and machine vision |
| CN113315803A (en) * | 2021-03-23 | 2021-08-27 | 中国地质大学（北京） | Short wave radiation mode data transmission performance optimization method based on nail-fixed memory |
| CN115408301B (en) * | 2022-10-31 | 2023-01-31 | 北京千尧新能源科技开发有限公司 | Test set construction method and system for fan simulation |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103308031B (en) * | 2013-05-23 | 2015-06-03 | 中国人民解放军理工大学 | Cloud top height retrieval method based on satellite tri-linear array CCD (charge coupled device) image |
| CN105046056B (en) * | 2015-06-24 | 2018-01-16 | 北京航天宏图信息技术股份有限公司 | The method of inspection and device of Cloud motion wind data |
-
2016
- 2016-08-08 CN CN201610642651.XA patent/CN106340004B/en not_active Expired - Fee Related
Also Published As
| Publication number | Publication date |
|---|---|
| CN106340004A (en) | 2017-01-18 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106340004B (en) | A kind of parallel Cloud-motion wind inversion method that cloud system is pre-processed based on fuzzy clustering | |
| CN109829399B (en) | Vehicle-mounted road scene point cloud automatic classification method based on deep learning | |
| CN103839261B (en) | SAR image segmentation method based on decomposition evolution multi-objective optimization and FCM | |
| CN104778476B (en) | A kind of image classification method | |
| CN102005034A (en) | Remote sensing image segmentation method based on region clustering | |
| CN101699514B (en) | Immune clone quantum clustering-based SAR image segmenting method | |
| Saleem et al. | Comparative analysis of recent architecture of convolutional neural network | |
| Chen et al. | The mixed kernel function SVM-based point cloud classification | |
| Hu et al. | A comparison and strategy of semantic segmentation on remote sensing images | |
| CN1636168A (en) | Clustering for data compression | |
| CN111797833A (en) | Automatic machine learning method and system oriented to remote sensing semantic segmentation | |
| CN115082790A (en) | Remote sensing image scene classification method based on continuous learning | |
| CN110647977B (en) | Method for optimizing Tiny-YOLO network for detecting ship target on satellite | |
| Marti Asenjo et al. | MRI brain tumor segmentation using a 2D-3D U-Net ensemble | |
| Alburshaid et al. | Palm trees detection using the integration between gis and deep learning | |
| CN111611960A (en) | Large-area ground surface coverage classification method based on multilayer perceptive neural network | |
| Hossam et al. | Accelerated hyperspectral image recursive hierarchical segmentation using GPUs, multicore CPUs, and hybrid CPU/GPU cluster | |
| Wang et al. | Cross-domain learning using optimized pseudo labels: toward adaptive car detection in different weather conditions and urban cities | |
| Chen et al. | Mapping urban form and land use with deep learning techniques: a case study of Dongguan City, China | |
| Song | A More Efficient Approach for Remote Sensing Image Classification. | |
| CN115797663A (en) | Space target material identification method under complex illumination condition | |
| CN106485686A (en) | One kind is based on gravitational spectral clustering image segmentation algorithm | |
| Shi et al. | Intelligent classification of land cover types in open-pit mine area using object-oriented method and multitask learning | |
| CN112991402B (en) | Wen Wudian cloud registration method and system based on improved differential evolution algorithm | |
| CN112308151A (en) | Weighting-based classification method for hyperspectral images of rotating forest |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20170901 Termination date: 20210808 |