CN104318550A - Eight-channel multi-spectral imaging data processing method - Google Patents

Abstract

The invention provides an eight-channel multi-spectral imaging data processing method. The eight-channel multi-spectral imaging data processing method comprises multi-spectral image preprocessing; multi-spectral image pre-rectification; wavelength calibration; radiation calibration; spectral curve building; multi-spectral image color calculation; multi-spectral image feature extraction; spectrum and color matching. According to the eight-channel multi-spectral imaging data processing method, the complete integrated data processing steps the detailed specific guidance are provided for the technical personnel, the corresponding method is selected according to the characteristics of eight-channel multi-spectral imaging data, and accordingly the good effect is achieved.

Description

Eight passage multispectral imaging data processing methods

Technical field

The present invention relates to a kind of multispectral imaging data processing method, particularly relate to a kind of eight passage multispectral imaging data processing methods, belong to multispectral imaging field.

Background technology

Multispectral imaging (Multispectral Imaging) is the abbreviation of multiband light spectral imaging technology, it is acquisition of information and the treatment technology of the unification of a kind of collection of illustrative plates, namely, while obtaining the space dimension information of target, the spectrum dimension information of this target is also obtained.According to the spectral band number of institute's obtaining information or the difference of spectral resolution size, current spectral imaging technology roughly can be divided into multispectral imaging (Multi-spectral Imaging), high light spectrum image-forming (Hyper-spectral Imaging) and Hyper spectral Imaging (Ultra-spectral Imaging).

Multispectral image data cube has the feature of multispectral section and high spatial resolution, therefore contains abundanter ground object target information than full-colour image or simple spectral measurement.According to spectral image data, people can analyze ground object target from spatial match and Spectral matching two aspects and identify.Therefore, spectral imaging technology is at palegeology, vegetation investigation, atmospheric exploration, ocean remote sensing, agricultural science and technology, and the aspects such as environmental monitoring, disaster reduction and prevention and military science are with a wide range of applications.Especially, in military surveillance, spectrum imaging system can judge the attribute of target according to the characteristic spectrum of the radiation of various weapon system-of-systems or reflection, find traditional not detectable military target of Visible Light Reconnaissance system; By the analysis to characteristic spectrum, kind and the model of weaponry can also be judged and biochemical poison gas identified and forecast.

From 20 century 70s, first multi-optical spectrum imaging technology is applied in airborne and spaceborne RS earth observation field, and is progressively played a role in fields such as agricultural, biomedicine, museum work collection, beauty treatment, high precision colour print, computer graphicss.Since the nineties in 20th century, the research of the multi-optical spectrum imaging technology of visible light wave range is paid attention to energetically, now become the study frontier of related science.International Commission on Illumination (CIE) lists the coloured image technical committee i.e. developing direction of the 8th branch (TC8-7) in November, 2002 in multi-optical spectrum imaging technology.

Domestic have certain basis in multispectral imaging and data processing, based on image procossing, multispectral imaging data are mainly concerned with the steps such as Image semantic classification, denoising, correction, registration, but often not there is when analyzing the image after above-mentioned treatment step comparatively science and comprehensive method, therefore needing to carry out systematicness, comprehensive and comprehensive analytical approach to whole multispectral imaging data processing.

Summary of the invention

In view of above-mentioned the deficiencies in the prior art part, the object of the present invention is to provide a kind of eight passage multispectral imaging data processing methods, the invention provides comprehensively complete multispectral imaging data processing step, for people's

Compared to prior art, eight passage multispectral imaging data processing methods provided by the invention comprise step one, multispectral image pre-service; The pre-registration of step 2, multispectral image; Step 3, wavelength scaling; Step 4, radiation calibration; Step 5, the curve of spectrum build; Step 6, multispectral image color calculate; Step 7, multispectral image feature extraction; Step 8, spectrum and color-match.

The present invention also provides a kind of eight passage multispectral imaging data processing methods, comprises the following steps:

Step one, multispectral image pre-service, comprise spectrum picture noise remove, CCD homogeneity correction and CCD gamma correction etc.; Described spectrum picture noise remove uses to face the territory method of average and median filtering method removal salt-pepper noise; Described CCD homogeneity correction is by carrying out imaging to the planar target with uniform luminance distribution, then calculate according to the response of each pixel of CCD and extract the homogeneity correction coefficient of each pixel of CCD surface, and these coefficients being saved in the calibration file of software systems; Described CCD gamma correction, by measuring the target with different brightness degree, then according to the response data of each brightness degree, adopts curve fit or other mathematical method to obtain non-linear correction factor;

The pre-registration of step 2, multispectral image, comprises Edge extraction and carries out image conversion according to the unique point after edge extracting; Described Edge extraction is carried out by the following method:

For two dimensional image signal, first come smoothing with following Gauss function:

$G(x,y,\mathrm{\σ})=\frac{1}{2\mathrm{\π}{\mathrm{\σ}}^2}\exp (-\frac{1}{{2\mathrm{\σ}}^2}(x^2+y^2))$

G (x, y, σ) is a function with circular symmetry, its level and smooth effect controls by σ, owing to carrying out linear smoothing to image, makes g (x, y) be level and smooth after image, obtain: g (x, y)=G (x, y, σ) * f (x, y), wherein f (x, y) is level and smooth front image;

Because marginal point is the place that in image, gray-value variation is violent, the sudden change of this image intensity will produce a peak in first order derivative, or be equivalent to a generation zero cross point in second derivative, and be nonlinear along the second derivative of gradient direction, so substitute with Laplacian operator, namely use:

${\▿}^{2}g(x,y)={\▿}^2(G(x,y)*f(x,y))=\left({\▿}^{2}G(x,y,\mathrm{\σ})\right)*f(x,y)$

Zero cross point as marginal point, in formula for LOG wave filter,

