CN106340004A - Fuzzy clustering preprocessing cloud system-based parallel cloud drift wind inversion method - Google Patents

Abstract

The invention discloses a fuzzy clustering preprocessing cloud system-based parallel cloud drift wind inversion method. According to the method, three adjacent equal-time interval cloud charts are selected; image extraction and preprocessing are carried out on the images of the cloud charts; each cloud chart is divided into a plurality of tracking blocks; the tracking blocks are assigned to a plurality of designated parallel computing processing units; and inversion calculation is performed on the wind vector of the corresponding tracking block is each independent computing processing unit through using a maximum correlation coefficient method, so that the wind speed and wind direction of the wind vector of the tracking block can be obtained. With the method of the invention adopted, the actual execution efficiency of the inversion of cloud drift wind can be improved effectively. The method is of great application value for the prediction and early warning of global wind fields at low, medium and high altitudes.

Description

A kind of parallel Cloud-motion wind inversion method based on fuzzy clustering pretreatment cloud system

Technical field

The present invention relates to satellite derived wind inverting field is and in particular to a kind of parallel based on fuzzy clustering pretreatment cloud system Cloud-motion wind inversion method.

Background technology

In the research of weather and damage forecasting, need substantial amounts of Wind Data, the observation of global atmosphere wind field serves Very positive effect.In recent years, with the constantly improve of land wind field observation network, research worker obtains land Wind Data 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 Very not enough, obtain these regional Wind Datas more difficult.By satellite remote dada inverting atmospheric wind, undoubtedly to acquisition Global atmosphere Wind Data serves the supplementary function being highly profitable.

