Abstract
Tree Planting Machine (TPM) is subject to a TreePlanting Robot (TPR) with desired tracking trajectory planning. In this topic, taking the TPR proposed as the analysis object, the positive and inverse solutions of the kinematics are analyzed to explore the optimal trajectory planning. An improved position/posture algorithm, based on the analytical solution of the inverse kinematics of the TPR, is proposed. The trajectory planning strategy for TPR in Cartesian coordinate system and Joint coordinate system is discussed, which is used for parabolic transition linear programming optimization, and the simulation model of TPR trajectory planning is constructed by MATLAB module. Numerical simulation results indicate that the deviation of the TPR trajectory from the expected value is significantly reduced. The proposed improved position/posture algorithm is verified by kinematic analysis, and the TPR followability and trajectory planning accuracy are greatly improved. Toward this goal, a variable trajectory planning can be effectively, and stability adjusted by predesigned TPM system in the field of ecological tree planting.
Highlights
 A TreePlanting Robot (TPR) with desired tracking trajectory planning is proposed.
 The positive and inverse solutions of the kinematics are analyzed to explore the optimal trajectory planning.
 An improved position/posture algorithm, based on the analytical solution of the inverse kinematics of the TPR, is proposed.
 The trajectory planning strategy for TPR in Cartesian coordinate system and Joint coordinate system is discussed.