${\▿}^{2}G(x,y,\mathrm{\σ})=\frac{{\∂}^{2}G}{\∂x^2}+\frac{{\∂}^{2}G}{\∂y^2}=\frac{1}{\mathrm{\π}{\mathrm{\σ}}^4}(\frac{x^2+y^2}{2{\mathrm{\σ}}^2}-1)\exp (-\frac{1}{2{\mathrm{\σ}}^2}(x^2+y^2))$

Above formula is exactly LOG edge detection operator, LOG operator is being similar to ganglia retinae receptive field spatial organization, can regard as and be made up of focus of excitation district and an inhibition surrounding zone, occur with the template form of 62 σ × 62 σ sizes, when σ gets different values, then the image border under available operators detection different scale, usually, we get σ >=1, carry out rim detection namely extract image border with above-mentioned LOG edge detection operator to image;

Described image conversion uses rigid transformation, and transformation for mula is as follows:

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\end{array}\right]=k\left[\begin{array}{cc}\mathrm{cod\θ}& \sin \mathrm{\θ}\\ -\sin \mathrm{\θ}& \cos \mathrm{\θ}\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]+\left[\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\end{array}\right]$

In formula: (x, y) is the point of piece image and benchmark image, after conversion in corresponding second width image and k, θ and Δ x and Δ y are the scale factor of the first width figure and the second width figure, twiddle factor and coordinate translation amount respectively, these parameters are obtained by manual operation;

Step 3, wavelength scaling, select cubic spline interpolation method, the SPL of respective numbers point known wavelength reflectivity data is set up according to the quantity of optical filter, i.e. spectral reflectivity function curve, according to obtained function curve, according to even step sizes, wavelength is split between every two known wavelength, just obtain the spectral reflectance data of other wavelength;

Step 4, radiation calibration, by various calibrated radiation source, quantitative relationship between different wave spectrum section sets up digital quantization value that the spectral radiance value at imaging spectrometer entrance pupil place and imaging spectrometer export, select the target light source of known spectra radiation characteristic as the measuring object of multi-spectral imager, then simulate duty and the environment of multi-spectral imager;

Step 5, the curve of spectrum build, and adopt cubic spline functions to set up spectral reflectance rate curve;

Step 6, multispectral image color calculate, and comprise computed image color tristimulus values, calculating coloured image RGB, computed image aberration and calculate near infrared luminance factor β value;

Step 7, multispectral image feature extraction, comprise feature selecting and feature extraction;

Step 8, spectrum and color-match, comprise the color-match based on tristimulus values, the color-match based on spectrum and Spectral matching.

Eight passage multispectral imaging data processing methods provided by the invention provide a kind of complete comprehensive data processing step, to the guidance that technician is concrete very in detail, and have selected corresponding method for the feature of eight passage multispectral imaging data, reach good effect.

Accompanying drawing explanation

Fig. 1 is that Laplacian operator commonly uses template schematic diagram;

Fig. 2 is LOG operator 5 × 5 rank template schematic diagram;

Fig. 3 is radiometric calibration device schematic diagram.

Embodiment

The invention provides a kind of eight passage multispectral imaging data processing methods, for making object of the present invention, technical scheme and effect clearly, clearly, developing simultaneously referring to accompanying drawing, the present invention is described in more detail for embodiment.Should be appreciated that specific embodiment described herein only in order to explain the present invention, be not intended to limit the present invention.

Embodiment

Eight passage multispectral imaging data processing methods disclosed in the present embodiment, its data acquisition system (DAS) comprises eight passages and is made up of the spike interference filter of the Visible-light CCD object lens of fixed focal length, 8 different wave lengths, ccd sensor, image pick-up card, industrial computer and data processing software system.Wherein, the spike interference filter covering visible light of 8 different wave lengths and near-infrared band; The optical lens being separately mounted on 8 passages is anterior; Centre wavelength and bandwidth (FWHM) parameter as shown in the table.

In order to improve actinometry precision and the accuracy of multi-optical spectrum imaging system, general needs carries out pre-service to the multispectral image data cube collected.Pre-service mainly comprises the contents such as spectrum picture noise remove, CCD homogeneity correction, CCD gamma correction.

Due to black, the white point shape random noise that exist in the imaging spectrometer product at different levels image that each side factors such as focal plane array sensor performance, target brightness level, circuit noise cause, be usually called as impulsive noise (ImPulse Noise) or salt-pepper noise (Salt and Pepper Noise).Salt-pepper noise significantly reduces picture quality, is therefore necessary to carry out detecting and processing.The feature of the images with salt and pepper noise is that noise spot is uniformly distributed in entire image, and the number percent that the noise number be distributed in image accounts for sum of all pixels is called noise rate.The detection method of salt-pepper noise has multiple.A kind of conventional method determines that whether the signal level of pixel to be measured is abnormal according to the level difference signal of pixel to be measured and surrounding pixel, if abnormal, belongs to salt-pepper noise point.

Traditional salt-pepper noise denoising method comprises faces the territory method of average and median filtering method.

(1) the territory method of average is faced

Facing the territory method of average is exactly replace noise value a (i, j) with the average a (i, j) ' of surrounding pixel DN value, such as, 8 point values before and after noise can be adopted to be averaged:

$a{(i,j)}^{\′}=\frac{1}{8}(\underset{j-4}{\overset{j-1}{\mathrm{\Σ}}}a(i,j)+\underset{j+1}{\overset{j+4}{\mathrm{\Σ}}}a(i,j))$

(2) intermediate value method of substitution

Intermediate value method of substitution adopts some sequences before and after noise, the intermediate value after sequence substituted noise spot, such as, chooses 9 points comprising noise spot:

a(i，j)′＝mid(a(i，x)，x∈(j-4，j+4)

For imaging-type radioactivity detection, identificationm, and computation, under ideal conditions, require that the whole image planes of imageing sensor have uniform radiometric response characteristic; Namely, when measured target has uniform luminance distribution, the response of the whole image planes of imageing sensor also should be uniform, otherwise measurement result accurately can not reflect the brightness relationship in tested scene between various target.But under many circumstances, preposition optical system produces gradual halation phenomena in focal plane and image sensor surface; In addition, likely there is the non-equal response characteristic on surface in CCD image sensor itself; Therefore need in these cases to carry out homogeneity correction to CCD image sensor.

In a lot of document, again the CCD homogeneity correction in spectrum imaging system is classified as " relative detector calibration " or " relative radiometric calibration ", their object is consistent.Also should be noted that another key concept in spectrum imaging system calibration work---" absolute radiation correction " or " absolute radiometric calibration ", this is one and " relative detector calibration " or " relative radiometric calibration " diverse concept." absolute radiometric calibration " refers to that the relation between the response of CCD and the radiancy value of measured target is calibrated.Generally, " absolute radiation correction " needs carry out after completing " relative detector calibration ".

One of basic skills of CCD homogeneity correction carries out imaging to the planar target with uniform luminance distribution, then calculate according to the response of each pixel of CCD and extract the homogeneity correction coefficient of each pixel of CCD surface, and these coefficients to be saved in the calibration file of software systems.CCD homogeneity correction can be combined with CCD gamma correction under many circumstances and jointly complete.

So-called CCD gamma correction refers to the correction carried out the non-linear response characteristic of multi-spectral imager CCD.Common CCD sensor has linear photoelectric response characteristic within the scope of certain intensity of illumination, but when the light signal strength that CCD receives is very little or very high, the response characteristic of CCD all can departs from linear.Therefore need under many circumstances to carry out gamma correction to CCD.

One of basic skills that CCD non-linearity corrects measures the target with different brightness degree, then according to the response data of each brightness degree, adopts curve fit or other mathematical method to obtain non-linear correction factor.

(1) under 10 radiancy grades, calibration experiment is done to multi-spectral imager CCD, with the even light of integrating sphere injection for light source, obtain the DN data of 10 768 × 576 pixels of CCD response.

(2) repetitive measurement CCD chip is under non-illuminated conditions, the output of each pixel, obtains the dark current data DN value of several 768 × 576 pixels.

Bearing calibration:

CCD homogeneity correction experimental data is utilized to calculate correction coefficient to the row, column nonlinear response of CCD chip,

Namely each pixel obtains one group of multinomial coefficient.Make a concrete analysis of trimming process below.

(1) each pixel dark current value in the CCD chip of repetitive measurement is averaged obtains Ai, j, then to full wafer chip is averaging and obtains A.

(2) the DN output valve of the i-th row jth row pixel of CCD chip under the n-th radiancy grade is denoted as (DN _n) _{i, j}, record the DN value of all pixels in the CCD chip under 10 radiancy grades, and calculate the mean value of all pixel data output DN of CCD chip under the n-th radiation levels, be denoted as DN _n.

(3) pixel dark current value A is utilized _{i, j}, CCD chip dark current mean value A, the response (DN of each pixel of CCD chip under each radiation levels _n) _{j, j}with mean value DN _n, calculate homogeneity correction coefficient a _n, i, j, computing formula is system of equations (1).

$\left\{\begin{array}{c}\stackrel{\‾}{{\mathrm{DN}}_1}-\stackrel{\‾}{A}=a_{1,i,j}{\left({\mathrm{DN}}_1\right)}_{i,j}-A_{i,j}\\ \stackrel{\‾}{{\mathrm{DN}}_2}-\stackrel{\‾}{A}=a_{2,i,j}{\left({\mathrm{DN}}_2\right)}_{i,j}-A_{i,j}\\ \·\·\·\\ \stackrel{\‾}{{\mathrm{DN}}_n}-\stackrel{\‾}{A}=a_{n,i,j}{\left({\mathrm{DN}}_n\right)}_{i,j}-A_{i,j}\end{array}\right.---\left(1\right)$

(4) solving equation group (1) can draw the multinomial coefficient ak corresponding to each pixel, k=1,2n.Complete to pixel corrections all on CCD can be the three-dimensional array form of a [k] [i] [j] by storage of array.For pixel same on CCD, different brightness degree can comprise a set of coefficient, carries out curve fitting to this set of coefficient, preserves fitting coefficient.

The registration of multispectral image is one of important step of whole spectrum imaging system data processing, and it is related to the accurate fusion problem of spectrum peacekeeping image dimension.In the data cube of multi-spectral imager, due to the motion of the target travel that exists in the rigging error of optical system and sampling process or instrument itself, cause the difference that there is the aspect such as geometric position and magnification between the multichannel image of Same Scene, therefore employing image registration link is necessary at the Data processing of multi-spectral imager, to make the image of data cube at space dimension accurate registration.

Cross-correlation method is the method based on gradation of image information, and it is suitable for the registration between the image of same sensor acquisition.The calculated amount that the method requires is large, is not suitable for process non-linear deformation and local deformation problems, need to improve.

Fourier transform is after carrying out fast fourier transform to image, and application phase such as to be correlated with at the image registration of technical finesse rotation, Pan and Zoom mismatch.But Fourier transform method can not process the registration of the problems such as nonlinear deformation and different gray scale attributed graph picture.

It is the method for registering the most often adopted when not knowing the mapping mode of two width images that point maps.But the positional precision of unique point easily by the impact of the subjective judgement of people, can not obtain precise and stable registration result, therefore, some mapping method also often utilizes the feedback between each stage to find optimal transformation.

Elastic model method is mainly used in the registration between medical image at present.Method for registering based on wavelet transformation causes the great attention of people in recent years, is characterized in calculated amount when can reduce image registration greatly.

Consider that 8 passage imaging spectrometers of the present invention have the factors such as number of active lanes is less, the processing time is unrestricted, therefore image registration of the present invention adopts a kind of manual intervention extraction image characteristic point to carry out the method for mating, and the method can classify as point mapping.

Point mapping needs to extract the edge details of each width image so that extraction to unique point or characteristic matching vector before carrying out registration.Method for detecting image edge at present existing ripe method can utilize, and can be divided into two large classes substantially at present, i.e. first order differential operator and Second Order Differential Operator by the edge detection operator that people adopt.First order differential operator mainly contains gradient operator, Roberts cross operator, Prewitt operator and Sobel operator; Second Order Differential Operator mainly refers to Laplacian (Laplce) operator.These operators are all in response to gray level change, or average gray level change.

If added Gauss (Gauss) conversion before second-order differential Laplacian operator, then form famous LOG operator and Gauss-Laplace operator.

Above-mentioned edge detection operator respectively has its feature.Adopt which kind of operator to carry out edge extracting in the present invention and then depend on analysis to practical application effect.Multiple arithmetic operators is applied in 8 channel image sequences of the present invention by we for this reason.

Although First-order Gradient operator and Roberts operator are not clear especially to Edge extraction but substantially remain the edge feature of original image, not fine to the little frame such as window extraction effect of image; Contrary Sobel operator and Prewitt operator not only can extract image border clearly, and the frame less for the relative overall profile in image also can be complete extracts, but these two kinds of operators easily produce dual edge when extracting image; And Laplacian operator is a kind of second derivative operator, so quite responsive to the noise in image.In addition for Laplacian operator, after second differential, the edge pixel extracted not on original border, namely because second differential operator does not have directivity, to location distortion also cannot revise, so actual edge extract in not with consideration.

It is the most clear that LOG operator extraction obtains contour edge, loses fringing and occur minimum, so the present invention finally selects LOG operator to complete the Edge extraction work of image registration.

The edge detection operator that the present invention adopts is LOG operator, and it on the basis of Laplacian operator, adds Gauss conversion and realizes, and is therefore necessary first to introduce Laplacian operator.

Laplacian operator is a kind of second derivative operator, and to 1 continuous function f (x, y), it is defined as follows in the Laplacian value of position (x, y):

${\mathrm{\Δ}}^{2}f=\frac{{\∂}^{2}f}{{\∂x}^2}+\frac{{\∂}^{2}f}{{\∂y}^2}$

$\frac{{\∂}^{2}f}{{\∂x}^2}=\frac{\∂\left(G_x\right)}{\∂x}=\frac{\∂\left(f\right[i,j+1]-f[i,j\left]\right)}{\∂x}=f[i,j+2]-2f[i,j+1]+f[i,j]$

$\frac{{\∂}^{2}f}{{\∂y}^2}=\frac{\∂\left(G_y\right)}{\∂y}=\frac{\∂\left(f\right[i+1,j]-f[i,j\left]\right)}{\∂y}=f[i+1,j+1]-2f[i,j+1]+f[i-1,j+1]$

Above approximate expression is centered by [i, j+1], replaces with j-1:

$\frac{{\∂}^{2}f}{\∂y^2}=f[i+1,j]-2f[i,j]+f[i-1,j]$

$\frac{{\∂}^{2}f}{\∂x^2}=f[i,j+1]-2f[i,j]+f[i,j-1]$

In the handling procedure of digital picture, the Laplacian value of computing function also can realize by various template.The coefficient being corresponding center pixel to the basic demand of template should be positive, and the coefficient of corresponding center adjacent pixel should be negative, and they and should be zero.Two kinds of conventional Laplacian operator templates as shown in Figure 1.

Laplacian operator is a kind of second derivative operator, so quite responsive to the noise in image.It often produces the wide edge of double image element in addition, and can not provide the information of edge direction.Due to above reason, Laplacian operator is seldom directly used in Edge detected, and determines after being mainly used in known edge pixel that this pixel is in the dark space of image or area pellucida.

For Laplacian operator, after second differential, the boundary pixel extracted is not on original border.Because second differential operator does not have directivity, also cannot revise location distortion.Therefore, locating distortion is Laplacian operator Problems existing.

First image registration of the present invention adopt LOG wave filter to carry out edge extracting.In the various method of Edge extraction, because the larger change of gray scale always corresponds to some larger derivatives, so the edge detection operator proposed the earliest is gradient operator and Laplacian operator.Gradient method computing is simple, software simulating is convenient, but produce wider response in its region near border, the result of gained usually needs refinement in addition, this not only have impact on the positioning precision on border, and the quality on border can be affected, and Laplacian operator is high frequency sensitivity, so larger by the impact of high frequency noise.In order to effectively suppress the impact of high frequency noise, a kind of improving one's methods first is carried out suitable level and smooth, and with restraint speckle, and then carry out asking micro-, namely LOG filter skirt detects.

For two dimensional image signal, first come smoothing with following Gauss function:

$G(x,y,\mathrm{\σ})=\frac{1}{2\mathrm{\π}{\mathrm{\σ}}^2}\exp (-\frac{1}{{2\mathrm{\σ}}^2}(x^2+y^2))$

G (x, y, σ) is a function with circular symmetry, its level and smooth effect controls by σ, owing to carrying out linear smoothing to image, is mathematically carry out convolution, make g (x, y) be level and smooth after image, obtain: g (x, y)=G (x, y, σ) * f (x, y), wherein f (x, y) is level and smooth front image.

Because marginal point is the place that in image, gray-value variation is violent, the sudden change of this image intensity will produce a peak in first order derivative, or be equivalent to a generation zero cross point in second derivative, and be nonlinear along the second derivative of gradient direction, calculate comparatively complicated, so substitute with Laplacian operator, namely use:

${\▿}^{2}g(x,y)={\▿}^2(G(x,y)*f(x,y))=\left({\▿}^{2}G(x,y,\mathrm{\σ})\right)*f(x,y)$

Zero cross point as marginal point, in formula for LOG wave filter.

${\▿}^{2}G(x,y,\mathrm{\σ})=\frac{{\∂}^{2}G}{\∂x^2}+\frac{{\∂}^{2}G}{\∂y^2}=\frac{1}{\mathrm{\π}{\mathrm{\σ}}^4}(\frac{x^2+y^2}{2{\mathrm{\σ}}^2}-1)\exp (1-\frac{1}{2{\mathrm{\σ}}^2}(x^2+y^2))$

Above formula is exactly LOG edge detection operator.LOG operator is Mexico's straw hat shape, being similar to ganglia retinae receptive field spatial organization, can regard as and be made up of focus of excitation district and an inhibition surrounding zone, usually occur with the template form of 62 σ × 62 σ sizes, when σ gets different values, then the image border under available operators detection different scale.Usually, we get σ >=1.

LOG wave filter has following two outstanding features:

(1) the Gauss function part G in this wave filter can image blur, effectively eliminates all yardsticks and changes much smaller than the image intensity of Gauss distribution space constant σ.Why selecting Gauss function to carry out blurred picture because of it is all level and smooth, localization in spatial domain and frequency domain inside, and the possibility therefore introducing any change do not occurred in original image is minimum.

(2) this wave filter adopts Laplacian operator to reduce calculated amount.If use picture or such first directional derivative, just must find out their peak, valley along each orientation; If use picture or such Second order directional just must detect their zero cross point.But all these operators have a common shortcoming: have directivity, and they are all relevant with orientation.In order to avoid the computation burden caused due to directivity, need to manage to select an operator had nothing to do with orientation.The isotropic differentiating operator of lowest-order is exactly just in time Laplacian operator

In actual program process of the present invention, LOG operator can realize by template.Fig. 2 is 5 × 5 templates of the LOG operator that the present invention adopts.

After the edge extracting completing image, namely can carry out image conversion according to unique point, thus realize registration.Aimed at another piece image by piece image, often need carry out a series of conversion to piece image, these conversion can be divided into rigid body translation, affined transformation, projective transformation and nonlinear transformation etc.

(1) rigid body translation

If the distance between 2 in the first sub-picture still remains unchanged through transforming to after in the second sub-picture, then this conversion is called rigid body translation.Rigid body translation can be analyzed to translation, rotation and reversion (mirror image), and in two-dimensional space, point (x, y) through rigid body translation to the transformation for mula of point (x ', y ') is:

Wherein for rotation angle, (t _x, t _y) be translation vector, k is zoom factor.

(2) affined transformation

It is still straight line that straight line after conversion on the first sub-picture is mapped to the second sub-picture, and keeping parallelism relation, such conversion is called affined transformation.Affined transformation can be decomposed into linearly (matrix) conversion and translation transformation.In 2D space, transformation for mula is:

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\end{array}\right]=\left[\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]+\left[\begin{array}{c}{t}_x\\ t_y\end{array}\right]$

Wherein, $\left[\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right]$ For real matrix

(3) projective transformation

It is still straight line that straight line after conversion on the first sub-picture is mapped on the second sub-picture, but parallel relation does not keep substantially, and such is transformed to projective transformation, and projective transformation can represent with linear (matrix) conversion on higher dimensional space.Transformation for mula is:

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\end{array}\right]=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\end{array}\right]\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]$

