Indoor nuclear environment multi-information double-layer hybrid map modeling method

Abstract

The invention discloses an indoor nuclear environment multi-information double-layer hybrid map modeling method. The method comprises the following steps: 1, acquiring a grid map containing grid elevation values; 2, establishing an IOP model to obtain a geometric grid map; 3, acquiring a radiation field grid map based on MCNP or a difference algorithm; 4, constructing a multi-information double-layer hybrid map: carrying out binary coding on the trafficability of each grid in the geometric grid map; coding the average radiation dose rate of each grid in the radiation field grid map by adopting a real number; and after the grid information is coded, the grid information is stored in a group number form to form a mapping relation. The invention provides a modeling method of a multi-information double-layer hybrid map for path planning of a mobile robot in a radioactive indoor environment, and the map is a grid map fusing trafficability information and environment radiation field information, and can describe multi-physical field coupling environment characteristics in an indoor nuclear environment. And a reliable environment map is provided for subsequent path planning of the nuclear robot.

Description

Technical Field

[0001] The invention relates to the field of map building algorithms, and in particular to a multi-information double-layer hybrid map modeling method for an indoor core environment. Background Art

[0002] The nuclear power plant indoor nuclear environment is a typical application scenario for mobile robots. In this scenario, mobile robots are often used for tasks such as component transportation, facility disassembly and assembly, and environmental detection. The prerequisite for performing the above tasks is to obtain the environmental map model, which is the basis for realizing the path planning and behavior decision-making of mobile robots.

[0003] When the nuclear environment in a nuclear power plant is unknown, mobile robots often use SLAM technology (based on lidar and cameras) to explore and model the working environment. The information in the map usually comes from obstacle data within the small-angle perception field of view. On the one hand, the 2D geometric grid map constructed in this way does not take into account the obstacle-crossing ability of the mobile robot, so it is easy to demarcate the traversable area as an inaccessible area, thereby reducing the drivable range of the mobile robot. On the other hand, the 2D geometric grid map constructed in this way does not integrate geometric information and radiation information, making it difficult to accurately describe the coupled environmental characteristics in the nuclear power plant, reducing the reliability of the mobile robot's path planning. Summary of the invention

[0004] The purpose of the present invention is to overcome the shortcomings of the prior art and provide a multi-information double-layer hybrid map modeling method for an indoor nuclear environment. It solves the problem that the 2D geometric grid map currently constructed based on SLAM technology does not combine the obstacle-crossing ability of the mobile robot itself to determine the passable area, and is easy to demarcate the passable area as an impassable area. It also solves the problem that the 2D geometric grid map currently constructed based on SLAM technology does not integrate geometric information and radiation information, and is difficult to accurately describe the coupled environmental characteristics within the nuclear power plant.

[0005] The technical solution of the present invention is: a method for modeling a multi-information double-layer hybrid map of an indoor nuclear environment, the steps of which are as follows:

[0006] S01, obtaining a grid map containing grid elevation values;

[0007] S02, establishing an IOP model to obtain a geometric grid map;

[0008] S03, obtaining a radiation field grid map based on MCNP or difference algorithm;

[0009] S04, multi-information two-layer hybrid map construction: the accessibility of each grid in the geometric grid map is binary-encoded; the average radiation dose rate of each grid in the radiation field grid map is encoded with real numbers; after the grid information is encoded, it is stored in the form of group numbers to form a mapping relationship.

[0010] A further technical solution of the present invention is: Step S01 is briefly described as follows:

[0011] A, obtain the 3D point cloud of the indoor nuclear environment and transform it into a rectangular coordinate system;

[0012] B. Eliminate low-precision data and invalid data in the 3D point cloud. Invalid data refers to the 3D point cloud where the laser beam is blocked by the mobile robot, and low-precision data refers to the 3D point cloud after the detection distance is 20m.

[0013] C, matching the virtual grid divided on the ground with the 3D point cloud after data removal to obtain a grid map containing grid elevation value information;

[0014] The S02 step is briefly described as follows:

[0015] A, calculate the maximum height of the steps that the mobile robot can climb;

[0016] B, calculate the maximum inclination angle that the mobile robot can pass;

[0017] C, calculate the height difference of the factors affecting traffic capacity;

[0018] D, calculate the longitudinal slope and transverse slope influencing factors of traffic capacity;

[0019] E, calculate the minimum ground clearance among the factors affecting traffic capacity;

[0020] F, IOP cost calculation of target grid: First, the passability of the target grid area is determined based on the values ​​of chassis ground clearance and minimum ground clearance; if it is passable, the IOP cost value of the target grid is calculated by combining the height range and the lateral slope and longitudinal slope; finally, the IOP threshold is set and the calculated IOP cost value is binarized to distinguish between impassable and passable areas;

[0021] The AB steps in this step are in no particular order, and the CE steps in this step are in no particular order;

[0022] The S03 step is briefly described as follows:

[0023] When the source information is known, the radiation field is calculated based on MCNP simulation, and the calculation results are fused with the grid map; when the source information is unknown, the radiation field is calculated using the difference algorithm, and the calculation results are fused with the grid map.

[0024] A further technical solution of the present invention is: the method is based on a mobile robot that uses double-sided tracks and is equipped with a laser radar, and a driven wheel and a driving wheel are respectively provided at the front and rear ends of the inner side of each track.

[0025] A further technical solution of the present invention is: Step S01 is described in detail as follows:

[0026] A. Obtain the 3D point cloud of the indoor nuclear environment and perform coordinate conversion: Control the mobile robot to move in the indoor nuclear environment and obtain complete environmental information through the laser radar; the environmental information is a frame-by-frame 3D point cloud in a spherical coordinate system with the laser radar as the origin; in the spherical coordinate system, the X axis is the forward direction of the mobile robot, the Y axis is the left and right direction of the mobile robot, the Z axis is the upward direction perpendicular to the XLY plane, and L is the origin of the laser radar; the target point P in the 3D point cloud is represented by P(D pl ,θ,φ), where D pl is the distance between the target point P and the laser radar origin L, θ is the scanning angle, is the rotation angle; refer to formula 1 to transform the 3D point cloud from the spherical coordinate system to the rectangular coordinate system;

[0027] Formula 1:

[0028] Where x, y, and z are the x-, y-, and z-axis coordinates of the target point P in the rectangular coordinate system, respectively;

[0029] B. Eliminate low-precision data and invalid data in the 3D point cloud: define the 3D point cloud after the detection distance R is 20m as low-precision data, and define the 3D point cloud where the laser beam is blocked by the mobile robot as invalid data; refer to formula 2 to eliminate low-precision data and invalid data from each frame of the 3D point cloud;

[0030] Formula 2:

[0031]

[0032] In the formula, st represents the constraint condition, B represents the width of the mobile robot, and L represents the length of the mobile robot;

[0033] C. Divide the ground into virtual grids and match the 3D point cloud with them: Divide the virtual grid using the ground as the reference plane, and set the side length of the grid to be consistent with the width B of the mobile robot; project the 3D point cloud obtained in step B into the virtual grid, and calculate the weighted average height of all discrete points in each grid as the elevation value of each grid, thus obtaining a grid map containing grid elevation value information.

[0034] A further technical solution of the present invention is: Step S02 is described in detail as follows:

[0035] A. Establish a climbing step model: The action of the mobile robot climbing the steps is divided into three processes: Ⅰ. The driven wheel of the mobile robot contacts the edge of the step. Driven by the mobile robot's own power system, the track of the mobile robot close to the step begins to rise to perform obstacle crossing; Ⅱ. The center of gravity G of the mobile robot moves forward and rises, and the angle α between the track and the ground gradually increases. When the center of gravity G of the mobile robot is directly above the edge of the step, it reaches the critical state of climbing the step; Ⅲ. The center of gravity G of the mobile robot passes the edge of the step, and the mobile robot returns to a horizontal state to complete the obstacle crossing;

