2004 3-D geometry enhancement by contour optimization in turntable sequences

合集下载

大地电磁全信息资料三维共轭梯度反演研究_英文_

大地电磁全信息资料三维共轭梯度反演研究_英文_

1Manuscript received by the Editor September 31, 2010; revised manuscript received January 11, 2011.*This work is jointly supported by the National Hi-tech Research and Development Program of China (863 Program) (No. 2007AA09Z310), National Natural Science Foundation of China (Grant No.40774029, 40374024), the Fundamental Research Funds for the Central Universities (Grant No. 2010ZY53), and the Program for New Century Excellent Talents in University (NCET).1. State Key Laboratory of Geological Processes and Mineral Resources, China University of Geosciences, Beijing, 100083, China.2. Key Laboratory of Geo-detection (China University of Geosciences), Ministry of Education, Beijing, 100083, China.3. School of Geophysics and Information Technology, China University of Geosciences, Beijing, 100083, China.Three-dimensional conjugate gradient inversion ofmagnetotelluric full information data*APPLIED GEOPHYSICS, Vol.8, No.1 (March 2011), P. 1 - 10, 3 Figures.DOI: 10.1007/s11770-010-0266-9Lin Chang-Hong 1,2,3, Tan Han-Dong 1,2,3, and Tong Tuo 1,2,3Abstract : Based on the analysis of impedance tensor data, tipper data, and the conjugate gradient algorithm, we develop a three-dimensional (3D) conjugate gradient algorithm for inverting magnetotelluric full information data determined from five electric and magnetic field components and discuss the method to use the full information data for quantitative interpretation of 3D inversion results. Results from the 3D inversion of synthetic data indicate that the results from inverting full information data which combine the impedance tensor and tipper data are better than results from inverting only the impedance tensor data (or tipper data) in improving resolution and reliability. The synthetic examples also demonstrate the validity and stability of this 3D inversion algorithm.Keywords : Magnetotelluric, full information data, 3D inversion, conjugate gradientIntroductionThe magnetotelluric (MT) method can provide more than 20 parameters for geological interpretation after the observation data are processed, such as apparent resistivity, phase, impedance, tipper, skew, and etc. Apparent resistivity, phase, and impedance are often applied for quantitative interpretation, while others are used for qualitative interpretation. Previously, the most frequently used parameters are the apparent resistivity and phase, especially the 2D apparent resistivity and phase. Some 2D inversion algorithms have been developed to invert 2D apparent resistivity and phase data for the E- and H-polarization modes (Jupp and Vozoff, 1977; DeGroot-Hedlin and Constable, 1990; Smithand Booker, 1991; Siripunvaraporn and Egbert, 2000; Rodi and Mackie, 2001). In recent years, 3D inversion has been a hot topic in electromagnetics. The apparent resistivity and phase data can be used for 3D inversion (Mackie and Madden, 1993; Spichak and Popova, 2000; Zhdanov et al., 2000; Tan et al., 2003b; Hu et al., 2006; Lin et al, 2009) and scholars also developed algorithms which use the impedance data for 3D inversion (Newman and Alumbaugh, 2000; Siripunvaraporn et al., 2005; Avdeev and Avdeeva, 2006; Lin et al., 2008).The magnetotelluric method usually measures five electric and magnetic field components: two horizontal electric components, two horizontal magnetic components, and one vertical magnetic component. Both the impedance tensor data and the apparent resistivity-phase data which can be calculated from impedance2Three-dimensional conjugate gradient inversionare determined from the relationship between the two electric and magnetic fields horizontal components,x x y y E H Z E H ªºªº «»«»¬¼¬¼, 21sZ U ZP , and Arg Z I . Thus, inverting only the impedance data or apparent resistivity-phase data means that only four electric and magneticfield components are used for interpretation. If we want to use the full information data (five electric and magnetic field components), another parameter which is related to the horizontal and vertical magnetic fields, tipper, should be considered in the inversion. Inmagnetotelluric theory, tipper is defined as x z y H H T H ªº«»¬¼.The tipper data were seldom obtained because past magnetotelluric field work usually did not acquire the vertical magnetic fields. It was considered that the tipper data signal-to-noise was low and it was difficult to apply the tipper data for quantitative interpretation. Therefore, the tipper data were only employed for recognition of geological structure strike. However, modern magnetotelluric field practice typically includes measurement of the vertical magnetic fields, particularly at long periods where a tri-axial magnetometer is used. Berdichevsky et al. (2003) studied the tipper data generated from 2D geological models and concluded that the tipper data is affected by magnetotelluric static shift smaller than the impedance tensor data (or the apparent resistivity-phase data). If the tipper data is included in the inversion, the reliability of geoelectrical models can be substantially improved. Some 2D inversion algorithms (Siripunvaraporn and Egbert, 2000; Rodi and Mackie, 2001) allowed inversion of the tipper data, which became included along with the TM and TE apparent resistivity-phase data in 2D profile interpretation. From the 2D inversion of field MT, Tuncer et al. (2006) concluded that the tipper data add useful information to the inversion and enhance its resolution. From the 2D inversion example of the synthetic data, Becken et al. (2008) found that the tipper data are most important for sensing deep anomalies.Based on the analysis of impedance tensor data, tipper data, and the conjugate gradient algorithm, we develop a 3D conjugate gradient inversion algorithm for inverting magnetotelluric full information data determined from five electric and magnetic field components and discuss the importance of full information data for 3D inversion. In this paper, we introduce a 3D conjugate gradient inversion algorithm for magnetotelluric full information data. Forward modeling of the impedance tensor andtipper response will be discussed first. Second, we will discuss the theory of the inversion, including the objective function and the pseudo-forward modeling problem. Finally, the results of 3D inversion with synthetic data will be shown.Forward problemBecause impedance and tipper data are used in 3D inversion, the 3D forward modeling problem needs to consider how to correctly simulate the impedancetensor response xx xy yx yy Z Z Z Z Z ªº«»¬¼ and the tipper response zx zy T T T ªº ¬¼ on a surface generated from complex realistic 3D geology.The magnetotelluric source S can always be divided into two orthogonal sources, SX and SY (Tan et al., 2004). Under the influence of SX, the electric field components on the surface can be denoted by E x 1 and E y 1 and the magnetic field components by H x 1, H y 1, and H z 1. The corresponding values under the influence of SY can be denoted by E x 2, E y 2, H x 2, H y 2, and H z 2. Then, the expressions used to compute the impedance tensor response are:12211221211212211221122121121221,,,.x y x y xx x y x y x x x x xyx y x y y y y y yx x y x y y x y x yy x y x y E H E H Z H H H H E H E H Z H H H H E H E H Z H H H H E H E H Z H H H H ­ ° °° ° °®° ° °°° ¯(1)The expressions used to compute the tipper responseare:2112122112211221,.y z y z zx x y x y x z x Z zy x y x y H H H H T H H H H H H H H T ­ ° °®° ° ¯(2)From equations (1) and (2), we know that the electricand magnetic fields on the surface are needed first if we want to obtain the impedance tensor or tipper responses. For a given polarization mode, the magnetic field H j at3Lin et al.the jth site can be expressed by a vector H (Newman and Alumbaugh, 2000):.h T j j H g H(3)where H j denotes the magnetic field at the center of jth site and H is a vector composed of the three magnetic field components at each grid subsurface cell, which can be transformed into each other using the interpolator vector, h g T j . For three directions x , y , and z , equation (3) can be correspondingly written as:,,.h Tjx jx h Tjy jy h Tjz jz H g H H g H H g H ­ °° ®° °¯ (4)Similarly, the electric field E j for a given polarization mode can also be expressed by a vector H :.e T j j E g H(5)where E j denotes the electric field at the center of jthsite and the vector H denotes the magnetic field onthe model boundaries calculated by the staggered-gridfinite difference method, which can be transformed into each other by using the interpolator vector, e g T j . For components at the two direction x and y of the electric field E j , equation (5) can be written as:,.e Tjx jx e T jy jy E g H E g H ­ °® °¯(6)From equations (4) and (6), we know that the vector H is needed first if we want to obtain the electric and magnetic field on the surface. Vector H can be calculated by the staggered-grid finite difference method (Newman and Alumbaugh, 2000; Tan et al., 2003a). With a staggered grid, a discretization of the integral form of the Maxwell equations: 0H l J S E S,E l H S.d d d d i d V P Z ­ °® °¯³³³³³³³³v v(7) (7)have the forward modeling equation of the magnetic field H at the sampling point of each grid cell as: ,KH s (8)where K is a sparse complex-symmetric matrix and sis the source vector related to the boundary conditions. Two forward modeling equations can be obtained, KH 1 = s 1 and KH 2 = s 2, if the magnetic field values to be solvedunder the action of the source SX and vectors at the right end of column of the equations are denoted by H 1 and s 1 and the magnetic field values to be solved under the action of the source SY and vectors at the right end of the equations are denoted by H 2 and s 2. These two forward modeling equations are solved to give magnetic fields H 1 and H 2 of the two polarization modes at sampling locations of each grid cell.Vector H is obtained by solving equation (8). Substituting the vector H into equations (3) and (5), the electric and magnetic fields on the surface are calculated. Then, the impedance tensor and tipper response of the 3D model is computed by equations (1) and (2).Inversion problemThe objective functionBased on the regularized solution method, the objective function of the 3D conjugate gradient inversionalgorithm for inverting magnetotelluric full informationdata can be defined as:d m \100obs T obs TT \\ 100(())(())()().obs T obs T T D F m V D F m m m L L m m O \\ 100(())(())()(obs T obs T T D F m V D F m m m L L m O (9)where 12,,...,T obs obs obs obsN D D D D ªº ¬¼ is a vector composed of the observed impedance or tipper data with D obsidenoting the real or imaginary part of impedance (or tipper) for a particular observation site and frequency, N is the total number of observations, F (m ) denotes the forward modeling function to calculate the impedance or tipper responses, V is the variance of error vector, λ is the regularization parameter, and L is the simple second-difference Laplacian. We have m = [m 1, m 2, …,m M ] where M is the total number of model blocks andm j is the logarithm of resistivity for a unique block (log ρ), and m 0 is the priori model. The objective functionhas two terms: the first term is data fit ψd and the second term is the model fit ψm being constraints for smoothingthe model. The regularization parameter λ balances the influence of ψd and ψm on the image solution. When boththe impedance data and the tipper data are used in theinversion, data fit is expressed as: .d z t \\Z\ (10)where ψz is the misfit of the impedance data, ψt is the4Three-dimensional conjugate gradient inversion misfit of the tipper data, and the weight factor ω is used to balance the ratio of ψz and ψt . Generally, the tipper data are small and easily influenced by noise in field work and it is difficult to acquire good-quality tipper data. Therefore, the value of ω can be decided by the data quality. The tipper data weights can increase (or decrease) according to the quality of the tipper data.The gradient of the objective function can be expressed as: 1022().T Tg A V e L L m m O (11)where A is the Jacobian matrix and the error vector ise = D obs-F (m ). Calculation of the Jacobian matrix – the pseudo-forward modeling problemThe most important part in the inversion process is to determine the Jacobian matrix. Once the value of theJacobian matrix is calculated, the model update at each iteration can be easily computed. The 3D conjugate gradient inversion algorithm flowchart for inverting magnetotelluric full information data is similar withthat for impedance data. The main computations in the complete inversion process are calculations of the valueof A (m ref )T V -1e and A (m ref )p . We will describe how to compute the values of A (m ref )T V -1e and A (m ref )p when the tipper data are inverted in this algorithm, which is quitesimilar with the process when the impedance data are inverted (Lin et al., 2008).Differentiating with respect to m k , equation (8) yields: 1/(/).k k H m K K m H w w w w (12)Substituting equation (4) into equation (12) and with ∂h g Tj /∂m = 0, we get:111/(/),/(/),/(/).h Tjx k jxk h T jy k jy k h T jz k jz k H m g K K m H H m g K K m HH m g K K m H ­w w w w °°w w w w ®°w w w w °¯ (13) Differentiating with respect to m k , equation (2) yields: k y z k z y y x y x k zx m H H m H H H H H H m T w w w w w w //)((/21121221)()//21121221z y z y k y z k z y H H H H m H H m H H w w k x y k y x k x y m H H m H H m H H w w w w w w u ///(2121122211221/)/(),x y k x y x y H H m H H H H º w w ¼>k x z k z x y x y x k xy m H H m H H H H H H m T w w //)((/12211221)()//12212112z x z x k x z k z x H H H H m H H m H H w w w w k x y k y x k x y m H H m H H m H H w w w w w w u ///(2121122211221/)/().x y k x y x y H H m H H H H º w w ¼ (14) >k y z k z y y x y x k zx m H H m H H H H H H m T w w w w w w //)((/21121221)()//21121221z y z y k y z k z y H H H H m H H m H H w w k x y k y x k x y m H H m H H m H H w w w w w w u ///(2121122211221/)/(),x y k x y x y m H H H H º w ¼>k x z k z x y x y x k xy m H H m H H H H H H m T w w //)((/12211221)()//12212112z x z x k x z k z x H H H H m H H m H H w w w w k x y k y x k x y m H H m H H m H H w w w w w w u ///(2121122211221/)/().x y k x y x y H H m H H H H º w w ¼ (14) >k y z k z y y x y x k zx m H H m H H H H H H m T w w w w w w //)((/21121221)()//21121221z y z y k y z k z y H H H H m H H m H H w w w w k x y k y x k x y m H H m H H m H H w w w w w w u ///(2121122211221/)/(),x y k x y x y H H m H H H H º w w ¼>k x z k z x y x y x k xy m H H m H H H H H H m T w w //)((/12211221)()//12212112z x z x k x z k z x H H H H m H H m H H w k x y k y x k x y m H H m H H m H H w w u ///(2121122211221/)/().x y k x y x y H H m H H H H º w w ¼ (14) k y z k z y m H H m H H w w //)(21121)()/211212z y z y k y z H H H H m H H k x y k y x H H m H H w /212121221/(),x y x y H H H H k x z k z x m H H m H H //)(12211)()/122121z x z x k x z H H H H m H H k x y k y x m H H m H H w w //2121k y z k z y y m H H m H H w w w //)(21121)()/211212z y z y k y z H H H H m H H k x y k y x m H H m H H w w w //212121221/(),x y x y H H H H º ¼k x z k z x y m H H m H H //)(12211)()122121z x z x k x z H H H H H k x y k y x m H H m H H w w w //2121ky z k z y y x y x k zx m H H m H H H H H H m T w //)((/21121221)()//21121221z y z y k y z k z y H H H H m H H m H H w w w k x y k y x k x y m H H m H H m H H w w w w w u ///(2121122211221/)/(),x y k x y x y H H m H H H º w w ¼>k x z k z x y x y x k xy m H H m H H H H H H m T w //)((/12211221)()//12212112z x z x k x z k z x H H H H m H H m H H w w w w k x y k y x k x y m H H m H H m H H w w w w w w u ///(2121122211221/)/().x y k x y x y H H m H H H º w w ¼ (14) k y z k z y y x y x m H H m H H H H H //)((21121221)()//21121221z y z y k y z k z y H H H m H H m H H w w w w k x y k y x k x y m H H m H H m H H w w w w w w ///(2121122211221/)/(),x y k x y x yH H m H H H H º w w ¼k x z k z x y x y x m H H m H H H H H H //)((12211221)()//12212112z x z x k x z k z x H H H H m H H m H w w w w k x y k y x k x y m H H m H H m H H w w w w w w ///(212112k y z k z y y x y x k zx m H H m H H H H H H m T w w //)((/21121221)()//21121221z y z y k y z k z y H H H H m H H m H H w w k x y k y x k x y m H H m H H m H H w w w w w w u ///(2121122211221/)/(),x y k x y x y H H m H H H H º w w ¼>k x z k z x y x y x k xy m H H m H H H H H H m T w w //)((/12211221)()//12212112z x z x k x z k z x H H H H m H m H H w w w w k x y k y x k x y m H m H H m H H w w w w w w u ///(212112 k y z k z y y x y x k zx m H H m H H H H H H m T w w //)((/21121221)()//21121221z y z y k y z k z y H H H H m H H m H H w w k x y k y x k x y m H H m H H m H H w w u ///(2121122211221/)/(),x y k x y x y H H m H H H H º w w ¼>k x z k z x y x y x k xy m H H m H H H H H H m T w w w w w w //)((/12211221)()//12212112z x z x k x z k z x H H H H m H H m H H w w w w k x y k y x k x y m H H m H H m H H w w w w w w u ///(2121122211221/)/().x y k x y x y H H m H H H H º w w ¼ (14) k y z k z y y x y x k zx m H H m H H H H H H m T w w //)((/21121221)()//21121221z y z y k y z k z y H H H H m H H m H H w w k x y k y x k x y m H H m H H m H H w w u ///(2121122211221/)/(),x y k x y x y H H m H H H H º w w ¼k x z k z x y x y x k xy m H H m H H H H H H m T w w w w //)((/12211221)()//12212112z x z x k x z k z x H H H H m H H m H H w w k x y k y x k x y m H H m H H m H H w w u ///(2121122211221/)/().x y k x y x y H H m H H H H º w w ¼ (14) k y z k z y y x y x k zx m H H m H H H H H H m T w w //)((/21121221)()//21121221z y z y k y z k z y H H H H m H H m H H w w k x y k y x k x y m H H m H H m H H w w u ///(2121122211221/)/(),x y k x y x y H H m H H H H º w w ¼k x z k z x y x y x k xy m H H m H H H H H H m T w w w w w //)((/12211221)()//12212112z x z x k x z k z x H H H H m H H m H H w w k x y k y x k x y m H H m H H m H H w w w w w u ///(2121122211221/)/().x y k x y x y H H m H H H H º w w ¼ (14) z k z y y x y x k zx H m H H H H H H m T w w /)((/121221()//21221y k y z k z y H m H H m H H w w y k y x k x y H m H H m H H w w u //(21122211221/)/(),x y k x y x y H H m H H H H º w w ¼z k z x y x y x k xy H m H H H H H H m T w w /)((/211221()//12112x k x z k z x H m H H m H H w w w w y k y x k x y H m H H m H H w w w w u //(21122211221/)/().x y k x y x y H H m H H H H º w w ¼ y z k z y y x y x k zx H H m H H H H H H m T w w //)((/21121221()//121221z y k y z k z y H H H m H H m H H w w x y k y x k x y H H m H H m H H w w u ///(2121122211221/)/(),x y k x y x y H H m H H H H º w w ¼>x z k z x y x y x k xy H H m H H H H H H m T w w w w w w //)((/12211221()//212112z x k x z k z x H H H m H H m H H w w w w x y k y x k x y H H m H H m H H w w w w w u ///(2121122211221/)/().x y k x y x y H H m H H H H º w w ¼ ky z k z y y x y x k zx m H H m H H H H H H m T w w //)((/21121221)()//21121221z y z y k y z k z y H H H H m H H m H H w w k x y k y x k x y m H H m H H m H H w w u ///(2121122211221/)/(),x y k x y x y H H m H H H H º w w ¼>k x z k z x y x y x k xy m H H m H H H H H H m T w w w w w w //)((/12211221)()//12212112z x z x k x z k z x H H H H m H H m H H w w w w k x y k y x k x y m H H m H H m H H w w w w w w u ///(2121122211221/)/().x y k x y x y H H m H H H H º w w ¼ (k y z k z y y x y x k zx m H H m H H H H H H m T w w //)((/21121221)()//21121221z y z y k y z k z y H H H H m H H m H H w w k x y k y x k x y m H H m H H m H H w w u ///(2121122211221/)/(),x y k x y x y H H m H H H H º w w ¼k x z k z x y x y x k xy m H H m H H H H H H m T w w //)((/12211221)()//12212112z x z x k x z k z x H H H H m H H m H H w w w w k x y k y x k x y m H H m H H m H H w w w w w w u ///(2121122211221/)/().x y k x y x y H H m H H H H º w w ¼ (14) k y z k z y y x y x k zx m H H m H H H H H H m T w w //)((/21121221)()//21121221z y z y k y z k z y H H H H m H H m H H w w k x y k y x k x y m H H m H H m H H w w u ///(2121122211221/)/(),x y k x y x y H H m H H H H º w w ¼k x z k z x y x y x k xy m H H m H H H H H H m T w w //)((/12211221)()//12212112z x z x k x z k z x H H H H m H H m H H w w k x y k y x k x y m H H m H H m H H w w w w w w u ///(2121122211221/)/().x y k x y x y H H m H H H H º w w ¼ (14) (14)Substituting equation (13) into equation (14) anddifferentiating each tipper element at the jth observationsite (∂T zxj /∂m k , ∂T zyj /∂m k ) with respect to the kth modelparameter m k yields:1112/(/)(/),T T zxj k jtx k jtx k T m g K K m H g K K m Hw w w w w w 2112/(/)(/),T T zxj k jtx k jtx k T m g K K m H g K K m H w w w w w w1112/(/)(/).TT zyj k jty k jty k T m g K K m H g K K m H w w w w w w 2112/(/)(/).T T zyj k jty k jty k T m g K K m H g K K m H w w w w w w >))((2212211Tjy h z T jz h y y x y x T jtx g H g H H H H H g 211222122()()/(h T h Ty z y z x jy y jx x y x H H H H H g H g H H H º ¼))((2212211T jy h z T jz h y y x y x T jtx g H g H H H H H g 22112221221()()/(),h T h T y z y z x jy y jx x y x y H H H H H g H g H H H H º ¼>))((1112212T jz h y T jy h z y x y x T jtx g H g H H H H H g 211211122()()/(h T h T y z y z y jx x jy x y x H H H H H g H g H H H º ¼>))((1112212Tjz h y T jy h z y x y x T jtx g H g H H H H H g 211211122()()/(h T h T y z y z y jx x jy x y x H H H H H g H g H H H º ¼))((1112212Tjz h y T jy h z y x y x T jtx g H g H H H H H g 22112111221()()/(),h T h T y z y z y jx x jy x y x y H H H H H g H g H H H H º ¼))((2212211T jzh x T jx h z y x y x T jty g H g H H H H H g 122122122()()/(h T h T x z x z x jy y jx x y x H H H H H g H g H H H º ¼>))((2212211Tjz h x T jx h z y x y x T jty g H g H H H H H g 122122122()()/(h T h T x z x z x jy y jx x y x H H H H H g H g H H H º ¼))((2212211T jz h x Tjxhz y x y x T jty g H g H H H H H g 21221221221()()/(),h T h T x z x z x jy y jx x y x y H H H H H g H g H H H H º ¼and>))((1112212T jx h z T jz h x y x y x T jty g H g H H H H H g 122111122()()/(h T h T x z x z y jx x jy x y x H H H H H g H g H H H Hº ¼>))((1112212Tjx h z T jz h x y x y x T jty g H g H H H H H g 122111122()()/(h T h T x z x z y jx x jy x y x H H H H H g H g H H H H º ¼))((1112212T jx h z T jz h x y x y x T jty g H g H H H H H g 21221111221()/().h T h T x z x z y jx x jy x y x y H H H H H g H g H H H H º ¼ , (15)where h g T jx is the interpolator vector between H jx and H , h g T jyis the interpolator vector between H jy and H , and h g T jz is the interpolator vector between H jz and H . H jx , H jy and H jz denote the x , y , and z components of the magnetic field at the jth observation site, respectively.5Lin et al.Since the tipper data are used in the inversion and the elements of the Jacobian matrix are the partial derivatives of the data with respect to the model parameters, ∂T zxj /∂m k and ∂T zyj /∂m k in equation (15) are the elements of the Jacobian matrix A . Then, A (m ref )T V -1e can be expressed as:111()(/)().NTT ref n k n n A m V e D m V e w w ¦,(16)where D n is the tipper data, N is the total number of tipperdata (number of observation sites × frequency × 4), and k = 1, 2…, M denotes model parameters.Substituting equation (15) into equation (16) and since matrix K is symmetric, we get:11111211()/()/().NNT T T ref k n n k n n n n A m V e H K m K g V e H K m K g V e w w w w ¦¦1211211/()/().NNT Tk n n k n n n n H K m K g V e H K m K g V ew w w w ¦¦ , (17)where 1g n and 2g n are described in equation (15). Forexample, 1g nT = 1g Tjty and 2g nT = 2g Tjty , if T zyj is the nth datavalue.If we set11111()Nn n n v K g V e ¦,(18)and12121(),Nn n n v K g V e ¦(19)then1111(),N n n n Kv g V e ¦(20) 2121(),N n n n Kv g V e ¦(21)11122()//.T T T ref k k A m V e H K m v H K m v w w w w (22)If we take v 1 as the vector to be solved and n Nn n e V g )(111 ¦as the vector on the right side of the equation, equation (20) is similar to the forward modeling equation (8). Therefore, we call equation (20) a pseudo-forward modeling equation. The value of v 1 can be obtained by solving the pseudo-forward modeling problem one time. Similarly, the value of v 2 can be acquired by solving equation (21). Finally, A (m ref )T V -1e in equation (22) can be obtained with the values of v 1 and v 2.Similarly, derivation of A (m ref )p can be expressed as:11111()(/)(/MMT T ref nk k nk k k A m p g K K m H p g K K m Hw w w w ¦¦211211()(/)(/).M MT Tref n k k n k k k k A m p g K K m H p g K K m H pww w w ¦¦ (23)If we set 1111(/),Mkkk t K K m H p w w ¦(24) 1221(/),Mkk k t K K m Hp w w ¦(25)Then111(/),M k k k Kt K m H p w w ¦ (26)and221(/),Mk k k Kt K m H p w w ¦(27) 1212().T Tref n n A m p g t g t(28)If we take t 1 as the vector to be solved and ¦ w w Mk kk p H m K 11)/(as the vector on the right side of the equation (26), equation (26) also is similar to the forward modeling equation (8). Therefore, we call equation (26) a pseudo-forward modeling equation. The value of t 1 can be obtained by solving the pseudo-forward modeling problem one time. Similarly, the value of t 2 can be acquired by solving equation (27). Finally, A (m ref )p in equation (28) can be obtained by substituting the values of t 1 and t 2 into equation (28).Synthetic exampleBased on the theory and equations, we have developed a procedure to realize the 3D conjugate gradient inversion algorithm for inverting full information data. To verify the correctness of this algorithm, we designed two test models composed of two prisms.Model 1This model is composed of two conductive rectangular prisms which are embedded in a 100 Ω•m half-space6Three-dimensional conjugate gradient inversion(see Figure 1). The 10 Ω•m prism (ρ1) has a dimensions of 1.25 km × 1.25 km × 0.4 km with the top located 0.2 km below the earth’s surface. The 1 Ω•m prism (ρ2) has a dimensions of 1.75 km × 1.75 km × 1.55 km with the top located 1.05 km below the earth’s surface. The size and position of the two prisms are shown in Table 1. The size of the 3D grid is 36 × 36 × 37 (including 10 air layers). The impedance tensor and the tipper data atfour frequencies (10 Hz, 3.3 Hz, 1 Hz, and 0.33 Hz) at the grid center on the earth’s surface (289 sites) were calculated by staggered-grid finite difference modeling. The impedance tensor and tipper data were inverted using the 3D conjugate gradient inversion algorithm after Gaussian random noise was added to the data. All results were computed on a PC with a 2.2 GHz Intel (R) core (TM) 2 Duo CPU and 2 GB of physical memory.Fig. 1 View of the two-prism model 1.(a) Side view along negative Y direction. (b) Top view from negative to positive Z direction.Table 2 Statistical data for three scenariosData used in the inversionIteration numbersRuntimes (minutes)Data fit Z xx + Z xy + Z yx + Z yy662212 1.047T zx + Z zy4755380.997Z xx + Z xy + Z yx + Z yy + T zx + Z zy8378011.071The 3D inversion was carried out in three scenarios: 1)inverting the impedance tensor data (Z xx , Z xy , Z yx and Z yy ); 2) inverting the tipper data (T zx and T zy ); and 3) jointly inverting both the impedance tensor and the tipper data (Z xx , Z xy , Z yx , Z yy , T zx and T zy ). One percent Gaussian random noise was added to the impedance tensor data and three percent Gaussian random noise was added to the tipper data. Table 2 shows the total number of iterations, the CPU runtime, and the data fit for the three scenarios.The inversion results are shown in Figure 2. The top row shows the test model, the second row shows the results from the inversion of the impedance data, the third row shows the results from the inversion of the tipper data, and the fourth row shows the results from the inversion of both the impedance and the tipper data. The black dashed lines show the prism margins. The first column shows the horizontal slices at 0.35 km depth, the second column shows the horizontal slices at 1.8 km depth, the third column shows the vertical slices at Y = -0.375 km along the X axis, and the fourth column shows the vertical slices at X = -0.375 km along the Y axis. From this figure we see that the results of inverting the impedance tensor data show the shallow prism better. For instance, the size, shape, and resistivity (the lowest resistivity value is the same as the real value, 10 Ω•m) of the shallow prism are well recovered. Only its center is a little bit deeper than the real model, the center of the deep prism is much shallower than the real model, and the lowest resistivity value is only recovered to 23 Ω•m. Relatively, the results of inverting the tipper data recover. the deep prism much better. The size and shape of theTable 1 The size and position of the two prismsPrism Resistivity In the x direction In the y direction In the z direction Depth of top ρ110 Ω•m -1.25 Km, 0 Km -1.25 Km, 0 Km -0.2 Km, -0.6 Km 0.2 Km ρ21 Ω•m-0.75 Km, 1 Km-0.75 Km, 1 Km-1.05 Km, -2.6 Km1.05 Kmρ1ρ21.75 km1.25 km1.75 km1.25 kmXYSurface0.2 km0.4 km1.05 km1.25 km 1.55 kmρ2ρ11.75 kmXZ(a) (b)。

多孔介质模型的三维重构方法_王波

多孔介质模型的三维重构方法_王波

0226 收稿日期: 2012“页岩气流动机理与产能预测模型研究 ” ( 编号: 311008 ) 基金项目: 教育部科学研究重大项目 ), mail: vongpo1986@ sina. com 作者简介: 王波( 1986男, 主要从事页岩气渗流机理和产能预测方面的研究 . E-
王波等: 多孔介质模型的三维重构方法
应 用 CT 扫 描 得 到 了 分 辨 率 小 于
[3 ]
1 μm的三维多孔介质, Arns 和 Knackstedt 等
利用
板, 将 MCMC 方法从二维扩展到三维, 选取了 9 种 分别重构出对应的三维多孔介质模 岩石样品切片, 还有一些组合方法: 过程法和模拟退火方 型. 此外, 法组合
[13 ] [14 ] , 高斯场方法和模拟退火方法组合 . 这
[6 ]
过程繁琐等因素, 实用性较差. CT 扫描法只 成本高、 可以得到足够精确的多孔介质模 要分辨率足够高, 型. 一般的 CT 扫描分辨率在微米级别, 虽然现在已 但是其最高分辨率只在 50nm 左 有 Nano - CT 技术, 右
[15 ]
, 因此应用 CT 扫描法容易忽略掉一些小尺度
在 1974 年提出. 该方
Comparison among digital reconstruction methods 适用性分析 仅适用于各向同性多孔介 重构的多孔介质连通性差, 质. 方法中可以考虑任意多的约束条件, 重构的多孔介质 连通性差, 仅适用于各向同性多孔介质 . 重构的多孔介质连通性差, 仅适用于各向同性多孔介 质. 可以建立各向异性的多孔介质, 重构的多孔介质连通 但过程复杂, 仅适用于成岩过程简单的岩石 . 性好, 可以建立各向异性的多孔介质, 重构的多孔介质连通 适用范围广泛. 但计算速度慢. 性好, 可以建立各向异性的多孔介质, 重构的多孔介质连通 性好, 计算速度快, 适用范围广泛 .

computational ˉuid dynamic models

computational ˉuid dynamic models

Calculating particle-to-wall distances in unstructured computational ¯uid dynamic modelsMauro Tambascoa,b ,David A.Steinman a,b,*a Imaging Research Labs,The John P.Robarts Research Institute,100Perth Dr.,P.O.Box 5015,London,Ont.,Canada N6A 5K8b Department of Medical Biophysics,University of Western Ontario,Medical Sciences Building,London,Canada N6A 5C1Received 27March 2000;received in revised form 31October 2000;accepted 4December 2000AbstractKnowledge of particle deposition is relevant in biomedical engineering situations such as computational modeling of aerosols in the lungs and blood particles in diseased arteries.To determine particle deposition distributions,one must track particles through the ¯ow ®eld,and compute each particle's distance to the wall as it approaches the geometric surface.For complex geometries,unstructured tetrahedral grids are a powerful tool for discretizing the model,but they complicate the particle-to-wall distance calculation,especially when non-linear mesh elements are used.In this paper,a general algorithm for ®nding minimum particle-to-wall distances in complex geometries constructed from unstructured tetrahedral grids will be presented.The algorithm is validated with a three-dimensional 90°bend geometry,and a comparison in accuracy is made between the use of linear and quadratic tetrahedral elements to calculate the minimum particle-to-wall distance.Ó2001Elsevier Science Inc.All rights reserved.Keywords:Particle deposition;Lagrangian particle tracking;Computation ¯uid dynamics;Finite element methods;Unstructured mesh generation;Boundary layer1.IntroductionPredictions of particle deposition are relevant in biomedical engineering situations involving laminar and turbulent ¯ow.Examples include deposition of aerosols in lung bifurcations [1,2],and deposition of blood particles in stenosed and/or branching arteries [3,4].In lung bifurcations,localized enhanced aerosol deposition has been related to clinically observed incidences of bronchogenic carcinomas [5].In stenosed arteries,locally enhanced platelet or monocyte depo-sition on the artery walls may lead to thrombosis and/or plaque development [6±8].To quantify the distribution of particle deposition using computational models,each point particle (with ``virtual radius''r p )must be tracked through the computational ¯ow ®eld,and deposition is as-sumed to occur if the particle's trajectory comes within a distance r p to the wall.One approach that may be used to determine whether a particle has intercepted the wall is to reduce the wall Applied Mathematical Modelling 25(2001)803±814*Corresponding author.Tel.:+1-519-685-8300;fax:+1-519-663-3900.E-mail address:steinman@irus.rri.on.ca (D.A.Steinman).0307-904X/01/$-see front matter Ó2001Elsevier Science Inc.All rights reserved.PII:S 0307-904X (01)00014-2distances by a particle diameter and check to see if the particle remains within the diameter-re-duced model.This approach was taken in work quantifying aerosol particle deposition in bronchial bifurcations [1±3,9,10].In these studies,the bronchial bifurcation geometries were constructed from a parent branch and two daughter branches that were approximated mathe-matically by three straight-tube sections joined by a central zone.Although this method works well in geometries de®ned mathematically by a set of parameters,it would be computationally intensive to carry out in complex geometries that cannot be divided into a small set of para-metrically de®ned regions.Another drawback of this method is that the particle-to-wall distances are not explicitly calculated.As a result,additional coordinate transformation computations are necessary if one wants to extend the simulation to more complex situations involving,for ex-ample,particles rolling along the model boundary.Another approach to determine whether a particle has intercepted the wall is to compute the distance of the particle to the wall as it approaches the surface geometry.However,due to computational discretization of the surface,one must address the problem of ®nding the surface face having a point on its surface that is closest to the particle.For a simple geometry represented by a structured mesh of hexahedral elements,this problem is trivial since the surface face that is closest to the particle will simply be a wall face of the element in which the particle resides.Moreover,if the structured mesh is composed of linear elements the surface geometry will be approximated by ¯at faces,and in this case a simple analytic formula for computing the minimum distance between the particle and the ¯at surface may be used [11].However,for a complex ge-ometry consisting of sharp bends and/or high curvatures,it may be necessary to compute the particle's distance to the wall in an unstructured mesh composed of (possibly superlinear)tet-rahedral elements.With advances in digital medical image processing and analysis techniques,computational models constructed for anatomically realistic geometries are becoming more common [12].These ``realistic''models usually have very complex geometries that cannot be represented by a set of mathematical parameters and are di cult to discretize with a structured mesh.In these situations,an unstructured tetrahedral element mesh permits an alternative approach that easily accom-modates the complex geometric domains.However,in using an unstructured ±and,in some cases,nonlinear ±grid [13]to subdivide the geometry,the particle-to-wall distance calculation is no longer trivial.In this paper,a method is presented for calculating a particle's distance to complex surfaces represented by linear and quadratic surface faces derived,respectively,from unstructured meshes composed of linear and quadratic tetrahedral ®nite elements.A general algorithm for determining the ®nite element surface face that contains the point closest to the particle position is also presented.To verify the algorithms and compare linear and quadratic surface approximations for unstructured meshes,particles are tracked through an analytical velocity ®eld de®ned on a three-dimensional 90°bend geometry.2.Minimum particle-to-wall distanceFor a given continuous parametric surface,S x * u ;v n x u ;v ;y u ;v ;z u ;vo de®ned on an open set D of R 2,the minimum distance between the surface and a point particle at position p *not lying on the surface may be found by minimizing804M.Tambasco,D.A.Steinman /Appl.Math.Modelling 25(2001)803±814d u ;v x * u ;v Àp * 2 1with respect to the parameters u and v .However,if the parametric surface is approximated by many unstructured piecewise continuous surface faces the problem of ®nding the minimum dis-tance of p *to the surface is no longer as straightforward as simply evaluating (1).For this discrete case,one must:1.Find the set of surface faces that form a surface patch that contains a point that is closest to the point particle.Such a surface patch will be thought of as ``covering''the point particle.(Fig.1).2.Calculate (using (1))the minimum distance between the point and each surface face that covers the point.3.Determine the smallest of all these distances.This value will be the closest distance between the point and the surface.The algorithm presented below for step (1)is a general method that may be applied to any 2-D surface face type.The minimization method presented for step (2)is also general,however the detailed calculations of the minimization depend on the nature of the surface faces used to dis-cretize the surface.An unstructured 3-D mesh is most commonly composed of tetrahedral ele-ments,and a subdivision of the geometry with these elements creates a surface consisting of contiguous triangular surface faces.Hence,as an example of the above calculation the case of computing the distance of a point to a surface that has been subdivided into contiguous triangular surface faces that bound an arbitrary 3-D tetrahedral mesh geometry will be considered.2.1.Patch ®nding algorithmLet S i denote a surface face on the wall of the mesh and let the set M f S 1;S 2;...;S m g denote the patch of surface faces that cover the point particle p *i x p ;y p ;z p within the mesh.In addi-tion,let q *i x q i ;y q i ;z q i denote the point on surface face S i that is closest to p *i .As mentioned above,the ®rst step to calculating the minimum distance of a point to a surface is to ®nd the patch M .If the point particle's ``virtual radius''r p is much smaller than the smallest distance between two connected nodes in the mesh,then the particle-to-wall distance calculation is only performed if the particle lies within a wall element (i.e.,a volume element that has at least one node attached to the wall).This is usually the case for particle sizes of interest in the bronchial or arterial models.Hence,for a ®nite element mesh composed of tetrahedral elements and triangular surface faces,M will be determined by applying the following algorithm:Case1:The wall element within which p*resides has only one node g1on the wall(Fig.2(a)). Carry out steps1±5below.1.Calculate the point q*i on each of the surface faces S i attached to the wall node.2.Check to see if each point q*i lies within its corresponding surface face S i that is,q*i P S i.The set of surface faces M f S1;S2;...;S m j q*i P S i g forms a patch covering the point.3.If M is an empty set,that is,if for each surface face S i;q*i P S i,then search the mesh for the surface corner wall node g c closest to q*i.4.If g c g1,then repeat steps2and3.5.If g c g1,or one still has that for each surface face S i,q*i P S i then one assumes the point par-ticle p*must lie at the base of a surface enclosing a concave volume(Fig.3).In this case thepoint on the surface that gives the minimum distance to p*will lie on the curve at the base of the surface of the concave volume,and this curve will form the patch covering the point.Proceed with steps6and7.6.Find nodes attached to g c that de®ne edges on the surface faces S i containing g c.7.Determine the closest distance of p*to each edge,and keep the smallest distance.Case2:The wall element within which p*resides has several nodes on the wall that de®ne a surface face S e(Fig.2(b)).1.Find the corner node of surface face S e that is closest to the point p*.2.Calculate the point q*i on each of the surface faces S i attached to this corner node.3.Repeat steps3±5in Case1.3.Minimum particle-to-wall distance calculations and surface face checks3.1.Determining whether a particle is within an elementThe global coordinates x* x;y;z de®ne a computational mesh that subdivides a geometry and forms the physical or global space.To®nd the location of a point particle p*that is being tracked through this global space one must®rst locate the element that contains the particle.To do this,the particle's global coordinates and the nodes of the element in question are mapped to a natural or local(non-dimensional)coordinate system a;b;c (Fig.4).Four nodal coordinates are used to de®ne a linear tetrahedral element,and the mapping from local to global coordinates is given by:x* a;b;c x*1 x*2Àx*1 a x*3Àx*1 b x*4Àx*1 c; 2 where x*i x i;y i;z i ;i 1;2;3;4,represent the nodal coordinates.Eq.(2)can be inverted ana-lytically because it does not have any non-linear terms.In the case of quadratic tetrahedral ele-ments,10nodal coordinates are used to specify an element,and the mapping from local to global coordinates is given by:x* a;b;c c*1 c*2a c*3b c*4c c*5ab c*6ac c*7bc c*8a2 c*9b2 c*10c2; 3wherec *1 x *10;c *2 4x *7Àx *1À3x *10;c *3 4x *8Àx *3À3x *10;c *4 4x *9Àx *5À3x *10;c *5 4x *2À4x *7À4x *8 4x *10;c *6 4x *6À4x *9À4x *7 4x *10;c *7 4x *4À4x *8À4x *9 4x *10;c *8 2x *1À4x *7 2x *10;c *9 2x *3À4x *8 2x *10;c *10 2x *5À4x *9 2x *10;4 and x *i x i ;y i ;z i ;i 1;...;10,are the nodal coordinates located on the corners and mid-sides of the elements (Fig.5).Eq.(3)is non-linear,and hence must be inverted numerically.Once Eq.(2)or (3)are inverted and the local coordinates a ;b ;c are evaluated,one may proceed to determine whether a particle lies within an element.In order for particle p *to lie within an element,its local coordinates must satisfy a set of conditions (which depend on the de®ned global to local coordinate mapping).For a tetrahedral mesh whose mapping is de®ned by Eq.(2)or (3),the following four conditions must be valid for p *to lie within the tetrahedron:06a 61:06b 61:06c 61:06k 1Àa Àb Àc 61:5 Condition (5)may be used to determine whether a particle has landed in a wall element.3.2.Minimum distance of a particle to a wall faceLet S e x *e u ;v n x e u ;v ;y e u ;v ;z e u ;v j 061Àu Àv 61o 6Fig.5.(a)A quadratic tetrahedral element shown in the local coordinate system.The corner and mid-side nodes de®ne the four surface faces as follows:Face (1):nodes 1,2,3,4,5,6;Face (2):nodes 1,7,10,8,3,2;Face (3):nodes 1,6,5,9,10,7;Face (4):nodes 3,8,10,9,5,4.808M.Tambasco,D.A.Steinman /Appl.Math.Modelling 25(2001)803±814represent a surface face parameterized by local coordinates u ;v ,then the minimum distance between the surface face and a point particle p *may be found by minimizingd e u ;v x *e u ;v Àp *2 7 with respect to u ;v .That is,the minimum of d e occurs when r d e 0,where r o =o u ;o =o v is the gradient operator.In the case of tetrahedral elements,each face may be parameterized in terms of two of the three local coordinates a ;b ;c (Fig.5).To illustrate this minimization procedure let us consider the case of quadratic tetrahedral el-ements.The parametric equation for a surface of a quadratic tetrahedral element may be obtained by collapsing the mapping described by Eq.(3)to one of the four faces.For example,for face 3in Fig.5,b 0,so Eq.(3)reducesx * a ;0;c x *e a ;c c *1 c *2a c *4c c *6ac c *8a 2 c *10c 2:8 To ®nd the minimum distance of p *to this face,Eq.(8)is substituted into (7)with u ;v a ;c ,and d e is minimized with respect to a ;c .In general,for a given tetrahedral face that lies on the surface of the geometry,Eq.(3)reduces to an equation of the formx *e u ;v c *a c *b u c *c v c *d uv c *e u 2 c *f v 2; 9 where u ;v represent two of the three local coordinates a ;b ;c that parametrize one of the four tetrahedral faces that may lie on the surface,and the c *'s represent the corresponding coe cients from Eq.(4).Next,Eq.(9)is substituted into (7),to getd e u ;v c *a c *b u c *c v c *d uv c *e u 2 c *f v 2Àp * 2:10 Taking the gradient of Eq.(10)to ®nd the minimum distance results in two coupled third-order polynomials of the formA 1u 3 A 2u 2 A 3u A 4 0;11 andB 1v 3 B 2v 2 B 3v B 4 0;12 where A i A i v ;c *a ;c *b ;c *c ;c *d ;c *e ;c *f ,B i B i u ;c *a ;c *b ;c *c ;c *d ;c *e ;c *f ,and i 1;2;3;4.Since Eqs.(11)and (12)cannot be solved analytically,Newton's method for nonlinear equations is used to ®nd the local coordinates u ;v .The initial local coordinate guess is taken to be the projection of the particle's local coordinates p *H (mapping of p *given by Eq.(3))onto the wall face.With this initial guess,convergence of the local coordinates occurs in one to ®ve iterations and typically in three.The global coordinates of the closest point x *e on the wall face to the particle p *are found by transforming u ;v to the global coordinate system using Eq.(9).In order for the global coordinate x *e to lie on the global surface face the corre-sponding local coordinates u ;v must satisfy condition (5)(that is,06u 61,06v 61,and 061Àu Àv 61).If this condition is satis®ed,the surface face lies above (covers)the point.Finally,the Euclidean distance from x *e to p *is calculated using the square root form of Eq.(7).M.Tambasco,D.A.Steinman /Appl.Math.Modelling 25(2001)803±8148094.Numerical validation and resultsTo verify the wall distance calculation codes,and to compare the accuracy of approximating a curved surface with quadratic as opposed to linear surface faces,the following 3-D geometry (Fig.6)was de®ned:16r 65;06h 6p =2;13 where r x 2 y 2p ,and h arctan y =x .On this 90°bend geometry an analytic velocity ®eld given byu *u x ;u y ;u z ÀÁ k sin h ; Àk cos h ;0 ; 14where k 0:003,was imposed.The above geometry (13),and velocity ®eld (14)were chosen so that a particle seeded within the geometry would be passively advected along a constant radial path (Fig.7).In this way,the particle's distance from the curved walls is known to be constant along its trajectory,and deviations in the particle wall distance calculations could be computed easily.In addition,the above geometry (13),allows one to compare the particle-to-wall distance calculation for two di erent surface curvatures (high curvature at r 1:0,lower curvature at r 5:0).The velocity ®eld was chosen to have a constant magnitude in order to yield identical particle trajectories for the linear and quadratic meshes.In doing this,the wall distance calcu-lations depend only on the order of surface face (quadratic or linear)used to approximate the surface geometry,and a meaningful comparison between the two surface face types may be made.To test the wall distance algorithm,quadratic and linear tetrahedral meshes for the 90°bend geometry were constructed.The linear meshes were derived from a subdivision of the quadratic meshes [14],and each linear mesh contained eight times the number of elements of the corre-sponding quadratic mesh from which it was derived.The total number of nodes in the meshes varied from 3972to 113702with corresponding nodal densities varying from 134nodes per unit volume to 3840nodes per unit volume (Fig.6).For each mesh,25particles were seeded equally spaced along the z -axis f z 0:15;0:20;...;1:35g above the lower wall at x * 0:0;1:003;z ,and below the upper wall atx* 0:0;4:997;z .The particle paths were calculated using the fourth-order Runge±Kutta method to solve the di erential equation:d x*d tu* r;h ; 15 where u*is given by Eq.(14).The time interval used for a given Runge±Kutta integration step was always chosen such that the particles traveled a®xed distance of3Â10À3units with each step along the velocity®eld.Also,since u*is only known at the nodes of a mesh,the velocity of a particle at each Runge±Kutta integration step was found by spatially interpolating the velocity ®eld from the known nodal velocities.Since the velocity®eld is de®ned such that a particle traces a constant radial path,it is known that particles seeded at radial distances of1.003and4.997units will trace paths that are0.003 units away from the lower and upper surfaces,respectively(Fig.7).In our numerical experiments, the algorithms described in this paper were used to compute the particle-to-wall distances for the lower and upper surfaces of the90°bend.To determine the accuracy of the wall distance cal-culation,the average percentage wall distance errorD errorD trueÀD compÀÁD trueÂ100% 16812M.Tambasco,D.A.Steinman/Appl.Math.Modelling25(2001)803±814was computed.In Eq.(16),D true represents the true distance of a particle to the surface(0.003 units in this case),and D comp represents the average computed distance of the particles to one of the surfaces.To compute D comp for a given surface(lower or upper),each particle's distance to the surface was calculated at equal intervals of0.03units as they traveled along their tra-jectories.Next,for each surface,the average of all the particle-to-wall distances was calculated. To illustrate the dependence of the wall distance error on mesh resolution,D error was plotted as a function of increasing average nodal density(average number of nodes per unit volume) (Fig.8).Due to surface discretization,one expects that the particle-to-wall distances will be overes-timated for particles traveling above the lower surface and underestimated for particles trav-eling below the upper surface resulting in positive and negative percentage wall distance errors, respectively.This is what was have obtained in our results(Fig.8).From the linear mesh results in Fig.8,one also observes that the percentage wall distance error is consistently larger for particles tracked just above the more highly curved lower surface,than for particles tracked just below the less curved upper surface.This is also expected since,for uniform element density,the discrete surface approximation will be better for lower curvatures,than for higher curvatures.In addition,it is observed,as expected,that in all the cases D error approaches zero as the average nodal density increases.Typical average nodal densities used in physiological models,for example the carotid artery bifurcation[13],vary from about80nodes per unit volume in the less re®ned regions to about600nodes per unit volume in the more re®ned regions.From Fig.8,one sees that in this range of nodal densities the absolute value of D error is between1%and2%for quadratic meshes and between15%and150%for linear meshes.HenceD error,as a function of the average nodal density.The error barsdeviations in D error.M.Tambasco,D.A.Steinman/Appl.Math.Modelling25(2001)803±814813the average particle-to-wall distance errors are signi®cantly smaller when quadratic surfaces are used.5.Discussion and conclusionMost®nite element and®nite volume CFD codes use linear elements,either tetrahedral or hexahedral,in which case,as mentioned above,it is trivial to compute the distance from a particle to the wall.FIDAP,a commercial®nite element CFD code that does support quadratic elements, is capable only of detecting whether a particle has left the computational domain.In this case the user may ignore the particle altogether,have the particle deposited at the wall,or have it ricochet back into the domain[15].In the present study,an e cient searching technique is presented for detecting not only when a particle has left the domain,but also when that particle is within a user-speci®c distance to the wall.This is particularly useful,for example,when tracking particles in the near wall region referred to as the``numerical boundary layer''[9,16],where even small pertur-bations to the velocity vector(i.e.,below convergence level or even within¯oating point error)can spuriously deposit a particle travelling along a near-wall ing the techniques de-scribed in this paper,such spurious depositions can be eliminated by convecting particles traveling within a®xed distance to the wall using only the component of the velocity vector instantaneously tangent to the wall.Admittedly,the quadratic particle-to-wall distance calculation involves more operations than the corresponding calculation for the linear case.However,one can signi®cantly reduce the time spent in calculating particle-to-wall distances for a given particle trajectory,by reducing the number of particle-to-wall distance checks.This may be done by allowing a particle within a wall element to travel®xed distances,D S i,between successive checks,where i 1to the total number of particle checks for a given trajectory.The distance D S i traveled between successive checks will depend on the particle-to-wall distance found at each check.For example,if D i is the particle's distance from the wall as determined from the i th check and r p is the particle radius,then a particle may safely travel a distance D S i D iÀr p in any direction before the next check must be made.In addition,it is important to note that the particle-to-wall distance calculation is per-formed for only a small percentage of the total number of particles seeded.For example,in our experience with particle deposition studies in the carotid artery bifurcation models it has been observed that for particle sizes of interest(0:5Â10À3to2:5Â10À3radial units)only3±4%of the total number of particles seeded come within a radial distance to the artery wall.Hence,this fact coupled with the e cient particle-to-wall distance checking scheme described above,leads to CPU times for particle deposition runs in quadratic meshes being comparable to those using linear meshes.This is because most of time is spent tracking the particles,and relatively little time is spent performing the particle-to-wall distance calculation.One may therefore conclude that for the large accuracy gained and only small computational e ciency lost,it is worth retaining quadratic surfaces to calculate particle-to-wall distances when quadratic®nite element meshes are employed.AcknowledgementsThe authors thank Jaques Milner for constructing the®nite element meshes used in this study. This work was supported by a Medical Research Council of Canada Group Grant(GR-14973). The®rst author was supported by the Ontario Graduate Scholarship in Science and Technology.814M.Tambasco,D.A.Steinman/Appl.Math.Modelling25(2001)803±814The second author was supported by a Research Scholarship from the Heart&Stroke Foun-dation of Canada.References[1]I.Bal a sh a zy,W.Hofmann,Particle deposition in airway bifurcations±I.Inspiratory¯ow,J.Aerosol Sci.24(6)(1993)745±772.[2]I.Bal a sh a zy,W.Hofmann,Particle deposition in airway bifurcations±II.Expiratory¯ow,J.Aerosol Sci.24(6)(1993)773±786.[3]I.Bal a sh a zy,W.Hofmann,T.Heistracher,Computation of local enhancement factors for the quanti®cation ofparticle deposition patterns in airway bifurcations,J.Aerosol Sci.30(2)(1999)185±203.[4]D.Bluestein et al.,Steady¯ow in an aneurysm model:correlation between¯uid dynamics and blood plateletdeposition,J.Biomech.Eng.118(1996)280±286.[5]R.B.Schlesinger,M.Lippmann,Selective particle deposition and bronchogenic carcinoma,Environ.Res.15(1978)424±431.[6]D.Bluestein et al.,Fluid mechanics of arterial stenosis:relationship to the development of mural thrombus,Ann.Biomed.Eng.5(1997)344±356.[7]W.F.Pritchard et al.,E ects of wall shear stress and¯uid recirculation on the localization of circulatingmonocytes in a three-dimensional¯ow model,J.Biomech.28(12)(1995)1459±1469.[8]R.T.Schoephoerster et al.,E ects of local geometry and¯uid dynamics on regional platelet deposition on arti®cialsurfaces,Arteriosclerosis and Thrombosis13(12)(1997)1806±1813.[9]I.Bal a sh a zy,Simulation of particle trajectories in bifurcating tubes,put.Phys.110(1994)11±22.[10]T.Heistracher,W.Hofmann,Physiologically realistic models of bronchial airway bifurcations,J.Aerosol Sci.26(3)(1995)497±509.[11]C.Green®eld,G.Quarini,A Lagrangian simulation of particle deposition in a turbulent boundary layer in thepresence of thermophoresis,Appl.Math.Model.22(1998)759±771.[12]Q.Long et al.,The combination of magnetic resonance angiography and computational¯uid dynamics:a criticalreview,Crit.Rev.Biomed.Eng.26(4)(1998)227±274.[13]ner et al.,Hemodynamics of human carotid artery bifurcations:computational studies with modelsreconstructed from magnetic resonance imaging of normal subjects,J.Vascular Surg.28(1)(1998)143±156. [14]A.Liu,B.Joe,Quality local re®nement of tetrahedral meshes based on8-subtetrahedron subdivision,Math.Comput.65(1996)1183±1200.[15]Fluent Inc.,Personal communication,2000.[16]K.A eld et al.,Fluid mechanics of the stagnation point¯ow chamber and its platelet deposition,Arti®cial Organs19(7)(1995)597±602.。

伊朗法尔斯地区第三系砂岩的古地磁

伊朗法尔斯地区第三系砂岩的古地磁

a,
Dominique
Frizon de Lamotte b
a ~
' , C h a r l e s A u b o u r g a,
Jamshid Hassanzadeh
" Universitg de Cergy-Pontoise, Dept. des Sciences de la Terre (CNRS, URAI759), F95011, Cergy-Pontoise Cedex, France I~Institute of Geophysics, Tehran UniversiO, PO. Box 14155-6466, Tehran, lran
Keywords: fold-thrust belt; magnetic fabric; sandstone; weak deformation; Arc of Fars; Zagros (Iran)
1. I n t r o d u c t i o n In sedimentary rocks undergoing horizontal shortening, the initial sedimentary fabric is progressively erased and replaced by a tectonic one (Ramsay and Huber, 1983). The analysis of these initial stages of deformation during which the inherited sedimentary fabric and the tectonic fabric interact, is generally not well documented due to the subtlety of the ini-

CdS薄膜的光学及其电学性能的研究

CdS薄膜的光学及其电学性能的研究

CdS薄膜的光学及其电学性能的研究何智兵 韩高荣3(浙江大学材料系,硅材料国家重点实验室 杭州 310027)Optical and E lectrical Properties of CdS FilmsHe Zhibing and Han G aorong3(Material Department,Zhejiang Univer sity,Hangzhou,310027,China) Abstract Optical and electrical properties of CdS films grown by vacuum deposition at different substrate temperatures,including its transmittance,its optical band gaps and its transversal and in2plane resistivites,were studied with conventional techniques.The results showed that the CdS film had a hexag onal structure with a preferential growth orientation of〈002〉.Large difference between transversal and in2plane resistivities was als o observed. K eyw ords CdS,Preferred orientation,Vacuum evaporation,Photosensitivity 摘要 本论文研究了真空热蒸发制备的不同衬底温度下CdS薄膜的光学及电学性能。

重点讨论了不同衬底温度下CdS 薄膜的光透过性能、光学带隙及其和薄膜结构的关系。

研究了CdS薄膜因高度定向生长而出现的薄膜平面与截面电阻率的差别。

Tunnelling and Underground Space Technology

Tunnelling and Underground Space Technology
Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction
Stress redistribution during tunneling should be three-dimensional (3D), with the exception of a two-dimensional (2D) plain strain condition. However, ideal assumptions of a circular-shaped tunnel and the 2D condition are typically invoked for the initial analytic solution. This stress distribution solution for a circular tunnel has been reported throughout the last century. Kirsch (1898) solved analytically the distribution of stress and displacement in an unsupported circular tunnel. His solution relied upon elasticity theory, using the plane stress condition with different K0 values. Bray (1967) proposed a theoretical model to permit analysis of the extent of failure, the plastic zone, based on Mohr–Coulomb failure criterion. Ladanyi (1974) discussed stress distribution around a circular opening in a hydrostatic stress field, and within annular failed rock generated in the excavation periphery, using Mohr–Coulomb elasto-plastic theory.

大口径KDP_晶体装配附加面形畸变抑制工艺优化

大口径KDP_晶体装配附加面形畸变抑制工艺优化

第 31 卷第 9 期2023 年 5 月Vol.31 No.9May 2023光学精密工程Optics and Precision Engineering大口径KDP晶体装配附加面形畸变抑制工艺优化全旭松1,独伟锋1,褚东亚1,2,周海1*,叶朗1(1.中国工程物理研究院激光聚变研究中心,四川绵阳 621900;2.清华大学机械工程系,北京100084)摘要:在高功率固体激光装置中,大口径KDP晶体的面形畸变控制是影响终端光学组件倍频转化效率的关键因素之一。

为了提高大口径KDP晶体的装配附加面形质量,提出了一种点支撑装配附加面形畸变抑制工艺方法。

首先,通过遗传算法对支撑点及其分布进行优化设计。

然后,采用有限元分析方法对KDP晶体的装配预紧工艺进行优化设计。

最后,开展优化后的装配工艺对KDP晶体装配附加面形畸变的抑制和倍频转换效率的实验验证。

实验结果表明:提出的工艺方法对KDP晶体装配附加面形畸变具有良好的抑制效果,实测面形PV值为6.51 μm,二倍频转化效率可达72.6%,且重复装配的一致性良好。

该方法大幅提升了晶体倍频效率和远场光斑质量,并在工程上得到应用与推广。

关键词:激光装置;KDP晶体;装配附加面形;点支撑;频率转换效率中图分类号:TN242;TH162 文献标识码:A doi:10.37188/OPE.20233109.1347Mounting optimization on large aperture KDP crystal tominimize assembling deformationQUAN Xusong1,DU Weifeng1,CHU Dongya1,2,ZHOU Hai1*,YE Lang1(1.Research Center of Laser Fusion, China Academy of Engineering Physics, Mianyang 621900, China;2.Mechanical Department, Tsinghua University, Beijing 100084, China)* Corresponding author, E-mail: a697097@Abstract:In the high-power laser facility,control of the surface deformation of the large-aperture KDP crystal is the key factor to reduce the frequency-conversion efficiency. To improve the assembling quality of the KDP crystal, a point-supporting process method is proposed for minimizing the assembly deforma⁃tion. First, a genetic algorithm is used to optimize the support points and their distribution scheme. Sec⁃ond, the finite-element method is used to optimize the assembling preload. Finally, mounting optimization design process experiments are conducted to evaluate the surface deformation and the frequency-doubling conversion efficiency. The experimental results indicate that the proposed method is effective for minimiz⁃ing the assembling deformation of the KDP crystal; the measured PV value is 6.51 μm, and the measured conversion efficiency of second-harmonic generation reaches 72.6% with excellent assembling repeatabili⁃ty. This result significantly improves the frequency-doubling efficiency and the quality of the far-field spot and has been widely used and promoted in engineering.文章编号1004-924X(2023)09-1347-10收稿日期:2022-11-04;修订日期:2022-12-17.基金项目:国家自然科学基金资助项目(No.51975322);北京市自然科学基金资助项目(No. 3212006)第 31 卷光学精密工程Key words: laser facility;KDP crystal;assembly deformation;point-supporting;frequency conversion efficiency1 引言终端光学组件作为高功率固体激光装置中末端的核心单元,其主要功能之一为将波长为1 053 nm的高能红外激光转化为波长为351 nm 的紫外激光。

THE INTERNATIONAL JOURNAL OF MEDICAL ROBOTICS AND COMPUTER ASSISTED SURGERY Int J Med Robot

THE INTERNATIONAL JOURNAL OF MEDICAL ROBOTICS AND COMPUTER ASSISTED SURGERY Int J Med Robot

Introduction
Computer-assisted surgery (CAS) is a methodology that translates into accurate and reliable image-to-surgical space guidance. Neurosurgery is a very complex procedure and the surgeon has to integrate multi-modal data to produce an optimal surgical plan. Often the lesion of interest is surrounded by vital structures, such as the motor cortex, temporal cortex, vision and audio sensors, etc., and has irregular configurations. Slight damage to such eloquent brain structures can severely impair the patient (1,2). CASMIL, an imageguided neurosurgery toolkit, is being developed to produce optimum plans resulting in minimally invasive surgeries. This system has many innovative features needed by neurosurgeons that are not available in other academic and commercial systems. CASMIL is an integration of various vital modules, such as rigid and non-rigid co-registration (image–image, image–atlas and

18 MeV质子辐照对TiNi形状记忆合金R相变的影响

18 MeV质子辐照对TiNi形状记忆合金R相变的影响

第15卷 第1期强激光与粒子束Vol.15,No.1 2003年1月HIGH POWER LASER AND PAR TICL E B EAMS Jan.,2003 文章编号:100124322(2003)012009720418Me V质子辐照对Ti Ni形状记忆合金R相变的影响Ξ王治国, 祖小涛, 封向东, 刘丽娟, 林理彬(四川大学物理系教育部辐射物理及技术重点实验室,四川成都610064) 摘 要: 研究了用HZ2B串列加速器的18MeV质子辐照对TiNi形状记忆合金R相变的影响,辐照在奥氏体母相状态下进行。

示差扫描量热法(DSC)表明,辐照后R相变开始温度T s R和逆马氏体相变结束温度T f A随辐照注量的增加而降低。

当注量为1.53×1014/cm2时,T s R和T f A分别下降6K和13K,辐照未引起R相变结束温度T f R和逆马氏体相变开始温度T s A的变化。

表明辐照后母相(奥氏体相)稳定。

透射电镜(TEM)分析表明辐照后没有引起合金可观察的微观组织变化。

辐照对R相变开始温度T s R和逆马氏体相变结束温度Af的影响可能是由于质子辐照后产生了孤立的缺陷团,形成了局部应力场,引起晶格有序度的下降所造成的。

关键词: TiNi形状记忆合金;质子辐照;R相变;示差扫描量热法;TEM 中图分类号:TG139.6 文献标识码:A TiNi形状记忆合金是目前应用最为广泛,也最成功的一种智能材料,集传感功能与驱动功能于一体,在核反应堆和太空等核辐射环境下用作传感与驱动元件已引起了关注[1,2]。

TiNi合金中R相变具有热滞后小,响应速度快的特点,在实际应用中得到了广泛的应用[3]。

在以前的研究中利用R相变得到了具有双向记忆效应的弹簧,伸缩率可达25%[4]。

由于核辐射会对形状记忆合金相变特性产生影响,因而研究其改变规律及机理对形状记忆合金在辐射环境下应用的可靠性和可行性是十分必要的。

2004) Three-Dimensional Dynamic Fracture Analysis

2004) Three-Dimensional Dynamic Fracture Analysis