Meanwhile, some land or water body (referred to as earth's surface) due to overhead cloudless or partly cloudy and direct " exposed " is in satellite cloud picture In.Research worker is mainly studied to the cloud in cloud atlas, and the such as identification of cloud classification, cloud, Cloud motion wind etc. therefore realize cloud sector The image segmentation of domain and earth surface area is requisite.Research work for satellite cloud picture image segmentation launches already, mesh The method of front use generally has the disadvantage that (1) is overly dependent upon cloud atlas image document over the years, needs to spend manpower and materials Collect and preserve；(2) realized using the image segmentation 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 comprises cloud and earth's surface, there are some and be in " transition rank also between cloud and earth's surface The cloud (cloud for example growing and dissipating) of section "；So, in conventional methods where, easily forced to be divided into a certain apoplexy due to endogenous wind, but this It is irrational, be also one of direction of research to the rationally segmentation of cloud atlas image.

In general, satellite cloud picture time interval 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 all proposes challenge to the efficiency of Cloud motion wind inversion algorithm. Although there being scholar to propose the simple algorithm of Cloud motion wind inversion algorithm, being also difficult to meet and increasingly require the wind field of real-time to calculate Efficiency requirements.Therefore, except the inversion algorithm studying relatively low calculation cost, make full use of the architecture of modern computer, grind Study carefully parallel Cloud motion wind inversion algorithm, be also an up one of effective way of inversion algorithm actual operating efficiency.

Gpu (graphic processing unit), graphic process unit is video card " heart ", is equivalent to cpu in microcomputer On effect.Gpu has very high transmission bandwidth and numerous computing units, and computing capability is very powerful.Continuous with gpu Update upgrading, the computing capability of gpu is more come outstanding, and this is that the inverting of Cloud-motion wind provides and a kind of parallel accelerates thinking；Gpu is in height Performance computing aspect has an advantage of three highly significants: first, the parallel processing capability of data is powerful, and nvidia company is current , more than 3tflops, this almost with small-sized cpu cluster is suitable for up-to-date gpu peak value floating-point operation ability；Second, gpu have Very high data throughput capabilities, the band block of k20 is more than 200gb/s；3rd, gpu have the DLL of close friend, nvidia company The cuda framework of exploitation supports multiple high-level programming languages, such as c, c++ and fortran language.

To sum up although in the prior art, there is the satellite Cloud motion wind parallel refutation algorithm based on multinuclear cpu for the article elaboration, But because cpu is good at the task with complex logic calculation procedure that processes, such as data compression, artificial intelligence, physical simulation etc., and Gpu is good at process computation-intensive and the better simply repetitive task of logic, such as image domains relevant treatment operations, and both are in thing Reason structure data is processed has very big difference, and on gpu, applied satellite Cloud-motion wind inversion method is to need to overcome a lot Difficulty and technology barrier, although gpu has hundreds and thousands of computing units, the computing capability of each computing unit is weak, difficult To execute 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 less, and it is real that this brings theory and technology to the transplanting of gpu parallel algorithm to the complicated algorithm on traditional cpu The difficulty now going up.How using the powerful computing capability of gpu, the satellite derived wind inverting to high density, high real-time is carried out simultaneously Row process is a challenge.

Content of the invention

The present invention has designed and developed a kind of parallel Cloud-motion wind inversion method based on fuzzy clustering pretreatment cloud system, the present invention One of purpose be Cloud-motion wind inverse time length, the reality solving to lead to because iterationses are more in serial Cloud-motion wind inversion method The problem of when property difference, by parallel refutation method by single iteration parallelization, reduces iterationses, and then significantly improves Cloud-motion wind The real-time of inverting.

The two of the goal of the invention of the present invention be solve regulation region of search in long to target cloud mass search time The problem of the Cloud-motion wind inverting poor real leading to, by optimizing the method for tracing of spike cloud mass, significantly simplifies 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 to the cloud atlas Image semantic classification needing to carry out inverting only to solve in prior art Using the hard dividing method of image, so exist to being difficult to during Image semantic classification obtain a suitable threshold value, and By such environmental effects, the problem that accuracy rate is difficult to ensure that.

The technical scheme that the present invention provides is:

A kind of parallel Cloud-motion wind inversion method based on fuzzy clustering pretreatment cloud system, comprising:

Cloud atlas image is divided into polylith to follow the trail of block, and is assigned to, by following the trail of block, the multiple gpu parallel computations specified Array, is mated to the target cloud mass in described tracking block by maximum correlation coefficient in each independent computing array After carry out Inversion Calculation, obtain described following the trail of the wind speed of block wind arrow and wind direction；

By described maximum correlation coefficient, the matching process of target cloud mass is being included: choose in described tracking block Spike cloud mass, in the region of search comprising described spike cloud mass, calculates the whole and described target cloud that step-size in search is more than 1 Block pixel is equal to the correlation coefficient of block, determines maximum correlation coefficient, as with described spike cloud mass in these correlation coefficienies Similarity highest block, is increased centered on this block in the region of pixel, 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 searching, as target cloud mass；And

Wherein, determine after separating with by fuzzy clustering method, original cloud atlas image being entered to rack and carry out Cloud-motion wind inverting Cloud atlas image, comprising:

The gray feature data of collection cloud atlas image mid-infrared passage, is processed to described gray feature data, obtains Sample gray feature data, described sample gray feature data is classified, and so that each class sample gray feature data is divided Not Dui Ying a kind of cloud atlas type, and then determine described cloud atlas image.

Preferably, also include, choose three width cloud atlas of adjacent constant duration, draw the wind of every width cloud atlas wind arrow respectively Speed and wind direction, by the wind speed of described wind arrow, wind direction and according to sequential by the wind arrow speed difference of cloud atlas adjacent in described cloud atlas, side It is compared to difference with threshold value, the following principle of execution:

If wind arrow speed is less than the minimum speed threshold value specified, reject this wind arrow；Or

If the difference of two wind arrow speed is more than the maximum speed discrepancy threshold value specified, reject this wind arrow；Or

If the difference in two wind arrow directions is more than the maximum direction difference limen value specified, reject this wind arrow；Or

If not having direction difference therewith less than the wind arrow of the maximum direction difference limen value specified, to pick in the peripheral region of two wind arrows Except this wind arrow.

Preferably, described cloud atlas view data is carried out with multi-channel data extraction, retains infrared channel data.

Preferably, also include: pretreatment is carried out to described cloud atlas image, it includes image being carried out wait longitude and latitude projection Convert and original disk is attempted to change and be changed to histogram.

Preferably, described cloud atlas image is carried out with pretreatment also to include carrying out grey level enhancement to described image.

Preferably, described gray feature data is carried out processing including:

The gray feature data input of cloud atlas image mid-infrared passage is obtained pixel grey scale value matrix, passes throughCalculate leading diagonal pixel grey scale average in a matrix, pass throughCalculate auxiliary in a matrix Diagonal pixels gray average, and then pass through m=(m₁+m₂Pixel grey scale average is obtained in)/2, i.e. sample gray feature data.

Preferably, the first LONG WAVE INFRARED, the second LONG WAVE INFRARED and short-wave infrared are chosen in described infrared channel Sample gray feature data, described sample gray feature data is set up 5 gray feature data statisticss amounts, Nephogram is entered Row image segmentation obtains described cloud atlas image.

Preferably, calculated according to the wide high of every width cloud atlas pixel and wind arrow density and determine described tracking block number, It is (width/ ρ) × (height/ ρ), and in formula, width is the width of described every width cloud atlas pixel, and height is described every width The height of cloud atlas pixel, ρ is wind arrow density.

Preferably, by described maximum correlation coefficient, the matching process of target cloud mass is being included: in described maximum In correlation coefficient process, choose spike cloud mass a in described tracking block, it is n × n-pixel image block, is comprising described spike In m × m pixel image block region of search of cloud mass a, calculate whole and described spike cloud mass a pixels that step-size in search is m and be equal to The correlation coefficient of block, determines maximum correlation coefficient, as with described spike cloud mass a similarity highest in these correlation coefficienies N × n-pixel image block a1, increases in described block a1Then step-length is kept to by individual pixelIndividual pixel, in Yi Qu Centered on block a1The correlation coefficient of each n × n-pixel image block is calculated, in these correlation coefficienies in the region of the length of side Middle determination maximum correlation coefficient so determine next with described spike cloud mass a similarity 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 searching, as target cloud mass, by institute The positional information calculation stating spike cloud mass and described target cloud mass obtains wind speed and the wind direction of wind arrow.

The present invention is had the advantages that compared with prior art

1st, the cascading algorithm in satellite derived wind inversion algorithm is changed to parallel algorithm, iterationses greatly reduces, By the computational methods of this single iteration parallelization, so that the wind arrow of whole Zhang Yun's in figure concurrently can execute in gpu, enter The wind arrow computing of the extensive quantity of row, in the case of not affecting the precision of Cloud motion wind product that is finally inversed by and performance, saves significantly Save the time needed for inverting and wind arrow verification, and then wind arrow computational efficiency has been greatly improved, significantly improve Cloud-motion wind The real-time of inverting；

2nd, in the present invention, by the way of variable step iterative calculation, 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, cloud atlas image segmentation pretreatment is applied to by fuzzy clustering method and can obtain more preferable cloud Ground separating effect, the method can effectively carry out the uncertain segmentation of cloud atlas and lifting segmentation accuracy rate, cloud separate after Carry out Cloud-motion wind inverting so that efficiency of inverse process is more accurate, simultaneously before carrying out fuzzy clustering pretreatment, to gray scale on image 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 producing in the calculating to Cloud-motion wind and identification work, further carries High efficiency of inverse process.

Brief description

Fig. 1 is that task in the satellite derived wind inversion method based on gpu platform for the cpu and gpu divides schematic diagram.

Fig. 2 is the fcm result figure of cloud atlas image.

Fig. 3 is the cloud atlas image segmentation flow chart based on fcm.

Fig. 4 basic flow sheet after parallel for Cloud-motion wind inversion algorithm.

Fig. 5 is the impact figure to efficiency for the stream handle quantity.

Fig. 6 is the schematic diagram based on gpu Cloud motion wind serial and concurrent inversion algorithm accuracy comparison 1.

Fig. 7 is the schematic diagram based on gpu Cloud motion wind serial and concurrent inversion algorithm accuracy comparison 2.

Fig. 8 is the Cloud motion wind product figure being finally inversed by based on gpu parallel Cloud motion wind algorithm.

Fig. 9 is the trendgram that during execution, personal attendant's thread block increases.

The trendgram that Figure 10 increases with thread block for speed-up ratio.

Figure 11 follows the trail of schematic diagram for spike cloud mass.

Figure 12 is wind arrow speed derivation schematic diagram.

Figure 13 is wind arrow direction derivation schematic diagram.

The iterative search schematic diagram that Figure 14 is 8 for step-length.

The iterative search schematic diagram that Figure 15 is 4 for step-length.

Specific 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 description literary composition 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 based on fuzzy clustering pretreatment cloud system Method, comprises the steps:

Step one, using adjacent constant duration three width infrared channel cloud atlas, set inverted parameters, including satellite remote sensing Read path on disk for the data, the density of Cloud-motion wind inverting, Cloud-motion wind inversion algorithm relevant parameter, gpu parallel parameter, cloud Mark wind verification dynamics parameter, Cloud-motion wind retrieval products storing path on disk etc.；

Step 2, using lossless decompression algorithm, the earth remote sensing data issuing from meteorological satellite load is decompressed Contracting；

Step 3, to satellite infrared image data extraction, satellite remote sensing date is decompressed with the satellite that flow process produces many Channel data carries out data extraction, retains infrared channel data, gives up other channel datas, and open up this high depositing of cloud atlas width Storage space, by the cloud atlas data duplication of the infrared channel extracting to this memory space；

Step 4, to extract after cloud atlas data carry out pretreatment, include image projection transformation, gradation of image enhancing and Determine cloud atlas image after fuzzy clustering method separates to the cloud atlas image extracting with entering to rack；

Step 5, performing environment is initialized, initial including the performing environment to inverted parameters and gpu equipment Change, comprising: calculating and need the wind arrow number of inverting according to the wide high of described every width cloud atlas pixel and wind arrow density, in gpu Video memory in open up Array for structural body for preserving Cloud-motion wind, open up in the video memory of gpu for preserving each wind arrow using Big correlation coefficient and the array of correlation coefficient value that produces, open up for storing 3 width infrared cloud image data in the video memory of gpu Array, cpu copies to three width cloud atlas data in hosting in three arrays that gpu video memory is opened up successively.

Step 6, carry out Cloud-motion wind parallel refutation calculating, series of parameters and multi-task planning and scheduling rule are carried by cpu Enter gpu equipment, and trigger the parallel refutation of the corresponding wind arrow of each thread of gpu to calculate, every width cloud atlas is divided into polylith to follow the trail of Block, and it is assigned to, by following the trail of block, the multiple parallel computation processing units specified, in each independent calculation processing unit By maximum correlation coefficient to described follow the trail of block in target cloud mass mate after carry out Inversion Calculation, obtain described in chase after The wind speed of track block wind arrow and wind direction, and be saved in the video card of gpu；

Wherein, the width of every width cloud atlas pixel, height and wind arrow density calculate and determine described tracking block number, and it is (width/ ρ) * (height/ ρ), in formula, width is the width of described every width cloud atlas pixel, and height is described every width cloud atlas The height of pixel, ρ is wind arrow density.

Step 7, quality indicator that Cloud-motion wind is carried out, 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 front two Zhang Yun's figures and the result that is finally inversed by of rear two Zhang Yun's figures, will The wind speed of wind arrow, wind direction and the wind arrow speed difference of adjacent cloud atlas in three width cloud atlas, direction difference are carried out with threshold value according to sequential Relatively, by the execution of following principle:

If wind arrow speed is less than the minimum speed threshold value specified, reject this wind arrow；Or

If the difference of two wind arrow speed is more than the maximum speed discrepancy threshold value specified, reject this wind arrow；Or

If the difference in two wind arrow directions is more than the maximum direction difference limen value specified, reject this wind arrow；Or

If not having direction difference therewith less than the wind arrow of the maximum direction difference limen value specified, to pick in the peripheral region of two wind arrows Except this wind arrow.

Step 8, data preservation is carried out to above-mentioned Cloud-motion wind product, obtain final retrieval products.

In another kind of embodiment, in step one, the remotely-sensed data read path in inverted parameters and product 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 comprised；Inverted parameters Include spike cloud mass and the number of pixels of trace regions ranks；Gpu parallel parameter is the net of the performance element of gpu kernel function Lattice division rule；Cloud-motion wind verification dynamics parameter is Cloud-motion wind credit rating, credit rating higher it is desirable to stricter, be finally inversed by Wind arrow with a high credibility, but number is fewer, otherwise quality is low, and number is many.

In another kind of embodiment, in step 2, satellite remote sensing date decompression algorithm can be using based on Run- Length Coding Decompression algorithm, the lzw decompression algorithm based on dictionary encoding technology, the decompression algorithm based on Huffman encoding principle or The decompression algorithm based on arithmetic coding for the person.

In another kind of embodiment, in step 4, to the longitude and latitude projective transformation such as image is carried out, it is transformed into Mercator and throws Shadow, to carry out wind arrow inverting, original for satellite disk being attempted to change and is changed to histogram, facilitate division and the execution of gpu task, to carry Rise inverting quality；

Image is carried out with grey level enhancement and processes and algorithm, piecewise linear gray transformation enhancing can be strengthened using linear gradation conversion Algorithm, histogram equalization strengthen algorithm or local histogram equalization strengthens algorithm and the contrast of image is strengthened；

In another kind of embodiment, in step 4, by Fuzzy c-means Clustering (fcm), original cloud atlas image is entered Described cloud atlas image is determined, it comprises the steps: after separating with racking

Input: satellite cloud picture image

Output: eliminate the cloud atlas image of earth surface area

Step a, first cloud atlas is carried out with pretreatment, the method for process is using median filtering method thus removing cloud atlas image In noise and interference.

Step b, the target data set of structure fcm.Choose infrared cloud image and visible cloud image builds ir-vis two-dimension spectrum Feature space, stores each pixel using a two-dimensional array data [256,2] in infrared, visible images respectively by row Gray feature value.

Step c, calculating sample general featuress value.If not having in data base to record 8 kinds of earth's surfaces, gray features of the varieties of clouds, Cloud atlas Sample Storehouse is analyzed, counts 8 kinds of earth's surfaces, the gray feature of the varieties of clouds and record in data base.If having, ignore This step, continues executing with.

Step d, initial cluster center is set to sample characteristics and carries out fcm.8 class earth's surfaces, cloud is obtained from data base The gray feature of class is as the initial cluster center of fcm, and is clustered in Infrared-Visible two-dimension spectrum feature space.

Step e, image segmentation is carried out to cloud atlas.Using the cluster block plan being obtained by step (4), all pixels are clicked through Row judges, if its two dimensional gray value is in earth's surface gray feature region, for earth's surface pixel, is labeled as black；If not existing, For cloud pixel, retain.

Finally obtain the cloud atlas after a rectangle, grey level enhancement and inactive area proposes.

In another kind of embodiment, in step 6, series of parameters and multi-task planning and scheduling rule are loaded into by cpu Gpu equipment, and trigger the parallel refutation of the corresponding wind arrow of each thread of gpu and calculate, all wind arrows of gpu executed in parallel anti- Drill, the inverting of each wind arrow executes in independent thread, the inverting of wind arrow uses maximum correlation coefficient to calculate target cloud The location of pixels that wind arrow is located, so that it is determined that wind arrow speed and direction, is finally converted to by block and the position relationship following the trail of cloud atlas Actual longitude and latitude position.

In another kind of embodiment, in step 6, Cloud motion wind inverting parallel algorithm adopts the fork- of process-thread Join model, after host process reads cloud atlas data and parameter and generates network subdivision relationships between nodes, is referred to by openmp compiling Lead sentence, the computational load of n wind arrow is dynamically assigning on m thread, using maximum correlation coefficient coupling target cloud mass, Calculate the wind speed and direction of wind arrow, the fixed step such as height, wind arrow quality indicator of wind arrow is independently calculated by each thread, reuses maximum phase Close Y-factor method Y coupling target cloud mass, in the overall wind arrow array step of write, each thread shares cloud atlas datarams region and wind arrow Array region of memory. all wind arrows calculate after finishing, and instruct sentence end lines journey by compiling, and host process is by the wind arrow calculating Data write disk preserves.

In another kind of embodiment, in step 6, maximum correlation coefficient refers to seek in satellite cloud picture trace regions Look for and follow the trail of cloud mass similarity highest cloud mass, the similarity between cloud mass then by correlation coefficient as index, chase after in cloud atlas Find and follow the trail of the maximum cloud mass of the correlation coefficient that calculates of cloud mass in track region, specifically include following step:

(1) cloud atlas 1 is divided into multiple tracking blocks, successively each tracking block of traversal cloud；

(2) in current gpu thread distribution wind arrow inverting block in centre of location position pixel region, in record Heart point position；

(3) calculate the gray variance in region and standard deviation in (2)；

(4) first pixel region in the upper left corner of the same block from cloud atlas 2 begins stepping through, until traversal knot Bundle, returns (1), start to 0, under from left to right travel through, until traversal terminates, calculate variance and the standard of this same block Difference；

(5) calculate the correlation coefficient of (3) and (4), and write the result into storage matrix；

(6) find out the maximum correlation coefficient in the storage matrix of (5), record its correspondence position；

(7) position coordinateses according to (2) and (6), calculate wind arrow distance and the angle of this region unit, are removed with the distance of wind arrow With time interval, then can obtain wind arrow speed；The form that the angle of the wind arrow being finally inversed by is converted into northern drift angle is wind arrow Direction.

Specifically, in step 6, by described maximum correlation coefficient, the matching process of target cloud mass is being included: In described maximum correlation coefficient, choose spike cloud mass a in described tracking block, it is n × n-pixel image block, in bag In m × m pixel image block region of search containing described spike cloud mass a, calculate the whole and described spike cloud mass that step-size in search is m A pixel is equal to the correlation coefficient of block, determines maximum correlation coefficient, as with described spike cloud mass a in these correlation coefficienies Similarity highest n × n-pixel image block a1, increases in described block a1Then step-length is kept to by individual pixelIndividual picture Vegetarian refreshments, centered on block a1The correlation coefficient of each n × n-pixel image block is calculated in the region of the length of side, Determine maximum correlation coefficient in these correlation coefficienies and then determine next and described spike cloud mass a similarity 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 searching, as target Cloud mass, obtains wind speed and the wind direction of wind arrow by the positional information calculation of described spike cloud mass and described target cloud mass.

In another kind of embodiment, in described step 8, Cloud-motion wind product storage method is cpu by the video memory of gpu Wind arrow data duplication in computer hosting, then output an external memory disk in.

In another kind of embodiment, the data object of inverting of the present invention is the cloud atlas data of fy-2e satellite infrared one passage, The size of this cloud atlas is 2581*1399 pixel, and that is, width is 2581 pixels, and height is 1399 pixels, chooses constant duration The cloud atlas of three adjacent fy-2e satellite infrared one passages is as the data processing object of the present embodiment.

In another kind of embodiment, inverted parameters in the present invention receive flow process, 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 the spike cloud mass in input reverse Long, that is, target cloud block size is p*p；Input Cloud-motion wind inverting grade d.

In another kind of embodiment, satellite remote sensing date in the present invention decompresses flow process, according to cloud provided above Figure read path path1, reads the remotely-sensed data of fy-2e satellite, using lossless decompression algorithm to fy-2e satellite remote sensing number 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 kind of embodiment, satellite infrared one image data extraction flow process in the present invention, decompress from above-mentioned To each passage cloud atlas data in, extract the satellite infrared one cloud atlas data of three Zhang Yun's figures respectively, copy to microcomputer host open In array data1, data2 and data3 warded off.

In another kind of embodiment, gradation of image enhancement process is contrasted to cloud atlas data using Nogata equalization algorithm Degree strengthens.

In another kind of embodiment, performing environment in the present invention initializes flow process, the execution of initialization cpu and gpu Environment and internal storage data.Calculated according to density and the wide high of cloud atlas, need Cloud-motion wind number n of inverting；Cpu calls cuda interface The n Array for structural body windrecords for preserving Cloud-motion wind is opened up in the shared drive of gpu；Cpu calls cuda interface Open up n size to produce using maximum correlation coefficient for preserving each wind arrow for (p+1) * (p+1) in the shared drive of gpu The array of the raw individual correlation coefficient value of (p+1) * (p+1)；Cpu calls cuda interface to open up three use in the shared drive of gpu In the array of three infrared cloud atlas data of storage, the size of array, according to claim 4, is the size of cloud atlas width*height；Cpu calls cuda interface that three width cloud atlas data in hosting are copied to gpu equipment successively in step 1 In three arrays that shared drive is opened up.

In another kind of 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 in a mission thread, all thread dividing are executed in grid to a gpu, gpu divides according to the task of thread lattice Join and dispatching principle, the inverting to all wind arrows concurrently executes, inverting determines the displacement of wind arrow using maximum correlation coefficient, enters And determine velocity magnitude and the direction of wind arrow.Finally the location of pixels that wind arrow is located is converted to actual longitude and latitude position, will be anti- Drill the wind arrow obtaining and write, by the index position that individual execution thread is distributed, the windrecords array opened up in gpu shared drive In.

In another kind of embodiment, preserve flow process in the Cloud-motion wind product of the present invention, the shared drive of gpu is fallen into a trap by cpu Cloud-motion wind data duplication after the quality indicator calculating in the main memory of host's microcomputer, then cpu by host in Cloud-motion wind number According to being preserved, storing path is that above-mentioned inverted parameters receive the Cloud-motion wind product storing path path2 that flow process receives, preservation Form is: wind arrow numbering, wind arrow longitude, wind arrow latitude, wind arrow speed and wind arrow direction.

In conjunction with specific embodiment, specific description is done to parallel refutation method and maximum correlation coefficient below.

In another kind of embodiment, as shown in Figure 2 and Figure 3, divide with density clustering method, cloud atlas being entered to rack From；According to " similarity " of gray feature between cloud atlas spectral digital image slices vegetarian refreshments pixel is carried out division 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 it is subject to such environmental effects, accuracy rate is difficult to ensure that.The method can effectively carry out the uncertain segmentation of cloud atlas and lifting point Cut accuracy rate, cloud separate after image on carry out Cloud-motion wind inverting so that efficiency of inverse process is more preferably accurate；In satellite digital figure During the generation of picture, process and projective transformation etc., cloud atlas generally all unavoidably produces some salt-pepper noises, and these noises are not But reduce picture quality, also make objective fuzzy, the calculating to Cloud-motion wind and identification work create very big interference, therefore exist It is necessary for removing noise and reduction interference in satellite cloud picture image.

Noise is removed by median filter method to gray scale characteristic in satellite cloud picture image and reduces at interference Reason obtains gray scale sample characteristics data, specific as follows:

Use box filter masterplate, define following two dimension median filter formula:

z_xy=med { g (i, j) }=med { g (x+r, y+s), (r, s) ∈ a (x, y) ∈ i²}

Wherein, med represents and takes median operation；During process, the gray value of p (x, y) is set to z_xy, as shown in table 1.

Table 13 × 3 medium filtering module

	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 gray value of each point in template is ranked up, according to row or column, the pixel in template is divided Three groups (every group of 3 pixels) are become to carry out median calculation, the intermediate value obtaining is sample gradation data feature, filters using Quick Median Wave method can reduce the number of comparisons between pixel, and can pass through Parallel Implementation, is greatly enhanced the place of image denoising Reason speed, Fast Median Filtering method specifically includes:

Input: the grey scale pixel value in 3 × 3 windows

Output: the gray scale intermediate value of 9 pixels

(1) to the every a line in template, ask for maximum, minima and intermediate value respectively, obtain data set a=[max0, Max1, max2], b=[min0, min1, min2], c=[med0, med1, med2], as follows.

A:max0=max [p₀, p₁, p₂], max1=max [p₃, p₄, p₅], max2=max [p₆, p₇, p₈]

B:min0=min [p₀, p₁, p₂], min1=min [p₃, p₄, p₅], min2=min [p₆, p₇, p₈]

C:med0=med [p₀, p₁, p₂], med1=med [p₃, p₄, p₅], med2=med [p₆, p₇, p₈]

(2) seek the minima in a: maxmin=min [max0, max1, max2]；

(3) seek the intermediate value in c: medmed=med [med0, med1, med2]；

(4) seek the maximum in b: minmax=max [min0, min1, min2]；

(5) med [maxmin, medmed, minmax], as required intermediate value are asked.

In theory, we only need 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, therefore by mistake 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 It is divided into 8 classes, i.e. land, water body, cumulonimbus, cumulus congestus, cirrus, medium cloud, cumulus humilis and low clouds, therefore herein fcm is clustered 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.Greatly reduce the number of times of circulation, so that the effect of algorithm is more notable.

In the present embodiment, choose the sample gradation data feature of infrared (ir) and steam (wv) two passage cloud atlas, Set up the gray feature statistic (averages of some samples) of this 8 class sample, as shown in table 2.

The gray feature statistics of table 2 cloud atlas sample

Determine that clusters number is 8, and the initial center that sample gray feature statistic is clustered as fcm, constraining and Instruct cluster process, in cluster iterative process, continuous adjusting and optimizing cluster centre and subordinated-degree matrix, finally obtain two dimensional character The cluster result in space.

Concrete fcm algorithm is as follows:

Input: clusters number c=8, data bulk n；Fixed value λ, 1≤λ≤∞；Threshold epsilon；

Output: subordinated-degree matrix u'=[u_ij]_n×n

(1) average of earth's surface and typical varieties of clouds characteristic quality of sample is assigned to v₀=[v_ij]_2×8；

(2) make cycle-index b=0, be subordinated-degree matrix u using formula in step (5)⁽⁰⁾Give initial value；

(3)b++；

(4) calculate c cluster centre v_i ^(b)(i=1,2 ..., c), formula such as (1-1).

$v_i=\left(\underset{j=1}{\overset{n}{\σ}}{u}_{ij}^{\mathrm{\λ}}{x}_j\right)/\left(\underset{j=1}{\overset{n}{\σ}}{u}_{ij}^{\mathrm{\λ}}\right)---(1-1)$

(5) for j=1～n, recalculate u^(b), and make u^(b+1)=u^(b)；

1. calculate i_jWith

i_j=i | 1≤i≤c；d_ij=| | x_j-v_i| |=0 } (1-2)