(4) nonlinear transformation

Nonlinear transformation can be transformed to curve straight line.In 2D space, can represent with once formula:

(x′，y′)＝F(x，y)

Wherein, F represents any one functional form the first sub-picture is mapped on the second sub-picture.Typical nonlinear transformation is as polynomial transformation, and in 2D space, polynomial function can be write as following form:

x′＝a ₀₀+a ₁₀x+a ₀₁y+a ₂₀x ²+a ₁₁xy+a ₀₂y ²+…

y′＝b ₀₀+b ₁₀x+b ₀₁y+b ₂₀x ²+b ₁₁xy+b ₀₂y ²+…

Nonlinear transformation is than being more suitable for the image registration problem that those have overall situation distortion, and global approximation rigid body but there is the registration situation of deformation local.

According to feature of image that the present invention gathers, mainly there is rotational differential and translational difference between each width image, but enlargement factor difference is very little, therefore mainly selects rigid body translation to process image.The rigid body translation formula adopted in the present invention is as follows:

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\end{array}\right]=k\left[\begin{array}{cc}\mathrm{cod\θ}& \sin \mathrm{\θ}\\ -\sin \mathrm{\θ}& \cos \mathrm{\θ}\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]+\left[\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\end{array}\right]$

In formula: (x, y) is the point of piece image and benchmark image, after conversion in corresponding second width image and k, θ and Δ x and Δ y are the scale factor of the first width figure and the second width figure, twiddle factor and coordinate translation amount respectively; These parameters are all obtained by manual operation in the image registration interface of data processing software of the present invention.