Abstract: This paper describes algorithms for threedimensional dynamic stress and fracture analysis using the material point method (MPM). By allowing dual velocity fields at background grid nodes, the method provides exact numerical implementation of explicit cracks in a predominantly meshless method. Crack contact schemes were included for automatically preventing crack surfaces from interpenetration. Crack-tip parameters, dynamic J -integral vector and mode I, II, and III stress intensity factor, were calculated from the dynamic stress analysis solution. Comparisons to finite difference method (FDM), finite element method (FEM), and boundary element method (BEM), as well as to static theories have shown that MPM can efficiently and accurately solve three-dimensional dynamic fracture problems. Since the crack description is independent of the object description, MPM could be useful for simulation of three-dimensional dynamic crack propagation in arbitrary directions. keyword: Material point method, MPM. threedimensional dynamic fracture, cracks, dynamic J-integral, dynamic stress intensity, contact, mode I, mode II, mode III. 1 Introduction

Dimensionless formulae for penetration depth of concrete target impacted by a non-deformable project

Dimensionless formulae for penetration depth of concrete target impacted by a non-deformable project

Q.M. Li, X.W. Chen / International Journal of Impact Engineering 28 (2003) 93–116
95
between 27 and 312 m/s and assessed published empirical formulae on concrete impact. Comparison between various empirical formulae and published test data was conducted by Williams [3]. Recently proposed formulae on penetration depth and perforation limit are discussed in Refs. [5,6]. Among the most commonly used formulae are the Petry formula, the Army Corps of Engineers formula (ACE), the UKAEA formula (i.e., the Barr formula) and the National Defence Research Committee (NDRC) formula [1–5]. These empirical formulae provide the most reliable, straightforward approach to design a concrete protective structure. However, some shortcomings, as shown below, in these empirical formulae confine their applications. First, most of the empirical formulae are not dimensionally homogeneous, leading to the disadvantage of unit dependency. This makes it difficult to identify important physical quantities in a given empirical formula. There are few exceptions, e.g., Chang formula [7], Haldar formula [8] and Hughes formula [9], which are unit independent. Second, definition of nose shape factor in the existing empirical formulae is not unique. For example, in NDRC and relevant formulae, the nose shape factor NN is defined as 0.72 for flat nose, 0.84 for blunt nose, 1.0 for spherical nose and 1.14 for sharp nose. In Hughes formula [9], the nose shape coefficient NH is chosen as 1.0, 1.12, 1.26 and 1.39 for a flat, blunt, spherical and sharp nose, respectively. Thus, it is necessary to define a unique nose shape factor for empirical penetration formulae. Thirdly, most of the published empirical formulae are valid for a small range of penetration depth in the low to medium impact velocity range. For example, it was shown in Refs. [2,3] that the modified NDRC formula and the Barr formula give good prediction of penetration depth only in the range of 0:6oX =d o2:0: When the impact velocity falls in the range of 500–1000 m/s, which may be associated with a KE weapon attack, a kinetic fragment or a tornado-generated projectile, X =d could be quite large. Nash et al. [10] reported penetration tests impacted by a solid steel projectile at velocities between 365 and 610 m/s. Iremonger et al. [11] conducted penetration tests using arms rounds at velocities between 800 and 924 m/s, which showed a relatively poor agreement with existing empirical formulae due to the high deformability of the projectile. Systematic studies on penetration of concrete with ogive-nose projectile have been conducted by Forrestal et al. [12,13] and Frew et al. [14], which covered a broad range of concrete strengths for striking velocities up to 1 km/s (X =d ranges from 6.21 to 92.8) until nose erosion becomes excessive. In the present paper, a dimensionless formula on penetration depth is proposed using two dimensionless numbers, i.e., the impact function I and the geometry function of projectile N ; obtained from a dimensional analysis and a penetration model based on dynamic cavity expansion theory. It is shown that these two dimensionless numbers are capable of representing test data on penetration depth in a broad range of impact velocity as long as the projectile deformation is negligible, which covers the applicable range of empirical formulae.