$\stackrel{\&overbar;}{i_j}=\{1,2,...,c\}-i_j---(1-3)$

2. for data x_j, update its degree of membership.IfThen to allMake u_ij=0,j =j+1；Otherwise

$u_{ij}=1/\underset{k=1}{\overset{c}{\σ}}{(d_{ij}/d_{kj})}^{2(\mathrm{\λ}-1)}---(1-4)$

(6) compare u^(b)And u^(b+1).IfExecution step (7)；Otherwise return (3)；

(7) make u '=u^(b+1)；Terminate.

Fy-2e meteorological satellite has visible ray, infrared and three spectrum cloud atlas of steam, each pixel therefore in cloud atlas There are visible ray, infrared and three kinds of gray feature attributes of steam；Fcm existing more intuitively effect in two dimensional surface space has again Higher execution efficiency, in the present embodiment, infrared and visible ray is made a concrete analysis of from three kind Sexual behavior mode, using infrared logical Road and visible channel build orthogonal two-dimension spectrum feature space, make point in pixel and two-dimensional feature space on image one by one Mapping, that is, cloud atlas pixel is r_i(m,n),r_v(m, n), two-dimension spectrum feature space point is s_iv(r_i,r_v), had by described above Body fcm algorithm carries out fcm cluster to two dimensional gray feature space.

