CN101606082A - Statistical tomographic reconstruction based on the charged particle measurement - Google Patents

Abstract

The system, device, computer program and the method that are used for charged particle detection comprise the statistics reconstruction of the object volume diffuse density section of charged particle chromatographic data, determine the probability distribution of charged particle scattering to use statistics multiple scattering model, and use expectation maximization (ML/EM) algorithm to determine the essence PRML estimation of object volume diffuse density, to rebuild the object volume diffuse density.Can from the volume scattering density profile of being rebuild, discern the existence and/or the type of the object that occupies interested volume.The charged particle chromatographic data can come the cosmic rays μ meson chromatographic data of the μ meson tracker of self-scanning parcel, container, vehicle or goods.Can use the computer program that to carry out on computers to realize such method.

Description

Statistical tomographic reconstruction based on the charged particle measurement

The cross reference of provisional application

It is 60/855 that this PCT application requires (1) application number, 064, title is that " SYSTEMS; METHODS AND APPARATUS FOR PARTICLE DETECTION ANDANALYSIS AND FIELD DEPLOYMENT OF THE SAME " and the U.S. Provisional Patent Application of submitting on October 27th, 2006 and (2) application number are 11/771,169, title be " RADIATION PORTAL MONITOR SYSTEM AND METHOD " and in the right of priority of the U.S. Provisional Patent Application of submission on June 29th, 2007.

By reference the disclosure of above two applications is incorporated into here.

Statement about federal right

The present invention carries out under the government that is authorized by USDOE, contract number is DE-AC52-06NA25396 supports.Government has certain right in the present invention.

Technical field

Embodiment relates to particle detection, analysis and control field, and more specifically and not exclusively say, relate to method and system, described method and system is used to analyze the data from the charged particle detection system with a plurality of Position-Sensitive Detectors, and be used to rebuild pass the charged particle detection system, such as the track of the charged particle of cosmic rays μ meson.

Background technology

Charged particle chromatography (charged particle tomography) is based on the scattering (scattering) of charged particle.A kind of form of charged particle chromatography is the cosmic rays chromatography that relies on the scattering of cosmic rays μ meson (cosmic raymuon).Stable particle major part from the outer space is a proton, and these particles bombard the earth continuously.Described particle and the atomic interaction in upper atmosphere comprise the shower (shower) of the particle of many short-life pi-meson (pion) with generation, and the μ meson of decay of described pi-meson and generation longer life.The μ meson is mainly by Coulomb force and matter interaction, and do not have nuclear interaction, and compares more difficult radiation with electronics.They only lose energy lentamente by electromagnetic interaction.Therefore, many μ mesons arrive earth surface as the charged radiation of high-penetrability.At the μ on sea level meson flux is every square centimeter of μ meson of about per minute.

When material was passed in the motion of μ meson, its track (trajectory) was disturbed in the Coulomb scattering of the electric charge of subatomic particle.Though different materials character is depended in total bias, dominant parameters is nuclear atomic number Z and density of material.

A kind of improved method and system need be provided, be used for rebuilding interested volume from the μ meson or other charged particles that pass volume (volume).

Summary of the invention

Following summary of the invention is provided to help to relating to some technical characterictics to technology, device and the system of the particle of surveying the such charged particle of all like μ mesons, and the understanding that the statistics of the object volume diffuse density section (object volume scattering density profile) of charged particle chromatographic data is rebuild, rather than for comprehensive description.Can be to each side overall understanding of the present invention by whole instructions, claim, accompanying drawing and summary are done as a whole the acquisition.

Above-mentioned aspect of the present invention and other targets and advantage can be reached as described herein now.

According to an aspect, describe a kind of detection system and be used for surveying object volume via the charged particle that passes object volume.Described system comprises: first group of Position-Sensitive Detector, and it is positioned at first side of object volume, measures position and angle to the charged particle of object volume incident; Second group of Position-Sensitive Detector, it is positioned at second side of the object volume relative with first side, measures the position and the angle of the charged particle that leaves the object volume outgoing; And signal processing unit, it receives from the measuring-signal of first group of Position-Sensitive Detector with from the data of the measuring-signal of second group of Position-Sensitive Detector.The data that signal processing unit processes received are to produce the statistics reconstruction that object volume inner volume diffuse density distributes.

Signal processing unit can be configured to: (a) obtain corresponding to the scattering angle of the charged particle that passes object volume and the charged particle chromatographic data of estimation momentum (estimated momenta); (b), provide the probability distribution of charged particle diffuse density based on statistics multiple scattering model (multiple scattering model); (c) use (ML/EM) algorithm of expectation maximization (expectation maximization), determine that the essence PRML of object volume diffuse density is estimated (substantially maximumlikelihood estimate); And, export the object volume diffuse density of reconstruction (d) based on the estimation of essence PRML.

According to another aspect, a kind of method that is used for detection object volume the charged particle chromatographic data that obtains from object volume comprises: (a) obtain corresponding to the scattering angle of the charged particle that passes object volume and the predetermined charged particle chromatographic data of estimation momentum; (b), provide the probability distribution of charged particle scattering based on statistics multiple scattering model; (c) use expectation maximization (ML/EM) algorithm, determine that the essence PRML of object volume diffuse density is estimated; (d) export the object volume diffuse density of being rebuild; And, if necessary, (e) make decision based on the object volume diffuse density of being rebuild.

This method allows the user to discern the existence and/or the type of the object that occupies (occupying) interested volume from the volume scattering density profile of being rebuild.Various application comprise the cosmic rays μ meson chromatography that is used for various Homeland Security inspections application, and in described application, vehicle or goods can scan by μ meson tracker.

The charged particle chromatographic data can comprise from the chromatographic data that is generated by other sources of cosmic rays or some, gather such as the charged particle of μ meson.

Based on the object volume diffuse density of being rebuild make decision (making a decision) can comprise: based on the object volume diffuse density of being rebuild, come to (1) of the target object in the object volume exist and (2) type at least one make decision.

The probability distribution that is provided for the charged particle scattering of expectation maximization (ML/EM) algorithm can comprise: the 2D probability distribution that (g) obtains charged particle based on the predefined diffuse density of homogeneous (homogenous) object; (h) obtain the 3D probability distribution of charged particle based on the 2D probability distribution; (i) obtain to pass probability distribution via the scattering of a plurality of charged particles of (characterized via basis functions) heterogeneous body (non-homogenous) object volume of basis function characterization; (j) based on the definition of multiple scattering and the scattering and the momentum survey of charged particle, determine the probability distribution of described multiple scattering (multiple scattering).

The 2D probability distribution that obtains charged particle based on the predefined diffuse density of homogeneous body can comprise: (k) determine the expected mean square scattering of the diffuse density of material as the charged particle of the unit depth (unit depth) of passing material; (l) be similar to the scatter distributions of (approximating) charged particle based on Gauss model; And the scattering and the displacement that (m) are similar to the charged particle ray based on (correlated) 2D Gaussian distribution that is associated.

The 3D probability distribution that obtains charged particle based on the 2D probability distribution can comprise: add coordinate and define three-dimensional path length; Calculate the 3D displacement; And definition 3D covariance matrix.

The probability distribution that acquisition is passed via the scattering of a plurality of charged particles of the heterogeneous body object volume of basis function characterization can comprise: the 3D grid (3Dgrid) of setting up the basis function of the discrete diffuse density estimation of representative; The covariance matrix of the scattering/displacement of the μ meson that each is independent is defined as the function of raypath and diffuse density estimation; And the product that the probability distribution of a plurality of charged particles is defined as the probability of independent charged particle.

Use expectation maximization (ML/EM) algorithm to determine that the essence PRML estimation of object volume diffuse density can comprise: to gather the scattering of each charged particle and the measurement of momentum; Estimate the interactional geometry (geometry ofinteraction) of each basis function of each charged particle and statistics scattering model; Right to each charged particle basis function, determine weight matrix: W _IjTo guess that (guess) comes the diffuse density in each basis function of initialization; And the approximate PRML of finding the solution (solving) object volume capacity (contents) is iteratively separated (solution); Wherein separate when becoming, stop iterative process less than predetermined tolerance value at the predetermined number place of iteration or when described.

Use expectation maximization (ML/EM) algorithm to determine that the essence PRML estimation of the diffuse density of object volume can comprise: to gather each charged particle i=1 to the scattering of M and the measurement of momentum (Δ θ _x, Δ θ _y, Δ x, Δ _y, p _r ²) _iEstimate each μ meson and each voxel (voxel) j=1 interactional geometry to N: (L, T) _IjRight to each charged particle voxel, Wij is calculated as with weight matrix $W_{\mathrm{ij}}\&equiv;\left[\begin{array}{cc}{L}_{\mathrm{<figure-callout\; id="2"\; label="ij\; L\; ij"\; filenames="A2007800487880033H2.png"\; state="\{\{state\}\}">ij</figure-callout>}}& L_{\mathrm{ij}}^{}{}_2^{}/2+L_{\mathrm{ij}}{T}_{\mathrm{<figure-callout\; id="2"\; label="ij\; L\; ij"\; filenames="A2007800487880033H2.png"\; state="\{\{state\}\}">ij</figure-callout>}}& \\ L_{\mathrm{ij}}^{}{}_2^{}/2+L_{\mathrm{ij}}{T}_{\mathrm{ij}}& L_{\mathrm{ij}}^3/3+L_{\mathrm{ij}}^2{T}_{\mathrm{ij}}+L_{\mathrm{<figure-callout\; id="2"\; label="ij\; T\; ij"\; filenames="A2007800487880033H2.png"\; state="\{\{state\}\}">ij</figure-callout>}}{T}_{\mathrm{ij}}^{}{}_2^{}\end{array}\right],$ The conjecture λ of initialization diffuse density in each voxel _{J, old}And use stops, and criterion process (a stopping criteria process) is next to be provided with λ to whole voxels _{J, old}=λ _{J, new}

Expectation maximization (ML/EM) algorithm can comprise average update rule or intermediate value update rule.