1. Introduction
Domestic forestry operations today still rely on human labor to get the job done, especially drilling deep afforestation, which is done by hand. This method is not only timeconsuming, but also laborintensive and inefficient, greatly reducing the speed of afforestation [13]. Forest tree planting is a complex undertaking, with labor costs generally account for more than 60 % to 80 % of the actual cost. Labor force is also facing a phase of transformation from manual labor to skilled labor, with the cost of manual labor increasing year on year [46]. Considering the increase in labor costs, available labor is limited, and production costs are increasing, therefore, TPRs are urgently needed to perform and simplify these field construction operations mechanically. In addition, as people’s living standards improve and their demand for timber increases, the rate of artificial afforestation is not yet as fast as the growth of people's demand for timber [7]. There has been some research and applications in addressing the inverse kinematics and path planning of robots. However, these studies have focused on industrial robots in unconstrained workspaces with homogeneous objects. TPRs integrate sensing tests, simulation calculations, manipulation and control systems to perform forest tree planting operations with artificial intelligence technology. Moreover, the successful adoption of these technologies requires a combination of the TPR's capabilities and the environment in which it operates in the forest. The environment in which trees are planted is known to be complex. TPR and end effectors attributed to them are likely to collide with adjacent branch obstacles, thus reducing the efficiency of tree planting. The application of TPR technology in forestry operating environments faces considerable challenges. In order to successfully adopt robots in forestry, it is necessary to develop an efficient motion planning algorithm that allows TPRs to complete tree planting tasks.
Forestry is a constrained and dynamic environment with working objects of varying shapes, sizes, positions and orientations. The target attitude of the TPR in each planting task is unknown and the path from the starting point to the target point needs to be replanned for each target point of the planting task. Furthermore, in order to meet the picking efficiency of TPR, the computational speed of the motion planning algorithm should be motion planning algorithm should be high. Therefore, in order to meet the above requirements, not only a motion planning algorithm is required, but also a motion planning algorithm. The motion planning algorithm for TPR needs to solve not only the inverse kinematic problem of obtaining picking poses for tree picking through a small number of iterations, but also to obtain picking poses through a small number of iterations based on information about the forest environment, and to determine the picking path quickly enough to avoid collisions between the TPR and branching obstacles.
This study is based on the above considerations, a motion path planning method for TPR based on the above research is proposed, the simulation results confirm that the proposed motion planning method enables the TPR to avoid obstacles in the workspace and to complete the tree planting task efficiently. The contradiction between the high cost of artificial afforestation, the shortage of labor, the slow pace of afforestation and the increasing demand for timber in people's daily lives is becoming increasingly apparent. The design of the robot has been specially developed by the team of coauthors (whose patent for the invention has been granted), as shown in Fig. 1.
This paper presents a novel treeplanting robot (TPR) with a 6degree of freedom (DOF) that grasps poplar saplings from a walking mechanism and places them in drilled holes. The most novel aspect of this topic is the optimization of a series of control strategy problems for the TPM. A more detailed kinematic analysis of the TPR is elaborated, which is a prerequisite for the TPM control strategy optimization scheme, which also plays an important foundation for the subsequent TPR dynamics, trajectory planning and offline programming.
Fig. 1Selfpropelled automatic poplar TPM
a) Axonometric view
b) Main view
2. DH model for kinematics solution
In this paper, a DenavitHartenberg (DH) coordinate system for TPR modelling is proposed based on the DH parametric approach, as shown in Fig. 2. The linkage parameters of the corresponding DH model are shown in Table 1, where, ${a}_{1}$ = 25 mm, ${a}_{2}$ = 560 mm, ${a}_{3}$ = 35 mm, ${d}_{4}$ = 515 mm. ${\theta}_{i}$ represents the joint angle variable of each joint of TPR, ${\alpha}_{i}$ indicates the joint torsion angle, ${a}_{i1}$ denotes the length of the link rod, ${d}_{i}$ is the offset of the link rod.
Using the modified DH model coordinate transformation [810], the transformation expression is written as Eq. (1):
$=\left[\begin{array}{cccc}\mathrm{cos}{\theta}_{i}& \mathrm{sin}{\theta}_{i}& 0& {\alpha}_{i1}\\ \mathrm{sin}{\theta}_{i}\mathrm{cos}{\alpha}_{i1}& \mathrm{cos}{\theta}_{i}\mathrm{cos}{\alpha}_{i1}& \mathrm{sin}{\alpha}_{i1}& {d}_{i}\mathrm{sin}{\alpha}_{i1}\\ \mathrm{sin}{\theta}_{i}\mathrm{cos}{\alpha}_{i1}& \mathrm{cos}{\theta}_{i}\mathrm{cos}{\alpha}_{i1}& \mathrm{cos}{\alpha}_{i1}& {d}_{i}\mathrm{cos}{\alpha}_{i1}\\ 0& 0& 0& 1\end{array}\right].$
Fig. 2DH coordinate system of a TPR
a) Axonometric view
b) DH coordinate system
Table 1DH model linkage parameters for TPR
Rod number  ${\theta}_{i}$_{}(°)  ${\alpha}_{i1}$ (°)  ${a}_{i1}$ (mm)  ${d}_{i}$ (mm)  Variable range (°) 
1  ${\theta}_{1}$  0  0  0  –170  +170 
2  ${\theta}_{2}$  –90°  ${a}_{1}$  0  –190  +45 
3  ${\theta}_{3}$  0  ${a}_{2}$  0  –120  +156 
4  ${\theta}_{4}$  –90°  ${a}_{3}$  ${d}_{4}$  –185  +185 
5  ${\theta}_{5}$  90°  0  0  –120  +120 
6  ${\theta}_{6}$  –90°  0  0  –350  +350 
The coordinate system is transformed according to the modified DH transformation order, and the connecting rod parameters in Table 1 are substituted into Eq. (1) to obtain the flush transformation matrix ${}_{1}{}^{0}\mathrm{}T$, ${}_{2}{}^{1}\mathrm{}T$, ${}_{3}{}^{2}\mathrm{}T$, ${}_{4}{}^{3}\mathrm{}T$, ${}_{5}{}^{4}\mathrm{}T$, ${}_{6}{}^{5}\mathrm{}T$. Multiplying the above matrices from left to right, the forward motion equation of the TPR is expressed as Eq. (2) [1113]:
where:
${c}_{i}=\mathrm{c}\mathrm{o}\mathrm{s}{\theta}_{i},\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{s}_{i}=\mathrm{s}\mathrm{i}\mathrm{n}{\theta}_{i},\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{c}_{ij}=\mathrm{c}\mathrm{o}\mathrm{s}\left({\theta}_{i}+{\theta}_{j}\right),{s}_{ij}=\mathrm{s}\mathrm{i}\mathrm{n}\left({\theta}_{i}+{\theta}_{j}\right).$
2.1. Reverse kinematic solutions for TPR position
This topic adopts the closed analysis method, and takes the KR10 R1100 robot, the main body of the TPR, as the target of the analysis of the kinematic mechanism. According to the judgment standard of the reverse solution of the articulated robot proposed by Pieper [1416], the 4th, 5th, and 6th axes of the TPR are preset to intersect at one point, which leads to a closed numerical solution of TPR.
Separate the part of Eq. (2) that contains from the righthand side of the equation ${\theta}_{1}$ to facilitate numerical solution of ${\theta}_{1}$. Rectify the two sides of Eq. (2) and multiply by the matrix ${}_{1}{}^{0}{T}^{1}$[17]:
Herein, the lefthand side of Eq. (3) is expanded and represented by the following matrix:
Herein, the righthand end of Eq. (3) is expanded and the following matrix is written:
Setting the elements of the variables on both sides of Eq. (3) equal, the following equation is listed:
Using the trigonometric constant conversion, the expression is written:
where, $\rho =\sqrt{{P}_{x}^{2}+{P}_{y}^{2}}$, $\phi =A\mathrm{tan}2({P}_{x},{P}_{y})$.
Substituting Eq. (5) into Eq. (4), the analytical solution for ${\theta}_{1}$ in terms of differential angular quantities is expressed as:
Assuming that the elements of the variables on both sides of Eq. (3) are equal [18], the following equation is written as:
By combining Eq. (7) and Eq. (4), the following equation is derived:
Using the trigonometric transformation, the analytic solution of the parameter ${\theta}_{3}$ is expressed as:
Rearranging from the Eq. (3), the following equation is described as:
Let, set the parameter elements on both sides of Eq. (10) to be equal [19], the following equation is shown as:
Modifying Eq. (11), the following equation is expressed as:
Simultaneous Eq. (9) and Eq. (12), the analytical solution of parameter ${\theta}_{2}$ is described as:
2.2. Reverse kinematic solutions for TPR posture
According to the above Eq. (10), the analytical solution of the parameter ${\theta}_{4}$ is written as [20]:
Reorganizing Eq. (10) further, the following equation is expressed as:
Let, set the parameter elements on both sides of Eq. (15) to be equal [21], the analytical solution of parameter ${\theta}_{5}$ is:
Rearrange the Eq. (15), set the parameter elements on both sides of the Eq. (15) to be equal, and the analytical solution of the parameter ${\theta}_{6}$ is shown as:
Throughout the above inverse solution analysis results, the inverse solution for the same posture of the TPR at the same point in space is not unique. The TPM is expressed as eight sets of inverse solutions for the same spatial posture, some of which are discarded due to limitations in the range of joint motion, and of the remaining inverse solutions, the one closest to the current TPR posture is preferred.
3. Trajectory planning analysis for TPM
In this paper, the linear interpolation of the parabolic transition is used to study the TPR trajectory planning of the TPM in the Cartesian coordinate system and the joint coordinate system, and through the numerical verification of MATLAB, the simulation results in the two coordinate systems are compared and analyzed, and the planning method is selected and evaluated.
Considering that the coordinate system 6 is located at the intersection of the 4th, 5th, and 6th axes of the manipulator, it is convenient to simplify, and the end effector of the manipulator is established in the coordinate system C, as shown in Fig. 3. The position transformation matrix of coordinate system C under coordinate system 6 is deduced.
Fig. 3Diagram of relationship between end effector in coordinate systems C and 6
In Fig. 3, the posture of the coordinate system C of the end effector of the TPM manipulator is set to be consistent with the posture of the coordinate system 6. The position moment of coordinate system C under coordinate system 6 is written as [22]:
where, $dc$ is the distance between coordinate system 6 and coordinate system C in the direction of the $Z$axis, the value of which is shown as $dc=$ 322.93 mm.
The starting point ${P}_{0}$, the ending point ${P}_{f}$ and the three intermediate points ${P}_{1}$, ${P}_{2}$ and ${P}_{3}$ of the TPR endeffector in the planned path are marked out, and the position and pose of each of these points and the moment of each point are planned under the base coordinate system 0 of TPR. The numerical simulation analysis data is shown in Table 2.
Table 2Position component of coordinate systems C relative to 0 at each point on the path
Track points  Time (s)  $X$ (mm)  $Y$ (mm)  $Z$ (mm)  ${\varphi}_{x}$ ()  ${\varphi}_{y}$ ()  ${\varphi}_{z}$ () 
${P}_{0}$  0  897.5  0  85  –90  0  –90 
${P}_{1}$  5  897.5  0  200  –90  0  –90 
${P}_{2}$  15  0  1238  500  0  0  0 
${P}_{3}$  20  0  1238  380  0  –90  0 
${P}_{f}$  30  0  1238  –370  0  –90  0 
3.1. Trajectory planning in Cartesian coordinate system
The position components of coordinate system 6 relative to coordinate system 0 at each point of the path were found using the “$XYZ$” transformation sequence according to Table 2 and ${}_{6}{}^{0}\mathrm{}{T}_{i}={}_{C}{}^{0}\mathrm{}{T}_{i}{}_{C}{}^{6}\mathrm{}{T}^{1}$$(i=0~f)$, as shown in Table 3. Setting the acceleration time of the parabolic transition section, solving for the velocity and acceleration of the attitude components at each point, the trajectory is simulated in MATLAB as shown in Fig. 4.
Table 3Positional components of coordinate systems 6 with respect to 0 at each point of the path
Track points  Time (s)  $X$ (mm)  $Y$ (mm)  $Z$ (mm)  ${\varphi}_{x}$ ()  ${\varphi}_{y}$ ()  ${\varphi}_{z}$ () 
${P}_{0}$  0  897.5  0  407.93  180  0  0 
${P}_{1}$  5  897.5  0  522.93  180  0  0 
${P}_{2}$  15  0  915.07  500  –90  0  –90 
${P}_{3}$  20  0  915.07  380  –90  0  –180 
${P}_{f}$  30  0  915.07  –370  –90  0  –180 
Fig. 4Parabolic linear interpolation of trajectory curves in Cartesian coordinate system
With the above planning method, the displacement and velocity profiles of the position component of coordinate system 6 with respect to coordinate system 0 and the angular displacement profiles of the six joints are shown in Fig. 5 and Fig. 6 respectively.
Fig. 5Displacement and velocity curves for displacement components
3.2. Trajectory planning in Joint coordinate system
According to the path points set in Table 2 for trajectory planning in the joint coordinate system, the kinematic inverse solution of the poses of the points in Table 2 is carried out to obtain the six joint angle values of the manipulator at each path point to, as shown in Table 4.
Setting the acceleration time ${t}_{a}=$ 1.5 s for the parabolic transition section, the velocities and accelerations of the six joint angles were solved for and the motion trajectory is simulated in MATLAB as shown in Fig. 7.
Fig. 6Angular displacement curve of six joints
Table 4Six joint angle values of the manipulator at each path point
Track points  Time(s)  ${\theta}_{1}$ ($\xb0$)  ${\theta}_{2}$ ($\xb0$)  ${\theta}_{3}$ ($\xb0$)  ${\theta}_{4}$ ($\xb0$)  ${\theta}_{5}$ ($\xb0$)  ${\theta}_{6}$ ($\xb0$) 
${P}_{0}$  0  0  –50.4138  –33.0731  0  83.4868  0 
${P}_{1}$  5  0  –49.2030  –47.9657  0  97.1687  0 
${P}_{2}$  15  90  –47.0087  –49.1878  0  6.1965  0 
${P}_{3}$  20  90  –47.9434  –34.1929  0  –7.8638  –90 
${P}_{f}$  30  90  –2.6970  –33.2546  0  –54.0484  –90 
Fig. 7Parabolic linear interpolation trajectory curves in the joint coordinate system
The angular displacement and angular velocity curves for the six joints under this planning method and the three displacement component curves for the origin of the coordinate system C under the Cartesian coordinate system are shown in Fig. 8 and Fig. 9 respectively.
By comparing Fig. 4 and Fig. 7, it can be seen that when the trajectory planning method with parabolic transition linear interpolation is used for trajectory planning in a rightangle coordinate system, the trajectory at the end of the robot arm basically consists of straightline segments and parabolic segments, and its trajectory matches the path planned at the beginning. When the trajectory planning method is used for trajectory planning in the joint coordinate system, the trajectory at the end of the robot arm basically consists of parabolic and some straightline segments, and its final trajectory is different from the path planned at the beginning.
Poplar saplings avoid the platen in the compaction mechanism during the planting process, so that the saplings can be put into the pit more smoothly through the reserved gap in the platen. TPM uses parabolic transition linear interpolation in the joint coordinate system for trajectory planning. This planning method significantly reduces the amount of manipulation of the manipulator during the tree planting process. Under the planning method, the angular displacement of the six joints of the manipulator varies linearly, which effectively avoids the creation of odd shapes of the manipulator.
Fig. 8Angular displacement and angular velocity curve of six joints
Fig. 9Displacement component curve of the endeffector
4. Conclusions
This paper proposes an improved DH parameter method, analyzes the mathematical model of the actuator TPR in the TPM, solves the forward kinematic equation of the manipulator in the TPR by matrix transformation, and deduces the joints of the manipulator by using the separation of variables method. Inverse numerical solution. A linear programming method based on parabolic transition is used to plan the motion trajectory in Cartesian coordinate system and joint coordinate system. By comparing the planning results in the two coordinate systems, considering the structural characteristics of TPM, an appropriate planning method is selected, which provides a theoretical basis for the subsequent tree planting process of TPM.
This paper addresses some of the research in inverse kinematics and path planning for TPRs, breaking through the limitations of industrial robots concentrated in unconstrained workspaces with homogeneous objects, and successfully combining the capabilities of TPRs with the environment in which it operates in forests. Tree planting environments are known to be complex.
References