Cloud atlas image is carried out with transition cloud system and the division of not distinguishable cloud system: the earth's surface of assuming in cloud atlas, the varieties of clouds can not be thin It is divided into c class: (x₁,x₂,...,x_c), s_iv(a_i,a_v) in each point degree of membership be u (u₁,u₂,...,u_c), then can be according to following Condition determines the classification in cloud atlas belonging to each pixel: if u_i=max (u) ＞ 0.5, then decision-pointCloud In figure pixel r_i(m,n)、r_v(m, n) belongs to classification i；IfThen this pixel, between two kinds of classifications, is subordinate to Genus degree is higher than bottom threshold and is less than upper threshold, can be classified as transition cloud system；If u_i≤ 0.3, then this pixel feature abnormalities mould Paste, it is difficult to distinguish, can be classified as not distinguishable cloud system.Threshold value bound in above-mentioned partiting step can adjust according to practical situation.

Introduce the concept of transition cloud system and not distinguishable cloud system, solve obscuring in image, uncertain pixel returning Genus problem, again carries out fcm cluster to cloud atlas, obtains Clustering Effect as shown in Figure 2.

It is seen that degree of membership be higher than upper threshold (0.5) pixel formed in two-dimensional feature space 8 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 represents transition cloud system and not distinguishable cloud system respectively；In conjunction with the spectroscopic data in table 2, according to pixel in two dimension It may be determined that a is land, b is water body, as shown in Figure 2 for the distribution of feature space.