三维编织复合材料几何建模及弹性常数预测

三维编织复合材料几何建模及弹性常数预测

c r e n tc l ,s prpo e o n r u i el i o s d.T e v ra in o e d n haa t rsis a d c o ss cin ld f r t n f r h a t fb n i g c r ce it n r s -e to a eo mai i o c o o
图 2 纤维 束 横 截 面 形 状

12 表面单 胞和 角单胞 .
UW, P U轴为对 角线 1 一1 W轴位 于单 胞 l一1对 角 ,
面内, 轴垂 直 于 u—W平 面 。纤维 束 1—1位 于坐
根 据 三维 四 向编 织 复合 材 料 的编 织 工艺 , 面 表
关键词 : 三维编织 复合材料 ; 单胞 ; 弹性性能 ; 工程弹性常数
中 图分 类 号 : 30,3 3 8 0 4 0 4 . 文献 标 识 码 : A
Ge m e r c m o e o h e - i e so a r i e o p st s o t i d lf r t r e d m n i n lb a d d c m o ie
c l u ae aue y t e p e e td mo e l a r ewi h e s rd v l e . ac ltd v l sb h r s n e d lwe l g e t te m a u e a u s Nume c lr s lsv rf h h i r a e u t e y t e i efci e e so h o e . f t n s ft e m d 1 e v K e o ds:h e — m e so a r i e o o ie untc l; lsi r p ry; l siiy c n t n yw r t r e di n i n lb a d d c mp st s; i e l e a t p o e c t ea tct o sa t