Z. Liqiu, A. Juan, Z. Ronghao, and M. Hairong, “Trajectory planning and simulation of industrial robot based on MATLAB and RobotStudio,” in 2021 IEEE 4th International Conference on Electronics Technology (ICET), pp. 910–914, May 2021, https://doi.org/10.1109/icet51757.2021.9451021

F. Sun, Y. Chen, Y. Wu, L. Li, and X. Ren, “Motion planning and cooperative manipulation for mobile robots with dual arms,” IEEE Transactions on Emerging Topics in Computational Intelligence, Vol. 6, No. 6, pp. 1345–1356, Dec. 2022, https://doi.org/10.1109/tetci.2022.3146387

Y. Bai, M. Luo, and F. Pang, “An algorithm for solving robot inverse kinematics based on FOA optimized BP neural network,” Applied Sciences, Vol. 11, No. 15, p. 7129, Aug. 2021, https://doi.org/10.3390/app11157129

W. Liu, C. Zhao, Y. Liu, H. Wang, W. Zhao, and H. Zhang, “Sim2real kinematics modeling of industrial robots based on FPGAacceleration,” Robotics and ComputerIntegrated Manufacturing, Vol. 77, p. 102350, Oct. 2022, https://doi.org/10.1016/j.rcim.2022.102350