If land, the cluster areas of water body are respectively s in gray scale two-dimensional feature space_iv(land), s_iv(water body), cloud atlas Middle pixel (x, y) is respectively p infrared with gray value in visible cloud image_i(x,y)、p_v(x, y), after image segmentation, point (x, y) exists Infrared and in visible cloud image gray value is g respectively_i(x,y)、g_v(x,y).Then combine the image segmentation thought based on threshold value, The criterion of the image segmentation using fcm method can be set up, as follows:

In another kind of embodiment, by diagonal filtering method, gray scale characteristic in satellite cloud picture image is gone Except noise and reduction interference process and obtain gray scale sample characteristics data, specific as follows:

Cloud atlas represents, using gray-scale level, the radiation that object sends, and the gray value that pixel has represents cloud or earth's surface surface Bright temperature value.Bright temperature value more low then object height above sea level is higher, such as altostratus；Bright temperature value more high then object height above sea level is lower, such as land, Ocean etc..During generation, digital processing and projective transformation etc., satellite cloud picture leads to unavoidably some noises of generation, noise Not only reduce picture quality, also make objective fuzzy, grind even more to follow-up image segmentation, feature extraction, Cloud-motion wind inverting etc. Study carefully work and cause inconvenience, it is necessary for removing noise and reduction interference therefore in satellite cloud picture image.

For the feature of Cloud-motion wind inverting main body, the globality especially cloud border of cloud affects relatively on Cloud-motion wind inversion result Greatly, therefore, it is proposed that not only having can overcome the disadvantages that cloud cluster noise but also diagonal (major-minor) the mean filter method on clear cloud border being retained, Algorithm is as follows:

Diagonal filtering method

Input: the pixel grey scale value matrix a=[a in n × n window_ij]_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=(m₁+m₂)/2；

By the purpose of cloud atlas image classification is had two, the first separates cloud and land and water body it is ensured that cloud mark The input of wind algorithm is purely " cloud ".It two is that cloud is divided into " high cloud ", " medium cloud ", " low clouds " three major types, anti-with Cloud-motion wind " upper-level winds ", " hollow wind " and " Low level wind " drilled is corresponding, 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 table 3 cloud atlas sample