The object of imaging spectrometer spectral calibration and wavelength scaling determines the centre wavelength of each spectral band of imaging spectrometer.Spectral calibration results is one of main performance index of imaging spectrometer, and its precision directly affects the reliability of measurement data.

The wavelength scaling method of imaging spectrometer can with reference to the calibrating method of similar spectrum measurement instruments.

Such as adopt the standard linear light source with known wavelength or known spectra reflectivity Characteristics target to carry out the demarcation of special wavelength point, all the other wavelength location adopt method that is non-linear or linear interpolation to determine.

Although the calibrating method of spectral instrument belongs to mature technology at present, the scheme that can use for reference is a lot, and the spectrum calibration method of 8 passage imaging spectrometers of the present invention has its singularity.Namely the beam splitting system of 8 road multi-spectral imagers of the present invention is made up of the narrow band pass filter of 8 fixed wave length, therefore measure raw data cube 8 the passage compositions obtained, and the wavelength of these 8 passages is known.The response data of all the other spectral bands is obtained by non-linear interpolation on the data basis of 8 passages, and therefore the concrete wavelength value of these interpolation wave bands is also obtained by non-linear interpolation.

The present invention selects cubic spline interpolation method, sets up the SPL of 8 known wavelength reflectivity datas, i.e. spectral reflectivity function curve.According to obtained function curve, according to even step sizes, wavelength is split between every two known wavelength, just obtain the spectral reflectance data of other wavelength.