According to another aspect, be used for comprising: (a) obtain corresponding to the scattering angle of the charged particle that passes object volume and the charged particle chromatographic data of estimation momentum in the computer-implemented method of the charged particle chromatographic data detection object volume that obtains from object volume; (b) be provided for the probability distribution of the charged particle diffuse density of expectation maximization (ML/EM) algorithm, and described probability distribution is based on statistics multiple scattering model; (c) use expectation maximization (ML/EM) algorithm to determine the essence PRML estimation of object volume diffuse density; And (d) the output object volume diffuse density of being rebuild.Can make decision based on the object scattering volume density of being rebuild.

According to another aspect, computer program comprises the data carrier storage instruction that computing machine can be used, and when carrying out by computing machine, the method that the statistics that described instruction causes computing machine to carry out the object volume density profile of charged particle chromatographic data is rebuild, described method comprises: (a) obtain corresponding to the scattering angle of the charged particle that passes object volume and the predetermined charged particle chromatographic data of estimation momentum; (b) be provided for the probability distribution of the charged particle scattering of expectation maximization (ML/EM) algorithm, and described probability distribution is based on statistics multiple scattering model; (c) use expectation maximization (ML/EM) algorithm to determine the essence PRML estimation of object volume density; And (d) the output object volume diffuse density of being rebuild.

Description of drawings

In the accompanying drawings, the same reference numerals that runs through the view of separation is censured identical or intimate element, and accompanying drawing is incorporated in the instructions, the part of formation instructions, and described accompanying drawing further illustrates the present invention, and, be used for explaining principle of the present invention with detailed description of the present invention.

Fig. 1 has illustrated an example of μ meson chromatography notion;

Fig. 2 has illustrated and has been used to determine the scattering of Coulomb scattering and the two-dimensional projection of displacement;

Fig. 3 has illustrated the parameter of the two-dimensional projection of the scattering that is used to adjust the 3D scattering model and displacement;

Fig. 4 has illustrated the scattering of passing multilayer material;

Fig. 5 to have illustrated and to have used closest approach (point of closest approach) in order calculating in the path of the projection shown in Fig. 3;

Fig. 6 has illustrated the automated system that the statistics of the volume scattering density profile that is used for μ meson chromatography is rebuild;

Fig. 7 has illustrated the skeleton view of the object of simulation;

Fig. 8 has illustrated the vertical view of the object of simulation;

Fig. 9 has illustrated to have two lines of hypothesis connect Gauss's scattering analogue of estimated path at PoCA point place reconstruction;

Figure 10 has illustrated reconstruction scattering, simulated data of using the afterbody (tails) that non-Gauss is arranged;

Figure 11 has illustrated reconstruction scattering, simulated data of using the afterbody that non-Gauss is arranged via median method;

Figure 12 has illustrated the main object in the passenger traffic railway carriage or compartment car (passenger van) of simulation;

Figure 13 has illustrated the reconstruction via one minute simulation μ meson exposure of the passenger traffic railway carriage or compartment car of Mean Method;

Figure 14 has illustrated the reconstruction via the scene of the passenger traffic railway carriage or compartment car of median method;

Figure 15 has illustrated the process flow diagram according to an embodiment, and this process flow diagram has usually been summarized the method for the statistics reconstruction of the volume scattering density profile that is used for the charged particle chromatography;

Figure 15 A has illustrated the process flow diagram according to an embodiment, and this process flow diagram has been summarized the example of process of probability distribution of the scattering of the μ meson that uses multiple statistics scattering analogue to estimate to pass object volume;

Figure 15 B has illustrated the process flow diagram according to an embodiment, and this process flow diagram has been summarized the example of estimating single μ meson process of the desired probability distribution of scattering in 2D based on the predefined diffuse density of object;

Figure 15 C has illustrated the process flow diagram according to an embodiment, and this process flow diagram has been summarized the example that statistical model is expanded to the process of 3D;

Figure 15 D has illustrated the process flow diagram according to an embodiment, and this process flow diagram has been summarized the example of process of the probability distribution of the scattering of a plurality of μ mesons of determining to pass heterogeneous material and displacement; And

Figure 15 E has illustrated the process flow diagram according to an embodiment, and this process flow diagram has been summarized the example of process of the likelihood of the estimation density profile that uses expectation-maximization algorithm to maximize object volume.

Embodiment

Particular value of being discussed in these unrestriced examples and configuration can change, and only are to propose at least one embodiment of the present invention is described, rather than in order to limit the scope of the invention.

Technical characterictic described in this application can be used for making up various particle detection system.For example, the particle detection system that is used to survey as the μ meson of charged particle can comprise: object storage district (object holding area) is used for placing with checked object; First group of position sensing muon detection device is positioned at first side in object storage district, measures position and angle to the μ of object storage district incident meson; Second group of position sensing muon detection device is positioned at second side relative with first side, the object storage district, measures the position and the angle of the μ meson that leaves the outgoing of object storage district; And signal processing unit, it can comprise, for example, microprocessor receives from the measuring-signal that enters the μ meson of first group of position sensing muon detection device with from the data of the measuring-signal of the μ meson of the outgoing of second group of position sensing muon detection device.As an example, in first group and the second group of particle detector each can be implemented as and comprise that drift tube (drift tubes), described drift tube are arranged to allow three charged particle position measurements and three the charged particle position measurements in being different from the second direction of first direction in first direction at least at least.Signal processing unit is configured to based on measured μ meson and enters position and angle with outgoing, analyze scattering behavior, to obtain the chromatography section or the space distribution of the scattering center (scattering centers) within the object storage district by the μ meson that scattering caused in the material of μ meson within the object storage district.The chromatography section of the scattering center that is obtained or space distribution can be used to be disclosed in the existence of the one or more objects in the object storage district or not exist, and described object comprises nuclear material or nuclear device such as the material with high atomic number.Each position sensing muon detection device can realize with various configurations, comprises drift unit (drift cells), as be full of can be by the drift tube of the gas of μ meson ionization (ionize).The cosmic rays μ meson that such system can be used for utilizing nature is as the source that is used for surveying at the μ meson of the one or more objects in object storage district.

In one embodiment, the method and system of rebuilding according to the statistics of the volume scattering density profile that is used for the charged particle chromatography of the embodiment of explanation provides a kind of scheme, and wherein the image or the model of object rebuild in the scattering of passing the cosmic rays charged particle of object based on motion.

Compare with the material (such as water, plastics, aluminium and steel) of forming more common object, special nuclear material (SNM) and the material (such as lead and tungsten) that can be used as good gamma ray shielding are understood influence track more consumingly.For the cosmic rays charged particle, cosmic rays μ meson particularly, each μ meson all carries the information of the object that penetrates about it, and by measuring the scattering of a plurality of μ mesons, the character that just can detect these objects.Specifically, people can survey high Z (high-Z) object in more typical low Z (low-Z) and medium Z (medium-Z) material.

In order to explain the various technical characterictics of rebuilding according to the statistics of the volume density section that is used for the charged particle chromatography of the embodiment of explanation, will be at first with reference to μ meson chromatography notion, the example of μ meson chromatography notion has been described in Fig. 1.

Position-Sensitive Detector group 10 be configured in by on (imaged) object volume 11 of imaging and under, with position and the angle (solid line by the band arrow illustrates) that all enter and charged particle trajectory outgoing (track) 12 is provided.Being arranged to be provided the position and the angle of the charged particle trajectory that enters by the two or more groups Position-Sensitive Detector 10 on the volume of imaging.These detectors are measured the position of charged particle in coordinate two quadratures or non-orthogonal.The position and the angle of the charged particle of another group Position-Sensitive Detector 10 record outgoing.Side detector (not shown) can be used for surveying more level towards the track of the charged particle of (orientate).From (coincident) of unanimity that enter with measurement outgoing calculate the scattering angle of each charged particle trajectory.From detector self or be placed in the slight scattering that takes place in the layer of scatterer (scatterers) of two groups of known properties between the Position-Sensitive Detector plane, estimate the charged particle momentum.

An example of position sensing charged particle detector is for being filled the drift tube of working gas (operating gas).Drift tube can be a cylindrical tube, and is filled the detector gas such as argon-isobutane (Argon-Isobutane), surveys the charged particle such as the cosmic rays of μ meson.With approximately+the positive HV (High Voltage, high voltage) of 2-3 kilovolt is applied to the center anode line that extends along the length of cylindrical tube and external surface grounded with pipe, so that there is the high-voltage electrostatic field.When charged particle and gas atom interact, from described atom, discharge many electronics along the straight line of the chord length of passing pipe (chord).Electrostatic field makes electronics " string " (string) to the anode line drift of positively charged, and anode line is read electronically by the TDCS (time-to-digital converter, time-digital quantizer) of data acquisition electronic device.Every group of detector can be a plurality of drift tubes, described a plurality of drift tube be arranged with allow in first direction at least three charged particle position measurements and be different from first direction, and can with at least three charged particle position measurements in the second direction of first direction quadrature.

Provide signal processing unit in the system in Fig. 1, for example, computing machine receives by at the μ meson that enters of the detector measurement on the object volume and by the data of the signal of the μ meson of the outgoing of the detector measurement under object volume.Described signal processing unit is configured to analyze scattering behavior by the μ meson that scattering caused in the volume based on enter and position outgoing of the μ meson of measuring and angle, to obtain the chromatography section or the space distribution of the scattering center within volume.The chromatography section of the scattering center within volume that is obtained or space distribution can disclose the existence of object in the volume or not exist.In some embodiments, can on the side of volume, realize additional drift tube detector, forming box or four sides structure (four sided structure), and wrap up, vehicle or cargo container can enter wherein, so that scan with system.Thereby, can use the multiple scattering of cosmic rays μ meson in the background of normal goods, optionally to survey high z material (z-material).Advantageously, present technique is passive, does not discharge any radiation dose on background, and to the fine and close object of high z for optionally.Can realize the chromatography processing section of signal processing unit at the local computer (on-premisecomputer) that is arranged in the position identical with detector 10.Replacedly, can in remote computer, realize the chromatography processing section of signal processing unit, and this remote computer is connected on the computer network such as dedicated network or common network (such as the internet).

In the illustrative embodiment of Fig. 1, charged particle is cosmic rays μ meson or other cosmic rays charged particles, and Position-Sensitive Detector 10 is for being full of the drift unit of the working gas that is used for sensing (sensing) charged particle.For example, can realize the unit that drifts about with drift tube, described drift tube has the center anode line along the longitudinal extension of each drift tube.Yet, can use the position sensitive detector that is different from the unit that drifts about to survey the charged particle that is different from the μ meson.In addition, can generate charged particle by the source that is different from cosmic rays.For example, the μ meson can generate as low-intensity beam (lowintensity beam) from accelerator.