Determine earth's surface, water body and cloud system total number be 5, using three-dimensional samples gray feature statistic as cluster initial in The heart, in cluster iterative process, continuous adjusting and optimizing cluster centre and subordinated-degree matrix, finally obtain the cluster in three-dimensional feature space As a result, cloud system clustering algorithm is as follows:

Cloud system clustering algorithm:

Input: clusters number c=5, Three-Dimensional Gray eigenmatrix value v₀=[v_ij]_3×5, data volume n；Algorithm terminates threshold epsilon

Output: subordinated-degree matrix u'=[u_ij]_n×n

Calculating process:

Step1: be assigned to v with 5 class three-dimensional samples characteristic quantity initial values₀, to subordinated-degree matrix u⁽⁰⁾Give as follows Initial value, b=0；

$u_{ij}=\underset{k=1}{\overset{c}{\σ}}{d}_{ij}/\underset{k=1}{\overset{c}{\σ}}{(d_{ij}/d_{kj})}^{3/\mathrm{\λ}},i=1,2,...,c;j=1,2,\mathrm{....}n;$

Step2:b++；Calculate new cluster centre v_i ^(b)(i=1,2 ..., c):

$v_i=\left(\underset{j=1}{\overset{n}{\σ}}{u}_{ij}^{2\mathrm{\λ}}{x}_j\right)/\left(\underset{j=1}{\overset{n}{\σ}}{u}_{ij}^{2\mathrm{\λ}}\right)$

Step3: for j=1～n, calculate u^(b), and make u^(b+1)=u^(b)；

1. calculate i_jWith

i_j=i | 1≤i≤c；d_ij=| | x_j-v_i| |=0 }

$\stackrel{\&overbar;}{i_j}=\{1,2,...,c\}-i_j$

2. for data x_j, update its degree of membership.

IfThen to allMake u_ij=0；

Otherwise:

$u_{ij}=\underset{k=1}{\overset{c}{\σ}}{d}_{ij}/\underset{k=1}{\overset{c}{\σ}}{(d_{ij}/d_{kj})}^{3/\mathrm{\λ}},i=1,2,...,c;j=1,2,\mathrm{....}n;$

Step4: if | | u^(b)-u^(b+1)| | ＜ ε, execution step step5；Otherwise return step3；

Step5:u '=u^(b+1)；Terminate

As shown in figure 4, because the calculating of wind arrow each 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) take the satellite infrared passage cloud atlas that three interval times are dt, be temporally sooner or later sequentially named as cloud atlas 1, cloud Fig. 2, cloud atlas 3；

(3) respectively three pictures are carried out grey level enhancement by gradation of image enhancing algorithm respectively, more accurately to track Target cloud mass；

(4) cloud atlas is divided into (width/ ρ) * (height/ ρ) block to follow the trail of block, in formula, width is described every width cloud atlas The width of pixel, height is the height of described every width cloud atlas pixel, and ρ is wind arrow density；

(5) calculation processing unit specified to all block distribution on cloud, uses certain to n parallel processing element Rule carries out task distribution.Parallel processing element herein includes the process of mpi, each stream of the thread in openmp and gpu Processor, and the complete wind arrow inverting task of the Charge-de-Mission in each parallel processing element one；

(6) start n parallel processing transaction and parallel computation is carried out to the process task distributing；

(7) data summarization.The wind arrow of all parallel processing element invertings is collected together, and by effective wind arrow data Write disk preserves；

(8) computing terminates.

In the present embodiment, as shown in figure 5, using the tall and handsome cuda framework reaching company, thread is stream handle in gpu The ultimate 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 refutation task, kernel function kernel_function<<<m, n>>>(arguments) of cuda of execution, Wherein m is the size of grid, and that is, a grid includes how many block；N is the size of block, has in each block How many thread, each thread distribute the calculating task of a wind arrow.

Using gpu concurrent technique, Cloud motion wind inversion algorithm is accelerated, mn wind arrow distribution of computation tasks is arrived gpu's Execute on thread, the experiment of the present invention for the sake of simplicity, thread lattice application is one-dimensional variable blockspergrid, each Threadsperblock thread, Thread Count on each block of the experiment of the present invention is opened up in thread lattice Threadsperblock is 256, the quantity of block in thread lattice

Blockspergrid=(imax+threadsperblock-1)/threadsperblock

Each wind arrow inverting task obtains wants the index of wind arrow of inverting to be

Windindex=blockdim.x*blockidx.x+threadidx.x

Finally call

Dev_makecloud<<<blockspergrid,threadsperblock>>>(arguments ...) is to each wind arrow Inverting task is allocated and executes, and the thread in each block is performed serially in a stream handle, and will be finally inversed by To in the corresponding position windindex of overall wind arrow array opening up.

Accuracy comparison is tested

(1) when inverted parameters are 22 × 22 pixels for target cloud block size, trace regions size is 44 × 44 pixels, wind arrow Density is 2/when spending, have chosen on equator continuous 10 wind arrows as comparison other, the wind of parallel and serial algorithm inverting Arrow accuracy comparison is as shown in table 4.(wherein s represents serial, and t-64 represents 64 thread of distribution, t-128 generation in a thread block In one thread block block of table distribute 128 thread thread, then the block number in thread lattice be general assignment quantity wn divided by Thread quantity in each block).

The serial and concurrent Cloud motion wind inversion algorithm accuracy comparison 1 based on gpu for the table 4

(2) inverted parameters are 32 × 32 pixels for target cloud block size, and trace regions size is 64 × 64, and wind arrow density is 4/degree, have chosen continuous 10 wind arrows of 32 degree of north latitude as comparison other, the wind arrow precision of parallel and serial algorithm inverting (wherein s represents serial, and t-64 represents 64 thread of distribution in a thread block, and t-128 represents one as shown in table 5 for contrast 128 thread thread are distributed, then the block number in thread lattice is general assignment quantity wn divided by each in thread block block Thread quantity in block).

The serial and concurrent Cloud motion wind inversion algorithm accuracy comparison 2 based on gpu for the table 5

Parallel Cloud motion wind inverting, serial inverting and the parallel refutation based on gpu is can be seen that from figure 6 above and Fig. 7 Obtain most of identical, numerical result that small part is similar, because the floating-point operation of gpu also has trueness error, therefore take Obtain this result ideal, thus may certify that the parallel method based on gpu is numerically credible and effective.

The Cloud motion wind product figure being finally inversed by based on the parallel refutation algorithm of gpu is as depicted in figure 8.

Performance comparison is tested

Experiment is 4 core 8 thread cpu using host machine, and the gpu of use is the tall and handsome nvidiageforce reaching 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, algorithm relative to speed-up ratio and parallel efficiency with inverted parameters and thread number of blocks change such as table 6 and table 7 Shown.

The parallel Cloud motion wind inversion algorithm speed-up ratio based on gpu for the table 6 and parallel efficiency 1

The parallel Cloud motion wind inversion algorithm speed-up ratio based on gpu for the table 7 and parallel efficiency 2