基于高斯映射的柱面与锥面点云拟合_李岸

基于高斯映射的柱面与锥面点云拟合_李岸

k的初值 ,
m1 |m 1
×m 2 ×m 2
|作为
n (即参数
和 θ)的初值 , 最小曲率方向 m1 作为 a, 可以确定参数
α的初值 , 参数 ρ的 初值设为零 。 迭代法的初值 S =
(ρ, , θ, k, σ, τ)的确定与圆柱面相类似 , 不同的只是
圆锥面旋转轴线的位置和方向需要计算出两个数据点
2 高斯映射
对曲面 S 上每点 P , 可作出它的单位法向量 n。 因 为 |n |=1, 所以把向量 n 的起点平行地移到原点 O 后 , n的终点就是以在 O 为球心的单位球面 S2 上的一点 P ′。 把这种点的映照 , 称为曲面 S 的高斯映照 (Gauss 映照 )。 在高斯映照下 , 曲面 S 的像是单位球面内的一 个点集 ΢, 这个点集可能是球面上的一个点 , 也可能是 一条球面曲线 。 点集 ΢一般叫做高斯映像[ 5] 。 高斯映 射在实际工程中已经有了广泛的应用 。 M. H evert和 J. Ponce利用曲面法矢在高斯球上的分布 , 通过 H ough 技术将曲面类型分为平面 、柱面 、锥面等 [ 6] 。 H. Pottm ann提出在曲面类型识别时 , 可以利用基于法矢的高 斯映射图像识别柱面 [ 7] 。 P. Benko″等在讨论特征识别 时 , 也提到了利用基于法矢的高斯球使曲面特征可视 化 , 从而辅助识别出拉伸曲面的算法 [ 8] 。
310023, Ch ina)
Ab stract:
-
-
K ey w ord s:
0 前言
大多实物特别是机械零件的表面常常是由平面 、 球面 、圆柱面 、圆锥面以及圆环面等二次曲面构成 , 或 者是其重要的组成部分 , 所以在反求工程中 , 研究基于 点云数据的二次曲面拟合具有重要的意义 。 许多研究 人员 对 这 一 课题 进 行 了 大 量 的 研究 工 作 , C hen 和 L iu[ 1] 提 出 了 一 种 基 于 遗 传 算 法 (GA , genetic algorithm s)的一般二次曲面提取算法 。 V aughan P ra tt[ 2] 提 出了一种 “准最小二 乘 (Quasi-Least-Squares)”的sina. com。