The μ meson ratio that penetrates fine and close object (black track) passes the strong more scattering (stronger) of μ meson of air (grey track).From a plurality of orbit measurements, can rebuild the object geometry and the electron density of all material.The μ meson that passes volume is to depend on the mode scattering of the material that the μ meson is passed.

The volume of being checked by the processing unit of system among Fig. 1 (for example, parcel, container or vehicle) in can comprise the processing of the measurement of cosmic rays μ meson: rebuild the track that the μ meson passes volume, the momentum of the μ meson that enters based on the signal measurement of the detector on each side of coming comfortable volume, and the space distribution of the diffuse density of definite volume.These and other results can be used for making up the various character of chromatography section and measurement volumes, as detection of a target object.

For example, the reconstruction of track of passing the charged particle of the detector 10 with one group of drift tube can comprise: (a) obtain representative by the identifier of the drift unit of charged particle bump and the bump signal of corresponding collision time; (b) will be identified as identical (in-time) a plurality of drifts unit bump marshalling (grouping) of time that is associated with the track of the specific charged particle that passes described detector; (c) estimate that initially described specific charged particle clashes into the time zero value (time zero value) in the moment of drift unit; (d) based on the estimation of time zero value, drift time translation data and bump time determine the drift radius; (e) with the track fitting (fitting) of linearity to drift radius corresponding to special time value at zero point; And the time zero value of (f) searching for and selecting to be associated, and the error in zero point computing time and the orbit parameter with the best track fitting of carrying out for specific charged particle.This linear track of rebuilding the reconstruction that the charged particle that passes charged particle detector is provided based on the track of time zero match, and needn't use other fast detector of fast detector (if any the photomultiplier of scintillater oar (scintillator paddles)) or some, described other fast detector detects nearest several nanoseconds with the μ meson to passing of this device, so that time zero is provided.

Again for example, owing to can comprise based on the processing of the momentum of the μ meson that enter or outgoing of the signal measurement from the detector among Fig. 1 10: (a) a plurality of Position-Sensitive Detectors of configuration come scattering to pass charged particle there; (b) measurement that scattering comprises at least three positions of the charged particle that obtains scattering is wherein measured in the scattering of charged particle in the sensing detector of measuring position; (c) determine at least one track of charged particle from the measurement of position; And at least one momentum survey of (d) determining charged particle from described at least one track.Present technique can be used for determining based on the track of charged particle the momentum of charged particle, and the track of charged particle is determined by the scattering of charged particle in the Position-Sensitive Detector self, and does not use sheet metal additional in the detector.

The example system that the statistics of the object volume diffuse density section of charged particle chromatographic data is rebuild and the details of method are provided below.

The example of the automated system of rebuilding according to the statistics of a volume scattering density profile embodiment, the charged particle chromatography has been described in the block diagram in Fig. 6.Automated system 50 has by controller 51 adaptive and that arrangement is analysed data 54 with receiving belt electrochondria sublayer.For example, the charged particle chromatographic data can be a μ meson chromatographic data, described μ meson chromatographic data uses the charged particle detector 1 of Fig. 1, or replacedly, use any have be configured to other charged particle detectors that enable the Position-Sensitive Detector of the tracking of the charged particle that passes volume, and from the measurement of μ meson, determine.As a result, μ meson or other charged particle chromatographic datas can be used to extract or determine to pass the scattering angle of the μ meson of object volume or other charged particles and estimate momentum.

Automated system 50 comprises the statistics reconstructor module 56 that is stored on the controller.Reconstructor module 56 is responsible for the volume scattering density profile that statistics is rebuild the chromatography of μ meson or other charged particles.Module 56 can be implemented as software module or hardware module.

In the illustrative embodiment of the automated system 50 of Fig. 6, use one or morely to form controller 51 based on the system of the computer processor unit (CPU) of link operationally or such as other systems based on microprocessor based on the system of digital signal processor such as computing machine (PC).Controller can be the single standard computing machine, but in order to reach real-time results, controller typically is included in enough to provide on the number and reaches a group of planes required processing power of real-time results, parallel processing computer (farm) (not shown).For example, controller can comprise, for example, and 20 CPU.The scan volume of muon detection device is big more, and desired resolution is good more, then needs a big more process computer group of planes (processing computer farm).

Operating system is moved on controller 51, and described operating system can be the operating system of commercially available or open source code, includes but not limited to from Apple, Windows, Linux or Unix or other operating system that may develop in the future.Because the instruction of operating system and application or program are stored in the memory device such as hard disk drive.And, in automated system 50, follow the tracks of reconstructor module 56 and be the software of the data carrier form that the computing machine of storage instruction can use, when being carried out by controller, described instruction makes controller carry out the method for rebuilding according to the statistics of volume scattering density profile embodiment, that be used for the charged particle chromatography that illustrates.Module can be installed in controller this locality as shown in Figure 6, or moves from remote location via the network that is couple to controller.Those skilled in the art are to be understood that the mode of this module of multiple realization.

If desired, automated system 50 also comprises display 58, and it operationally is couple to controller 51, because to image or the data of user's demonstration by the object density section of system reconstructing.If desired, the user interface (not shown) can be operatively attached to disposal system, comes the operational processes system to allow operating personnel.

The example of the embodiment of automated system 50 is only represented in the explanation that it will be understood by those skilled in the art that Fig. 6, and embodiment is not limited thereto.For example, in the reconstructor functions of modules some or all can be embodied as hardware, and not use microprocessor such as the analog or digital circuit.

With reference to Figure 15, its illustrated according to an embodiment, general introduction usually is used for the process flow diagram of the method that the statistics of the volume density section of charged particle chromatography rebuilds.Shown in treatment step 101, method 100 is by obtaining to come initialization corresponding to the predetermined charged particle chromatographic data of the charged particle trajectory position of passing object volume, scattering angle and estimation momentum.For example, can obtain predetermined charged particle chromatographic data from the detector of Fig. 1.Afterwards, shown in treatment step 102, based on multiple statistics scattering model, provide by the space distribution of diffuse density (will to give a definition) representative, a plurality of charged particles pass the probability distribution of the scattering of object volume.Then, shown in treatment step 103, use expectation-maximization algorithm to determine the PRML estimation of object volume diffuse density section.Then, the volume scattering density profile of being rebuild is exported with make decision (treatment step 104).Shown in treatment step 105, the process of making decision is optionally, and can be to be used to discern the existence of the object that occupies volume and/or the process of type.The process of making decision can relate to the human interpretation of image of density profile of the reconstruction of representing object volume and/or the robotization decision by additional algorithm.

The method of embodiment and automated system allow to carry out based on the data that provided by many charged particles the discrete tomographic reconstruction of interested volume.The example of the expectation maximization of iteration (EM) algorithm is used to seek the PRML of the density profile of object and estimates.The method and system of embodiment allows user's identification from the volume density section of rebuilding to occupy the existence and/or the type of the object of interested volume.Various application comprise the cosmic rays μ meson chromatography that is used for various Homeland Security inspections application, wherein can be by μ meson tracker (tracker) scanning vehicle or goods.The method of operation instruction embodiment and automated system, consequent μ meson chromatographic data can be used for rebuilding and showing the density profile of vehicle or goods, to allow any object that constitutes a threat to of identification.

PRML is being used for the medical image reconstruction, especially, in the time of for PET and SPECT reconstruction, several important differences hinder the use for the standard method of those application and developments.At first, measuring-signal-scattering angle-be at random, and have null average (mean) and by the standard deviation (deviation) of the property definition of the material that penetrates.Secondly, cosmic rays μ meson is not from discrete direction or the angle that provides, and also almost extends to the broadness of local horizon angle, continuous angular distribution but have around the zenith angle.At last, the track of μ meson is not straight; The bending of track makes us can find the rough position of strong scatterer.The EM algorithm is efficient with calculating flexibly, and its application to the geometry of complexity can be described.

To describe treatment step 102 to 104 according to an embodiment now, wherein data are the cosmic rays μ meson chromatographic datas that obtain from the detector of the measuring μ meson, Fig. 1 that passes volume.

In the process flow diagram of Figure 15 A, summarized according to process (treatment step a 102) embodiment, that use multiple statistics scattering model to provide the estimated probability of the scattering that the μ meson passes object volume to distribute.Shown in treatment step 110 to 113, described process has four main ingredients.At first, estimation is based on the 2D probability distribution predefined diffuse density, single μ meson (treatment step 110) of homogeneous body.Then, the 2D distributed expansion is arrived 3D (treatment step 111).Next, in treatment step 112, use the voxel basis function to express the heterogeneous body object volume, and express probability distribution given voxelization (voxelized) scatterings diffuse density, a plurality of μ mesons.At last, probability distribution expression formula is expanded to the μ meson scattering and the momentum survey (treatment step 113) of limited precision.

Use the multiple scattering statistical model to realize treatment step 110 to 113, and will be at first with reference to the scattering in the individual layer homogeneous material, describe described model with reference to the scattering in the heterogeneous material then.

As illustrated in fig. 2, the cosmic rays μ meson that passes material experiences multiple Coulomb scattering (multiple Coulomb scattering), and Fig. 2 explanation is used to describe the two-dimensional projection of the scattering and the displacement of multiple Coulomb scattering.In this and other accompanying drawing, for illustrative purposes, the amplitude of greatly exaggerative scattering.With the μ meson of incident towards and position versus, can come the μ meson track of characterization outgoing by scattering angle and displacement.Typical scattering angle is tens milliradians (milliradian) (1 milliradian ≈, 0.06 degree), and the scattering angle bigger than the several years is very rare.The distribution of the central authorities 98% of scattering angle can be approximately zero-mean Gaussian (zero-mean Gaussian).

$f_{\mathrm{\&Delta;\&theta;}}\left(\mathrm{\&Delta;\&theta;}\right)\&cong;\frac{1}{\sqrt{2\mathrm{\&pi;}}{\mathrm{\&sigma;}}_{\mathrm{\&theta;}}}\exp (-\frac{{\mathrm{\&Delta;\&theta;}}^2}{{2\mathrm{\&sigma;}}_{\mathrm{\&theta;}}^2}),$ Formula (1)