H. Ren and P. BenTzvi, “Learning inverse kinematics and dynamics of a robotic manipulator using generative adversarial networks,” Robotics and Autonomous Systems, Vol. 124, p. 103386, Feb. 2020, https://doi.org/10.1016/j.robot.2019.103386

S. Dereli and R. Köker, “Calculation of the inverse kinematics solution of the 7DOF redundant robot manipulator by the firefly algorithm and statistical analysis of the results in terms of speed and accuracy,” Inverse Problems in Science and Engineering, Vol. 28, No. 5, pp. 601–613, May 2020, https://doi.org/10.1080/17415977.2019.1602124

L. Wu, Q. Liu, X. Tian, J. Zhang, and W. Xiao, “A new improved fruit fly optimization algorithm IAFOA and its application to solve engineering optimization problems,” KnowledgeBased Systems, Vol. 144, pp. 153–173, Mar. 2018, https://doi.org/10.1016/j.knosys.2017.12.031

C. LopezFranco, J. HernandezBarragan, A. Y. Alanis, and N. AranaDaniel, “A soft computing approach for inverse kinematics of robot manipulators,” Engineering Applications of Artificial Intelligence, Vol. 74, pp. 104–120, Sep. 2018, https://doi.org/10.1016/j.engappai.2018.06.001