超振荡透镜的三维全矢量电磁分析方法

超振荡透镜的三维全矢量电磁分析方法
L I U Ta o, YANG S h u mi n g,J I ANG Zh u a n g d e
( S t a t e Ke y La b o r a t o r y f o r Ma nu f a c t u r i ng Sy s t e ms En g i n e e r i n g,Xi ’ a r t J i a o t o n g Un i v e r s i t y,Xi ’ a n 71 0 0 4 9,Ch i n a )
Abs t r a c t : A t h r e e — di me n s i on a 1 c o m pl e t e ve c t or i a l e l e c t r o ma gn e t i c a n a l ys i s i s c o nd uc t e d t o
维 全 矢 量 电磁 分 析 方 法 , 包 括 矢 量 电 磁 聚 焦 和 矢 量 电磁 成 像 两 个 连 续 物 理 过 程 。 采 用 矢 量 角谱 理
论, 计 算预 测超振 荡透 镜后 多衍 射光 束精 细干 涉形 成 的矢量 光场 分布 , 利 用等价磁 偶板 子 矢量成像
理论, 定量 计算探 测 面 内像 场 分布 。研 究结 果表 明 : 三 维 时域有 限差 分法 建模 的严格 求解 结果与 矢 量 角谱 结 果基本 吻合 ; 高倍 、 大数值 孔径 显微 成像 系统具 有 选择 性偏 振 滤波 成像 机 理 , 特 别是 可 对
d e mo ns t r a t e t he s u bwa ve l e n gt h f o c u s i ng a n d hi gh— — r e s o l ut i o n i ma gi n g m e c h a n i s ms o f s u pe r — — o s c i l l a t or y l e ns ( SOL), wh i c h i nc l u d e s t wo s u c c e s s i v e p hy s i c a l p r o c e s s e s, e l e c t r o ma g ne t i c f o c us i n g a nd e l e c t r o ma gn e t i c i ma g i ng . The v e c t o r i a l f i e l d di s t r i bu t i o n f or me d b y d e l i c a t e i nt e r f e r e nc e o f m ul t i — d i f f r a c t i on be a ms f r o m SOL i s t he o r e t i c a l l y p r e d i c t e d a nd t he n v a l i d a t e d by a r i go r ou s r e s ul t c a l c ul a t e d b y t h e t hr e e — d i me ns i o na l f i ni t e — d i f f e r e n c e t i me — d o ma i n me t h o d. Th e e qu i v a l e nt m a g ne t i c — d i po l e v e c t o r i a l i ma gi n g t h e or y i s a p pl i e d t o qu a n t i t a t i v e l y e v a l ua t e t he i ma g i n g f i e l d i n t he de t e c t i o n pl a ne,a nd t he po l a r i z a t i o n f i l t e r i n g t r a n s mi s s i o n me c ha n i s m o f a hi gh - nu me r i c a l — a p e r t ur e mi c r os c o pi c i ma gi ng s ys t e m wi t h h i gh m a gn i f i c a t i on i s r e v e a l e d. Thi s a ppr o a c h pr o v i de s a t he o r e t i c a l b a s i s f or a p pl y i ng S OL t o t he s e f i e l d s s u c h a s n a no l i t ho g r a phy,