By being analyzed to Fig. 9 and Figure 10 as can be seen that in the case of inverted parameters identical, substantially following thread block Number is more, the bigger rule of speed-up ratio lifting.When thread block, < when 64, parallel speedup ratio becomes 2 times of increase trend substantially, when being more than When 64, parallel speedup ratio is relevant with the division of thread block, 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 block) is comprising the mesh of a Slip in mark field of search s (m × m pixel image block), all tracking cloud masses in calculating s and the correlation coefficient of spike cloud mass a, Big correlation coefficient corresponding tracking cloud mass b is target cloud mass.

The computing formula of maximum correlation coefficient r (δ x, δ y) is as shown in (2-1):

$r(\mathrm{\&delta;}x,\mathrm{\&delta;}y)=\frac{{\mathrm{\&sigma;}}_i^n{\mathrm{\&sigma;}}_j^n\&lsqb;a(i,j)-\stackrel{\&overbar;}{a}\&rsqb;\&lsqb;b(\mathrm{\&delta;}x+i,\mathrm{\&delta;}y+j)-\stackrel{\&overbar;}{b}\&rsqb;}{\sqrt{{\mathrm{\&sigma;}}_i^n{\mathrm{\&sigma;}}_j^n{\&lsqb;a(i,j)-\stackrel{\&overbar;}{a}\&rsqb;}^2}\sqrt{{\mathrm{\&sigma;}}_i^n{\mathrm{\&sigma;}}_j^n{\&lsqb;b(\mathrm{\&delta;}y+i,\mathrm{\&delta;}x+j)-\stackrel{\&overbar;}{b}\&rsqb;}^2}}---(2-1)$

Wherein,It is respectively and follow the trail of block deviation 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 spike cloud mass and follow the trail of the picture of block Plain gray value,The average gray being respectively spike cloud mass and following the trail of block.

After obtaining maximum correlation coefficient, the tracking cloud mass b in target search area is with respect to spike cloud mass a on x, y direction The distance of movement is respectively δ x and δ y, if the longitude and latitude of the central point of spike cloud mass a is (x₁,y₁), follow the trail of the center of cloud mass b The longitude and latitude of point is (x₂,y₂), set up earth model (as shown in figure 12) and target cloud mass b can be calculated with respect to spike cloud mass a Apart from d, just can get the speed of wind arrow with d divided by the observation interval t of two width images；Derivation is as follows:

The speed of wind arrow

$v=\frac{d}{t}---(2-2)$

The distance (re be earth radius) of wherein a, b point-to-point transmission

D=re × γ (2-3)

The centre of sphere angle of 2 points of a, b and centre of sphere formation can be calculated with space vector scalar product

$cos\left(\mathrm{\&gamma;}\right)=\frac{\stackrel{\&rightarrow;}{oa}\&centerdot;\stackrel{\&rightarrow;}{ob}}{\left|\stackrel{\&rightarrow;}{oa}\right|\&centerdot;\left|\stackrel{\&rightarrow;}{ob}\right|}---(2-4)$

If vectorialFor (a₁,b₁,c₁), vectorFor (a₂,b₂,c₂), then

a₁=re × cos (x₁)×cos(y₁) (2-5)

b₁=re × sin (x₁)×cos(y₁) (2-6)

c₁=re × sin (y₁) (2-7)

In the same manner

a₂=re × cos (x₂)×cos(y₂) (2-8)

b₂=re × sin (x₂)×cos(y₂) (2-9)

c₂=re × sin (y₂) (2-10)

Formula (2-5)～(2-10) is entered the size that (2-4) can obtain centre of sphere angle γ, then the γ obtaining is substituted into (2-3) Can obtain a, b point-to-point transmission apart from d, finally d is substituted into speed υ that (2-2) just can obtain 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 point a, and a, b are respectively the beginning and end of wind arrow.

Triangle ode is had with the triangle cosine law

$cos\left(a\right)=\frac{{\mathrm{oe}}^2+{\mathrm{od}}^2+{\mathrm{de}}^2}{2oe\&centerdot;od}---(2-10)$

Because triangle oae and triangle oad is right angled triangle, therefore have

oe²=oa²+ae²(2-11)

od²=oa²+ad²(2-12)

Have by (2-12) (2-13) substitution (2-11) and with the triangle cosine law

Cos (α)=cos (β) × cos (γ)+sin (β) × sin (γ) × cos (θ) (2-13)

Wherein

α=x₂-x₁(2-14)

β=y₂-y₁(2-15)

The γ obtaining and above (2-15) (2-16) substitutes into (2-14), can obtain the size of θ using arcsin function. If a point northern area under the line, wind arrow deflection is γ, if x₂＜ x₁, then wind arrow deflection is 2 π-θ；If a point under the line with South area, then wind arrow deflection is π-θ, if x₂＜ x₁, then wind arrow deflection is π+θ.

The height of wind arrow can be specified by the bright temperature value of infrared the one of wind arrow position and vapor channel.

Wind arrow calculates to be needed after finishing to carry out wind speed and direction verification, to reject wind speed and direction in primary Calculation result The larger wind arrow of deviation.Defining the cloud atlas observed three continuous times is c₁、c₂、c₃.If utilizing c₁And c₂The wind arrow cv calculating₁'s Wind speed and direction is respectively υ₁And θ₁, using c₂And c₃Calculate wind arrow cv₂Wind speed and direction be respectively υ₂And θ₂, then two wind arrows Wind direction difference δ θ=| θ₁-θ₂|, the relative mistake δ υ of wind speed=| υ₁-υ₂|, when the wind direction difference δ θ or wind speed difference δ υ of wind arrow are more than During given threshold value, then reject this wind arrow.

In the present embodiment, the scene that maximum correlation coefficient uses is in the wind arrow place block of pixels found out and want inverting The maximum 32*32 region of heart 32*32 pixel region similarity, similarity is monotonically increasing with the relation of correlation coefficient, so The position of target block can be found out using this algorithm, thus calculating the correlation values of wind arrow, schematic diagram is as shown in figure 11.

In the algorithm, the correlation coefficient time calculating between two blocks is less, and the calculating of correlation coefficient array takes Longer, if the size of block is n*n, the trace regions that block is located are m*m, and such as each correlation coefficient calculates, then need Calculate the secondary correlation coefficient of (m-n+1) * (m-n+1), the time is calculate a correlation coefficient time more than 1,000 times, therefore, this Bright this amount of calculation will be optimized.

In the present embodiment, use " variable step iterative calculation ", variable step-size search method refers in trace regions s, spike The moving step length of cloud mass is changed to the strategy of variable step movement by the original strategy being always 1, and the initial step length of search is using relatively Big be worth, progressively reduce afterwards step-length reduce simultaneously hunting zone until step-length be 1, in the case of every kind of step-length all finds out current step Optimum target coupling cloud mass, step-length is just to have found approximately overall optimal coupling cloud mass b when 1.

In this embodiment it is assumed that initial step length is 8, then there is a following derivation:

First pass moving step length is 8, calculates the individual correlation coefficient value of (m/8) * (m/8), finds out in these correlation coefficient value Maximum value, thus finding similarity highest n*n block a1, method is as shown in figure 14.

