When using quantitative photoacoustic tomography (q-PAT) reconstruction to recover the optical absorption coefficients of tissue, the commonly used diffusion equation has several limitations in the case of the objects that have small geometries and high-absorption or low-scattering areas. Furthermore, the conventional perturbation reconstruction strategy is unsatisfactory when the target tissue containing large heterogeneous features. We herein present a modified q-PAT implementation that employs the higher-order photon migration model achieving the tradeoff between mathematical rigidity and computational efficiency. Besides, a nonlinear iterative method is proposed to obtain the perturbations of optical absorption considering the updating of the sensitivity matrix in calculating the fluence perturbations. Consequently, the distribution of tissue optical properties can be recovered in a robust way even if the targets with high absorption are included. The proposed approach has been validated by simulation, phantom and in vivo experiments, exhibiting promising performances in image fidelity and quantitative feasibility for practical applications.