三维时域有限差分法对金sers基底电磁场增强的初步模拟(理学)

三维时域有限差分法对金sers基底电磁场增强的初步模拟(理学)
3.企业期间费用包括制造费用、销售费用、管理费用和财务费用。()
4.会计科目是对会计对象的具体内容进行分类的,它既有分类的名称,又有一定的格式。()
5.一项经济业务的发生,没有影响资产总额,是因为资产与负债项目增加或减少同一金额。
()
6.一般地讲,总分类账户要根据明细分类账户来登记,这样才有利于保证总分类账户与明细分类账户的一致性。()
20.根据我国企业会计制度的规定,企业的应收账款应按净价法确认。()
四、计算分析题
1.资料:湖南某公司2013 年3 月上旬发生下列部分经济业务:
(1)5 日财务科王明因业务需要到西安出差,经领导同意预借其差旅费,开出现金支票5 000 元付讫。
(2)11 日财务科王明出差回来,报销差旅费4 800 元,经审核无误同意报销,同时退回现金
D.各项财产物资账面余额与实有数额之间的核对
11.关于记账凭证账务处理程序,下列说法不正确的是()。
A.根据记账凭证逐笔登记总分类账,是最基本的账务处理程序
B.简单明了,易于理解,总分类账可以较详细地反映经济业务的发生情况
C.登记总分类账的工作量较大
D.适用于规模较大、经济业务量较多的单位
12.在科目汇总表账务处理程序下,记账凭证不可以用来()。
A.单位所有的会计档案均不得销毁
B.会计档案在保管期满后,可以直接销毁
C.对于保管期满但尚未结清的债权债务的原始凭证不得销毁
D.国家财政部门销毁会计档案可以不派人监销19.下列各项中,不属于增值税一般纳税人存货成本的是()。
A.商品的买价
B.商品的增值税(取得专用发票)
C.商品的消费税
D.商品的运输费
A.进行备忘登记,但不做任何账务处理