After finding a1, for avoiding error, the a1 block length of side is increased step pixel, then step-length is reduced to 8/2 I.e. 4 pixels, that is, step is 4, to the n*n region that correlation coefficient is maximum in the region of the n+step length of side centered on a1 Scan for, calculate the correlation coefficient of each n*n block and b block respectively, the coefficient matrix size calculating is (n+step)/4* (n+step)/4, i.e. (n+step)²/ 16, this just represents these step numbers of iterative search and just can find the maximum that step is in the case of 4 Correlation coefficient 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 searching when step step-length is 1 Correlation coefficient block block is and center block b similarity highest block, thus the positional information calculation according to two blocks Go out speed and the angle of wind arrow.

The former method that target cloud is followed the trail of finds Best similarity block needs our calculating of number of times of iterative search to understand,

Stepcount1=(m-n+1)²

Finding Best similarity block after algorithm optimization needs the number of times of iterative search

Stepcount2=(m/step)²+step/2+step/4+…+1

Step-length be revised as m be 64, n be in the case that 32, step is 8, former method needs 1089 steps search, after optimization Method needs 71 step search, greatly simplifies the time that search needs, can accelerate more than 15 times in theory, and the precision in later stage Optimization can do the adjustment related to m size for the scope searching again for after variable step.

Although embodiment of the present invention is disclosed as above, it is not restricted to listed in description and embodiment With, it can be applied to various suitable the field of the invention completely, for those skilled in the art, can be easily Realize other modification, therefore under the general concept being limited without departing substantially from claim and equivalency range, the present invention does not limit In specific details with shown here as the legend with description.

Claims (9)

1. a kind of parallel Cloud-motion wind inversion method based on fuzzy clustering pretreatment cloud system is it is characterised in that include:

Cloud atlas image is divided into polylith to follow the trail of block, and is assigned to, by following the trail of block, the multiple gpu parallel computation arrays specified, By maximum correlation coefficient, the described target cloud mass followed the trail of in block is carried out mating laggard in each independent computing array Row Inversion Calculation, obtains the described wind speed following the trail of block wind arrow and wind direction；

By described maximum correlation coefficient, the matching process of target cloud mass is being included: in described tracking block, choose spike Cloud mass, in the region of search comprising described spike cloud mass, calculates the whole and described spike cloud mass picture that step-size in search is more than 1 The correlation coefficient of the equivalent block of element, determines maximum correlation coefficient in these correlation coefficienies, as similar to described spike cloud mass Degree highest block, is increased centered on this block in the region of pixel, step-length is progressively halved and searches again for out similarity Highest block, until searching out the block that step-length is 1 maximum correlation coefficient searching, as target cloud mass；And

Wherein, determine, after fuzzy clustering method separates with original cloud atlas image being entered rack, the cloud carrying out Cloud-motion wind inverting Figure image, comprising:

The gray feature data of collection cloud atlas image mid-infrared passage, is processed to described gray feature data, obtains sample Gray feature data, described sample gray feature data is classified, and makes each class sample gray feature data right respectively Answer a kind of cloud atlas type, and then determine described cloud atlas image.

2. the parallel Cloud-motion wind inversion method based on fuzzy clustering pretreatment cloud system as claimed in claim 1 it is characterised in that Also include, choose three width cloud atlas of adjacent constant duration, draw wind speed and the wind direction of every width cloud atlas wind arrow respectively, by described wind The wind speed of arrow, the wind direction and wind arrow speed difference of cloud atlas adjacent in described cloud atlas, direction difference and threshold value being compared according to sequential Relatively, execute following principle:

If wind arrow speed is less than the minimum speed threshold value specified, reject this wind arrow；Or

If the difference of two wind arrow speed is more than the maximum speed discrepancy threshold value specified, reject this wind arrow；Or

If the difference in two wind arrow directions is more than the maximum direction difference limen value specified, reject this wind arrow；Or

If not having direction difference therewith less than the wind arrow of the maximum direction difference limen value specified, to reject this in the peripheral region of two wind arrows Wind arrow.

3. the parallel Cloud-motion wind inversion method based on fuzzy clustering pretreatment cloud system as claimed in claim 1 or 2, its feature exists In described cloud atlas view data being carried out with multi-channel data extraction, retains infrared channel data.

4. the parallel Cloud-motion wind inversion method based on fuzzy clustering pretreatment cloud system as claimed in claim 3 it is characterised in that Also include: pretreatment is carried out to described cloud atlas image, it is included to the longitude and latitude projective transformation such as image is carried out with by original disk Attempt to change and be changed to histogram.

5. the parallel Cloud-motion wind inversion method based on fuzzy clustering pretreatment cloud system as claimed in claim 4 it is characterised in that Described cloud atlas image is carried out with pretreatment also include carrying out grey level enhancement to described image.

6. the parallel Cloud-motion wind inversion method based on fuzzy clustering pretreatment cloud system as claimed in claim 1 it is characterised in that Described gray feature data is carried out process and includes:

The gray feature data input of cloud atlas image mid-infrared passage is obtained pixel grey scale value matrix, passes through Calculate leading diagonal pixel grey scale average in a matrix, pass throughCalculate auxiliary diagonal pixel in a matrix Gray average, and then pass through m=(m₁+m₂Pixel grey scale average is obtained in)/2, i.e. sample gray feature data.

7. the parallel Cloud-motion wind inversion method based on fuzzy clustering pretreatment cloud system as claimed in claim 6 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 is chosen in described infrared channel, Described sample gray feature data is set up 5 gray feature data statisticss amounts, image segmentation is carried out to Nephogram and obtains institute State cloud atlas image.

8. the parallel Cloud-motion wind inverting side based on fuzzy clustering pretreatment cloud system as any one of claim 1,2,4-7 Method determines described tracking block number it is characterised in that calculating according to the wide high of every width cloud atlas pixel and wind arrow density, and it is (width/ ρ) × (height/ ρ), in formula, width is the width of described every width cloud atlas pixel, and height is described every width cloud atlas The height of pixel, ρ is wind arrow density.

9. the parallel Cloud-motion wind inversion method based on fuzzy clustering pretreatment cloud system as claimed in claim 8 it is characterised in that By described maximum correlation coefficient, the matching process of target cloud mass is being included: in described maximum correlation coefficient, in institute State in tracking block and choose spike cloud mass a, it is n × n-pixel image block, in the m × m pixel map comprising described spike cloud mass a It is equal to the correlation coefficient of block as in block search region, calculating whole and described spike cloud mass a pixels that step-size in search is m, Maximum correlation coefficient is determined, as with described spike cloud mass a similarity highest n × n-pixel image block in these correlation coefficienies A1, increases in described block a1Then step-length is kept to by individual pixelIndividual pixel, centered on block a1Calculate the correlation coefficient of each n × n-pixel image block in the region of the length of side, these correlation coefficienies determine maximum phase Close coefficient and then determine next and described spike cloud mass a similarity highest n × n-pixel image block a2, step-length is progressively subtracted Half until searching out the block that step-length is 1 maximum correlation coefficient searching, as target cloud mass, by described spike cloud mass with And the positional information calculation of described target cloud mass obtains wind speed and the wind direction of wind arrow.