[0036] The critical state analysis of the mobile robot during the climbing process is carried out and formula 3 is established;

[0037] Formula 3: h0=(l1-htanα)sinα

[0038] Where h0 is the maximum height of the step that the mobile robot can climb, l1 is the distance between the center of gravity of the mobile robot and the vertical plane where the center line of the driving wheel is located, h is the height from the center of gravity of the mobile robot to the ground below the step, and α is the angle between the crawler of the mobile robot and the ground below the step;

[0039] B. Establish a slope driving model: Slope driving includes three situations: longitudinal uphill, longitudinal downhill and lateral driving. As the slope angle increases, the risk of the mobile robot tipping over backwards gradually increases when it is longitudinally uphill, the risk of the mobile robot tipping over forward gradually increases when it is longitudinally downhill, and the risk of the mobile robot tipping over sideways gradually increases when it is lateral driving. The maximum slope angle θ of the mobile robot longitudinal uphill is obtained by analysis. max1 As shown in Formula 4, the maximum slope angle θ of the longitudinal downhill slope max2 As shown in Formula 5, the maximum roll angle θ for lateral travel max3 As shown in formula 6;

[0040] Formula 4:

[0041] Formula 5:

[0042] Formula 6:

[0043] Where l2 is the distance between the center of the mobile robot and the vertical plane where the center line of the driven wheel of the mobile robot is located; B is the width of the mobile robot, and b is the track width of the mobile robot;

[0044] Combining formulas 4-6, we can get the maximum inclination angle θ that the mobile robot can pass. max As shown in formula 7;

[0045] Formula 7: θ max =min{θmax1 ,θ max2 ,θ max3};

[0046] In the formula, θ max is the maximum inclination angle that the mobile robot can pass;

[0047] C, calculate the height difference among the factors affecting traffic capacity:

[0048] The center of gravity of the mobile robot is projected onto the grid below. The grid into which the projection is projected is the central grid. The center of the central grid is the center of the circle. The maximum height change value of a total of 9 grids in a circular area with a grid side length as the radius is the height range, which is calculated according to Formula 8.

[0049] Formula 8:

[0050] Among them, δH is the height difference value, H x,y Represents the elevation value of the grid at the xth row and yth column, x = (i-1, i, i+1), y = (j-1, j, j+1);

[0051] D. Calculate the longitudinal slope and transverse slope influencing factors of traffic capacity: define the central grid and the nearest 8 grids around it as the analysis area, use the data of the nearest 8 grids around the central grid where the mobile robot is located to estimate the slope value of the central grid. The global coordinate system coordinates of the grid in the analysis area are (x i ,y i , z i ), the subscript i is the grid number, z i is the elevation value of the grid numbered i; based on the least squares method, the fitting surface of the analysis area is obtained as shown in formula 9;

[0052] Formula 9: Z = Ix i +Jy i +K

[0053] Among them, Z is the fitting surface function of the analysis area, and I, J, and K are all surface constant coefficients to be determined;

[0054] The point cloud coordinates in the analysis area are transformed from the rectangular coordinate system to the global coordinate system, and the global coordinate system coordinates of the point cloud in the analysis area or the global coordinate system coordinates of all grids in the analysis area are substituted into Formula 9 to obtain the cumulative error formula, as shown in Formula 10;

[0055] Formula 10:

[0056] Wherein, the superscript n is the maximum number of grids or the maximum number of point clouds in the analysis area, and Ω is the cumulative error;

[0057] Taking the minimum cumulative error Ω as the expectation, we take the partial derivative of I, J, and K in formula 10 to obtain formula 11;

[0058] Formula 11:

[0059] Solve formula 11 to obtain the values ​​of the surface constant coefficients I, J, and K, and then obtain the normal vector of the fitting surface in the analysis area In addition, the normal vector of the horizontal plane is Therefore, the slope of the i-th grid in the analysis area is the angle between the two normal vectors, as shown in formula 12. Formula 12 is then split into longitudinal slope and transverse slope to obtain formula 13;

[0060] Formula 12:

[0061] Formula 13:

[0062] in, is the set of longitudinal slope and transverse slope of the i-th grid in the analysis area, is the longitudinal slope of the i-th grid in the analysis area, is the transverse slope of the i-th grid in the analysis area;

[0063] E, calculate the minimum ground clearance among the factors affecting the traffic capacity: define the minimum vertical distance between the horizontal plane where the chassis of the mobile robot is located and the horizontal plane where the elevation value of the target grid is located when the mobile robot successfully passes through the target grid area, that is, the minimum ground clearance; set the center of gravity G of the mobile robot to be vertically opposite to the center point O of the target grid, and calculate the horizontal position coordinates of the contact points of the left driven wheel T1, the left driving wheel T2, the right driven wheel T3 and the right driving wheel T4 with the ground according to the size parameters of the mobile robot and the slope of the target grid, as shown in Formula 14;

[0064] Formula 14:

[0065] In the formula, the coordinates of the four wheels are all in the global coordinate system, B is the width of the mobile robot, b is the track width of the mobile robot, θ is the forward direction of the mobile robot in the global coordinate system, and θ v is the longitudinal slope of the target grid, θ c is the lateral slope of the target grid, l1 is the distance between the center of gravity of the mobile robot and the vertical plane where the center line of the driving wheel is located, l2 is the distance between the center of gravity of the mobile robot and the vertical plane where the center line of the driven wheel is located, x i ,y i is the horizontal position coordinate of the center point of the target grid in the global coordinate system;

[0066] Any three of the left driven wheel T1, the left driving wheel T2, the right driven wheel T3 and the right driving wheel T4 constitute a contact plane, and a plane set Ψ and a plane set Ω parallel to the plane set Ψ are constructed based on this; the plane set Ψ includes four planes formed by assembling any three points of T1, T2, T3 and T4, wherein the horizontal position coordinates x, y of the three points in each plane are obtained by referring to formula 14, wherein the height coordinates z of the three points in each plane are calculated by using the IDW interpolation algorithm; correspondingly, the plane set Ω also includes four planes corresponding to the four planes in the plane set Ψ; finally, the minimum ground clearance calculation formula is constructed in combination with the plane distance formula, as shown in formula 15;

[0067] Formula 15:

[0068] In the formula, H grad is the minimum ground clearance, R i is the constant term coefficient of the plane set Ψ, S i is the constant term coefficient of the plane set Ω, and the subscript i is the common number of a plane in the plane set Ψ and the corresponding plane in the plane set Ω;

[0069] F, IOP cost calculation of the target grid: Assuming that the mobile robot is in the target grid area, before doing IOP cost calculation, based on H grad With H stand The ratio of is used for pre-pass determination, refer to formula 16;

[0070] Formula 16:

[0071] In the formula, H grad is the chassis ground clearance, Θ(H grad ) is the passability determination factor. If the value is less than 1, the target grid area is determined to be a passable area. If the value is greater than 1, the target grid area is determined to be impassable and the target grid is assigned a value of 1.

[0072] If the result of the above traffic determination is a passable area, the height difference, transverse slope and longitudinal slope in the traffic capacity influencing factors are used as the IOP cost for calculation, and the piecewise function is defined as shown in Formula 17;

[0073] Formula 17:

[0074] In the formula, represents one of the three factors affecting traffic capacity; when i=1, Represents the height difference; when i=2, represents the longitudinal slope. When i=3, represents the transverse slope; maxRepresents the limit value of the factors affecting traffic capacity. When i=1, is the maximum height h0 of the mobile robot climbing the steps, when i = 2, is the maximum slope angle θ of the mobile robot going uphill longitudinally max1 , when i=3, is the maximum roll angle θ of the mobile robot when traveling laterally max3 , is the weight factor of the corresponding capacity influencing factor; when i=1, Represents the weight factor of the height range. When i=2, Represents the weight factor of the longitudinal slope. When i=3, represents the weight factor of the transverse slope;