三维时域有限差分法对金SERS基底电磁场增强的初步模拟的开题报告

三维时域有限差分法对金SERS基底电磁场增强的初步模拟的开题报告

三维时域有限差分法对金SERS基底电磁场增强的初步模拟的开题报告一、研究背景和意义金表面增强拉曼光谱(SERS)是近年来快速发展的光学技术,用于检测和识别微量分子。

金基底作为最重要的组成部分之一,对于SERS效应的产生和提高有着重要的影响。

其中,电磁场是SERS效应产生的关键因素之一,能够引起分子的局部电场增强,从而大大提高其敏感度。

因此,研究金基底电磁场增强机理及方式,探究其内部电磁场分布图案对SERS响应产生的影响是十分必要的。

而有限差分法是计算电磁场分布的常用方法之一,它可以通过数值计算得到近似的电磁场分布,从而有效地研究金基底的电磁场增强现象,是该研究的选择方法。

二、研究内容和方法本研究计划通过三维时域有限差分法(FDTD)模拟金基底中的电磁场分布,以探究金基底电磁场增强的机理和规律。

具体而言,本研究将采取以下步骤:1. 建立金基底FDTD模型根据实际金基底的几何尺寸、电磁参数等物理参数建立三维时域有限差分法模型,利用MATLAB等软件对模型进行搭建。

2. 模拟金基底中的电磁场分布通过对电磁场的初始化,定义入射光的频率、极化方向、入射角度和偏振状态,对金基底的电磁场进行模拟计算,得到金基底内电磁场分布。

3. 分析金基底电磁场增强机理和规律通过对金基底内的电磁场分布进行分析,探究不同几何形状、电磁参数和光学周期性结构的金基底电磁场增强机理和规律,为进一步的理论分析和实验研究提供科学依据。

三、预期成果本研究通过FDTD模拟金基底电磁场分布,有望得到金基底中电磁场强度和分布等方面的定量分析结果,探究金基底电磁场增强机理和规律,达到以下预期成果:1. 构建金基底FDTD模型,并得到金基底内电磁场分布;2. 探究不同金基底电磁参数、几何形状和周期性结构对其电磁场强度和分布的影响;3. 为进一步实验和理论分析提供基础和指导,提升金基底电磁场增强现象的理论解释和实验验证能力。

四、进度计划本研究计划在以下时间节点完成相关任务:1. 第一阶段(1个月):文献综述和研究框架的构建;2. 第二阶段(2个月):FDTD模拟模型的建立和参数调整;3. 第三阶段(2个月):FDTD模拟电磁场分布的计算和分析;4. 第四阶段(1个月):结果总结、论文撰写和答辩准备。

一个不确定度模型的三维褶皱与断层重建(中文)

一个不确定度模型的三维褶皱与断层重建(中文)

一个不确定度模型的三维褶皱与断层重建--以高山隧道为例为了改善奥地利与意大利之间的铁路连接,一座隧道,从福尔因斯布鲁克(57公里),正在研究。

设计走廊横切大和强烈的构造部分东阿尔卑斯山脉,由复杂的变质岩和火成岩岩性和多相结构下开发的韧脆性变形条件。

为了建模的子表面地质的地区,表面和子表面地质数据已被集成在一个空间数据库。

利用2种方法对走廊的意大利部分的三维地质模型进行了构建。

第一种是一种比较传统的方法,涉及多个平行和交叉的横截面的重建。

它已与定制开发的脚本,使一个自动工程结构数据使用arcgiss软件实现,在表面和沿钻孔采集到的横截面。

的投影方向可以被控制,并从一个详细的统计分析的方向数据的基础上获得的结构趋势。

其他arcgiss脚本使横切剖面连接网络和帮助保护他们的一致性。

第二方法涉及到一个真正的三维地质模型gocads汇编。

就时间效率和可视化而言,二是更强大的方法。

的基本结构地质假设,然而,在第一种方法适用于类似。

除了3D模型,编写脚本(arcgiss和gocads)已被开发,它允许在表面或沿钻孔观测结构的深度外推的不确定性估计。

这些脚本允许分配的每个投影的结构元素(即,地质边界,断层和剪切带)的参数估计的可靠性。

基本区别的'data-driven”插值和外推'knowledge-based”地质特征的深度进行了讨论和后果的不确定性估计三维地质模型进行评估。

1、介绍和动机三维地质模型是基于现场数据的目的,在预测地质条件的深度,但强烈影响不同来源的不确定性。

这些来源包括整体结构和地层现场的复杂性,需要预测,地形的深度、地质构造区域的态度和岩层的连续性。

传统上,子表面模型是基于网络的横截面,有时使用合理的定量结构地质方法,但更多的往往是设计的手使用地质洞察力和毫米纸。

这种地质模型是很难实现的,因此,在复杂的多聚变形地区,因此而不可靠的。

在这样的网站,没有区域的结构趋势可以被定义为东方的横截面与适当的角关系的基本结构(见拉姆齐和胡贝尔,1987,366,讨论如何东方的横截面与结构趋势)。

基于图形电磁学雷达散射截面计算方法之改进

基于图形电磁学雷达散射截面计算方法之改进

基于图形电磁学雷达散射截面计算方法之改进崔俊伟;杨飏【期刊名称】《电子学报》【年(卷),期】2014(000)012【摘要】图形电磁学(GRaphical Electromagnetic COmputing,GRECO)利用图形加速卡和Z-Buffer技术可较为有效地解决传统电磁计算方法中存在的消隐困难和非可视化难题,是求解高频电大尺寸目标特性最有效的方法之一。

但传统GRECO算法存在着无法精确提取目标法矢信息、计算精度依赖屏幕分辨率和多次反射计算困难等缺点,限制了这种方法的使用。

本文针对GRECO方法就精确提取像素几何信息方法进行了简要改进,将其与基于帧缓存对象(Frame-Buffer Object,FBO)的离屏渲染技术相结合,提出了改进的GRECO算法,克服了传统GRECO算法无法精确提取像素法矢信息和计算精度依赖屏幕分辨率的缺点。

进而,采用AP/PO(Area Projection/Physical Optics)法,并对传统的多次散射面元对判别方法进行了适当改进,实现了对产生多次反射目标的雷达截面计算。

%Based on 3-D graphics hardware accelerator and Z-Buffer technique,graphic electromagnetic computing (GRE-CO)algorithm can efficiently resolve blanking difficulties and visualization problems of traditional electromagnetic calculation pro-cedures .Therefore,GRECO algorithm is considered as one of the most efficient methods to acquire characteristics of high-frequency and electricity large-sizedtarget .However,there are disadvantages of traditional GRECO algorithm,as follows:normal vector of tar-get cannot be extracted accurately,calculationaccuracy is affected by screen resolution greatly and multiple reflections cannot be calculated directly .As a result,traditional GRECO algorithm is limited for this reason and cannot be used in some region widely . The traditional GRECO is improved in this paper,so that the geometric information of pixel can be extracted accurately .Technique of off-screen rendering based on frame buffer object (FBO)is used for improving the algorithm .Then the normal vector of target can be obtained precisely and effectively .Traditional discriminated method of facet pairs is improved by using area projection/physi-cal optics to adapt the computation of RCS multiple scattering .【总页数】6页(P2409-2414)【作者】崔俊伟;杨飏【作者单位】大连理工大学船舶工程学院,辽宁大连 116024;大连理工大学船舶工程学院,辽宁大连 116024【正文语种】中文【中图分类】TN911.23【相关文献】1.图形学在进气道雷达散射截面计算中的应用 [J], 张丽星;吴强;官火梁2.雷达散射截面图形算法加速技巧 [J], 李建周;毛继志;许家栋3.基于图形电磁学的雷达散射截面计算方法改进 [J], 刘立国;张国军;莫锦军;袁乃昌4.基于改进的图形电磁学法的近场RCS计算 [J], 王盟;李燕;程志华;牛立强5.基于改进型CLEAN算法三维成像的雷达散射截面积反演 [J], 任浩田;廖可非因版权原因,仅展示原文概要,查看原文内容请购买。

柳纱科技 AAAnalyst 400 原子吸收光谱仪器说明书

柳纱科技 AAAnalyst 400 原子吸收光谱仪器说明书

Atomic Absorption AAnalyst 400 AASpectrometerControl and Data SystemUser Interface Complete PC control of all functions of the AAnalyst™ 400 using WinLab32™ for AA software. WinLab32 for AA includes an innovative user interface that makes the software easy to learn and use, including a clear graphical design, task-oriented organization of the windows, an understandablevocabulary, extensive tool tips in multiple languages, simple data displays and Wizards for the simplification of many tasks.WinLab32 is fully multitasking, allowing the analyst to report analytical results, view data or add priority samples without interrupting the analysis in progress.Using WinLab32 software, setup is flexible and easy. Standard operating conditions for flame, graphite furnace and FIAS techniques are included.Auto Analysis Control links methods for each technique with a sample-information file. The sample list can be created by third-party software orLIMS and downloaded to the system.WinLab32 software provides many tools to increase lab productivity. With WinLab32 Offline, method and sample-information files can be createdand data reviewed or reprocessed without interrupting the current analysis.WinLab32 software provides extensive QC protocols to meet internal and regulatory requirements. Data Reprocessing allows changes to manymethod and sample-information parameters after data collection and the recalculation of the results using the new parameters. With DataReprocessing, the raw data are never altered, thereby ensuring data integrity is maintained.WinLab32 Reporter provides the ability to generate post-run reports in a variety of formats. The Export feature of Data Manager can be used toexport results as comma-delimited ASCII files for compatibility with commercial third-party programs such as Microsoft® Excel®, Access® and Word.21 CFR Part 11 An optional WinLab32 Enhanced Security™ package is available for labs needing to be compliant with 21 CFR Part 11 regulations.HardwareSystem True double-beam echelle optical system. Front surfaced, reflecting optics with protective coating. Deuterium background corrector and two built-in EDL power supplies.Optical System Echelle monochromator. Focal length: 300 mm. Grating: 36 x 185 mm area, 79 lines/mm, blaze angle 76˚. Fused-quartz prism: 95 x 40 mm, 60°.Wavelength: 189-900 nm. Spectral bandpass: 0.15 nm at 200 nm. Reciprocal linear dispersion: 2.4 nm/mm. The photometer optics are covered toprotect against dust and corrosive vapors. For maximum protection, the optical system can be purged with an inert gas.Detector High-efficiency, segmented solid-state detector.Light Sources Hollow cathode or electrodeless discharge lamps (EDLs). EDLs provide much higher light output and longer lifetime when compared to conventional hollow cathode lamps. Lamp elements, recommended operating currents and slit selection are automatically recognized and set when usingPerkinElmer® Lumina™ series AA lamps. Lamp alignment is completely automatic with the four-lamp turret.E-box All electronics are located in a single user-replaceable module that the operator can easily replace without requiring a service visit.For a complete listing of our global offices, visit /ContactUsCopyright ©2004-2010, PerkinElmer, Inc. All rights reserved. PerkinElmer ® is a registered trademark of PerkinElmer, Inc. All other trademarks are the property of their respective owners. PerkinElmer reserves the right to change this document at any time without notice and disclaims liability for editorial, pictorial or typographical errors. 006675E_01PerkinElmer, Inc. 940 Winter Street Waltham, MA 02451 USA P: (800) 762-4000 or (+1) Gas Controls and Burner SystemFlame GasFully automated gas box with computer-controlled oxidant selection, automatic gas sequencing, oxidant and fuel monitoring and control.ControlSoftware-actuated ignition with air/acetylene. Acetylene flow is automatically adjusted when switching to or from nitrous-oxide/acetylene operation.Flame SafetyFully interlocked operation prevents ignition if the proper burner head, the nebulizer, end cap or burner drain system are not correctly Featuresinstalled, the level of the liquid in the drain vessel is incorrect, or gas pressures are too low. Interlocks will automatically shut d own thegases if a flame is not detected. The flame is automatically and safely extinguished in the event of a power failure or when theemergency flame-off button is used.Burner SystemAn inert-polymer mixing chamber provides superior analysis of corrosive and high-solid matrices. The spray chamber is manufactured froma high-strength composite, eliminating the need for pressure-relief devices. The high-precision inert nebulizer maximizes stability andsensitivity. A 10-cm single-slot solid titanium burner head for air/acetylene operation is included. Optional burner heads include: 5-cmnitrous-oxide/acetylene, 10-cm three-slot air/acetylene and 5-cm single-slot air/acetylene.Sample Area 25 cm wide x 25 cm deep sample compartment for easy access to burner components.Accessories for the AAnalyst 400AutosamplersFlame autosamplers automate standard and sample introductions for instrument calibration and sample analysis, extending thespectrometer’s capabilities to those of a fully automated analytical workstation. Sample Dilution The AutoPrep 50 sample-dilution system provides an optimized tool for truly automated flame AA. With automatic, intelligent on-line dilution capabilities, the AutoPrep 50 eliminates the time-consuming, manual, error-prone portion of your flame AA analyses.Mercury/Hydride For the analysis of mercury or hydride-forming elements, an optional automated flow injection system or a manual mercury/hydride System system can be added. Flow Injection Atomic Spectroscopy (FIAS) combines the advantages of mercury/hydride AA with those of the flow injection, enabling mercury/hydride AA procedures to be truly automated.System SpecificationsDimensions70 x 65 (0.46 m 2) x 65 cm (W x D x H)Weight49 kg Power100-230 V (±10%), 50/60 Hz (±1%), 300 VA (maximum)TechnicalClassified as a laboratory instrument. Complies with the applicable European Union directives and standards for safety and electro-magnetic compatibility for CE Marking, the safety requirements for Canada and the United States for CSA/NRTL certification and the FCC requirements for radio-frequency emissions. The instrument was developed and produced in compliance with ISO 9001.Environmental Dust-free, free of vibrations, ambient temperatures: +15 ˚C to +35 ˚C with a change rate of a maximum 3 ˚C per hour. Relative humidity: 20% to 80% non-condensing.。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