I. Pikalov, E. Spirin, M. Saramud, and M. Kubrikov, “Vector model for solving the inverse kinematics problem in the system of external adaptive control of robotic manipulators,” Mechanism and Machine Theory, Vol. 174, p. 104912, Aug. 2022, https://doi.org/10.1016/j.mechmachtheory.2022.104912

L. Ye, J. Duan, Z. Yang, X. Zou, M. Chen, and S. Zhang, “Collisionfree motion planning for the litchipicking robot,” Computers and Electronics in Agriculture, Vol. 185, p. 106151, Jun. 2021, https://doi.org/10.1016/j.compag.2021.106151

Y. Wang et al., “Rapid citrus harvesting motion planning with preharvesting point and quadtree,” Computers and Electronics in Agriculture, Vol. 202, p. 107348, Nov. 2022, https://doi.org/10.1016/j.compag.2022.107348

H. Ye, D. Wang, J. Wu, Y. Yue, and Y. Zhou, “Forward and inverse kinematics of a 5DOF hybrid robot for composite material machining,” Robotics and ComputerIntegrated Manufacturing, Vol. 65, p. 101961, Oct. 2020, https://doi.org/10.1016/j.rcim.2020.101961

P. Sudhakara, V. Ganapathy, B. Priyadharshini, and K. Sundaran, “Obstacle avoidance and navigation planning of a wheeled mobile robot using amended artificial potential field method,” Procedia Computer Science, Vol. 133, pp. 998–1004, 2018, https://doi.org/10.1016/j.procs.2018.07.076