Compare with Gaussian though distribute really and to have heavier or bigger afterbody.Can come the width of expression and distribution approx according to the character of material.As people such as S.Eidelman, Phys.Lett, vol.B592, p.1, that is commented in 2004 " Review of particle physics " is such, many researchers have provided the empirical representation that is used for scattering as the function of various material characters, by reference its disclosure are incorporated into here.A simple especially form is

${\mathrm{\&sigma;}}_{\mathrm{\&theta;}}\&cong;\frac{15\mathrm{MeV}}{\mathrm{\&beta;cp}}\sqrt{\frac{H}{L_{\mathrm{rad}}}}.$ Formula (2)

Here, p is the particle momentum in MeV/c, and H is a depth of material, and L _RadBe the radiation length of material, β c is speed (c is the light velocity), and uses the approximate of i, β=1.Radiation length reduces with the increase of atomic number and density of material.We set up nominal (nominal) μ meson momentum, p ₀, and will have radiation length L _RadThe diffuse density of material be defined as:

$\mathrm{\&lambda;}\left(L_{\mathrm{rad}}\right)\&equiv;{\left(\frac{15}{p_0}\right)}^2\frac{1}{L_{\mathrm{rad}}}.$ Formula (3)

Thereby the diffuse density λ of material represents all square scattering angle of the μ meson of momentum unit depth, that have nominal that passes that material.For example, for some material, in the value of the diffuse density of every centimetre of square milliradian be: aluminium is approximately 3, iron is approximately 14, and uranium is approximately 78.Thereby the variance of passing the scattering of μ meson material, that have momentum p with diffuse density λ and depth H is

${\mathrm{\&sigma;}}_{\mathrm{\&theta;}}^2=\mathrm{\&lambda;H}{\left(\frac{p_0}{p}\right)}^2.$ Formula (4)

Order

${\mathrm{<figure-callout\; id="2"\; label="p\; r"\; filenames="A2007800487880033H2.png"\; state="\{\{state\}\}">p</figure-callout>}}_r^{}={(p_0/p)}^2,$ Formula (5)

Draw

${\mathrm{\&sigma;}}_{\mathrm{\&theta;}}^2=\mathrm{\&lambda;}{\mathrm{<figure-callout\; id="2"\; label="Hp\; r"\; filenames="A2007800487880033H2.png"\; state="\{\{state\}\}">Hp</figure-callout>}}_r^{}.$ Formula (6)

With displacement x be associated with scattering angle Δ θ (correlated).Consider that in the lump scattering angle and displacement provide the information of the local contribution of scatters person's (contributors) who points out in larger volume position, and be (kinks) pointed as " bending " in path among Fig. 1.The distribution character of scattering angle and displacement can be turned to associating Gaussian (jointly Gaussian) with zero-mean, and

${\mathrm{\&sigma;}}_{\mathrm{\&Delta;x}}=\frac{H}{\sqrt{3}}{\mathrm{\&sigma;}}_{\mathrm{\&Delta;\&theta;}},$ Formula (7)

${\mathrm{\&rho;}}_{\mathrm{\&Delta;\&theta;\&Delta;x}}=\frac{\sqrt{3}}{2}.$ Formula (8)

We can be expressed as covariance matrix

$\mathrm{\&Sigma;}\&equiv;\left[\begin{array}{cc}{\mathrm{\&sigma;}}_{\mathrm{\&Delta;\&theta;}}^2& {\mathrm{\&sigma;}}_{\mathrm{\&Delta;\&theta;\&Delta;x}}\\ {\mathrm{\&sigma;}}_{\mathrm{\&Delta;\&theta;\&Delta;x}}& {\mathrm{\&sigma;}}_{\mathrm{\&Delta;x}}^2\end{array}\right]=\mathrm{\&lambda;}\left[\begin{array}{cc}\mathrm{<figure-callout\; id="2"\; label="H\; H"\; filenames="A2007800487880033H2.png"\; state="\{\{state\}\}">H</figure-callout>}& H^{}{}^2/2\\ {\mathrm{<figure-callout\; id="2"\; label="H"\; filenames="A2007800487880033H2.png"\; state="\{\{state\}\}">H</figure-callout>}}^2/2& H^3/3\end{array}\right]{\mathrm{<figure-callout\; id="2"\; label="p\; r"\; filenames="A2007800487880033H2.png"\; state="\{\{state\}\}">p</figure-callout>}}_r^{}.$ Formula (9)

Order

$A\&equiv;\left[\begin{array}{cc}\mathrm{<figure-callout\; id="2"\; label="H\; H"\; filenames="A2007800487880033H2.png"\; state="\{\{state\}\}">H</figure-callout>}& H^{}{}^2/2\\ {\mathrm{<figure-callout\; id="2"\; label="H"\; filenames="A2007800487880033H2.png"\; state="\{\{state\}\}">H</figure-callout>}}^2/2& H^3/3\end{array}\right],$ Formula (10)

Draw

$\mathrm{\&Sigma;}=\mathrm{\&lambda;A}{\mathrm{<figure-callout\; id="2"\; label="p\; r"\; filenames="A2007800487880033H2.png"\; state="\{\{state\}\}">p</figure-callout>}}_r^{}.$ Formula (11)

Notice above-mentioned, can be as summarizing in the process flow diagram of Figure 15 B, describing the probability distribution (treatment step 110) that obtains scattering single μ meson scattering, 2D according to an embodiment.Shown in treatment step 150,, the diffuse density of material is defined as p according to formula (3) ₀The μ meson of=3GeV/c passes the expected mean square scattering of the unit depth of described material.Then, shown in treatment step 151, formula (1,5,6), Gaussian approximation is carried out in the RMS scattering.At last, shown in treatment step 152, summed up scattering and Displacements Distribution that the 2D Gaussian distribution that is associated by zero-mean is come approximatelyc ray via formula (10,11).

In three-dimensional, consider that the y coordinate that is orthogonal to x comes the characterization scattering, and with reference to scattering angle Δ θ _xWith Δ θ _y, and displacement x and Δ y.To the bias on x and y plane for independently and same (the seeing Eidelman etc.) that distributes.Above development (development) based on coordinate system towards the direction of the μ meson that is orthogonal to incident.In the 3-D model, we must consider the 3-D path, and displacement measurement are projected to the plane in the μ meson path that is orthogonal to incident.Be used for adjusting Fig. 3 of the parameter of 3-D scattering model in explanation, illustrated to become θ with vertical direction _{X, 0}The μ meson of projection angle incident.

In order to help to understand this 3-D scattering, the projection angle θ that the imagination is associated from page quadrature y coordinate outwardly _{Y, 0}Be useful.Pass (the non-scattering) point (x of described layer to institute's projection _p, y _p) the straight line in μ meson path extend (that is 3-D path) and be

$L=\mathrm{<figure-callout\; id="1"\; label="H"\; filenames="A2007800487880032H1.png,A2007800487880036H1.png"\; state="\{\{state\}\}">H</figure-callout>}\sqrt{1+{\tan }^2{\mathrm{\&theta;}}_{x,0}+{\tan }^2{\mathrm{\&theta;}}_{y,0}}\&equiv;{\mathrm{HL}}_{\mathrm{xy}}.$ Formula (12)

With the x position and the viewpoint definition of the μ meson of outgoing is (x ₁, θ _{X, 1}), order then

Δ θ _x=θ _{X, 1}-θ _{X, 0}. formula (13)

Though the x displacement of measuring may be calculated x _m=x ₁-x _pBut we must rotate to this measurement the plane that is orthogonal to raypath, and adjust for the 3-D path.Displacement is defined as

${\mathrm{\&Delta;}}_x=(x_1-x_p)\cos \left({\mathrm{\&theta;}}_{x,0}\right)L_{\mathrm{xy}}\frac{\cos (\mathrm{\&Delta;}{\mathrm{\&theta;}}_x+{\mathrm{\&theta;}}_{x,0})}{\cos \left(\mathrm{\&Delta;}{\mathrm{\&theta;}}_x\right)},$ Formula (14)

Wherein, consider the 3-D paths for middle two, and last with measurement project to suitable towards.

At last, the covariance weight is redefined into

$A\&equiv;\left[\begin{array}{cc}\mathrm{<figure-callout\; id="2"\; label="L\; L"\; filenames="A2007800487880033H2.png"\; state="\{\{state\}\}">L</figure-callout>}& L^{}{}^2/2\\ {\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="A2007800487880033H2.png"\; state="\{\{state\}\}">L</figure-callout>}}^2/2& L^3/3\end{array}\right].$ Formula (15)

Then, need proceed in a similar manner, and formula (11) definition is to the covariance matrix of the two and coordinate scattering to scattering and displacement.Two quadratures, independently carry out scatterometry in the horizontal coordinate.In order to simplify mark, we develop (develop) to the only analysis of a coordinate.After a while with the incompatible information of discussion group from two coordinates.We must be noted that this model is effective (valid) to " little " scattering angle and displacement.The second order term that model is ignored in deriving may become remarkable for large angle scattering.

Notice above-mentionedly,,, obtain to extend to the statistical model of 3D as summarizing in the flow process of Figure 15 C according to an embodiment.At first, add the y coordinate, and definition three-dimensional path length (treatment step 160, formula (12)).Next, in treatment step 161, (13-14) calculates the 3D displacement according to formula.At last, express 3D covariance matrix (treatment step 162) according to formula (15).