3-D GEOMETRY ENHANCEMENT BY CONTOUR OPTIMIZATIONIN TURNTABLE SEQUENCESPeter EisertFraunhofer Institute for Telecommunications,Heinrich-Hertz-InstituteEinsteinufer37,D-10587Berlin,GermanyEmail:eisert@hhi.fhg.deABSTRACTA method for the enhancement of geometry accuracy in shape-from-shading frameworks is presented.For the particular case of turntable scenarios,an optimization scheme is presented that minimizes silhouette deviations which correspond to shape errors. Only three unknown parameters have to be optimized leading to a robust and relatively fast framework.In spite of the small number of parameters,experiments have shown that the silhouette error can be reduced by a factor of more than10even after an already quite accurate camera calibration step.The quality of an additional texture mapping can also be drastically improved making the pro-posed scheme applicable as a preprocessing step in many different 3-D multimedia applications.1.INTRODUCTIONPhotorealistic3-D computer models of real objects and scenes are key components in many3-D multimedia systems.The quality of these models often has a large impact on the acceptability of applications like virtual walk-throughs(e.g.,city guides or vir-tual museums),caves,computer games,product presentations in e-commerce,or other virtual reality systems.Although the ren-dering of3-D computer scenes can be performed in real-time even on hand-held devices,the creation and design of3-D models with high quality is still time consuming and thus expensive.This has motivated the investigation of a large number of methods for the automatic acquisition of textured3-D geometry models from mul-tiple views of the object.One class of approaches often used due to its robustness and simplicity is called shape-from-silhouette or shape-from-contours [1,2,3].For shape estimation,multiple views of an object are captured and in the images,the object is segmented from the back-ground.The contour in the image plane forms–together with the focal point of the camera–a3-D viewing cone that contains the entire object.Intersection of all viewing cones leads to an esti-mate for the object geometry.In order to circumvent the limitation of reconstructing the visual hull only,many extensions have been proposed that incorporate also color information as,e.g.,voxel col-oring or space carving[4,5].Once the images are reliably segmented,shape-from-silhouette methods work very robust,are hardly affected by il-lumination variations,and are computationally quite efficient. With today’s computers,real-time implementations now become available,enabling new applications like3-D communication or 3-D video recording[6,7,8,9].In these scenarios,data is usually captured with multiple cameras in a cave,where the background can easily be controlled for simplesegmentation.Fig.1.Turntable setup.However,the accuracy of the reconstructed geometry is con-siderably affected by the knowledge of the true camera parameters which requires an accurate camera calibration.Deviations from the correct values can lead to smaller3-D models,since valid ob-ject parts might be cut off during viewing cone intersection.As a result,the silhouette of the reconstructed object is always smaller or equal to the true contour.However,this unwanted effect can be used to refine the camera parameters again.In[10],Grattarola uses pairs of images and projects the viewing cone of one image con-tour into the other image plane,computing silhouette mismatches. Minimization of these mismatches over all camera parameters op-timizes the calibration.Similarly,Niem[11]minimizes the devi-ation of the back-projected silhouette of the reconstructed object to the true silhouette in the camera images.Both approaches have in common,that a non-linear optimization in a high dimensional space has to be performed.In this paper,we focus on the turntable3-D acquisition sce-nario as shown in Fig.1.Instead of placing multiple cameras around the object,a single camera is used capturing images while the object slowly rotates.The rotation angle between two shots can usually be accurately controlled,whereas the position of the camera relative to the turntable is in general unknown and requires camera calibration.In the following,we present a method for opti-mizing the camera’s position from silhouette mismatches.In con-trast to the work in[11],the circular motion adds severe constraints which increase robustness and result in a very low dimensional pa-rameter space that has to be searched.Such constraints are also utilized in[12],where extrinsic and intrinsic camera parameters are derived from feature point correspondences instead of silhou-view i object contourof camera imagecontour ofdifference area Abetween real andreconstructedsegmentation maskFig.2.Contours of the real and reconstructed object for an incor-rectly calibrated camera.ette.Experiments show that the proposed silhouette optimization converges even for large deviations making camera calibration less important or even unnecessary.Due to the improved camera pa-rameters,the accuracy of the reconstructed3-D models is highly increased by this preprocessing step.2.CONTOUR OPTIMIZATION FOR TURNTABLESETUPSWith perfectly calibrated cameras,shape-from-silhouette recon-structs the visual hull of an object by volume intersection of the viewing frustrums defined by the segmented contours.If the pa-rameters of the real cameras are not accurately known,an erro-neous object shape is reconstructed.Due to the intersection pro-cess,the reconstructed object contour always lies within the seg-mentation mask of the camera image as shown in Fig.2.The area difference A i between correct and reconstructed ob-ject masks for view i is directly related to the amount of camera parameter deviation.Therefore,it can be used to optimize the M extrinsic and intrinsic camera parameters p j by minimizing the area difference.Following[11],we define the cost function to be minimized as the sum of all area deviations A i over all available N viewsc(p0,...,p M−1)=1dRxdRz e r r o rFig.3.Silhouette error as a function of deviations in the rotational angles R x and R z .Fig.4.Geometry enhancement for a Moche vase.Left:two frames of the original sequence;middle:contour error after camera cali-bration;right:contour error after proposed optimization.iments are performed on natural image sequences.For all se-quences,every 10degrees of turntable rotation an image is cap-tured leading to 36views with a resolution of 432x576pixels.Based on background color information,the object is segmented from the background leading to 36binary object masks.The inter-nal camera parameters are determined once for the camera with a model-based camera calibration technique [13].With these initial values,the optimization is started.Fig.4shows the results for two frames of the sequence Moche contain-ing a vase from the Peruvian Moche culture.The left images de-pict the original frames and the middle ones the silhouette error A ifinal error Moche calib.271Moche 264Tree1467∆t x ∆R z 1.5mm 0.05o −3.4mm 0.95o −5.7mm1.42oTable 2.Estimated parameter changes for the unknown parame-ters R x ,R z ,and t x .for the corresponding view after reconstruction with the parame-ters from the camera calibration step.The right images show the deviations between real and reconstructed contours which are sig-nificantly reduced by the proposed algorithm.The cost function (1)describing parameter inaccuracy is reduced by a factor of more than 10as shown in the first row of Table 1.The resulting camera parameter updates,however,are rela-tively small as depicted in the first row of Table 2.The initial position error of the turntable axis is about 1.5mm and their ori-entation is determined with an accuracy of less than 1degree.This indicates that even small parameter deviations from the camera calibration can lead to rather large geometry errors and optimizing these values can significantly enhance the accuracy of the recon-structed geometry.In a second experiment,the initial values for R x ,R z ,and t x are set to zero and no information from the previous camera cali-bration is exploited.The algorithm,however,still converges as in-dicated in Fig.5.Although the initial silhouettes does not match,the final silhouette error is in the same range as indicated in the middle row of Table 1.The correction of the rotation angle R x is more than 20degrees (see Table 2).Similar experiments are performed for a second sequence Tree that has a much more sophisticated contour.The algorithm still converges even for unknown initial values of the camera position and orientation.Similarly,silhouette errors befor and after opti-mization are depicted in Fig.6.The improvement of geometry accuracy resulting from the proposed method also affects the quality of the texture map esti-mated from the views.Incorrect shape prevents an accurate stitch-ing of multiple texture parts from different views leading to blurred images.The left side of Fig.7shows the quality of the final tex-tured 3-D model of the vase after reconstruction with parameters from the initial camera calibration.Although the refined values are quite similar,the sharpness of the texture map after the pro-posed geometry enhancement is drastically improved.This indi-cates that an additional refinement of camera parameters for shape-from-silhouette reconstruction can help to improve the quality of the results.4.CONCLUSIONSIn this paper,we have presented a method for the refinement of camera parameters for turntable scenarios.The contour devia-tion between reconstructed object and segmentation mask of the camera frames is minimized in a shape-from-silhouette framework leading to an optimal set of camera parameters.For the particular case of turntable acquisition,we have shown that only 3parame-ters are sufficient to calibrate the external camera parameters.In-Fig.5.Silhouette optimization with no initial pose information from calibration.Middle:Large initial silhouette error;right:sig-nificantly reduced deviations afteroptimization.Fig.6.Geometry deviations for the sequence Tree .Left:original camera frames;middle:initial error;right:silhouette error after refinement.ternal camera parameter optimization can be added similarly with-out changing the algorithm.The low-dimensional parameter space ensures high robustness and reduced complexity.Experiments have shown that the quality of 3-D models can be drastically im-proved with the proposed optimizationtechnique.Fig.7.Magnification of a reconstructed object without (left)and with (right)camera position refinement.The improved geometry results in a much sharper texture.5.ACKNOWLEDGMENTThe work presented in this paper has been developed with the sup-port of the European Network of Excellence VISNET (IST Con-tract 506946).6.REFERENCES[1]R.Szeliski,“Rapid octree construction from image sequences,”inComputer Vision,Graphics and Image Processing ,Jul.1993,pp.23–32.[2] urentini,“The visual hull concept for silhouette-based imageunderstanding,”IEEE Transactions on Pattern Analysis and Machine Intelligence ,vol.16,no.2,pp.150–162,Feb.1994.[3]W.Niem and J.Wingberm¨u hle,“Automatic reconstruction of 3Dobjects using a mobile monoscopic camera,”in Proceedings of the International Conference on Recent Advances in 3D Imaging and Modelling ,Ottawa,Canada,May 1997.[4]S.M.Seitz and C.R.Dyer,“Photorealistic scene reconstruction byvoxel coloring,”in puter Vision and Pattern Recognition ,Puerto Rico,1997,pp.1067–1073.[5]P.Eisert,E.Steinbach,and B.Girod,“Multi-hypothesis,volumet-ric reconstruction of 3-D objects from multiple calibrated camera views,”in Proc.International Conference on Acoustics,Speech,and Signal Processing (ICASSP),Phoenix,Mar.1999,pp.3509–3512.[6]X.Wu and T.Matsuyama,“Real-time active 3D shape reconstructionfor 3D video,”in Proc.of 3rd International Symposium on Image and Signal Processing and Analysis ,Rome,Italy,Sep.2003,pp.186–191.[7] B.Goldl¨u cke and M.Magnor,“Real-time,free-viewpoint video ren-dering from volumetric geometry,”in Proc.Visual Computation andImage Processing (VCIP),Lugano,Switzerland,Jul.2003.[8]O.Grau,“Studio production system for dynamic 3D content,”inProc.Visual Computation and Image Processing (VCIP),Lugano,Switzerland,Jul.2003.[9]M.Gross,“blue-c:A spatially immersive display and 3D video portalfor telepresence,”in puter Graphics (SIGGRAPH),San Diego,USA,Jul.2003,pp.819–827.[10] A.A.Grattarola,“V olumetric reconstruction from object silhouettes:A regularization procedure,”Signal Processing ,vol.27,pp.27–35,1992.[11]W.Niem,Automatische Rekonstruktion starrer dreidimensionalerObjekte aus Kamerabildern ,Ph.D.thesis,University of Hanover,Hanover,Germany,1998.[12] A.W.Fitzgibbon,G.Cross,and A.Zisserman,“Automatic 3D modelconstruction for turn-table sequences,”in Proc.ECCV 98Workshop on 3D Structure from Multiple Images in Large-Scale Environments ,Freiburg,June 1998,pp.154–170.[13]P.Eisert,“Model-based camera calibration using analysis by syn-thesis techniques,”in Proc.Vision,Modeling,and Visualization VMV’02,Erlangen,Germany,Nov.2002,pp.307–314.。

相关文档
最新文档