It can thus be appreciated that wavelength scaling precision of the present invention depends primarily on the nonlinearities change degree of the curve of spectrum that the accuracy of 8 fixed wave length, the interval width of every two wavelength points and interpolation obtain.

Radiometric calibration refers to " absolute radiometric calibration " of multi-spectral imager, also known as work " absolute radiation correction " in some document, its task to set up the quantitative relationship in the digital quantization value (DN) of image device detecting element output of each spectrum channel of multi-spectral imager (the present invention has 8 passages) and visual field corresponding to it between output radiation brightness value.

When multi-spectral imager is through radiation calibration and after obtaining radiant correction coefficient, its measurement data can be just the spectral radiance degrees of data of standard according to radiometric calibration transformation of coefficient, as spectrum width brightness curve, thus the spectral characteristics of radiation of correct reflection measured target.Also the parameters such as the spectral reflectivity of target can be calculated further by the spectral characteristic of known spectra width brightness curve and lighting source.

As previously mentioned, radiation calibration is by various calibrated radiation source, the quantitative relationship between different wave spectrum section sets up digital quantization value that the spectral radiance value at imaging spectrometer entrance pupil place and imaging spectrometer export.Therefore be all generally select the target light source of known spectra radiation characteristic as the measuring object of multi-spectral imager, then simulate duty and the environment of multi-spectral imager.

Consider that luminous energy that multi-spectral imager receives mainly comes from the solar radiation of earth reflection, the spectral distribution of solar simulating radiation of should trying one's best for the light source of radiation calibration.But in Laboratory Calibration equipment, the spectral distribution of solar simulating radiation is very difficult accurately, therefore generally adopt the Halogen lamp LED that colour temperature is higher as luminophor.In addition, for the calibration of the focus planardetector of multi-spectral imager, the outlet of the integrating sphere of larger caliber should be adopted as light source, and irradiate the whole visual field of multi-spectral imager sensor.In addition, also need the spectral radiant emittance meter of use standard to the spectral radiance of light source distribute (such as spectral radiance) measure, and in this, as spectral radiant emittance calibration benchmark.Robot scaling equipment as shown in Figure 3.

In order to set up the analytical expression of the radiometric calibration of each spectrum channel, if the spectral radiance obtained in a certain channel sensor standard field of view of multi-spectral imager is Y (can be obtained by apparatus measures), the DN value that sensor exports is X, in radiation calibration coefficient, slope is the A (enlargement factor of instrument, be assumed to be linear), intercept is B (DC component of circuit or end level), then the pass that spectral radiance and image export between DN value is:

Y＝B+AX

In formula: the unit of Y is Wcm ^-2sr ^-1nm ^-1.

In the dynamic range of multi-spectral imager, regulate the radiance of light source, the spectral radiance value of one group of measured target and the relational expression (radiometric calibration formula) of multi-spectral imager output DN value can be set up:

L _j(λ _i)＝a _iDN _(j，i)+b _i

Wherein: L _j(λ _i) for radiance grade is j, spectrum segment is the standard sources radiance of i, known quantity; DN _{(j, i)}for radiance grade is j, spectrum segment is the output DN value of the multi-spectral imager of i, known quantity.A _i, b _ifor the radiation calibration coefficient that spectrum segment is i, undetermined.Work as a _i, b _iafter determining, above formula just can as the radiometric calibration formula of multi-spectral imager.Namely each spectrum channel has respective radiometric calibration formula.

Due to some spectral measurement passages (i.e. CCD passage) of the corresponding multi-spectral imager of wavelength i difference of radiometric calibration formula, therefore a _i, b _idetermination complete in this passage.Notice: by L _j(λ _i) and DN _{(j, i)}data group, therefore can according to these L _j(λ _i), DN _{(j, i)}the coefficient a of data group determination formula (4-19) _i, b _i.Method that linear fit fits such as can be adopted to determine a _i, b _i.

In some cases, there is nonlinear situation in the photoelectric response of CCD, now then can directly by L _j(λ _i) and DN _{(j, i)}data group fits out quafric curve as calibration formula:

Y＝P1 _ji+P2 _jiX+P3 _jiX ²

Be described above the basic skills of multi-spectral imager radiometric calibration.Such as, but the optical system parameter that the method is generally applicable to imaging spectrometer is fixing situation, the radiometric calibration of satellite remote sensing imaging spectrometer.But for multi-spectral imager of the present invention, due to the size of measured target with apart from very unfixing, therefore above-mentioned calibration scheme is not best selection.

Consider that multi-spectral imager is mainly used in spectral reflectance factor or the spectral reflectance (i.e. reflecting body) of measurement target, therefore can adopt the calibrating method being similar to reflectance spectrophotomete.

The radiometric calibration method of reflectance spectrophotomete is otherwise known as the calibration of spectral reflectivity, and its ultimate principle is similar with " comparing colour examining " method in spectrophotometry.

What is called compares the monochromatic radiation power that colorimetry method reflects on Same Wavelength by comparing some known spectral characteristic " standard " (object of reference) and sample quantitatively, thus measures spectral reflectance or the spectral radiance factor of sample.

When measuring reflected sample, be generally finished total reflection diffuser as reference standard.The reflectance of complete reflected diffusion body is 1 on each wavelength.But do not have a kind of real material to have such characteristic, the material close with its Nature comparison can only be selected as working stamndard.The present invention adopts lacklustre plane white model as object of reference, and its spectral reflectivity can be obtained by standard spectrophotometers measurement.

Relatively the colour examining principle of colorimetry method is summarized as follows:

If the instrumental response value that W (λ) is standard white plate; W ₀(λ) spectral reflectivity of standard white plate for demarcating in advance; The instrumental response value that S (λ) is sample (target); The spectral reflectivity that R (λ) is sample to be calculated;

