In this paper, we study the polynomial optimization problem of multi-forms over the intersection of the multi-spheres and the nonnegative orthants. This class of problems is NP-hard in general, and includes the problem of finding the best nonnegative rank-one approximation of a given tensor. A Positivstellensatz is given for this class of polynomial optimization problems, based on which a globally convergent hierarchy of doubly nonnegative (DNN) relaxations is proposed. A (zero-th order) DNN relaxation method is applied to solve these problems, resulting in linear matrix optimization problems under both the positive semidefinite and nonnegative conic constraints. A worst case approximation bound is given for this relaxation method. Then, the recent solver SDPNAL+ is adopted to solve this class of matrix optimization problems. Typically, the DNN relaxations are tight, and hence the best nonnegative rank-one approximation of a tensor can be revealed frequently. Extensive numerical experiments show that this approach is quite promising.2010 Mathematics Subject Classification. 15A18; 15A42; 15A69; 90C22. Key words and phrases. Tensor, nonnegative rank-1 approximation, polynomial, multi-forms, doubly nonnegative semidefinite program, doubly nonnegative relaxation method.