[0075] Combining Formula 16 and Formula 17, the IOP cost calculation function can be obtained, as shown in Formula 18;

[0076] Formula 18:

[0077]

[0078] In the formula, represents the IOP cost calculation function, st represents the constraint condition, is the weight factor of the height extreme, is the weight factor of the longitudinal slope, is the weight factor of the transverse slope, IOP1(h) is the IOP cost calculation value of the height extreme difference, IOP2(θ1) is the IOP cost calculation value of the longitudinal slope, and IOP3(θ2) is the IOP cost calculation value of the transverse slope;

[0079] The IOP threshold is set, and the IOP cost value calculated by formula 18 is binarized using formula 19 to realize the binarization of the target grid; the binarization result of the target grid is 0, indicating a passable area, and the binarization result of the target grid is 1, indicating an impassable area;

[0080] Formula 19:

[0081] Where Gmap(i) represents the geometric cost value of the grid numbered i. To customize the IOP threshold.

[0082] A further technical solution of the present invention is: Step S03 is described in detail as follows:

[0083] When the source information is known, the specific operations are as follows:

[0084] A. Construction of 3D model of indoor nuclear environment: Establish a 3D model of the indoor nuclear environment corresponding to the grid map, and determine the material properties of each geometric body in the 3D model according to the actual situation on site, which will be used for material definition in the INP file in the MCNP calculation code;

[0085] B. MCNP simulation calculation of radiation dose rate: According to the source information and environmental characteristics, write MCNP calculation code to perform simulation calculation to obtain the radiation dose rate at various locations in the 3D model;

[0086] C. Count the average radiation dose rate of each grid: Based on the calculated radiation dose rate, count the average radiation dose rate of all grids in the grid map within the range of height 0-a, so as to obtain the radiation field grid map;

[0087] The specific operations when the source item information is unknown are as follows:

[0088] A. Construction of 3D model of indoor nuclear environment: Establish a 3D model of the indoor nuclear environment corresponding to the grid map, and determine the material properties of each geometric body in the 3D model according to the actual situation on site, which will be used for material definition in the INP file in the MCNP calculation code;

[0089] B. Obtaining the reference radiation field: Enter the above-mentioned indoor nuclear environment to perform random multi-point radiation data sampling. Based on the sampled data and environmental characteristics, use MCNP to simulate and calculate the radiation field to obtain the reference radiation field for reconstructing the radiation field;

[0090] C. Reconstruct the radiation field based on the difference algorithm: extract some grids in the reference radiation field as sampling points, re-enter the above-mentioned indoor nuclear environment to sample radiation data at the extracted grid positions, use the sampled data as the initial data set, program the interpolation function, and reconstruct the radiation field through the interpolation program to obtain the radiation dose rate at each location in the 3D model;

[0091] D. Count the average radiation dose rate of each grid: Based on the reconstructed radiation dose rate, count the average radiation dose rate of all grids in the grid map within the height range of 0-a to obtain the radiation field grid map.

[0092] A further technical solution of the present invention is: Step S04 is described in detail as follows:

[0093] A. Encode each grid information in the grid map: Binary encode the accessibility information of the grid obtained in step S02, where the grid is not accessible and is coded as 1, and the grid is accessible and is coded as 0; the average radiation dose rate information of the grid obtained in step S03 is encoded using real numbers, that is, the average radiation dose rate corresponding to each grid is recorded in the grid; after the grid information is encoded, it is stored in the form of groups to form a mapping relationship, thereby obtaining a multi-information double-layer hybrid map.

[0094] A further technical solution of the present invention is that: step S04 further includes a sub-step B following sub-step A;

[0095] B. Perform grid identification on the multi-information double-layer hybrid map: Use a hybrid identification method that combines the rectangular coordinate system method and the serial number method to equivalently represent the grid; Let the rectangular coordinate of any grid k be (i, j), where 0 < i < M and 0 < j < N, M is the number of grid columns in the x-axis direction of the rectangular coordinate system, and N is the number of grid rows in the y-axis direction of the rectangular coordinate system. Then the grid k is represented as k(i, j). If the corresponding number of grid k is n, then the grid k is equivalently represented in the form of k(n); Define the grid numbering to start from the first grid in the upper left corner of the rectangular coordinate system, and determine the serial number in the order from left to right and top to bottom. The serial number starts from 1; According to the above definition, the mapping relationship from coordinate (i, j) to number n is shown in formula 20;

[0096] Formula 20:

[0097] In the formula, mod is the remainder function, that is, the remainder after the division operation of two numerical expressions; The ceil function is used for rounding up, that is, to obtain an integer value greater than or equal to the parameter.

[0098] A further technical solution of the present invention is that: step S03 further includes a sub-step D following sub-step C of the known part of the source item information;

[0099] D. Visualize the dose rate in the grid map: Divide the dose rates in each grid of the grid map into low radiation dose rate values, medium radiation dose rate values, and high radiation dose rate values from high to low, and define the low radiation dose rate values, medium radiation dose rate values, and high radiation dose rate values with the limit value of the radiation dose rate tolerated by the mobile robot; Specify different colors for the low radiation dose rate values, medium radiation dose rate values, and high radiation dose rate values, and finally fill the corresponding specified colors in each grid of the grid map to achieve the visualization of the low radiation dose rate values, medium radiation dose rate values, and high radiation dose rate values.

[0100] A further technical solution of the present invention is that: step S03 further includes a sub-step E following sub-step D of the unknown part of the source item information;

[0101] E. Visualize the dose rate in the grid map: divide the dose rate in each grid in the grid map into low radiation dose rate value, medium radiation dose rate value and high radiation dose rate value from high to low, and define the low radiation dose rate value, medium radiation dose rate value and high radiation dose rate value according to the limit value of the radiation dose rate that the mobile robot can tolerate; assign different colors to the low radiation dose rate value, medium radiation dose rate value and high radiation dose rate value, and finally fill the corresponding assigned colors in each grid of the grid map to realize the visual display of the low radiation dose rate value, medium radiation dose rate value and high radiation dose rate.

[0102] Compared with the prior art, the present invention has the following advantages:

[0103] 1. A modeling method for a multi-information double-layer hybrid map for mobile robot path planning in a radioactive indoor environment is proposed. The multi-information double-layer hybrid map is a raster map that integrates passability information (i.e., environmental obstacle information) and environmental radiation field information. It can describe the multi-physical field coupling environmental characteristics in an indoor nuclear environment, and provides a reliable environmental map that can perform multi-information value indexing for subsequent nuclear robot path planning.

[0104] 2. In the process of establishing the IOP model, the obstacle-crossing ability of the mobile robot itself and several factors affecting its trafficability are fully considered, so that the trafficability of the target area can be accurately determined.

[0105] The present invention is further described below in conjunction with the figures and embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0106] Figure 1 It is the technical roadmap of the present invention;

[0107] Figure 2 is a schematic diagram of the rectangular coordinate system in step S01;

[0108] Figure 3 This is a schematic diagram of 3D point cloud data removal in step S01;

[0109] Figure 4 It is a schematic diagram of the step climbing model process I in step S02;

[0110] Figure 5 It is a schematic diagram of the step climbing model process II in step S02;

[0111] Figure 6 It is a schematic diagram of the step climbing model process III in step S02;

[0112] Figure 7 It is the longitudinal uphill state diagram of the slope form model in step S02;

[0113] Figure 8 It is the longitudinal downhill state diagram of the slope driving model in step S02;

[0114] Fig. 9 is the lateral driving state diagram of the slope driving model in step S02;