φ ₀(λ) be incident radiation flux; φ _s(λ) and φ _w(λ) reflected flux received when being respectively test sample product and when surveying blank;

K is the conversion coefficient of instrument.According to definition, reflectivity is reflected flux and the ratio of incident flux, instrumental response value expression when can write out test sample product respectively and when surveying blank:

S(λ)＝kΦ _S(λ)＝kΦ ₀(λ)R(λ)

W(λ)＝kΦ _W(λ)＝kΦ ₀(λ)W ₀(λ)

Upper two formulas are compared the spectral reflectivity calculating formula obtaining sample:

$R\left(\mathrm{\λ}\right)=\frac{S\left(\mathrm{\λ}\right)}{W\left(\mathrm{\λ}\right)}\×W_0\left(\mathrm{\λ}\right)$

From above formula, in the spectral analysis that measurement image is carried out, we only need to obtain the instrumental response value W (λ) of on-gauge plate and the instrumental response value S (λ) of measured target, just can calculate the spectral reflectivity R (λ) of this target.

From analyzing above, in the actual measurement work of multi-spectral imager, in order to obtain spectral calibration data, require that target scene internal memory is at the standard white plate of a known spectra reflectivity or approximate standard white plate, the method is all in most of the cases feasible.

As previously mentioned, 8 passage multi-spectral imagers of the present invention utilize the raw data of 8 narrow-cut filters acquisitions for Wavelength distribution is in the spectral response value of eight wave bands of (420-940nm) from visible ray to infrared radiation.According to 8 original values of gained, we can only be linked to be a spectrum broken line, this broken line truly can not reflect the spectral characteristic of object, and want to carry out further spectral analysis and color synthesis to object scene, just must obtain a continuous print curve of spectrum, therefore we need to adopt interpolation algorithm to come original the smooth operation of spectrum broken line containing 8 data points to obtain other band spectrum response with this.

Light is mapped on object, may reflect, absorb and transmission.Wherein reflected flux is called reflectance with the ratio of incident flux.According to the reflection case of body surface to light, can be divided into

(1) reflected light observes reflection law, and from the injection of minute surface transmit direction, this part light is normal reflection ratio.

(2) radiation flux of incidence nondestructively all reflects away by persect reflecting diffuser, and reflectance equals 1, and has identical brightness in all directions.

(3) normal reflection and irreflexive combination, in certain solid angle, the size of reflected flux is namely relevant with incident direction, relevant with measurement direction again, i.e. spectral reflectance factor.Spectral reflectance factor is under specific illumination condition, on the direction that regulation solid angle limits, from the spectral radiant flux φ of the wavelength X of object _s(λ) with the spectral radiant flux φ of wavelength X that reflects from persect reflecting diffuser under the same terms _n(λ) ratio.When given solid angle ^ωbe exactly spectral reflectance ρ (λ) close to the spectral reflectance factor measured during 2 π.If what measure is spectral radiance, then that obtain is spectral radiance factor beta (λ).

Spectral reflectance, spectral reflection factor, spectral radiance factor can reflect that object selects the characteristic of reflection to incident light spectrum, and geometric condition when just measuring is different.The color characteristics of reflecting surface can have its spectral reflection characteristic to calculate, so the measurement of reflections off objects characteristic is of great significance in colorimetry.

When body surface reflectance is not with body surface variation in thickness, reflectance is also known as reflectivity.Reflectivity is relevant with chemical composition with the physical arrangement of body surface.Strict differentiation is not carried out to reflectance and reflectivity herein, be all called reflectivity.

No matter utilize human eye to measure or apparatus measures, light source, sample, detector (or human eye) are corresponding three key elements of output finally obtaining sample.Relation between this three, is called illumination and accepts geometric condition, is called for short geometric condition.The inconsistent meeting of geometric condition causes the difference of measurement result.For the unified measurement result CIE 15:2004 third edition recommends 10 kinds of geometric conditions to reflected sample:

(1) .di:8 ° (diffuse illumination, 8 degree of directions accept, and comprise specular reflection component)

(2) .de:8 ° (diffuse illumination, 8 degree of directions accept, and get rid of specular reflection component)

(3) .8 °: di (8 degree of directional lightings, diffuse reflection accepts, and comprises specular reflection component)

(4) .8 °: de (8 degree of directional lightings, diffuse reflection accepts, and gets rid of specular reflection component)

(5) .d:d (diffuse illumination, diffuse reflection accepts)

(6) .d:0 ° (diffuse illumination, 0 degree of direction accepts, and gets rid of specular reflection component)

(7) .45 ° a:0 ° (45 degree of ring illuminations, 0 degree of direction accepts)

(8) .0 °: 45 ° of a (0 degree of directional lighting, 45 degree of annulars accept)

(9) .45 ° x:0 ° (45 degree of directional lightings, 0 degree of direction accepts)

(10) .0 °: 45 ° of x (0 degree of directional lighting, 45 degree of directed acceptance)

Known according to introducing above, when multi-spectral imager work of the present invention, its optical system receives the optical radiation of distant object, and aperture angle very little (target is to the long angle of camera lens), therefore geometric condition is normal incidence, lighting condition is sunshine diffusion, namely close to the d:0 geometric condition that CIE recommends.Its measurement result strictly goes up and should be called spectral radiance factor.But in some document, be generally referred to as spectral reflectivity.

Usually there is such problem in real work: provide a collection of discrete sampling point, require to make a smooth curve by these points, to meet design requirement or to process.This problem is concluded in mathematical method: the value of known function on some points, seek its analytical expression.Above problem is exactly the interpolation problem that multi-spectral imager of the present invention needs the multispectral datacube solved.