To the heterogeneous body volume of material, for the purpose of rebuilding, according to usage factor { v ₁..., v _j, v _N3-D basis function { φ ₁..., φ _j, φ _NLinear combination represent density profile, that is,

$\mathrm{\&lambda;}(x,y,z)=\underset{j}{\mathrm{\&Sigma;}}{v}_j{\mathrm{\&phi;}}_j(x,y,z).$ Formula (16)

Though have many selections for basis function, our notice is the 3D voxel of rectangle herein.λ _jBe used to represent the coefficient of j basis function, that is, and the diffuse density in j voxel.Consider Fig. 4, three layers (or voxel) are shown, pass described stack (stack), pass on (delivering) the information Δ θ and the Δ x of observation with ray." hide " scattering in j voxel and displacement is expressed as Δ θ respectively _jWith Δ x _jOnce more, the amount of scattering is exaggerated in the drawings.We can connect the data of being observed with hiding by following formula

Δ θ=Δ θ ₁+ Δ θ ₂+ Δ θ ₃, formula (17)

Δx＝Δx ₁+L ₂tan(Δθ ₁)+Δx ₂+

L ₃tan(Δθ ₁+Δθ ₂)+Δx ₃

≈ Δ x ₁+ Δ x ₂+ Δ x ₃+ T ₁Δ θ ₁+ T ₂Δ θ ₂. formula (18)

Herein, we rely on the hypothesis of small angle scattering in second formula, and with T _jBe defined as from the exit point of j voxel to 3-D raypath length from the exit point of reconstruct volume.More generally say, for passing one group of voxel

Ray,

$\mathrm{\&Delta;\&theta;}=\underset{j\&Element;\&aleph;}{\mathrm{\&Sigma;}}{\mathrm{\&Delta;\&theta;}}_j,$ Formula (19)

$\mathrm{\&Delta;x}=\underset{j\&Element;\&aleph;}{\mathrm{\&Sigma;}}({\mathrm{\&Delta;x}}_j+T_j\mathrm{\&Delta;}{\mathrm{\&theta;}}_j).$ Formula (20)

At last, we can express the covariance to scattering/displacement i ray, that add up to, by at first noticing, for j voxel

${\mathrm{\&Sigma;}}_{\mathrm{ij}}={\mathrm{\&lambda;}}_j{A}_{\mathrm{ij}}{p}_{r,\mathrm{<figure-callout\; id="2"\; label="i"\; filenames="A2007800487880033H2.png"\; state="\{\{state\}\}">i</figure-callout>}}^2,$ Formula (21)

Wherein

$A_{\mathrm{ij}}\&equiv;\left[\begin{array}{cc}{L}_{\mathrm{ij}}& L_{\mathrm{ij}}^2/2\\ {\mathrm{<figure-callout\; id="2"\; label="L\; ij"\; filenames="A2007800487880033H2.png"\; state="\{\{state\}\}">L</figure-callout>}}_{\mathrm{ij}}^{}/2& L_{\mathrm{ij}}^3/3\end{array}\right]$ Formula (22)

And, L _IjBe i the path that ray passes j voxel, to not being defined as zero by the voxel of ray " bump ".Combinatorial formula (19) is to (22), and we can write out

${\mathrm{\&Sigma;}}_i=p_{r,\mathrm{<figure-callout\; id="2"\; label="i"\; filenames="A2007800487880033H2.png"\; state="\{\{state\}\}">i</figure-callout>}}^2\underset{j\&le;N}{\mathrm{\&Sigma;}}{\mathrm{\&lambda;}}_j{W}_{\mathrm{ij}}.$ Formula (23)

Herein, N is the voxel sum, and we are defined as weight matrix based on the simple but tediously long calculating to key element (elements)

$W_{\mathrm{ij}}\&equiv;\left[\begin{array}{cc}{L}_{\mathrm{ij}}& L_{\mathrm{ij}}^2/2+L_{\mathrm{ij}}{T}_{\mathrm{ij}}\\ L_{\mathrm{ij}}^2/2+L_{\mathrm{ij}}{T}_{\mathrm{ij}}& L_{\mathrm{ij}}^3/3+L_{\mathrm{ij}}^2{T}_{\mathrm{ij}}+L_{\mathrm{ij}}{\mathrm{<figure-callout\; id="2"\; label="T\; ij"\; filenames="A2007800487880033H2.png"\; state="\{\{state\}\}">T</figure-callout>}}_{\mathrm{ij}}^{}\end{array}\right],$ Formula (24)

Some hypothesis are carried out in unknown μ meson path, so that the raypath length that voxel is passed in estimation.With reference to figure 5, closest approach (point of closestapproach, the PoCA) (x of described approximate that enter with calculating and track outgoing _Ca, y _Ca) beginning.Then, the inlet (entry) of PoCA is connected to exit point, so that estimate the voxel path.

At last, definition of data vector

$D_i\&equiv;\left[\begin{array}{c}\mathrm{\&Delta;}{\mathrm{\&theta;}}_i\\ {\mathrm{\&Delta;x}}_i\end{array}\right]$ Formula (25)

And make D represent whole measurements from M μ meson.We are written as the data likelihood of given density profile λ

$P\left(D\right|\mathrm{\&lambda;})=\underset{i\&le;M}{\mathrm{\&Pi;}}P\left(D_i\right|\mathrm{\&lambda;})$ Formula (26)

Wherein utilize the factor

$P\left(D_i\right|\mathrm{\&lambda;})=\frac{1}{2\mathrm{\&pi;}{\left|{\mathrm{\&Sigma;}}_i\right|}^{1/2}}\exp (-\frac{1}{2}{D}_i^T{\mathrm{\&Sigma;}}_i^{-1}{D}_i).$ Formula (27)

Notice above-mentionedly, as summarizing among Figure 15 D,, can obtain the probability distribution that a plurality of μ mesons pass the scattering and the displacement of heterogeneous material according to an embodiment.At first, set up the 3D grid (perhaps other basis functions) (treatment step 170) of voxel.Then, in treatment step 171, calculate the covariance matrix of the scattering/displacement of each μ meson according to formula (23,24).At last, shown in treatment step 172, according to formula (25-27), given raypath and voxel diffuse density are calculated the overall probability distribution of whole μ mesons.

The statistical model of multiple scattering has been described, referring now to extension (treatment step 113) to the model of experiment effect.Real muon detection device presents limited position resolution.Come the characterization track with μ meson outgoing that enter by angle of deriving and position from the track that is fitted to a plurality of position measurements.Thereby measuring error propagates into the scattering angle and the displacement measurement of the data set (dataset) that constitutes μ meson chromatography.By the RMS error e _pThe precision of the given detector of characterization.For the specific arrangements of detector, can how to propagate based on error and define error matrix

$E\&equiv;\left[\begin{array}{cc}{e}_{\mathrm{\&Delta;\&theta;}}^2& e_{\mathrm{\&Delta;\&theta;\&Delta;x}}\\ e_{\mathrm{\&Delta;\&theta;\&Delta;x}}& e_{\mathrm{\&Delta;x}}^2\end{array}\right]$ Formula (28).

In the method for reconstructing of iteration, such error is relatively easy the processing.Under our situation, we can consider the detector error by the covariance matrix that replenishes formula (23)

${\mathrm{\&Sigma;}}_i=E+p_{r,\mathrm{<figure-callout\; id="2"\; label="i"\; filenames="A2007800487880033H2.png"\; state="\{\{state\}\}">i</figure-callout>}}^2\underset{j\&le;N}{\mathrm{\&Sigma;}}\left({\mathrm{\&lambda;}}_j{W}_{\mathrm{ij}}\right).$ Formula (29)

By this way, reduce noise, otherwise described noise will appear in the reconstruction owing to the detector error.Precise analytic model more to the detector error should be considered the momentum dependence, because a source of tracking error (tracking error) is the scattering in the detector self, and described scattering increases along with particle momentum and reduces.If the estimation of independent (individual) μ meson momentum

Be available, so can be to each ray evaluated error matrix

As obviously finding out from formula (2), the width of multiple Coulomb scattering depends on particle momentum.By in formula (5), introducing factor p _r ²Consider different μ meson momentum.In the practice, μ meson momentum is not accurately known, but can estimate the estimation of the momentum of independent μ meson from the scatterometry the known scatterer (as the known spectra of cosmic rays μ meson).Suppose that we have the good estimation to each μ meson herein,

Can use expectation-maximization algorithm to determine the PRML estimation (treatment step 103 of method 100) of object volume density.The EM algorithm relies on the likelihood of expressing " imperfect " data according to " complete " data (that is, observation data adds hiding data).In our application, observation data D={D _i: 1≤i≤M} is the scattering of measuring.Hiding data H={H _Ij: 1≤i≤M﹠amp; 1≤j≤N} is the scattering angle and the displacement of passing i μ meson of j voxel.At A.Dempster, " Maximum likelihood from incomplete data via the EM algorithm " (J.Roy.Statist.Soc.B of N.Laird and D.Rubin, vol.39, pp.1-78,1977) (by reference its disclosure is incorporated into here) in, is described algorithm according to following auxiliary function:

$Q_{\mathrm{DLR}}=E_{H|D,{\mathrm{\&lambda;}}^{\left(n\right)}}[log\left(P(D,H|\mathrm{\&lambda;})\right)].$ Formula (30)

Given and the parameter vector λ for H ⁽ⁿ⁾Condition distribute (conditional distribution), given parameter vector λ, this function be observation and the expectation values of the log-likelihood (loglikelihood) of the data of observation not.Each iteration of forming algorithm by two following steps.

E step: estimate or characterization P (H|D, λ ⁽ⁿ⁾), the conditional probability of hiding data.

The M step: maximization auxiliary function Q, this auxiliary function Q is the expectation value about the distribution of characterization in the E step.

Under our situation, because the definite uniquely data of being observed of hiding data, wherein by using simpler auxiliary function

$Q(\mathrm{\&lambda;};{\mathrm{\&lambda;}}^{\left(n\right)})=E_{H|D,{\mathrm{\&lambda;}}^{\left(n\right)}}\left[\log \left(P\left(H\right|\mathrm{\&lambda;})\right)\right]$ Formula (31)

We obtain can be by using Q _DLRThe identical sequence estimation that obtains

From parameter estimation λ ⁽ⁿ⁾In, by following formula, the iteration of described algorithm produces new estimation λ ⁽ⁿ⁺¹⁾

λ ⁽ⁿ⁺¹⁾=arg max _λ(Q (λ; λ ⁽ⁿ⁾). formula (32)

We notice that at first the probability distribution that single μ meson passes the scattering of single voxel draws from the statistical model of individual layer simply.

$P\left(H_{\mathrm{ij}}\right|\mathrm{\&lambda;})=\frac{1}{2\mathrm{\&pi;}{\left|{\mathrm{\&Sigma;}}_{\mathrm{ij}}\right|}^{1/2}}\exp (-\frac{1}{2}{H}_{\mathrm{ij}}^T{\mathrm{\&Sigma;}}_{\mathrm{ij}}^{-1}{H}_{\mathrm{ij}}),$ Formula (33)

Wherein, ${\mathrm{\&Sigma;}}_{\mathrm{ij}}={\mathrm{\&lambda;}}_j{A}_{\mathrm{ij}}{p}_{r,\mathrm{<figure-callout\; id="2"\; label="i"\; filenames="A2007800487880033H2.png"\; state="\{\{state\}\}">i</figure-callout>}}^2,$ Be defined in the formula (21).Because the unconditional distribution of the scattering in each voxel is independent of the scattering in other voxels, therefore the probability of the hiding data collection that adds up to is the product of the probability of each key element.Thereby, log-likelihood can be written as

$\log \left(P\left(H\right|\mathrm{\&lambda;})\right)=\underset{j\&le;N}{\mathrm{\&Sigma;}}\underset{i:L_{\mathrm{ij}}\&NotEqual;0}{\mathrm{\&Sigma;}}(-\log {\mathrm{\&lambda;}}_j-\frac{H_{\mathrm{ij}}^T{A_{\mathrm{ij}}}^{-1}{\mathrm{<figure-callout\; id="2"\; label="H\; ij"\; filenames="A2007800487880033H2.png"\; state="\{\{state\}\}">H</figure-callout>}}_{\mathrm{ij}}}{{2\mathrm{\&lambda;}}_j{p}_{r,i}^2})+C,$ Formula (34)

Wherein the C representative does not comprise the item of λ.Consider conditional expectation, we are written as the Q function

$Q(\mathrm{\&lambda;};{\mathrm{\&lambda;}}^{\left(n\right)})=C+\underset{j\&le;N}{\mathrm{\&Sigma;}}{Q}_j({\mathrm{\&lambda;}}_j;{\mathrm{\&lambda;}}_j^{\left(n\right)})$ Formula (35)

Wherein utilize sum term (summands)

$Q_j({\mathrm{\&lambda;}}_j;{\mathrm{\&lambda;}}_j^{\left(n\right)})=-M_j\log {\mathrm{\&lambda;}}_j-\frac{1}{2{\mathrm{\&lambda;}}_j}\underset{i:L_{\mathrm{ij}}\&NotEqual;0}{\mathrm{\&Sigma;}}{S}_{\mathrm{ij}}^{\left(n\right)}.$ Formula (36)

Here, M _jBe its L _IjThe number of ≠ 0 ray (that is, clashing into the number of the ray of j voxel), and S _Ij ⁽ⁿ⁾Be defined as

$S_{\mathrm{ij}}^{\left(n\right)}\&equiv;E_{H|D,{\mathrm{\&lambda;}}^{\left(n\right)}}\left[p_{r,i}^{-2}{H}_{\mathrm{ij}}^T{A}_{\mathrm{ij}}^{-1}{H}_{\mathrm{ij}}\right].$ Formula (37)

With formula (36) with respect to λ _jDerivative (derivative) be set to zero, we find the iterative formula of following maximization auxiliary function (M step)

${\mathrm{\&lambda;}}_j^{(n+1)}=\frac{2}{{2M}_j}\underset{i:L_{\mathrm{ij}}\&NotEqual;0}{\mathrm{\&Sigma;}}{S}_{\mathrm{ij}}^{\left(n\right)}.$ Formula (38)

S _IjThe quadratic form form guarantee λ ⁽ⁿ⁺¹⁾Positivity (positivity).It keeps with design conditions expectation S _IjMake X represent stochastic variable H _Ij| D _iQuadratic form X ^TA ^-1The expectation value of X is

$E\left[X^T{A}^{-1}X\right]=\mathrm{Tr}\left(A^{-1}{\mathrm{\&Sigma;}}_X\right)+{\mathrm{\&mu;}}_X^T{A}^{-1}{\mathrm{\&mu;}}_X,$ Formula (39)

μ wherein _XAnd ∑ _XBe respectively average and the covariance of X.Because D _iDepend on H linearly _Ij, so they are the associating Gaussian.Given D _iH _IjCondition to distribute also be Gaussian, from the result of multivariate (multivariate) distribution theory.Use described theory and H _IjWith D _iEach all has the fact of zero-mean, and we find

${\mathrm{\&mu;}}_X={\mathrm{\&Sigma;}}_{D_i{H}_{\mathrm{ij}}}^T{\mathrm{\&Sigma;}}_{D_i}^{-1}{D}_i,$ Formula (40)

${\mathrm{\&Sigma;}}_X={\mathrm{\&Sigma;}}_{H_{\mathrm{ij}}}-{\mathrm{\&Sigma;}}_{D_i{H}_{\mathrm{ij}}}^T{\mathrm{\&Sigma;}}_{D_i}^{-1}{\mathrm{\&Sigma;}}_{D_i{H}_{\mathrm{ij}}.}$ Formula (41)

Here, by the covariance of the given observation data of formula (29)

And can express the covariance of hiding data key element via formula (21) Matrix computations that we can carry out simply (though for tediously long) (but not is write out the hiding data of observation clearly to illustrate

Covariance)

${\mathrm{\&Sigma;}}_{D_i{H}_{\mathrm{ij}}}{A}_{\mathrm{ij}}^{-1}{\mathrm{\&Sigma;}}_{D_i{H}_{\mathrm{ij}}}^T=W_{\mathrm{ij}}{\left(p_{r,\mathrm{<figure-callout\; id="2"\; label="i"\; filenames="A2007800487880033H2.png"\; state="\{\{state\}\}">i</figure-callout>}}^2{\mathrm{\&lambda;}}_j\right)}^2.$ Formula (42)

To be substituted in the formula (37) from the result of formula (39) to (42), we find

$S_{\mathrm{ij}}^{\left(n\right)}=p_{r,i}^{-2}\mathrm{Tr}(A_{\mathrm{ij}}^{-1}{\mathrm{\&Sigma;}}_{H_{\mathrm{ij}}}-A_{\mathrm{ij}}^{-1}{\mathrm{\&Sigma;}}_{D_i{H}_{\mathrm{ij}}}^T{\mathrm{\&Sigma;}}_{D_i}^{-1}{\mathrm{\&Sigma;}}_{D_i{H}_{\mathrm{ij}}})+p_{r,i}^{-2}{D}_i^T{\mathrm{\&Sigma;}}_{D_i}^{-1}{W}_{\mathrm{ij}}{\mathrm{\&Sigma;}}_{D_i}^{-1}{D}_i{\left(p_{r,\mathrm{<figure-callout\; id="2"\; label="i"\; filenames="A2007800487880033H2.png"\; state="\{\{state\}\}">i</figure-callout>}}^2{\mathrm{\&lambda;}}_j^{\left(n\right)}\right)}^2$

$=2{\mathrm{\&lambda;}}_j^{\left(n\right)}+(D_i^T{\mathrm{\&Sigma;}}_{D_i}^{-1}{W}_{\mathrm{ij}}{\mathrm{\&Sigma;}}_{D_i}^{-1}{D}_i-\mathrm{Tr}\left({\mathrm{\&Sigma;}}_{D_i}^{-1}{W}_{\mathrm{ij}}\right))p_{r,i}^2{\left({\mathrm{\&lambda;}}_j^{\left(n\right)}\right)}^2,$ Formula (43)

Wherein, we have used Tr (AB)=Tr (BA) in a last step.

At last, in order to incorporate x and y coordinate scattering data into, we use average simply in the formula (38) that upgrades

$S_{\mathrm{ij}}^{\left(n\right)}=\frac{S_{\mathrm{ij},x}^{\left(n\right)}+S_{\mathrm{ij},y}^{\left(n\right)}}{2},$ Formula (44)

What notice is the average that the ray of voxel is being clashed in formula (38) representative.After this, will be referred to as Mean Method to the use of this formula.To illustrate below, this more replaceable form of new formula reducing owing to be useful in the noise that outlier (outlier) μ meson data cause.Define the median method of algorithm by the more new formula that is changed:

${\mathrm{\&lambda;}}_j^{(n+1)}=\frac{1}{2}\mathrm{media}{n}_{i:L_{\mathrm{ij}}\&NotEqual;0}\left(S_{\mathrm{ij}}^{\left(n\right)}\right),$ Formula (45)

Notice above-mentionedly, in Figure 15 E, summarized the process (treatment step 103 of method 100) of the maximization likelihood that uses estimation density profile expectation-maximization algorithm, object volume.Shown in treatment step 180, each μ meson i=1 is gathered the measurement of scattering and momentum to M: (Δ θ _x, Δ θ _y, Δ x, Δ _y, p _r ²) _iEstimate the interactional geometry of each μ meson and each voxel j=1 to N: (L, T) _Ij(treatment step 181).Right for each μ meson voxel, use formula (24) to calculate weight matrix: W _Ij(treatment step 182).With conjecture λ _{J, old}Come the diffuse density (treatment step 183) in each voxel of initialization.As following, indication stops criterion in treatment step 184.For each μ meson, use formula (29) and get inverse matrix (inverse taken) and calculate Right to each μ meson voxel, use formula (43,44) design conditions expectation item: S _Ij, use formula (38) or formula (45) to calculate λ _{J, old}, and to whole voxels setting λ _{J, old}=λ _{J, new}

Referring now to numerical example, so that further specify method 100.Simulation and the similar setting shown in Fig. 1.As initial validity test, use the simple analog that is designed to mate closely multiple statistics scattering model.There is 2 * 2 square metres size on single detector plane (rather than go out as shown in FIG. 3), and the vertical range between the detector of top and bottom is 1.1 meters.These detectors write down the position and the angle of μ meson ideally.Use the μ meson spectrum of simplifying, wherein the momentum of μ meson is uniformly distributed in 500 to 10000MeV/c.Particle is at the top detector plane, enter volume with the vertical projection angle of even leap (spanning).

Multiple scattering and displacement according to treatment step 110 to 113 simulation μ mesons.As manifesting in Fig. 7 and 8, object is placed in 1.1 * 1.1 * 1.1 cubic metres of parts in center of volume.Simulate three 10 * 10 * 10 cubic centimetres cubical material of tungsten (W), iron (Fe) and aluminium (Al), have diffuse density 71.5,14.2 and 2.8mrad respectively ²/ cm.Simulation hypothesis corresponding to 400000 μ mesons of about 10 minutes exposure in top detector stack incident.Do not run into (missed) bottom detector plane for about 160000 in the described μ meson, stay 240000 and be used for rebuilding.Suppose to know ideally the momentum of each μ meson, 5 * 5 * 5 cubic centimetres voxel size is used to rebuild, and realizes Mean Method described above.Simulation begins with the volume that is filled air, and with 100 iteration of algorithm operation (convergence of enough block features (block features)).The result appears among Fig. 9.Corresponding in three objects each, the reconstructed value of 8 voxels average, for (Al) piece is respectively (74.0,14.7,2.7) for W, Fe.The mark distribution (fractional spread) of forming each cubical 8 voxel is (12.6%, 13.2%, 12.1%) (rms/mean).Give the coupling that fixes between simulation and the inverse model (inversionmodel), this result makes inversion algorithm and implements effectively.

It is identical with the object scene that reconstruction shows as.Corresponding in three objects each, the reconstructed value of 8 voxels average, for (Al) piece is respectively (74.0,14.7,2.7) for W, Fe.The mark of forming each cubical 8 voxel scatters (rms/mean) and is (12.6%, 13.2%, 12.1%).Give the coupling that fixes between simulation and the inverse model, this result makes inversion algorithm and implements effectively.

Secondly, use GEANT4 Monte Carlo routine package to simulate identical scene again.Can be at J.Allison " Geant4 developments and applications " (IEEE Trans.Nucl.Sci., vol.53, no.1, pp.270-278, Feb.2006) find the details of GEANT4 in the publication, by reference its disclosure is incorporated into here.GEANT4 has realized more complete, the accurate and effective model to multiple scattering.This model comprises calculating, the realization of heavy-tailed (heavy tails) and the simulation of the energy loss of μ meson when passing material more accurately of the center Gauss width partly of scatter distributions.Also use μ meson event generator, it duplicates the sea level angle and the momentum spectrum of cosmic rays μ meson.Detector is desirable in the hypothetical simulation, knows the momentum of each μ meson ideally, does not comprise cosmic-ray electron or track secondary particle.The result appears among Figure 10.Corresponding to (W, Fe, Al) on average being respectively of the voxel value of piece (674.4,63.4,5.4).

Voxel value is too high, and the several mis-classification (missclassification) in middle and the low Z zone is obvious.By whole voxel values are approximately rebuild with the correct average voxel value normalization (normalizing) of Z voxel in the middle of producing divided by 4, do not produce correct value to high and low Z voxel, do not eliminate whole mis-classifications yet.To be sub-fraction μ meson fail the mode scattering of fine description with Gauss model to the reason of this effect.Require the projection angle center of distribution 98% of scattering to be approximately Gaussian well.All about 2% of the μ meson scatter to the angle very big, that is, occur much larger than Gaussian distribution with respect to statistical model described herein.Because square definite match of scattering angle, so effect can be significant.The μ meson scattering that drops on described afterbody produces too big diffuse density and estimates.

And, can be recorded as very large scattered through angles incident (though these sources do not occur in our simulation) mistakenly such as other processes or the significant detector error of the decay of the μ meson in the instrument of Fig. 1.This can be in volume generation Anywhere, and tend to produce single voxel with big unreasonably diffuse density.Such incident should be eliminated, because they provide false positive (false-positive) indication of SNM.

In order to make the EM algorithm tolerate the data of non-Gaussian, can replace average update rule formula (38) by formula (45), that is, use median method.

Figure 11 illustrates the result who uses median method.To (Al) Qu Yu voxel on average is respectively (79.2,14.2,2.1) for W, Fe, and the mark distribution is respectively (21.5%, 26.3%, 23.2%).Very clear, use the intermediate value update rule to improve the robustness of inversion algorithm.

Figure 12 illustrates the real example more of rebuilding density profile, it has illustrated the detailed GEANT4 simulation of passenger traffic railway carriage or compartment car.The explanation of removing the chief component of railway carriage or compartment car car body appears among Figure 12.Represent 10 * 10 * 10 cubic centimetres solid piece of tungsten at the piece of redness at the center of explanation, as the representative (proxy) of high Z threat objects.In this case, we use the detector plane four long limits, simulation be positioned at scene, so as to use more multilevel towards the μ meson.According to Mean Method and median method, simulate one minute cosmic rays μ meson exposure, and carry out the reconstruction of the data of using 5 * 5 * 5 cubic centimetres of big or small voxels.Figure 13 and 14 illustrates visual (visualization) that uses the reconstruction that average EM method and intermediate value EM method carry out respectively.The effect of non-Gaussian data is significantly in the Mean Method of this scene is rebuild, and shows as the darker point that is scattered on the image.In the reconstruction of median method, these illusions all disappear, and the fine and close more ingredient (engine, battery, kinematic train) of railway carriage or compartment car is depicted as (low Z) or (middle Z), and threat objects darker highlightedly (high Z).The use of median method has produced tends to resist the result that false positive that scatter distributions and other anomalous events by non-Gauss cause has robustness.

Here embodiment of Chan Shuing and example are used for explaining the present invention and practical application thereof best, and make those skilled in the art can make and utilize the present invention thus.Yet those skilled in the art will be appreciated that aforementioned description and example only are for the purpose with example is described.

Other changes of the present invention and modification (as the adjustment of process of reconstruction) will be clearly to those skilled in the art, and contain such change and modification is the intention of appended claims.

The description of being set forth do not plan to be exhaustive (exhaustive), or in order to limit the scope of the invention.According to top instruction, many modifications and change all are possible, and do not deviate from the scope of following claims.Can be contemplated that use of the present invention can relate to the ingredient with different characteristics.

Claims (25)

1. As given a definition and wherein asked for protection the inventive embodiment of exclusive property or right.Having invention such description, that ask for protection is:

   1. method of from the charged particle chromatographic data, surveying object volume, described charged particle chromatographic data obtains from described object volume, and described method comprises:

   (a) obtain and pass object volume charged particle scattering angle and estimate the corresponding predetermined charged particle chromatographic data of momentum;

   (b) be provided for the probability distribution of the charged particle scattering of expectation maximization (ML/EM) algorithm, described probability distribution is based on statistics multiple scattering model;

   (c) use described expectation maximization (ML/EM) algorithm to determine the essence PRML estimation of object volume diffuse density; And

   (d) estimate the object volume diffuse density that output is rebuild based on described essence PRML.
2. 2. the method for claim 1 comprises:

   Based on the object volume diffuse density of described reconstruction, (1) existence of the target object in the described object volume and at least one in (2) type are made decision.
3. 3. the probability distribution that the method for claim 1, wherein is provided for the charged particle scattering of expectation maximization (ML/EM) algorithm comprises:

   (g) obtain the 2D probability distribution of charged particle based on the predefined diffuse density of homogeneous body;

   (h) obtain the 3D probability distribution of described charged particle based on described 2D probability distribution;

   (i) obtain to pass probability distribution via the scattering of a plurality of charged particles of the heterogeneous body object volume of basis function characterization; And

   (j), obtain the probability distribution of multiple scattering based on the definition of multiple scattering and the scattering and the momentum survey of described charged particle.
4. 4. method as claimed in claim 3, wherein, the 2D probability distribution that obtains charged particle based on the predefined diffuse density of homogeneous body comprises:

   (k) determine of the expected mean square scattering of the diffuse density of material as the described charged particle of the unit depth of passing described material;

   (l) be similar to the charged particle scatter distributions based on Gauss model; And

   (m) based on the 2D that is associated

   Gaussian distribution is similar to charged particle ray scattering and displacement.
5. 5. method as claimed in claim 3, wherein, the 3D probability distribution that obtains described charged particle based on described 2D probability distribution comprises:

   Add coordinate and define three-dimensional path length;

   Calculate the 3D displacement; And

   Definition 3D covariance matrix.
6. 6. method as claimed in claim 5, wherein, the probability distribution that acquisition is passed via the scattering of a plurality of charged particles of the heterogeneous body object volume of basis function characterization comprises:

   Set up the 3D grid of the basis function of the discrete diffuse density estimation of representative;

   Scattering/displacement the covariance matrix of the μ meson that each is independent is defined as the function of raypath and diffuse density estimation; And

   The probability distribution of a plurality of charged particles is defined as the product of independent charged particle probability.
7. 7. the method for claim 1, wherein use described expectation maximization (ML/EM) algorithm to determine that the essence PRML estimation of described object volume diffuse density comprises:

   Collection is to the scattering of each charged particle and the measurement of momentum;

   Estimate the interactional geometry of each basis function of each charged particle and described statistics scattering model;

   Right to each charged particle basis function, determine weight matrix: W _Ij

   Utilize the diffuse density in each basis function of conjecture initialization; And

   Finding the solution the approximate PRML of object volume capacity iteratively separates;

   Wherein, separate when becoming, stop described iterative process less than predetermined tolerance value at the iteration place of predetermined number or when described.
8. 8. the method for claim 1, wherein use described expectation maximization (ML/EM) algorithm to determine that the essence PRML estimation of described object volume diffuse density comprises:

   Collection to each charged particle i=1 to the scattering of M and the measurement of momentum: (Δ θ _x, Δ θ _y, Δ x, Δ _y, p _r ²) _i

   Estimate each μ meson and each voxel j=1 interactional geometry to N: (L, T) _Ij

   Right to each charged particle voxel, with described weight matrix W _IjBe calculated as

   $W_{\mathrm{ij}}\&equiv;\left[\begin{array}{cc}{L}_{\mathrm{ij}}& L_{\mathrm{ij}}^2/2+L_{\mathrm{ij}}{T}_{\mathrm{ij}}\\ L_{\mathrm{ij}}^2/2+L_{\mathrm{ij}}{T}_{\mathrm{ij}}& L_{\mathrm{ij}}^3/3+L_{\mathrm{ij}}^2{T}_{\mathrm{ij}}+L_{\mathrm{ij}}{T}_{\mathrm{ij}}^2\end{array}\right],$

   The conjecture λ of initialization diffuse density in each voxel _{J, old}And

   Whole voxels uses are stopped the criterion process λ is set _{J, old}=λ _{J, new}
9. 9. the method for claim 1, wherein described expectation maximization (ML/EM) algorithm comprises average update rule or intermediate value update rule.
10. 10. the method for claim 1, wherein described charged particle chromatographic data comprises cosmic rays μ meson chromatographic data.
11. 11. survey method object volume, computer-implemented for one kind from the charged particle chromatographic data, described charged particle chromatographic data obtains from described object volume, described method comprises:

    (a) obtain and pass object volume charged particle scattering angle and estimate the corresponding charged particle chromatographic data of momentum;

    (b) be provided for the probability distribution of the charged particle scattering of expectation maximization (ML/EM) algorithm, described probability distribution is based on statistics multiple scattering model;

    (c) use described expectation maximization (ML/EM) algorithm to determine the essence PRML estimation of object volume diffuse density; And

    (d) the object volume diffuse density of output reconstruction.
12. 12. method as claimed in claim 11 further comprises based on the object scattering volume density of described reconstruction and makes decision.
13. 13. method as claimed in claim 11, wherein, the probability distribution that is provided for the charged particle scattering of expectation maximization (ML/EM) algorithm comprises:

    (g) obtain the 2D probability distribution of charged particle based on the predefined diffuse density of homogeneous body;

    (h) obtain the 3D probability distribution of described charged particle based on described 2D probability distribution;

    (i) obtain to pass probability distribution via the scattering of a plurality of charged particles of the heterogeneous body object volume of basis function characterization; And

    (j), obtain the probability distribution of multiple scattering based on the definition of multiple scattering and the scattering and the momentum survey of described charged particle.
14. 14. method as claimed in claim 13, wherein, described charged particle chromatographic data comprises cosmic rays charged particle chromatographic data.
15. 15. method as claimed in claim 14, wherein, the 2D probability distribution that obtains charged particle based on the predefined diffuse density of homogeneous body comprises:

    The diffuse density of determining material has a nominal momentum p as the unit depth of passing this material ₀The expected mean square scattering of the cosmic rays charged particle of=3GeV;

    Use Gauss model to be similar to cosmic rays charged particle scatter distributions; And

    The 2D Gaussian distribution is similar to cosmic rays scattering and displacement by being associated;

    Wherein, the 3D probability distribution that obtains described charged particle based on described 2D probability distribution comprises:

    Add coordinate and define three-dimensional path length;

    Calculate the 3D displacement; And

    Definition 3D covariance matrix, and

    Wherein, the probability distribution that obtains to pass via the scattering of a plurality of charged particles of the heterogeneous body object volume of basis function characterization comprises:

    Set up the 3D grid of the basis function of the discrete diffuse density estimation of representative;

    Scattering/displacement the covariance matrix of the cosmic rays charged particle that each is independent is defined as the function of raypath and diffuse density estimation; And

    The probability distribution of a plurality of cosmic rays charged particles is defined as the product of the probability of independent charged particle.
16. 16. method as claimed in claim 14 wherein, uses described expectation maximization (ML/EM) algorithm to determine that the essence PRML estimation of object volume diffuse density comprises:

    Gather the scattering of each cosmic rays charged particle and the measurement of momentum;

    Estimate the interactional geometry of each basis function of each charged particle and described statistics multiple scattering model;

    Right to each charged particle basis function, determine weight matrix: W _Ij

    Utilize the diffuse density in each basis function of conjecture initialization; And

    Finding the solution the approximate PRML of object volume capacity iteratively separates;

    Wherein, separate when becoming, stop described iterative process less than predetermined tolerance value at the iteration place of predetermined number or when described.
17. 17. method as claimed in claim 10 wherein, uses described expectation maximization (ML/EM) algorithm to determine that the essence PRML estimation of object volume diffuse density comprises:

    Gather each charged particle i=1 to the scattering of M and the measurement of momentum: (Δ θ _x, Δ θ _y, Δ x, Δ _y, p _r ²) _i

    Estimate each charged particle and each voxel j=1 interactional geometry to N: (L, T) _Ij

    Right to each charged particle voxel, with described weight matrix W _IjBe calculated as

    $W_{\mathrm{ij}}\&equiv;\left[\begin{array}{cc}{L}_{\mathrm{ij}}& L_{\mathrm{ij}}^2/2+L_{\mathrm{ij}}{T}_{\mathrm{ij}}\\ L_{\mathrm{ij}}^2/2+L_{\mathrm{ij}}{T}_{\mathrm{ij}}& L_{\mathrm{ij}}^3/3+L_{\mathrm{ij}}^2{T}_{\mathrm{ij}}+L_{\mathrm{ij}}{T}_{\mathrm{ij}}^2\end{array}\right],$

    The conjecture λ of initialization diffuse density in each voxel _{J, old}And

    Whole voxels uses are stopped the criterion process λ is set _{J, old}=λ _{J, new}
18. 18. method as claimed in claim 17, wherein, the described criterion process that stops to comprise: to each charged particle, use $W_{\mathrm{ij}}\&equiv;\left[\begin{array}{cc}{L}_{\mathrm{ij}}& L_{\mathrm{ij}}^2/2+L_{\mathrm{ij}}{T}_{\mathrm{ij}}\\ L_{\mathrm{ij}}^2/2+L_{\mathrm{ij}}{T}_{\mathrm{ij}}& L_{\mathrm{ij}}^3/3+L_{\mathrm{ij}}^2{T}_{\mathrm{ij}}+L_{\mathrm{ij}}{T}_{\mathrm{ij}}^2\end{array}\right],$ Calculate

    

    , and get inverse matrix,

    Right to each charged particle voxel, use

    $S_{\mathrm{ij}}^{\left(n\right)}=p_{r,i}^{-2}\mathrm{Tr}(A_{\mathrm{ij}}^{-1}{\mathrm{\&Sigma;}}_{H_{\mathrm{ij}}}-A_{\mathrm{ij}}^{-1}{\mathrm{\&Sigma;}}_{D_i{H}_{\mathrm{ij}}}^T{\mathrm{\&Sigma;}}_{D_i}^{-1}{\mathrm{\&Sigma;}}_{D_i{H}_{\mathrm{ij}}})+p_{r,i}^{-2}{D}_i^T{\mathrm{\&Sigma;}}_{D_i}^{-1}{W}_{\mathrm{ij}}{\mathrm{\&Sigma;}}_{D_i}^{-1}{D}_i{\left(p_{r,i}^2{\mathrm{\&lambda;}}_j^{\left(n\right)}\right)}^2$

    $=2{\mathrm{\&lambda;}}_j^{\left(n\right)}+(D_i^T{\mathrm{\&Sigma;}}_{D_i}^{-1}{W}_{\mathrm{ij}}{\mathrm{\&Sigma;}}_{D_i}^{-1}{D}_i-\mathrm{Tr}\left({\mathrm{\&Sigma;}}_{D_i}^{-1}{W}_{\mathrm{ij}}\right))p_{r,i}^2{\left({\mathrm{\&lambda;}}_j^{\left(n\right)}\right)}^2,$

    Come design conditions expectation item: S _Ij,

    Wherein, Tr (AB)=Tr (BA) in the step in the end;

    And, in update rule, use

    $S_{\mathrm{ij}}^{\left(n\right)}=\frac{S_{\mathrm{ij},x}^{\left(n\right)}+S_{\mathrm{ij},y}^{\left(n\right)}}{2},$

    Merge x and y coordinate scattering data.
19. 19. method as claimed in claim 18, wherein, described charged particle comprises the μ meson.
20. 20. method as claimed in claim 18, wherein, described ML/EM algorithm comprises and being defined as ${\mathrm{\&lambda;}}_j^{(n+1)}=\frac{1}{{2M}_j}\underset{i:L_{\mathrm{ij}}\&NotEqual;0}{\mathrm{\&Sigma;}}{S}_{\mathrm{ij}}^{\left(n\right)}.$ The average update rule, or be defined as ${\text{\&lambda;}}_j^{(n+1)}=\frac{1}{2}{\mathrm{median}}_{i:L_{\mathrm{ij}}\&NotEqual;0}\left(S_{\mathrm{ij}}^{\left(n\right)}\right),$ The intermediate value update rule.
21. 21. computer program, comprise: the data carrier that the computing machine of storage instruction can be used, when being carried out by computing machine, the method that the statistics that described instruction makes described computing machine carry out the object volume density profile that is used for the charged particle chromatographic data is rebuild, described method comprises:

    (a) obtain and pass object volume charged particle scattering angle and estimate the corresponding predetermined charged particle chromatographic data of momentum;

    (b) be provided for the probability distribution of the charged particle scattering of expectation maximization (ML/EM) algorithm, described probability distribution is based on statistics multiple scattering model;

    (c) use described expectation maximization (ML/EM) algorithm to determine the essence PRML estimation of object volume diffuse density; And

    (d) the object volume diffuse density of output reconstruction.
22. 22. a detection system that is used for via the charged particle detection object volume that passes object volume comprises:

    First group of Position-Sensitive Detector, it is positioned at first side of object volume, to measure position and the angle to the charged particle of described object volume incident;

    Second group of Position-Sensitive Detector, it is positioned at second side relative with described first side, described object volume, leaves the position and the angle of the charged particle of described object volume outgoing with measurement; With

    Signal processing unit, it receives from the measuring-signal of described first group of Position-Sensitive Detector with from the data of the measuring-signal of described second group of Position-Sensitive Detector, wherein, the data that described signal processing unit processes received are so that produce the statistics reconstruction that described object volume inner volume diffuse density distributes.
23. 23. the system as claimed in claim 22, wherein, described signal processing unit is configured to:

    (a) obtain and pass object volume charged particle scattering angle and estimate the corresponding charged particle chromatographic data of momentum;

    (b) provide the probability distribution of charged particle diffuse density based on adding up the multiple scattering model;

    (c) use described expectation maximization (ML/EM) algorithm to determine the essence PRML estimation of object volume diffuse density; And

    (d) estimate the object volume diffuse density that output is rebuild based on described essence PRML.
24. 24. the system as claimed in claim 22, wherein:

    In described first group and the second group of particle detector each comprises drift tube, and this drift tube is arranged to allow three charged particle position measurements and three the charged particle position measurements in being different from the second direction of first direction in first direction at least at least.
25. 25. the system as claimed in claim 22, wherein, described charged particle is the natural cosmic rays μ meson that incides described object volume, and whether described signal processing unit is configured to rebuild the indicating target object to be present in the described object volume based on the statistics that described object volume inner volume diffuse density distributes.

Images