[0115] Fig.10 This is a schematic diagram of the height range calculation in step S02;

[0116] Fig.11 A schematic diagram of a geometric grid map obtained in step S02;

[0117] Fig.12 A schematic diagram of the radiation field grid map obtained in step S03;

[0118] Fig.13 A schematic diagram of a multi-information double-layer hybrid map obtained in step S04;

[0119] Fig.14 The technology roadmap for step S02;

[0120] Fig.15 A flow chart for obtaining a radiation field grid map when source information is known in step S03;

[0121] Fig.16 This is a flow chart for obtaining a radiation field grid map when the source item information in step S03 is unknown. DETAILED DESCRIPTION

[0122] Embodiment 1:

[0123] like Figure 1-16 As shown in the figure, the multi-information double-layer hybrid map modeling method for indoor nuclear environment is based on a mobile robot with double-sided tracks and equipped with a laser radar, and a driven wheel and a driving wheel are respectively provided at the front and rear ends of the inner side of each track.

[0124] Here are the steps:

[0125] S01, get the grid map containing the grid elevation values:

[0126] A. Obtain the 3D point cloud of the indoor nuclear environment and perform coordinate conversion: Control the mobile robot to move in the indoor nuclear environment and obtain complete environmental information through the laser radar; the environmental information is a frame-by-frame 3D point cloud in a spherical coordinate system with the laser radar as the origin; in the spherical coordinate system, the X axis is the forward direction of the mobile robot, the Y axis is the left and right direction of the mobile robot, the Z axis is the upward direction perpendicular to the XLY plane, and L is the origin of the laser radar; the target point P in the 3D point cloud is represented by P(D pl ,θ,φ), where D plis the distance between the target point P and the laser radar origin L, θ is the scanning angle (i.e., line segment D pl The angle between the plane XLY and the is the rotation angle (i.e. line segment D pl The angle between the projection line on the plane XLY and the X-axis in the spherical coordinate system); refer to formula 1 to transform the 3D point cloud from the spherical coordinate system to the rectangular coordinate system;

[0127] Formula 1:

[0128] Where x, y, and z are the x-, y-, and z-axis coordinates of the target point P in the rectangular coordinate system, respectively.

[0129] B. Eliminate low-precision data and invalid data in the 3D point cloud: within the vertical viewing angle of the laser radar, the smaller the scanning angle θ, the longer the detection distance R, the greater the vertical indirection between adjacent laser beams, the sparser the 3D point cloud, and the lower the measurement accuracy; define the 3D point cloud after the detection distance R is 20m as low-precision data, and define the 3D point cloud whose laser beam is blocked by the mobile robot as invalid data; refer to formula 2 to eliminate low-precision data and invalid data from each frame of the 3D point cloud;

[0130] Formula 2:

[0131]

[0132] Where st represents the constraint condition, B represents the width of the mobile robot, and L represents the length of the mobile robot.

[0133] C. Divide the ground into a virtual grid and match the 3D point cloud with it: Use the ground as the reference plane to divide the virtual grid (the grid is square by default unless otherwise specified), and set the side length of the grid to be consistent with the width B of the mobile robot; project the 3D point cloud obtained in step B into the virtual grid, and calculate the weighted average height of all discrete points in each grid as the elevation value of each grid, thus obtaining a grid map containing grid elevation value information.

[0134] S02, establish IOP model to obtain geometric grid map:

[0135] A. Establish a climbing step model: The action of the mobile robot climbing the steps is divided into three processes: Ⅰ. The driven wheel of the mobile robot contacts the edge of the step. Driven by the mobile robot's own power system, the track of the mobile robot close to the step begins to rise to perform obstacle crossing; Ⅱ. The center of gravity G of the mobile robot moves forward and rises, and the angle α between the track and the ground gradually increases. When the center of gravity G of the mobile robot is directly above the edge of the step, it reaches the critical state of climbing the step; Ⅲ. The center of gravity G of the mobile robot passes the edge of the step, and the mobile robot returns to a horizontal state to complete the obstacle crossing;

[0136] The critical state analysis of the mobile robot during the climbing process is carried out and formula 3 is established;

[0137] Formula 3: h0=(l1-htanα)sinα

[0138] Where h0 is the maximum height of the steps that the mobile robot can climb, l1 is the distance between the center of gravity of the mobile robot and the vertical plane where the center line of the driving wheel is located, h is the height from the center of gravity of the mobile robot to the ground below the steps, and α is the angle between the crawler of the mobile robot and the ground below the steps.

[0139] B. Establish a slope driving model: Slope driving includes three situations: longitudinal uphill, longitudinal downhill and lateral driving. As the slope angle increases, the risk of the mobile robot tipping over backwards gradually increases when it is longitudinally uphill, the risk of the mobile robot tipping over forward gradually increases when it is longitudinally downhill, and the risk of the mobile robot tipping over sideways gradually increases when it is lateral driving. The maximum slope angle θ of the mobile robot longitudinal uphill is obtained by analysis. max1 As shown in Formula 4, the maximum slope angle θ of the longitudinal downhill slope max2 As shown in Formula 5, the maximum roll angle θ for lateral travel max3 As shown in formula 6;

[0140] Formula 4:

[0141] Formula 5:

[0142] Formula 6:

[0143] Where l2 is the distance between the center of the mobile robot and the vertical plane where the center line of the driven wheel of the mobile robot is located; B is the width of the mobile robot, and b is the track width of the mobile robot;

[0144] In order to ensure the driving safety of the mobile robot, the maximum inclination angle θ that the mobile robot can pass is obtained by combining formulas 4-6: max As shown in formula 7;

[0145] Formula 7: θ max =min{θ max1 ,θ max2 ,θ max3}

[0146] In the formula, θ max is the maximum inclination angle that the mobile robot can pass.

[0147] C, calculate the height difference among the factors affecting traffic capacity:

[0148] The center of gravity of the mobile robot is projected onto the grid below. The grid into which the projection is projected is the central grid. The center of the central grid is the center of the circle. The maximum height change value of a total of 9 grids in a circular area with a grid side length as the radius is the height range, which is calculated according to Formula 8.

[0149] Formula 8:

[0150] Among them, δH is the height difference value, H x,y Represents the elevation value of the grid in the x-th row and y-th column, x = (i-1, i, i+1), y = (j-1, j, j+1).

[0151] D. Calculate the longitudinal slope and transverse slope influencing factors of traffic capacity: define the central grid and the nearest 8 grids around it as the analysis area, use the data of the nearest 8 grids around the central grid where the mobile robot is located to estimate the slope value of the central grid. The global coordinate system coordinates of the grid in the analysis area are (x i ,y i , z i ), the subscript i is the grid number, z i is the elevation value of the grid numbered i; based on the least squares method, the fitting surface of the analysis area is obtained as shown in formula 9;

[0152] Formula 9: Z = Ix i +Jy i +K

[0153] Among them, Z is the fitting surface function of the analysis area, and I, J, and K are all surface constant coefficients to be determined;

[0154] The point cloud coordinates in the analysis area are transformed from the rectangular coordinate system to the global coordinate system, and the global coordinate system coordinates of the point cloud in the analysis area or the global coordinate system coordinates of all grids in the analysis area are substituted into Formula 9 to obtain the cumulative error formula, as shown in Formula 10;

[0155] Formula 10:

[0156] Wherein, the superscript n is the maximum number of grids or the maximum number of point clouds in the analysis area, and Ω is the cumulative error;

[0157] Taking the minimum cumulative error Ω as the expectation, we take the partial derivative of I, J, and K in formula 10 to obtain formula 11;

[0158] Formula 11:

[0159] Solve formula 11 to obtain the values ​​of the surface constant coefficients I, J, and K, and then obtain the normal vector of the fitting surface in the analysis area In addition, the normal vector of the horizontal plane is known to be Therefore, the slope of the i-th grid in the analysis area is the angle between the two normal vectors, as shown in formula 12. Formula 12 is then split into longitudinal slope and transverse slope to obtain formula 13;

[0160] Formula 12:

[0161] Formula 13:

[0162] in, is the set of longitudinal slope and transverse slope of the i-th grid in the analysis area, is the longitudinal slope of the i-th grid in the analysis area, is the transverse slope of the i-th grid in the analysis area.

[0163] E, calculate the minimum ground clearance among the factors affecting the traffic capacity: define the minimum vertical distance between the horizontal plane where the chassis of the mobile robot is located and the horizontal plane where the elevation value of the target grid is located when the mobile robot successfully passes through the target grid area, that is, the minimum ground clearance; set the center of gravity G of the mobile robot and the center point O of the target grid to be directly opposite to each other (that is, on the same vertical line), and calculate the horizontal position coordinates of the contact points of the left driven wheel T1, the left driving wheel T2, the right driven wheel T3 and the right driving wheel T4 with the ground according to the size parameters of the mobile robot and the slope of the target grid, as shown in Formula 14;

[0164] Formula 14:

[0165] In the formula, the coordinates of the four wheels are all in the global coordinate system, B is the width of the mobile robot, b is the track width of the mobile robot, θ is the forward direction of the mobile robot in the global coordinate system, and θ v is the longitudinal slope of the target grid, θ c is the lateral slope of the target grid, l1 is the distance between the center of gravity of the mobile robot and the vertical plane where the center line of the driving wheel is located, l2 is the distance between the center of gravity of the mobile robot and the vertical plane where the center line of the driven wheel is located, x i ,y i is the horizontal position coordinate of the center point of the target grid in the global coordinate system;

[0166] Considering that the crawler of the mobile robot is a flexible body, it is approximately considered that any three of the left driven wheel T1, the left driving wheel T2, the right driven wheel T3 and the right driving wheel T4 constitute a contact plane, and a plane set Ψ and a plane set Ω parallel to the plane set ψ are constructed accordingly; the plane set ψ includes four planes formed by assembling any three points of T1, T2, T3 and T4, wherein the horizontal position coordinates x, y of the three points in each plane are obtained by referring to formula 14, wherein the height coordinates z of the three points in each plane are calculated by using the IDW interpolation algorithm; correspondingly, the plane set Ω also includes four planes corresponding to the four planes in the plane set Ψ; finally, the minimum ground clearance calculation formula is constructed in combination with the plane distance formula, as shown in formula 15;

[0167] Formula 15:

[0168] In the formula, H grad is the minimum ground clearance, R i is the constant term coefficient of the plane set Ψ, S i is the constant coefficient of the plane set Ω, and the subscript i is the common number of a plane in the plane set Ψ and the corresponding plane in the plane set Ω.

[0169] F, IOP cost calculation of the target grid: Assuming that the mobile robot is in the target grid area, before doing IOP cost calculation, based on H grad With H stand The ratio of is used for pre-pass determination, refer to formula 16;

[0170] Formula 16:

[0171] In the formula, H grad is the chassis ground clearance (an inherent parameter of the mobile robot), Θ(H grad ) is the passability determination factor. If the value is less than 1 (i.e. the minimum ground clearance H grad Less than chassis ground clearance H stand ), the target grid area is determined to be a passable area. If the value is greater than 1 (i.e. the minimum ground clearance H grad Greater than chassis ground clearance H stand ), determine that the target grid area is impassable and assign the target grid a value of 1;

[0172] If the result of the above traffic determination is a passable area, the height difference, transverse slope and longitudinal slope in the traffic capacity influencing factors are used as the IOP cost for calculation, and the piecewise function is defined as shown in Formula 17;

[0173] Formula 17:

[0174] In the formula, represents one of the three factors affecting traffic capacity; when i=1, represents the height difference (calculated by formula 8); when i = 2, represents the longitudinal slope (obtained from formula 13), when i = 3, represents the transverse slope (obtained by formula 13); χ max Represents the limit value of the factors affecting traffic capacity. When i=1, is the maximum height h0 of the mobile robot climbing the steps (obtained from Formula 3), when i = 2, is the maximum slope angle θ of the mobile robot going uphill longitudinally max1 (obtained from formula 4), when i=3, is the maximum roll angle θ of the mobile robot when traveling laterally max3 (obtained from Formula 6), is the weight factor of the corresponding capacity influencing factor; when i=1, Represents the weight factor of the height range. When i=2, Represents the weight factor of the longitudinal slope. When i=3, represents the weight factor of the transverse slope;

[0175] Combining Formula 16 and Formula 17, the IOP cost calculation function can be obtained, as shown in Formula 18;

[0176] Formula 18:

[0177]

[0178] In the formula, represents the IOP cost calculation function, st represents the constraint condition, is the weight factor of the height extreme, is the weight factor of the longitudinal slope, is the weight factor of the transverse slope, IOP1(h) is the IOP cost calculation value of the height extreme (calculated by formula 17), IOP2(θ1) is the IOP cost calculation value of the longitudinal slope (calculated by formula 17), and IOP3(θ2) is the IOP cost calculation value of the transverse slope (calculated by formula 17);

[0179] Formula 18 can be used to assign the target grid an IOP cost value. However, in view of the fact that the environment involved is an indoor core environment where the terrain change trend is not obvious, in order to simplify the IOP cost calculation, the IOP threshold is set again, and the IOP cost value calculated by Formula 18 is binarized using Formula 19 to achieve the binarization of the target grid. The target grid binarization result of 0 indicates a passable area, and the target grid binarization result of 1 indicates an impassable area. Thus, the geometric grid map is successfully constructed.

[0180] Formula 19:

[0181] Where Gmap(i) represents the geometric cost value of the grid numbered i. To customize the IOP threshold.

[0182] S03, obtain the radiation field grid map based on MCNP or difference algorithm:

[0183] When the source information is known, the radiation field is calculated based on MCNP simulation (specifically, the radiation dose rate is calculated), and the calculation results are merged with the grid map (here refers to the grid map in step S01); the specific operations are as follows:

[0184] A. Construction of 3D model of indoor nuclear environment: Establish a 3D model of the indoor nuclear environment corresponding to the grid map, and determine the material properties of each geometric body in the 3D model according to the actual situation on site, which will be used for material definition in the INP file in the MCNP calculation code;

[0185] B. MCNP simulation calculation of radiation dose rate: According to the source information (including the number, type, location and intensity of radiation sources) and environmental characteristics (reference the established 3D model), write MCNP calculation code to perform simulation calculation to obtain the radiation dose rate at each location in the 3D model;

[0186] C. Statistical average radiation dose rate of each grid: Based on the calculated radiation dose rate, the average radiation dose rate of all grids in the grid map within the range of height 0-a is statistically calculated (the purpose is to facilitate the subsequent calculation of the cumulative radiation dose rate when the nuclear robot does path planning), so as to obtain the radiation field grid map (the corresponding average radiation dose rate is recorded in each grid);

[0187] D. Visualize the dose rate in the grid map: divide the average radiation dose rate of each grid in the radiation field grid map into low radiation dose rate value, medium radiation dose rate value and high radiation dose rate value from high to low, and define the low radiation dose rate value, medium radiation dose rate value and high radiation dose rate value according to the limit value of the radiation dose rate that the mobile robot can tolerate; assign different colors to the low radiation dose rate value, medium radiation dose rate value and high radiation dose rate value, and finally fill the corresponding assigned colors in each grid of the radiation field grid map to realize the visual display of the low radiation dose rate value, medium radiation dose rate value and high radiation dose rate value.

[0188] When the source information is unknown, the difference algorithm is used to calculate the radiation field (specifically calculate the radiation source position, radiation source intensity and radiation dose rate), and the calculation results are merged with the grid map (here refers to the grid map in step S01). The specific operations are as follows:

[0189] A. Construction of 3D model of indoor nuclear environment: Establish a 3D model of the indoor nuclear environment corresponding to the grid map, and determine the material properties of each geometric body in the 3D model according to the actual situation on site, which will be used for material definition in the INP file in the MCNP calculation code;

[0190] B. Obtaining the reference radiation field: Enter the above-mentioned indoor nuclear environment to perform random multi-point radiation data sampling. Based on the sampled data and environmental characteristics, use MCNP to simulate and calculate the radiation field to obtain the reference radiation field for reconstructing the radiation field;

[0191] C. Reconstruct the radiation field based on the difference algorithm: extract some grids in the reference radiation field as sampling points, re-enter the above-mentioned indoor nuclear environment to sample radiation data at the extracted grid positions, use the sampled data as the initial data set, program the interpolation function, and reconstruct the radiation field through the interpolation program to obtain the radiation dose rate at each location in the 3D model;

[0192] D. Count the average radiation dose rate of each grid: Based on the reconstructed radiation dose rate, count the average radiation dose rate of all grids in the grid map within the height range of 0-a (the purpose is to facilitate the calculation of the cumulative radiation dose rate during the subsequent path planning of the nuclear robot), so as to obtain the radiation field grid map (the corresponding average radiation dose rate is recorded in each grid).

[0193] E. Visualize the dose rate in the grid map: divide the average radiation dose rate of each grid in the radiation field grid map into low radiation dose rate value, medium radiation dose rate value and high radiation dose rate value from high to low, and define the low radiation dose rate value, medium radiation dose rate value and high radiation dose rate value according to the limit value of the radiation dose rate that the mobile robot can tolerate; assign different colors to the low radiation dose rate value, medium radiation dose rate value and high radiation dose rate value, and finally fill the corresponding assigned colors in each grid of the radiation field grid map to realize the visual display of the low radiation dose rate value, medium radiation dose rate value and high radiation dose rate value.

[0194] In this step, 95% and above of the limit value of the radiation dose rate that the mobile robot can tolerate is taken as the high radiation dose rate value, 40%-95% of the limit value of the radiation dose rate that the mobile robot can tolerate is taken as the medium radiation dose rate value, and 40% and below of the limit value of the radiation dose rate that the mobile robot can tolerate is taken as the low radiation dose rate value.

[0195] In this step, the difference function is any one of cubic spline interpolation, linear interpolation, bilinear interpolation, k-nearest neighbor interpolation and Gaussian regression interpolation functions.

[0196] S04, multi-information dual-layer hybrid map construction:

[0197] In order to provide more accurate environmental data input for the path planning of mobile robots, after the construction of the geometric grid map and the radiation field grid map is completed, it is necessary to fuse them to form a two-layer hybrid map (the hybrid map is in the form of a grid map) containing passability information (equivalent to obstacle information) and radiation dose rate information, so as to realize the expression of multi-physical field coupling information;

[0198] A. Encode each grid information in the grid map: For the geometric grid map, the passability of each grid is encoded in binary, with non-passable encoded as 1 and passable encoded as 0; for the radiation field grid map, the average radiation dose rate of each grid is encoded in real numbers, that is, the average radiation dose rate corresponding to each grid is recorded in that grid. After the grid information encoding is completed, it is stored in the form of groups to form a mapping relationship, that is, a multi-information two-layer hybrid map is obtained.

[0199] B. Identify the grids in the multi-information two-layer hybrid map: Use a hybrid identification method that combines the rectangular coordinate system method and the serial number method to equivalently represent the grids; assume that the rectangular coordinate of any grid k is (i, j), 0 < i < M, 0 < j < N, where M is the number of grid columns in the x-axis direction of the rectangular coordinate system and N is the number of grid rows in the y-axis direction of the rectangular coordinate system, then the grid k is represented as k(i, j). If the corresponding number of the grid k is n, then the grid k is equivalently represented in the form of k(n); define the grid numbering starting from the first grid in the upper left corner of the rectangular coordinate system, and determine the serial number in the order from left to right and top to bottom, with the serial number starting from 1. According to the above definition, the mapping relationship from the coordinates (i, j) to the number n is shown in Equation 20;

[0200] Equation 20:

[0201] In the formula, mod is the remainder function, that is, the remainder after the division operation of two numerical expressions; the ceil function is used for rounding up, that is, obtaining an integer value greater than or equal to the parameter.

Claims

1. A multi-information double-layer hybrid map modeling method for indoor nuclear environment, characterized by the following steps: as follows: S01, obtaining a grid map containing grid elevation values; S02, establishing an IOP model to obtain a geometric grid map; S03, obtaining a radiation field grid map based on MCNP or difference algorithm; S04, multi-information dual-layer hybrid map construction: binary encoding of the accessibility of each grid in the geometric grid map; The average radiation dose rate of each grid in the radiation field grid map is encoded with real numbers; after the grid information is encoded, it is stored in the form of group numbers to form a mapping relationship.

2. The method for modeling a multi-information double-layer hybrid map of an indoor nuclear environment as claimed in claim 1, characterized in that: The S01 step is briefly described as follows: A, obtain the 3D point cloud of the indoor nuclear environment and transform it into a rectangular coordinate system; B. Eliminate low-precision data and invalid data in the 3D point cloud. Invalid data refers to the 3D point cloud where the laser beam is blocked by the mobile robot, and low-precision data refers to the 3D point cloud after the detection distance is 20m. C, matching the virtual grid divided on the ground with the 3D point cloud after data removal to obtain a grid map containing grid elevation value information; The S02 step is briefly described as follows: A, calculate the maximum height of the steps that the mobile robot can climb; B, calculate the maximum inclination angle that the mobile robot can pass; C, calculate the height difference of the factors affecting traffic capacity; D, calculate the longitudinal slope and transverse slope influencing factors of traffic capacity; E, calculate the minimum ground clearance among the factors affecting traffic capacity; F, IOP cost calculation of the target grid: First, the passability of the target grid area is determined based on the values ​​of chassis ground clearance and minimum ground clearance; If it is passable, the IOP cost value of the target grid is calculated by combining the height range, the transverse slope and the longitudinal slope; finally, the IOP threshold is set and the calculated IOP cost value is binarized to distinguish between the impassable and passable areas; The AB steps in this step are in no particular order, and the CE steps in this step are in no particular order; The S03 step is briefly described as follows: When the source information is known, the radiation field is calculated based on MCNP simulation, and the calculation results are fused with the grid map; when the source information is unknown, the radiation field is calculated using the difference algorithm, and the calculation results are fused with the grid map.

3. The method for modeling a multi-information double-layer hybrid map of an indoor nuclear environment as claimed in claim 2, characterized in that: The method is based on a mobile robot with double-sided tracks and equipped with a laser radar, wherein a driven wheel and a driving wheel are respectively arranged at the front and rear ends of the inner side of each track.

4. The method for modeling a multi-information double-layer hybrid map of an indoor nuclear environment as claimed in claim 3 is characterized by: The S01 step is described in detail as follows: A. Obtain the 3D point cloud of the indoor nuclear environment and perform coordinate conversion: Control the mobile robot to move in the indoor nuclear environment and obtain complete environmental information through the laser radar; the environmental information is a frame-by-frame 3D point cloud in a spherical coordinate system with the laser radar as the origin; in the spherical coordinate system, the X axis is the forward direction of the mobile robot, the Y axis is the left and right direction of the mobile robot, the Z axis is the upward direction perpendicular to the XLY plane, and L is the origin of the laser radar; the target point P in the 3D point cloud is represented by P(D pl ,θ,φ), where D pl is the distance between the target point P and the laser radar origin L, θ is the scanning angle, is the rotation angle; refer to formula 1 to transform the 3D point cloud from the spherical coordinate system to the rectangular coordinate system; Formula 1: Where x, y, and z are the x-, y-, and z-axis coordinates of the target point P in the rectangular coordinate system, respectively; B. Eliminate low-precision data and invalid data in the 3D point cloud: define the 3D point cloud after the detection distance R is 20m as low-precision data, and define the 3D point cloud where the laser beam is blocked by the mobile robot as invalid data; refer to formula 2 to eliminate low-precision data and invalid data from each frame of the 3D point cloud; Formula 2: In the formula, st represents the constraint condition, B represents the width of the mobile robot, and L represents the length of the mobile robot; C. Divide the ground into virtual grids and match the 3D point cloud with them: Divide the virtual grid using the ground as the reference plane, and set the side length of the grid to be consistent with the width B of the mobile robot; project the 3D point cloud obtained in step B into the virtual grid, and calculate the weighted average height of all discrete points in each grid as the elevation value of each grid, thus obtaining a grid map containing grid elevation value information.

5. The method for modeling a multi-information double-layer hybrid map of an indoor nuclear environment as claimed in claim 4, characterized in that: The S02 step is described in detail as follows: A. Establish a climbing step model: The action of the mobile robot climbing the steps is divided into three processes: Ⅰ. The driven wheel of the mobile robot contacts the edge of the step. Driven by the mobile robot's own power system, the track of the mobile robot close to the step begins to rise to perform obstacle crossing; Ⅱ. The center of gravity G of the mobile robot moves forward and rises, and the angle α between the track and the ground gradually increases. When the center of gravity G of the mobile robot is directly above the edge of the step, it reaches the critical state of climbing the step; Ⅲ. The center of gravity G of the mobile robot passes the edge of the step, and the mobile robot returns to a horizontal state to complete the obstacle crossing; The critical state analysis of the mobile robot during the climbing process is carried out and formula 3 is established; Formula 3: h0=(l1-htanα)sinα Where h0 is the maximum height of the step that the mobile robot can climb, l1 is the distance between the center of gravity of the mobile robot and the vertical plane where the center line of the driving wheel is located, h is the height from the center of gravity of the mobile robot to the ground below the step, and α is the angle between the crawler of the mobile robot and the ground below the step; B. Establish a slope driving model: Slope driving includes three situations: longitudinal uphill, longitudinal downhill and lateral driving. As the slope angle increases, the risk of the mobile robot tipping over backwards gradually increases when it is longitudinally uphill, the risk of the mobile robot tipping over forward gradually increases when it is longitudinally downhill, and the risk of the mobile robot tipping over sideways gradually increases when it is lateral driving. The maximum slope angle θ of the mobile robot longitudinal uphill is obtained by analysis. max1 As shown in Formula 4, the maximum slope angle θ of the longitudinal downhill slope max2 As shown in Formula 5, the maximum roll angle θ for lateral travel max3 As shown in formula 6; Formula 4: Formula 5: Formula 6: Where l2 is the distance between the center of the mobile robot and the vertical plane where the center line of the driven wheel of the mobile robot is located; B is the width of the mobile robot, and b is the track width of the mobile robot; Combining formulas 4-6, we can get the maximum inclination angle θ that the mobile robot can pass. max As shown in formula 7; Formula 7: θ max = min{θ max1 , θ max2 , θ max3}; In the formula, θ max is the maximum inclination angle that the mobile robot can pass; C, calculate the height difference among the factors affecting traffic capacity: The center of gravity of the mobile robot is projected onto the grid below. The grid into which the projection is projected is the central grid. The center of the central grid is the center of the circle. The maximum height change value of a total of 9 grids in a circular area with a grid side length as the radius is the height range, which is calculated according to Formula 8. Among them, δH is the height difference value, H x,y Represents the elevation value of the grid at the xth row and yth column, x = (i-1, i, i+1), y = (j-1, j, j+1); D. Calculate the longitudinal slope and transverse slope influencing factors of traffic capacity: define the central grid and the nearest 8 grids around it as the analysis area, use the data of the nearest 8 grids around the central grid where the mobile robot is located to estimate the slope value of the central grid. The global coordinate system coordinates of the grid in the analysis area are (x i ,y i , z i ), the subscript i is the grid number, z i is the elevation value of the grid numbered i; based on the least squares method, the fitting surface of the analysis area is obtained as shown in formula 9; Formula 9: Z = Ix i + Jy i + K Among them, Z is the fitting surface function of the analysis area, and I, J, and K are all surface constant coefficients to be determined; The point cloud coordinates in the analysis area are transformed from the rectangular coordinate system to the global coordinate system, and the global coordinate system coordinates of the point cloud in the analysis area or the global coordinate system coordinates of all grids in the analysis area are substituted into Formula 9 to obtain the cumulative error formula, as shown in Formula 10; Formula 10: Wherein, the superscript n is the maximum number of grids or point clouds in the analysis area, and Ω is the cumulative error; Taking the minimum cumulative error Ω as the expectation, we take the partial derivative of I, J, and K in formula 10 to obtain formula 11; Formula 11: Solve formula 11 to obtain the values ​​of the surface constant coefficients I, J, and K, and then obtain the normal vector of the fitting surface in the analysis area In addition, the normal vector of the horizontal plane is known to be Therefore, the slope of the i-th grid in the analysis area is the angle between the two normal vectors, as shown in formula 12. Formula 12 is then split into longitudinal slope and transverse slope to obtain formula 13; Formula 12: Formula 13: in, is the set of longitudinal slope and transverse slope of the i-th grid in the analysis area, is the longitudinal slope of the i-th grid in the analysis area, is the transverse slope of the i-th grid in the analysis area; E, calculate the minimum ground clearance among the factors affecting the traffic capacity: define the minimum vertical distance between the horizontal plane where the chassis of the mobile robot is located and the horizontal plane where the elevation value of the target grid is located when the mobile robot successfully passes through the target grid area, that is, the minimum ground clearance; set the center of gravity G of the mobile robot to be vertically opposite to the center point O of the target grid, and calculate the horizontal position coordinates of the contact points of the left driven wheel T1, the left driving wheel T2, the right driven wheel T3 and the right driving wheel T4 with the ground according to the size parameters of the mobile robot and the slope of the target grid, as shown in Formula 14; Formula 14: In the formula, the coordinates of the four wheels are all in the global coordinate system, B is the width of the mobile robot, b is the track width of the mobile robot, θ is the forward direction of the mobile robot in the global coordinate system, and θ v is the longitudinal slope of the target grid, θ c is the lateral slope of the target grid, l1 is the distance between the center of gravity of the mobile robot and the vertical plane where the center line of the driving wheel is located, l2 is the distance between the center of gravity of the mobile robot and the vertical plane where the center line of the driven wheel is located, x i ,y i is the horizontal position coordinate of the center point of the target grid in the global coordinate system; Any three of the left driven wheel T1, the left driving wheel T2, the right driven wheel T3 and the right driving wheel T4 constitute a contact plane, and a plane set Ψ and a plane set Ω parallel to the plane set Ψ are constructed based on this; the plane set Ψ includes four planes formed by assembling any three points of T1, T2, T3 and T4, wherein the horizontal position coordinates x, y of the three points in each plane are obtained by referring to formula 14, wherein the height coordinates z of the three points in each plane are calculated by using the IDW interpolation algorithm; correspondingly, the plane set Ω also includes four planes corresponding to the four planes in the plane set Ψ; finally, the minimum ground clearance calculation formula is constructed in combination with the plane distance formula, as shown in formula 15; Formula 15: In the formula, H grad is the minimum ground clearance, R i is the constant term coefficient of the plane set Ψ, S i is the constant term coefficient of the plane set Ω, and the subscript i is the common number of a plane in the plane set Ψ and the corresponding plane in the plane set Ω; F, IOP cost calculation of the target grid: Assuming that the mobile robot is in the target grid area, before doing IOP cost calculation, based on H grad With H stand The ratio of is used for pre-pass determination, refer to formula 16; Formula 16: In the formula, H grad is the chassis ground clearance, Θ(H grad ) is the passability determination factor. If the value is less than 1, the target grid area is determined to be a passable area. If the value is greater than 1, the target grid area is determined to be impassable and the target grid is assigned a value of 1. If the result of the above traffic determination is a passable area, the height difference, transverse slope and longitudinal slope in the traffic capacity influencing factors are used as the IOP cost for calculation, and the piecewise function is defined as shown in Formula 17; Formula 17: In the formula, χ represents one of the three factors affecting traffic capacity; when i = 1, χ represents the height difference; when i = 2, χ represents the longitudinal slope; when i = 3, χ represents the transverse slope; max represents the limit value of the factors affecting traffic capacity. When i=1, χ max is the maximum height h0 of the mobile robot climbing the stairs. When i = 2, X max is the maximum slope angle θ of the mobile robot going uphill longitudinally max1 , when i = 3, χ max is the maximum roll angle θ of the mobile robot when traveling laterally max3 , is the weight factor of the corresponding capacity influencing factor; when i=1, Represents the weight factor of the height range. When i=2, Represents the weight factor of the longitudinal slope. When i=3, represents the weight factor of the transverse slope; Combining Formula 16 and Formula 17, the IOP cost calculation function can be obtained, as shown in Formula 18; Formula 18: In the formula, represents the IOP cost calculation function, st represents the constraint condition, is the weight factor of the height extreme, is the weight factor of the longitudinal slope, is the weight factor of the transverse slope, IOP1(h) is the IOP cost calculation value of the height extreme difference, IOP2(θ1) is the IOP cost calculation value of the longitudinal slope, and IOP3(θ2) is the IOP cost calculation value of the transverse slope; The IOP threshold is set, and the IOP cost value calculated by formula 18 is binarized using formula 19 to realize the binarization of the target grid; the binarization result of the target grid is 0, indicating a passable area, and the binarization result of the target grid is 1, indicating an impassable area; Formula 19: Where Gmap(i) represents the geometric cost value of the grid numbered i. To customize the IOP threshold.

6. The method for modeling a multi-information double-layer hybrid map of an indoor nuclear environment as claimed in claim 5, characterized in that: The S03 step is described in detail as follows: When the source information is known, the specific operations are as follows: A. Construction of 3D model of indoor nuclear environment: Establish a 3D model of the indoor nuclear environment corresponding to the grid map, and determine the material properties of each geometric body in the 3D model according to the actual situation on site, which will be used for material definition in the INP file in the MCNP calculation code; B. MCNP simulation calculation of radiation dose rate: According to the source information and environmental characteristics, write MCNP calculation code to perform simulation calculation to obtain the radiation dose rate at various locations in the 3D model; C. Count the average radiation dose rate of each grid: Based on the calculated radiation dose rate, count the average radiation dose rate of all grids in the grid map within the range of height 0-a, so as to obtain the radiation field grid map; The specific operations when the source item information is unknown are as follows: A. Construction of 3D model of indoor nuclear environment: Establish a 3D model of the indoor nuclear environment corresponding to the grid map, and determine the material properties of each geometric body in the 3D model according to the actual situation on site, which will be used for material definition in the INP file in the MCNP calculation code; B. Obtain the reference radiation field: Enter the above-mentioned indoor nuclear environment to perform random multi-point radiation data sampling. Based on the sampling data and environmental characteristics, use MCNP to simulate and calculate the radiation field to obtain the reference radiation field for reconstructing the radiation field. C. Reconstruct the radiation field based on the difference algorithm: Extract some grids as sampling points in the reference radiation field, and enter the above-mentioned indoor nuclear environment again to sample the radiation data at the positions of the extracted grids. Using the sampling data as the initial data set, program the interpolation function, and realize the reconstruction of the radiation field through the interpolation program to obtain the radiation dose rate at each location in the 3D model. D. Statistically calculate the average radiation dose rate of each grid: Based on the reconstructed radiation dose rate, statistically calculate the average radiation dose rate of all grids in the grid map within the height range of 0 - a, so as to obtain the radiation field grid map.

7. The indoor nuclear environment multi-information double-layer hybrid map modeling method according to claim 6, characterized in that: Step S04 is described in detail as follows: A. Encode each grid information in the grid map: Perform binary encoding on the passability information of the grids obtained in step S02. The grid is encoded as 1 if it is impassable and 0 if it is passable. Adopt real number encoding for the average radiation dose rate information of the grids obtained in step S03, that is, record the average radiation dose rate corresponding to each grid in the grid. After the grid information encoding is completed, store it in the form of groups to form a mapping relationship, that is, obtain the multi-information double-layer hybrid map.

8. The method for modeling a multi-information double-layer hybrid map of an indoor nuclear environment as claimed in claim 7, characterized in that: Step S04 further includes sub-step B that follows sub-step A. B. Perform grid identification on the multi-information double-layer hybrid map: Use a hybrid identification method that combines the rectangular coordinate system method and the serial number method to equivalently represent the grids. Let the rectangular coordinate system coordinates of any grid k be (i, j), 0 < i < M, 0 < j < N, where M is the number of grid columns in the x-axis direction of the rectangular coordinate system and N is the number of grid rows in the y-axis direction of the rectangular coordinate system. Then the grid k is represented as k(i, j). If the corresponding number of the grid k is n, then the grid k is equivalently represented in the form of k(n). Define the grid numbering to start from the first grid in the upper left corner of the rectangular coordinate system, and determine the serial number in the order from left to right and top to bottom. The serial number starts from 1. According to the above definition, the mapping relationship from the coordinates (i, j) to the number n is shown in formula 20. Formula 20: In the formula, mod is the remainder function, that is, the remainder after the division operation of two numerical expressions; the ceil function is used for rounding up, that is, obtaining an integer value greater than or equal to the parameter.

9. The method for modeling a multi-information double-layer hybrid map of an indoor nuclear environment as claimed in claim 8, characterized in that: Step S03 further includes sub-step D that follows sub-step C in the part where the source term information is known. D. Visualize the dose rate in the grid map: divide the average radiation dose rate of each grid in the radiation field grid map into low radiation dose rate value, medium radiation dose rate value and high radiation dose rate value from high to low, and define the low radiation dose rate value, medium radiation dose rate value and high radiation dose rate value according to the limit value of the radiation dose rate that the mobile robot can tolerate; assign different colors to the low radiation dose rate value, medium radiation dose rate value and high radiation dose rate value, and finally fill the corresponding assigned colors in each grid of the radiation field grid map to realize the visual display of the low radiation dose rate value, medium radiation dose rate value and high radiation dose rate value.

10. The method for modeling a multi-information double-layer hybrid map of an indoor nuclear environment as claimed in claim 9, characterized in that: Step S03 also includes sub-step E following sub-step D of the unknown part of the source item information; E. Visualize the dose rate in the grid map: divide the average radiation dose rate of each grid in the radiation field grid map into low radiation dose rate value, medium radiation dose rate value and high radiation dose rate value from high to low, and define the low radiation dose rate value, medium radiation dose rate value and high radiation dose rate value according to the limit value of the radiation dose rate that the mobile robot can tolerate; assign different colors to the low radiation dose rate value, medium radiation dose rate value and high radiation dose rate value, and finally fill the corresponding assigned colors in each grid of the radiation field grid map to realize the visual display of the low radiation dose rate value, medium radiation dose rate value and high radiation dose rate value.

Images