M. Kang et al., “Divisionmerge based inverse kinematics for multiDOFs humanoid robots in unstructured environments,” Computers and Electronics in Agriculture, Vol. 198, p. 107090, Jul. 2022, https://doi.org/10.1016/j.compag.2022.107090

H. Wang, Y. Lin, X. Xu, Z. Chen, Z. Wu, and Y. Tang, “A study on longclose distance coordination control strategy for litchi picking,” Agronomy, Vol. 12, No. 7, p. 1520, Jun. 2022, https://doi.org/10.3390/agronomy12071520

A. Hentout, A. Maoudj, and M. Aouache, “A review of the literature on fuzzylogic approaches for collisionfree path planning of manipulator robots,” Artificial Intelligence Review, Sep. 2022, https://doi.org/10.1007/s10462022102577

M. Chen et al., “Threedimensional perception of orchard banana central stock enhanced by adaptive multivision technology,” Computers and Electronics in Agriculture, Vol. 174, p. 105508, Jul. 2020, https://doi.org/10.1016/j.compag.2020.105508

X. Cao, X. Zou, C. Jia, M. Chen, and Z. Zeng, “Rrtbased path planning for an intelligent litchipicking manipulator,” Computers and Electronics in Agriculture, Vol. 156, pp. 105–118, Jan. 2019, https://doi.org/10.1016/j.compag.2018.10.031

J. Wang, B. Li, and M. Q.H. Meng, “Kinematic constrained bidirectional RRT with efficient branch pruning for robot path planning,” Expert Systems with Applications, Vol. 170, p. 114541, May 2021, https://doi.org/10.1016/j.eswa.2020.114541

G. Lin, L. Zhu, J. Li, X. Zou, and Y. Tang, “Collisionfree path planning for a guavaharvesting robot based on recurrent deep reinforcement learning,” Computers and Electronics in Agriculture, Vol. 188, p. 106350, Sep. 2021, https://doi.org/10.1016/j.compag.2021.106350

H. Wang, C. J. Hohimer, S. Bhusal, M. Karkee, C. Mo, and J. H. Miller, “Simulation as a tool in designing and evaluating a robotic apple harvesting system,” IFACPapersOnLine, Vol. 51, No. 17, pp. 135–140, 2018, https://doi.org/10.1016/j.ifacol.2018.08.076

Y. Li, W. Wei, Y. Gao, D. Wang, and Z. Fan, “Pqrrt*: an improved path planning algorithm for mobile robots,” Expert Systems with Applications, Vol. 152, p. 113425, Aug. 2020, https://doi.org/10.1016/j.eswa.2020.113425
About this article
The authors would like to thank the Northeast Forestry University (NEFU), Heilongjiang Institute of Technology (HLJIT), and Huaqiao University for their support. The research topic was supported by the 2022 Heilongjiang Provincial Key Research and Development Plan Key Area Special Project (SC2022ZX01A0052, Xigui Wang, NEFU), the Innovation Foundation for Doctoral Program of Forestry Engineering of Northeast Forestry University (Grant No. 000/4111410203, Jie Tang, NEFU), the Doctoral Research Startup Foundation Project of Heilongjiang Institute of Technology (Grant No. 2020BJ06, Yongmei Wang, HLJIT), the Natural Science Foundation Project of Heilongjiang Province (Grant No. LH2019E114, Baixue Fu, HLJIT), the Basic Scientific Research Business Expenses (Innovation Team Category) Project of Heilongjiang Institute of Engineering (Grant No. 2020CX02, Baixue Fu, HLJIT), the Special Project for Double FirstClassCultivation of Innovative Talents (Grant No. 000/41113102, Jiafu Ruan, NEFU), the Special Scientific Research Funds for Forest Nonprofit Industry (Grant No.201504508), the Youth Science Fund of Heilongjiang Institute of Technology (Grant No. 2015QJ02), and the Fundamental Research Funds for the Central Universities (Grant No. 2572016CB15).
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.
The authors declare that they have no conflict of interest.