In existing numerical analysis method, the approach solved the problem mainly contains two classes: a class is some sample value providing function f (x), select a functional form being convenient to calculate, as polynomial expression, fractional linear function and trigonometric polynomial etc., require that it is by known sampling point, determine function (x) being similar to as f (x) thus, so-called method of interpolation that Here it is.Another kind of method, after the form of selected approximate function, does not require that approximate function crosses known sampling point, and to only require that under certain meaning its total departure on these aspects is minimum, these class methods are called curve (data) fitting process.Interpolation algorithm conventional in current numerical analysis has Lagrange's interpolation, Newton interpolation, piecewise linear interpolation, Hermite interpolation and spline interpolation (Spline) etc.Spline interpolation wherein can obtain level and smooth curve, meets the actual conditions of multi-spectral imager, is therefore the interpolation algorithm selected by the present invention.

Eight passage multispectral imaging data processing methods provided by the invention provide a kind of complete comprehensive data processing step, to the guidance that technician is concrete very in detail, and have selected corresponding method for the feature of eight passage multispectral imaging data, reach good effect.

Be understandable that, for those of ordinary skills, can be equal to according to technical scheme of the present invention and inventive concept thereof and replace or change, and all these change or replace the protection domain that all should belong to the claim appended by the present invention.

Claims (2)

1. eight passage multispectral imaging data processing methods, is characterized in that comprising the following steps:

Step one, multispectral image pre-service;

The pre-registration of step 2, multispectral image;

Step 3, wavelength scaling;

Step 4, radiation calibration;

Step 5, the curve of spectrum build;

Step 6, multispectral image color calculate;

Step 7, multispectral image feature extraction;

Step 8, spectrum and color-match.

2. eight passage multispectral imaging data processing methods, is characterized in that comprising the following steps:

Step one, multispectral image pre-service, comprise spectrum picture noise remove, CCD homogeneity correction and CCD gamma correction etc.;

Described spectrum picture noise remove uses to face the territory method of average and median filtering method removal salt-pepper noise;

Described CCD homogeneity correction is by carrying out imaging to the planar target with uniform luminance distribution, then calculate according to the response of each pixel of CCD and extract the homogeneity correction coefficient of each pixel of CCD surface, and these coefficients being saved in the calibration file of software systems;

Described CCD gamma correction, by measuring the target with different brightness degree, then according to the response data of each brightness degree, adopts curve fit or other mathematical method to obtain non-linear correction factor;

The pre-registration of step 2, multispectral image, comprises Edge extraction and carries out image conversion according to the unique point after edge extracting;

Described Edge extraction is carried out by the following method:

For two dimensional image signal, first come smoothing with following Gauss function:

$G(x,y,\mathrm{\σ})=\frac{1}{2\mathrm{\π}{\mathrm{\σ}}^2}\exp (-\frac{1}{2{\mathrm{\σ}}^2}(x^2+y^2))$

G (x, y, σ) is a function with circular symmetry, its level and smooth effect controls by σ, owing to carrying out linear smoothing to image, makes g (x, y) be level and smooth after image, obtain: g (x, y)=G (x, y, σ) * f (x, y), wherein f (x, y) is level and smooth front image;

Because marginal point is the place that in image, gray-value variation is violent, the sudden change of this image intensity will produce a peak in first order derivative, or be equivalent to a generation zero cross point in second derivative, and be nonlinear along the second derivative of gradient direction, so substitute with Laplacian operator, namely use:

${\▿}^{2}g(x,y)={\▿}^2(G(x,y)*f(x,y))=({\▿}^{2}G(x,y,\mathrm{\σ}))*f(x,y)$

Zero cross point as marginal point, in formula for LOG wave filter,

${\▿}^{2}G(x,y,\mathrm{\σ})=\frac{{\∂}^{2}G}{\∂x^2}+\frac{{\∂}^{2}G}{\∂y^2}=\frac{1}{\mathrm{\π}{\mathrm{\σ}}^4}(\frac{x^2+y^2}{2{\mathrm{\σ}}^2}-1)\exp (-\frac{1}{2{\mathrm{\σ}}^2}(x^2+y^2))$

Above formula is exactly LOG edge detection operator, LOG operator is being similar to ganglia retinae receptive field spatial organization, can regard as and be made up of focus of excitation district and an inhibition surrounding zone, occur with the template form of 62 σ × 62 σ sizes, when σ gets different values, then the image border under available operators detection different scale, usually, we get σ >=1, carry out rim detection namely extract image border with above-mentioned LOG edge detection operator to image;

Described image conversion uses rigid transformation, and transformation for mula is as follows:

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\end{array}\right]=k\left[\begin{array}{cc}\cos \mathrm{\θ}& \sin \mathrm{\θ}\\ -\sin \mathrm{\θ}& \cos \mathrm{\θ}\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]+\left[\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\end{array}\right]$

In formula: (x, y) be the point of piece image and benchmark image, after conversion in corresponding second width image (x ', y '), and k, θ and Δ x and Δ y are the scale factor of the first width figure and the second width figure, twiddle factor and coordinate translation amount respectively, these parameters are obtained by manual operation;

Step 3, wavelength scaling, select cubic spline interpolation method, the SPL of respective numbers point known wavelength reflectivity data is set up according to the quantity of optical filter, i.e. spectral reflectivity function curve, according to obtained function curve, according to even step sizes, wavelength is split between every two known wavelength, just obtain the spectral reflectance data of other wavelength;

Step 4, radiation calibration, by various calibrated radiation source, quantitative relationship between different wave spectrum section sets up digital quantization value that the spectral radiance value at imaging spectrometer entrance pupil place and imaging spectrometer export, select the target light source of known spectra radiation characteristic as the measuring object of multi-spectral imager, then simulate duty and the environment of multi-spectral imager;

Step 5, the curve of spectrum build, and adopt cubic spline functions to set up spectral reflectance rate curve;

Step 6, multispectral image color calculate, and comprise computed image color tristimulus values, calculating coloured image RGB, computed image aberration and calculate near infrared luminance factor β value;

Step 7, multispectral image feature extraction, comprise feature selecting and feature extraction;

Step 8, spectrum and color-match, comprise the color-match based on tristimulus values, the color-match based on spectrum and Spectral matching.